Start supporting flexible 2D volume geometry in Cuda FBP too

3 files changed, 61 insertions, 39 deletions

diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index 087905d..be273eb 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -68,6 +68,7 @@ public:
 	SDimensions dims;
 	float* angles;
 	float* TOffsets;
+	astraCUDA::SFanProjection* fanProjections;
 
 	float fOriginSourceDistance;
 	float fOriginDetectorDistance;
@@ -94,6 +95,8 @@ AstraFBP::AstraFBP()
 	pData = new AstraFBP_internal();
 
 	pData->angles = 0;
+	pData->fanProjections = 0;
+	pData->TOffsets = 0;
 	pData->D_sinoData = 0;
 	pData->D_volumeData = 0;
 
@@ -117,6 +120,9 @@ AstraFBP::~AstraFBP()
 	delete[] pData->TOffsets;
 	pData->TOffsets = 0;
 
+	delete[] pData->fanProjections;
+	pData->fanProjections = 0;
+
 	cudaFree(pData->D_sinoData);
 	pData->D_sinoData = 0;
 
@@ -173,6 +179,7 @@ bool AstraFBP::setProjectionGeometry(unsigned int iProjAngles,
 
 bool AstraFBP::setFanGeometry(unsigned int iProjAngles,
                               unsigned int iProjDets,
+                              const astraCUDA::SFanProjection *fanProjs,
                               const float* pfAngles,
                               float fOriginSourceDistance,
                               float fOriginDetectorDistance,
@@ -186,6 +193,9 @@ bool AstraFBP::setFanGeometry(unsigned int iProjAngles,
 	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;
 
@@ -314,7 +324,7 @@ bool AstraFBP::run()
 		// Call FDK_PreWeight to handle fan beam geometry. We treat
 		// this as a cone beam setup of a single slice:
 
-		// TODO: TOffsets...
+		// TODO: TOffsets affects this preweighting...
 
 		// We create a fake cudaPitchedPtr
 		cudaPitchedPtr tmp;
@@ -358,29 +368,7 @@ bool AstraFBP::run()
 	}
 
 	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;
+		ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->fanProjections);
 
 	} else {
 		ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets);
diff --git a/cuda/2d/astra.h b/cuda/2d/astra.h
index d2d7059..efbac1d 100644
--- a/cuda/2d/astra.h
+++ b/cuda/2d/astra.h
@@ -78,9 +78,10 @@ public:
 	// fDetSize indicates the size of a detector pixel compared to a
 	// volume pixel edge.
 	//
-	// pfAngles will only be read from during this call.
+	// pfAngles, fanProjs will only be read from during this call.
 	bool setFanGeometry(unsigned int iProjAngles,
 	                    unsigned int iProjDets,
+	                    const astraCUDA::SFanProjection *fanProjs,
 	                    const float *pfAngles,
 	                    float fOriginSourceDistance,
 	                    float fOriginDetectorDistance,
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 39740c8..63f4f40 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -32,6 +32,7 @@ $Id$
 #include <cstring>
 
 #include "astra/AstraObjectManager.h"
+#include "../cuda/2d/astra.h"
 
 #include <astra/Logger.h>
 
@@ -241,28 +242,65 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
 
 		bool ok = true;
 
-		// TODO: off-center geometry, non-square pixels
+		// TODO: non-square pixels?
 		ok &= m_pFBP->setReconstructionGeometry(volgeom.getGridColCount(),
 		                                         volgeom.getGridRowCount(),
 		                                         volgeom.getPixelLengthX());
-		// TODO: off-center geometry
 		int iDetectorCount;
 		if (parprojgeom) {
+
+			float *offsets, *angles, detSize, outputScale;
+
+			ok = convertAstraGeometry(&volgeom, parprojgeom, offsets, angles, detSize, outputScale);
+
+
 			ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(),
 			                                     parprojgeom->getDetectorCount(),
-			                                     parprojgeom->getProjectionAngles(),
+			                                     angles,
 			                                     parprojgeom->getDetectorWidth());
 			iDetectorCount = parprojgeom->getDetectorCount();
+
+			// TODO: Are detSize and outputScale handled correctly?
+
+			if (offsets)
+				ok &= m_pFBP->setTOffsets(offsets);
+			ASTRA_ASSERT(ok);
+
+			delete[] offsets;
+			delete[] angles;
+
 		} else if (fanprojgeom) {
+
+			astraCUDA::SFanProjection* projs;
+			float outputScale;
+
+			// FIXME: Implement this, and clean up the interface to AstraFBP.
+			if (abs(volgeom.getWindowMinX() + volgeom.getWindowMaxX()) > 0.00001 * volgeom.getPixelLengthX()) {
+				// Off-center volume geometry isn't supported yet
+				ASTRA_ASSERT(false);
+			}
+			if (abs(volgeom.getWindowMinY() + volgeom.getWindowMaxY()) > 0.00001 * volgeom.getPixelLengthY()) {
+				// Off-center volume geometry isn't supported yet
+				ASTRA_ASSERT(false);
+			}
+
+			ok = convertAstraGeometry(&volgeom, fanprojgeom, projs, outputScale);
+
+			// CHECKME: outputScale?
+
 			ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(),
-			                                     fanprojgeom->getDetectorCount(),
-			                                     fanprojgeom->getProjectionAngles(),
-			                                     fanprojgeom->getOriginSourceDistance(),
-			                                     fanprojgeom->getOriginDetectorDistance(),
-			                                     fanprojgeom->getDetectorWidth(),
-			                                     m_bShortScan);
+			                             fanprojgeom->getDetectorCount(),
+			                             projs,
+			                             fanprojgeom->getProjectionAngles(),
+			                             fanprojgeom->getOriginSourceDistance(),
+			                             fanprojgeom->getOriginDetectorDistance(),
+
+		                                     fanprojgeom->getDetectorWidth(),
+			                             m_bShortScan);
 
 			iDetectorCount = fanprojgeom->getDetectorCount();
+
+			delete[] projs;
 		} else {
 			assert(false);
 		}
@@ -271,11 +309,6 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
 
 		ASTRA_ASSERT(ok);
 
-		const float *pfTOffsets = m_pSinogram->getGeometry()->getExtraDetectorOffset();
-		if (pfTOffsets)
-			ok &= m_pFBP->setTOffsets(pfTOffsets);
-		ASTRA_ASSERT(ok);
-
 		ok &= m_pFBP->init(m_iGPUIndex);
 		ASTRA_ASSERT(ok);