Add fan beam support to FBP_CUDA

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 <astra/CudaFilteredBackProjectionAlgorithm.h>
+#include <astra/FanFlatProjectionGeometry2D.h>
 #include <boost/lexical_cast.hpp>
 #include <cstring>
 
@@ -226,7 +227,8 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
 	if (!m_bAstraFBPInit) {
 
 		const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
-		const CParallelProjectionGeometry2D& projgeom = *dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
+		const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
+		const CFanFlatProjectionGeometry2D* fanprojgeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(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);