From 0c77eee16e2f4161c1ebc110b15ce6563d4a96c2 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 16 Apr 2014 11:13:22 +0000 Subject: Add fan beam support to FBP_CUDA --- cuda/2d/astra.cu | 91 ++++++++++++++++++++++++++++- cuda/2d/astra.h | 13 +++++ cuda/3d/fdk.cu | 23 +++++++- cuda/3d/fdk.h | 8 +++ src/CudaFilteredBackProjectionAlgorithm.cpp | 30 ++++++++-- 5 files changed, 155 insertions(+), 10 deletions(-) diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index 1c2e623..e317f6c 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -33,6 +33,7 @@ $Id$ #include "par_fp.h" #include "fan_fp.h" #include "par_bp.h" +#include "fan_bp.h" #include "arith.h" #include "astra.h" @@ -43,6 +44,9 @@ $Id$ #include "../../include/astra/Logger.h" +// For fan beam FBP weighting +#include "../3d/fdk.h" + using namespace astraCUDA; using namespace std; @@ -61,8 +65,14 @@ public: float* angles; float* TOffsets; + float fOriginSourceDistance; + float fOriginDetectorDistance; + float fPixelSize; + bool bFanBeam; + bool bShortScan; + bool initialized; bool setStartReconstruction; @@ -152,9 +162,33 @@ bool AstraFBP::setProjectionGeometry(unsigned int iProjAngles, 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 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->bFanBeam = true; + pData->bShortScan = bShortScan; + return true; } + bool AstraFBP::setPixelSuperSampling(unsigned int iPixelSuperSampling) { if (pData->initialized) @@ -274,6 +308,34 @@ bool AstraFBP::run() 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... + + // 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); @@ -293,7 +355,34 @@ bool AstraFBP::run() } - ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets); + if (pData->bFanBeam) { + // TODO: TOffsets? + // TODO: Remove this code duplication with CudaReconstructionAlgorithm + SFanProjection* projs; + projs = new SFanProjection[pData->dims.iProjAngles]; + + float fSrcX0 = 0.0f; + float fSrcY0 = -pData->fOriginSourceDistance / pData->fPixelSize; + float fDetUX0 = pData->dims.fDetScale; + float fDetUY0 = 0.0f; + float fDetSX0 = pData->dims.iProjDets * fDetUX0 / -2.0f; + float fDetSY0 = pData->fOriginDetectorDistance / pData->fPixelSize; + +#define ROTATE0(name,i,alpha) do { projs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); projs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } while(0) + for (unsigned int i = 0; i < pData->dims.iProjAngles; ++i) { + ROTATE0(Src, i, pData->angles[i]); + ROTATE0(DetS, i, pData->angles[i]); + ROTATE0(DetU, i, pData->angles[i]); + } + +#undef ROTATE0 + ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, projs); + + delete[] projs; + + } else { + ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets); + } if(!ok) { return false; diff --git a/cuda/2d/astra.h b/cuda/2d/astra.h index 9e58301..42b15e8 100644 --- a/cuda/2d/astra.h +++ b/cuda/2d/astra.h @@ -68,6 +68,19 @@ public: unsigned int iProjDets, const float *pfAngles, float fDetSize = 1.0f); + // Set the projection angles and number of detector pixels per angle. + // pfAngles must be a float array of length iProjAngles. + // fDetSize indicates the size of a detector pixel compared to a + // volume pixel edge. + // + // pfAngles will only be read from during this call. + bool setFanGeometry(unsigned int iProjAngles, + unsigned int iProjDets, + const float *pfAngles, + float fOriginSourceDistance, + float fOriginDetectorDistance, + float fDetSize = 1.0f, + bool bShortScan = false); // Set linear supersampling factor for the BP. // (The number of rays is the square of this) diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 883187d..fd5a8a2 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -298,8 +298,7 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un // Perform the FDK pre-weighting and filtering -bool FDK_Filter(cudaPitchedPtr D_projData, - cufftComplex * D_filter, +bool FDK_PreWeight(cudaPitchedPtr D_projData, float fSrcOrigin, float fDetOrigin, float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize, bool bShortScan, @@ -335,13 +334,22 @@ bool FDK_Filter(cudaPitchedPtr D_projData, } cudaTextForceKernelsCompletion(); + return true; +} +bool FDK_Filter(cudaPitchedPtr D_projData, + cufftComplex * D_filter, + float fSrcOrigin, float fDetOrigin, + float fSrcZ, float fDetZ, + float fDetUSize, float fDetVSize, bool bShortScan, + const SDimensions3D& dims, const float* angles) +{ // The filtering is a regular ramp filter per detector line. int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU); int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount); - + int projPitch = D_projData.pitch/sizeof(float); // We process one sinogram at a time. @@ -390,11 +398,18 @@ bool FDK(cudaPitchedPtr D_volumeData, int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU); int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount); + ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin, + fSrcZ, fDetZ, fDetUSize, fDetVSize, + bShortScan, dims, angles); + if (!ok) + return false; + cufftComplex *pHostFilter = new cufftComplex[dims.iProjAngles * iHalfFFTSize]; memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize); genFilter(FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); + allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_filter); uploadComplexArrayToDevice(dims.iProjAngles, iHalfFFTSize, pHostFilter, D_filter); @@ -403,6 +418,8 @@ bool FDK(cudaPitchedPtr D_volumeData, // Perform filtering + + ok = FDK_Filter(D_projData, D_filter, fSrcOrigin, fDetOrigin, fSrcZ, fDetZ, fDetUSize, fDetVSize, bShortScan, dims, angles); diff --git a/cuda/3d/fdk.h b/cuda/3d/fdk.h index 5443b19..c9123b9 100644 --- a/cuda/3d/fdk.h +++ b/cuda/3d/fdk.h @@ -29,8 +29,16 @@ $Id$ #ifndef _CUDA_FDK_H #define _CUDA_FDK_H +#include "dims3d.h" + namespace astraCUDA3d { +bool FDK_PreWeight(cudaPitchedPtr D_projData, + float fSrcOrigin, float fDetOrigin, + float fSrcZ, float fDetZ, + float fDetUSize, float fDetVSize, bool bShortScan, + const SDimensions3D& dims, const float* angles); + bool FDK(cudaPitchedPtr D_volumeData, cudaPitchedPtr D_projData, float fSrcOrigin, float fDetOrigin, diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 3656f3f..c53ac2d 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -27,6 +27,7 @@ $Id$ */ #include +#include #include #include @@ -226,7 +227,8 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) if (!m_bAstraFBPInit) { const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - const CParallelProjectionGeometry2D& projgeom = *dynamic_cast(m_pSinogram->getGeometry()); + const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast(m_pSinogram->getGeometry()); + const CFanFlatProjectionGeometry2D* fanprojgeom = dynamic_cast(m_pSinogram->getGeometry()); bool ok = true; @@ -235,10 +237,26 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) volgeom.getGridRowCount(), volgeom.getPixelLengthX()); // TODO: off-center geometry - ok &= m_pFBP->setProjectionGeometry(projgeom.getProjectionAngleCount(), - projgeom.getDetectorCount(), - projgeom.getProjectionAngles(), - projgeom.getDetectorWidth()); + int iDetectorCount; + if (parprojgeom) { + ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(), + parprojgeom->getDetectorCount(), + parprojgeom->getProjectionAngles(), + parprojgeom->getDetectorWidth()); + iDetectorCount = parprojgeom->getDetectorCount(); + } else if (fanprojgeom) { + ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(), + fanprojgeom->getDetectorCount(), + fanprojgeom->getProjectionAngles(), + fanprojgeom->getOriginSourceDistance(), + fanprojgeom->getOriginDetectorDistance(), + fanprojgeom->getDetectorWidth(), + false); // TODO: Support short-scan + + iDetectorCount = fanprojgeom->getDetectorCount(); + } else { + assert(false); + } ok &= m_pFBP->setPixelSuperSampling(m_iPixelSuperSampling); @@ -252,7 +270,7 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) ok &= m_pFBP->init(m_iGPUIndex); ASTRA_ASSERT(ok); - ok &= m_pFBP->setSinogram(m_pSinogram->getDataConst(), projgeom.getDetectorCount()); + ok &= m_pFBP->setSinogram(m_pSinogram->getDataConst(), iDetectorCount); ASTRA_ASSERT(ok); ok &= m_pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter); -- cgit v1.2.3