diff options
Diffstat (limited to 'cuda/2d/astra.cu')
| -rw-r--r-- | cuda/2d/astra.cu | 726 |
1 files changed, 116 insertions, 610 deletions
diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index c1e6566..b39e0a9 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -42,8 +42,10 @@ $Id$ #include <fstream> #include <cuda.h> +#include "../../include/astra/GeometryUtil2D.h" #include "../../include/astra/VolumeGeometry2D.h" #include "../../include/astra/ParallelProjectionGeometry2D.h" +#include "../../include/astra/ParallelVecProjectionGeometry2D.h" #include "../../include/astra/FanFlatProjectionGeometry2D.h" #include "../../include/astra/FanFlatVecProjectionGeometry2D.h" @@ -64,515 +66,6 @@ enum CUDAProjectionType { }; -class AstraFBP_internal { -public: - SDimensions dims; - float* angles; - float* TOffsets; - astraCUDA::SFanProjection* fanProjections; - - float fOriginSourceDistance; - float fOriginDetectorDistance; - - float fPixelSize; - - bool bFanBeam; - bool bShortScan; - - bool initialized; - bool setStartReconstruction; - - float* D_sinoData; - unsigned int sinoPitch; - - float* D_volumeData; - unsigned int volumePitch; - - cufftComplex * m_pDevFilter; -}; - -AstraFBP::AstraFBP() -{ - pData = new AstraFBP_internal(); - - pData->angles = 0; - pData->fanProjections = 0; - pData->TOffsets = 0; - pData->D_sinoData = 0; - pData->D_volumeData = 0; - - pData->dims.iVolWidth = 0; - pData->dims.iProjAngles = 0; - pData->dims.fDetScale = 1.0f; - pData->dims.iRaysPerDet = 1; - pData->dims.iRaysPerPixelDim = 1; - - pData->initialized = false; - pData->setStartReconstruction = false; - - pData->m_pDevFilter = NULL; -} - -AstraFBP::~AstraFBP() -{ - delete[] pData->angles; - pData->angles = 0; - - delete[] pData->TOffsets; - pData->TOffsets = 0; - - delete[] pData->fanProjections; - pData->fanProjections = 0; - - cudaFree(pData->D_sinoData); - pData->D_sinoData = 0; - - cudaFree(pData->D_volumeData); - pData->D_volumeData = 0; - - if(pData->m_pDevFilter != NULL) - { - freeComplexOnDevice(pData->m_pDevFilter); - pData->m_pDevFilter = NULL; - } - - delete pData; - pData = 0; -} - -bool AstraFBP::setReconstructionGeometry(unsigned int iVolWidth, - unsigned int iVolHeight, - float fPixelSize) -{ - if (pData->initialized) - return false; - - pData->dims.iVolWidth = iVolWidth; - pData->dims.iVolHeight = iVolHeight; - - pData->fPixelSize = fPixelSize; - - return (iVolWidth > 0 && iVolHeight > 0 && fPixelSize > 0.0f); -} - -bool AstraFBP::setProjectionGeometry(unsigned int iProjAngles, - unsigned int iProjDets, - const float* pfAngles, - float fDetSize) -{ - if (pData->initialized) - return false; - - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjDets = iProjDets; - pData->dims.fDetScale = fDetSize / pData->fPixelSize; - - if (iProjAngles == 0 || iProjDets == 0 || pfAngles == 0) - return false; - - pData->angles = new float[iProjAngles]; - memcpy(pData->angles, pfAngles, iProjAngles * sizeof(pfAngles[0])); - - pData->bFanBeam = false; - - return true; -} - -bool AstraFBP::setFanGeometry(unsigned int iProjAngles, - unsigned int iProjDets, - const astraCUDA::SFanProjection *fanProjs, - const float* pfAngles, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetSize, - bool bShortScan) -{ - // Slightly abusing setProjectionGeometry for this... - if (!setProjectionGeometry(iProjAngles, iProjDets, pfAngles, fDetSize)) - return false; - - pData->fOriginSourceDistance = fOriginSourceDistance; - pData->fOriginDetectorDistance = fOriginDetectorDistance; - - pData->fanProjections = new astraCUDA::SFanProjection[iProjAngles]; - memcpy(pData->fanProjections, fanProjs, iProjAngles * sizeof(fanProjs[0])); - - pData->bFanBeam = true; - pData->bShortScan = bShortScan; - - return true; -} - - -bool AstraFBP::setPixelSuperSampling(unsigned int iPixelSuperSampling) -{ - if (pData->initialized) - return false; - - if (iPixelSuperSampling == 0) - return false; - - pData->dims.iRaysPerPixelDim = iPixelSuperSampling; - - return true; -} - - -bool AstraFBP::setTOffsets(const float* pfTOffsets) -{ - if (pData->initialized) - return false; - - if (pfTOffsets == 0) - return false; - - pData->TOffsets = new float[pData->dims.iProjAngles]; - memcpy(pData->TOffsets, pfTOffsets, pData->dims.iProjAngles * sizeof(pfTOffsets[0])); - - return true; -} - -bool AstraFBP::init(int iGPUIndex) -{ - if (pData->initialized) - { - return false; - } - - if (pData->dims.iProjAngles == 0 || pData->dims.iVolWidth == 0) - { - return false; - } - - if (iGPUIndex != -1) { - cudaSetDevice(iGPUIndex); - cudaError_t err = cudaGetLastError(); - - // Ignore errors caused by calling cudaSetDevice multiple times - if (err != cudaSuccess && err != cudaErrorSetOnActiveProcess) - { - return false; - } - } - - bool ok = allocateVolumeData(pData->D_volumeData, pData->volumePitch, pData->dims); - if (!ok) - { - return false; - } - - ok = allocateProjectionData(pData->D_sinoData, pData->sinoPitch, pData->dims); - if (!ok) - { - cudaFree(pData->D_volumeData); - pData->D_volumeData = 0; - return false; - } - - pData->initialized = true; - - return true; -} - -bool AstraFBP::setSinogram(const float* pfSinogram, - unsigned int iSinogramPitch) -{ - if (!pData->initialized) - return false; - if (!pfSinogram) - return false; - - bool ok = copySinogramToDevice(pfSinogram, iSinogramPitch, - pData->dims, - pData->D_sinoData, pData->sinoPitch); - if (!ok) - return false; - - // rescale sinogram to adjust for pixel size - processSino<opMul>(pData->D_sinoData, - 1.0f/(pData->fPixelSize*pData->fPixelSize), - pData->sinoPitch, pData->dims); - - pData->setStartReconstruction = false; - - return true; -} - -static int calcNextPowerOfTwo(int _iValue) -{ - int iOutput = 1; - - while(iOutput < _iValue) - { - iOutput *= 2; - } - - return iOutput; -} - -bool AstraFBP::run() -{ - if (!pData->initialized) - { - return false; - } - - zeroVolumeData(pData->D_volumeData, pData->volumePitch, pData->dims); - - bool ok = false; - - if (pData->bFanBeam) { - // Call FDK_PreWeight to handle fan beam geometry. We treat - // this as a cone beam setup of a single slice: - - // TODO: TOffsets affects this preweighting... - - // We create a fake cudaPitchedPtr - cudaPitchedPtr tmp; - tmp.ptr = pData->D_sinoData; - tmp.pitch = pData->sinoPitch * sizeof(float); - tmp.xsize = pData->dims.iProjDets; - tmp.ysize = pData->dims.iProjAngles; - // and a fake Dimensions3D - astraCUDA3d::SDimensions3D dims3d; - dims3d.iVolX = pData->dims.iVolWidth; - dims3d.iVolY = pData->dims.iVolHeight; - dims3d.iVolZ = 1; - dims3d.iProjAngles = pData->dims.iProjAngles; - dims3d.iProjU = pData->dims.iProjDets; - dims3d.iProjV = 1; - dims3d.iRaysPerDetDim = dims3d.iRaysPerVoxelDim = 1; - - astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance, - pData->fOriginDetectorDistance, 0.0f, 0.0f, - pData->dims.fDetScale, 1.0f, // TODO: Are these correct? - pData->bShortScan, dims3d, pData->angles); - } - - if (pData->m_pDevFilter) { - - int iFFTRealDetCount = calcNextPowerOfTwo(2 * pData->dims.iProjDets); - int iFFTFourDetCount = calcFFTFourSize(iFFTRealDetCount); - - cufftComplex * pDevComplexSinogram = NULL; - - allocateComplexOnDevice(pData->dims.iProjAngles, iFFTFourDetCount, &pDevComplexSinogram); - - runCudaFFT(pData->dims.iProjAngles, pData->D_sinoData, pData->sinoPitch, pData->dims.iProjDets, iFFTRealDetCount, iFFTFourDetCount, pDevComplexSinogram); - - applyFilter(pData->dims.iProjAngles, iFFTFourDetCount, pDevComplexSinogram, pData->m_pDevFilter); - - runCudaIFFT(pData->dims.iProjAngles, pDevComplexSinogram, pData->D_sinoData, pData->sinoPitch, pData->dims.iProjDets, iFFTRealDetCount, iFFTFourDetCount); - - freeComplexOnDevice(pDevComplexSinogram); - - } - - float fOutputScale = (M_PI / 2.0f) / (float)pData->dims.iProjAngles; - - if (pData->bFanBeam) { - ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->fanProjections, fOutputScale); - - } else { - ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets, fOutputScale); - } - if(!ok) - { - return false; - } - - return true; -} - -bool AstraFBP::getReconstruction(float* pfReconstruction, unsigned int iReconstructionPitch) const -{ - if (!pData->initialized) - return false; - - bool ok = copyVolumeFromDevice(pfReconstruction, iReconstructionPitch, - pData->dims, - pData->D_volumeData, pData->volumePitch); - if (!ok) - return false; - - return true; -} - -int AstraFBP::calcFourierFilterSize(int _iDetectorCount) -{ - int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount); - int iFreqBinCount = calcFFTFourSize(iFFTRealDetCount); - - // CHECKME: Matlab makes this at least 64. Do we also need to? - return iFreqBinCount; -} - -bool AstraFBP::setFilter(E_FBPFILTER _eFilter, const float * _pfHostFilter /* = NULL */, int _iFilterWidth /* = 0 */, float _fD /* = 1.0f */, float _fFilterParameter /* = -1.0f */) -{ - if(pData->m_pDevFilter != 0) - { - freeComplexOnDevice(pData->m_pDevFilter); - pData->m_pDevFilter = 0; - } - - if (_eFilter == FILTER_NONE) - return true; // leave pData->m_pDevFilter set to 0 - - - int iFFTRealDetCount = calcNextPowerOfTwo(2 * pData->dims.iProjDets); - int iFreqBinCount = calcFFTFourSize(iFFTRealDetCount); - - cufftComplex * pHostFilter = new cufftComplex[pData->dims.iProjAngles * iFreqBinCount]; - memset(pHostFilter, 0, sizeof(cufftComplex) * pData->dims.iProjAngles * iFreqBinCount); - - allocateComplexOnDevice(pData->dims.iProjAngles, iFreqBinCount, &(pData->m_pDevFilter)); - - switch(_eFilter) - { - case FILTER_NONE: - // handled above - break; - case FILTER_RAMLAK: - case FILTER_SHEPPLOGAN: - case FILTER_COSINE: - case FILTER_HAMMING: - case FILTER_HANN: - case FILTER_TUKEY: - case FILTER_LANCZOS: - case FILTER_TRIANGULAR: - case FILTER_GAUSSIAN: - case FILTER_BARTLETTHANN: - case FILTER_BLACKMAN: - case FILTER_NUTTALL: - case FILTER_BLACKMANHARRIS: - case FILTER_BLACKMANNUTTALL: - case FILTER_FLATTOP: - case FILTER_PARZEN: - { - genFilter(_eFilter, _fD, pData->dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); - uploadComplexArrayToDevice(pData->dims.iProjAngles, iFreqBinCount, pHostFilter, pData->m_pDevFilter); - - break; - } - case FILTER_PROJECTION: - { - // make sure the offered filter has the correct size - assert(_iFilterWidth == iFreqBinCount); - - for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) - { - float fValue = _pfHostFilter[iFreqBinIndex]; - - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; - } - } - uploadComplexArrayToDevice(pData->dims.iProjAngles, iFreqBinCount, pHostFilter, pData->m_pDevFilter); - break; - } - case FILTER_SINOGRAM: - { - // make sure the offered filter has the correct size - assert(_iFilterWidth == iFreqBinCount); - - for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) - { - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - float fValue = _pfHostFilter[iFreqBinIndex + iProjectionIndex * _iFilterWidth]; - - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; - } - } - uploadComplexArrayToDevice(pData->dims.iProjAngles, iFreqBinCount, pHostFilter, pData->m_pDevFilter); - break; - } - case FILTER_RPROJECTION: - { - int iProjectionCount = pData->dims.iProjAngles; - int iRealFilterElementCount = iProjectionCount * iFFTRealDetCount; - float * pfHostRealFilter = new float[iRealFilterElementCount]; - memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); - - int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); - int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; - int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; - - int iFilterShiftSize = _iFilterWidth / 2; - - for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) - { - int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; - float fValue = _pfHostFilter[iDetectorIndex]; - - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; - } - } - - float* pfDevRealFilter = NULL; - cudaMalloc((void **)&pfDevRealFilter, sizeof(float) * iRealFilterElementCount); // TODO: check for errors - cudaMemcpy(pfDevRealFilter, pfHostRealFilter, sizeof(float) * iRealFilterElementCount, cudaMemcpyHostToDevice); - delete[] pfHostRealFilter; - - runCudaFFT(iProjectionCount, pfDevRealFilter, iFFTRealDetCount, iFFTRealDetCount, iFFTRealDetCount, iFreqBinCount, pData->m_pDevFilter); - - cudaFree(pfDevRealFilter); - - break; - } - case FILTER_RSINOGRAM: - { - int iProjectionCount = pData->dims.iProjAngles; - int iRealFilterElementCount = iProjectionCount * iFFTRealDetCount; - float* pfHostRealFilter = new float[iRealFilterElementCount]; - memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); - - int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); - int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; - int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; - - int iFilterShiftSize = _iFilterWidth / 2; - - for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) - { - int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; - - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - float fValue = _pfHostFilter[iDetectorIndex + iProjectionIndex * _iFilterWidth]; - pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; - } - } - - float* pfDevRealFilter = NULL; - cudaMalloc((void **)&pfDevRealFilter, sizeof(float) * iRealFilterElementCount); // TODO: check for errors - cudaMemcpy(pfDevRealFilter, pfHostRealFilter, sizeof(float) * iRealFilterElementCount, cudaMemcpyHostToDevice); - delete[] pfHostRealFilter; - - runCudaFFT(iProjectionCount, pfDevRealFilter, iFFTRealDetCount, iFFTRealDetCount, iFFTRealDetCount, iFreqBinCount, pData->m_pDevFilter); - - cudaFree(pfDevRealFilter); - - break; - } - default: - { - ASTRA_ERROR("AstraFBP::setFilter: Unknown filter type requested"); - delete [] pHostFilter; - return false; - } - } - - delete [] pHostFilter; - - return true; -} - BPalgo::BPalgo() { @@ -617,18 +110,17 @@ float BPalgo::computeDiffNorm() bool astraCudaFP(const float* pfVolume, float* pfSinogram, unsigned int iVolWidth, unsigned int iVolHeight, unsigned int iProjAngles, unsigned int iProjDets, - const float *pfAngles, const float *pfOffsets, - float fDetSize, unsigned int iDetSuperSampling, + const SParProjection *pAngles, + unsigned int iDetSuperSampling, float fOutputScale, int iGPUIndex) { SDimensions dims; - if (iProjAngles == 0 || iProjDets == 0 || pfAngles == 0) + if (iProjAngles == 0 || iProjDets == 0 || pAngles == 0) return false; dims.iProjAngles = iProjAngles; dims.iProjDets = iProjDets; - dims.fDetScale = fDetSize; if (iDetSuperSampling == 0) return false; @@ -678,7 +170,7 @@ bool astraCudaFP(const float* pfVolume, float* pfSinogram, } zeroProjectionData(D_sinoData, sinoPitch, dims); - ok = FP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, pfAngles, pfOffsets, fOutputScale); + ok = FP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, pAngles, fOutputScale); if (!ok) { cudaFree(D_volumeData); cudaFree(D_sinoData); @@ -713,7 +205,6 @@ bool astraCudaFanFP(const float* pfVolume, float* pfSinogram, dims.iProjAngles = iProjAngles; dims.iProjDets = iProjDets; - dims.fDetScale = 1.0f; // TODO? if (iDetSuperSampling == 0) return false; @@ -789,118 +280,99 @@ bool astraCudaFanFP(const float* pfVolume, float* pfSinogram, } -bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, - const CParallelProjectionGeometry2D* pProjGeom, - float*& detectorOffsets, float*& projectionAngles, - float& detSize, float& outputScale) +// adjust pProjs to normalize volume geometry +template<typename ProjectionT> +static bool convertAstraGeometry_internal(const CVolumeGeometry2D* pVolGeom, + unsigned int iProjectionAngleCount, + ProjectionT*& pProjs, + float& fOutputScale) { - assert(pVolGeom); - assert(pProjGeom); - assert(pProjGeom->getProjectionAngles()); - + // TODO: Make EPS relative const float EPS = 0.00001f; - int nth = pProjGeom->getProjectionAngleCount(); - // Check if pixels are square if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthY()) > EPS) return false; + float dx = -(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; + float dy = -(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; - // Scale volume pixels to 1x1 - detSize = pProjGeom->getDetectorWidth() / pVolGeom->getPixelLengthX(); + float factor = 1.0f / pVolGeom->getPixelLengthX(); - // Copy angles - float *angles = new float[nth]; - for (int i = 0; i < nth; ++i) - angles[i] = pProjGeom->getProjectionAngles()[i]; - projectionAngles = angles; - - // Check if we need to translate - bool offCenter = false; - if (abs(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) > EPS || - abs(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) > EPS) - { - offCenter = true; + for (int i = 0; i < iProjectionAngleCount; ++i) { + // CHECKME: Order of scaling and translation + pProjs[i].translate(dx, dy); + pProjs[i].scale(factor); } + // CHECKME: Check factor + fOutputScale *= pVolGeom->getPixelLengthX() * pVolGeom->getPixelLengthY(); - // If there are existing detector offsets, or if we need to translate, - // we need to return offsets - if (offCenter) - { - float* offset = new float[nth]; + return true; +} - for (int i = 0; i < nth; ++i) - offset[i] = 0.0f; - if (offCenter) { - float dx = (pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; - float dy = (pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; - // CHECKME: Is d in pixels or in units? +bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CParallelProjectionGeometry2D* pProjGeom, + SParProjection*& pProjs, + float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionAngles()); - for (int i = 0; i < nth; ++i) { - float d = dx * cos(angles[i]) + dy * sin(angles[i]); - offset[i] += d; - } - } + int nth = pProjGeom->getProjectionAngleCount(); - // CHECKME: Order of scaling and translation + pProjs = genParProjections(nth, + pProjGeom->getDetectorCount(), + pProjGeom->getDetectorWidth(), + pProjGeom->getProjectionAngles(), 0); - // Scale volume pixels to 1x1 - for (int i = 0; i < nth; ++i) { - //offset[i] /= pVolGeom->getPixelLengthX(); - //offset[i] *= detSize; - } + bool ok; + fOutputScale = 1.0f; + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); - detectorOffsets = offset; - } else { - detectorOffsets = 0; + if (!ok) { + delete[] pProjs; + pProjs = 0; } - outputScale = pVolGeom->getPixelLengthX(); - outputScale *= outputScale; - - return true; + return ok; } -static void convertAstraGeometry_internal(const CVolumeGeometry2D* pVolGeom, - unsigned int iProjectionAngleCount, - astraCUDA::SFanProjection*& pProjs, - float& outputScale) +bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CParallelVecProjectionGeometry2D* pProjGeom, + SParProjection*& pProjs, + float& fOutputScale) { - // Translate - float dx = (pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; - float dy = (pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionVectors()); - for (int i = 0; i < iProjectionAngleCount; ++i) { - pProjs[i].fSrcX -= dx; - pProjs[i].fSrcY -= dy; - pProjs[i].fDetSX -= dx; - pProjs[i].fDetSY -= dy; + int nth = pProjGeom->getProjectionAngleCount(); + + pProjs = new SParProjection[nth]; + + for (int i = 0; i < nth; ++i) { + pProjs[i] = pProjGeom->getProjectionVectors()[i]; } - // CHECKME: Order of scaling and translation + bool ok; + fOutputScale = 1.0f; - // Scale - float factor = 1.0f / pVolGeom->getPixelLengthX(); - for (int i = 0; i < iProjectionAngleCount; ++i) { - pProjs[i].fSrcX *= factor; - pProjs[i].fSrcY *= factor; - pProjs[i].fDetSX *= factor; - pProjs[i].fDetSY *= factor; - pProjs[i].fDetUX *= factor; - pProjs[i].fDetUY *= factor; + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + if (!ok) { + delete[] pProjs; + pProjs = 0; } - // CHECKME: Check factor - outputScale = pVolGeom->getPixelLengthX(); -// outputScale *= outputScale; + return ok; } + bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, const CFanFlatProjectionGeometry2D* pProjGeom, astraCUDA::SFanProjection*& pProjs, @@ -910,6 +382,7 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, assert(pProjGeom); assert(pProjGeom->getProjectionAngles()); + // TODO: Make EPS relative const float EPS = 0.00001f; int nth = pProjGeom->getProjectionAngleCount(); @@ -928,23 +401,9 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, float fDetSize = pProjGeom->getDetectorWidth(); const float *pfAngles = pProjGeom->getProjectionAngles(); - pProjs = new SFanProjection[nth]; - - float fSrcX0 = 0.0f; - float fSrcY0 = -fOriginSourceDistance; - float fDetUX0 = fDetSize; - float fDetUY0 = 0.0f; - float fDetSX0 = pProjGeom->getDetectorCount() * fDetUX0 / -2.0f; - float fDetSY0 = fOriginDetectorDistance; - -#define ROTATE0(name,i,alpha) do { pProjs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); pProjs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } while(0) - for (int i = 0; i < nth; ++i) { - ROTATE0(Src, i, pfAngles[i]); - ROTATE0(DetS, i, pfAngles[i]); - ROTATE0(DetU, i, pfAngles[i]); - } - -#undef ROTATE0 + pProjs = genFanProjections(nth, pProjGeom->getDetectorCount(), + fOriginSourceDistance, fOriginDetectorDistance, + fDetSize, pfAngles); convertAstraGeometry_internal(pVolGeom, nth, pProjs, outputScale); @@ -961,6 +420,7 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, assert(pProjGeom); assert(pProjGeom->getProjectionVectors()); + // TODO: Make EPS relative const float EPS = 0.00001f; int nx = pVolGeom->getGridColCount(); @@ -983,6 +443,52 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, } +bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CProjectionGeometry2D* pProjGeom, + astraCUDA::SParProjection*& pParProjs, + astraCUDA::SFanProjection*& pFanProjs, + float& outputScale) +{ + const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<const CParallelProjectionGeometry2D*>(pProjGeom); + const CParallelVecProjectionGeometry2D* parVecProjGeom = dynamic_cast<const CParallelVecProjectionGeometry2D*>(pProjGeom); + const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<const CFanFlatProjectionGeometry2D*>(pProjGeom); + const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<const CFanFlatVecProjectionGeometry2D*>(pProjGeom); + + bool ok = false; + + if (parProjGeom) { + ok = convertAstraGeometry(pVolGeom, parProjGeom, pParProjs, outputScale); + } else if (parVecProjGeom) { + ok = convertAstraGeometry(pVolGeom, parVecProjGeom, pParProjs, outputScale); + } else if (fanProjGeom) { + ok = convertAstraGeometry(pVolGeom, fanProjGeom, pFanProjs, outputScale); + } else if (fanVecProjGeom) { + ok = convertAstraGeometry(pVolGeom, fanVecProjGeom, pFanProjs, outputScale); + } else { + ok = false; + } + + return ok; +} + +bool convertAstraGeometry_dims(const CVolumeGeometry2D* pVolGeom, + const CProjectionGeometry2D* pProjGeom, + SDimensions& dims) +{ + dims.iVolWidth = pVolGeom->getGridColCount(); + dims.iVolHeight = pVolGeom->getGridRowCount(); + + dims.iProjAngles = pProjGeom->getProjectionAngleCount(); + dims.iProjDets = pProjGeom->getDetectorCount(); + + dims.iRaysPerDet = 1; + dims.iRaysPerPixelDim = 1; + + return true; +} + + + } |