Adjust FBP to line integral scaling

8 files changed, 31 insertions, 19 deletions

diff --git a/cuda/2d/algo.cu b/cuda/2d/algo.cu
index 11422ff..59ea984 100644
--- a/cuda/2d/algo.cu
+++ b/cuda/2d/algo.cu
@@ -153,6 +153,12 @@ bool ReconAlgo::setSuperSampling(int raysPerDet, int raysPerPixelDim)
 	return true;
 }
 
+bool ReconAlgo::setReconstructionScale(float fScale)
+{
+	fOutputScale *= fScale;
+	return true;
+}
+
 bool ReconAlgo::setVolumeMask(float* _D_maskData, unsigned int _maskPitch)
 {
 	assert(useVolumeMask);
@@ -242,7 +248,7 @@ bool ReconAlgo::allocateBuffers()
 	return true;
 }
 
-bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
                               const float* pfReconstruction, unsigned int iReconstructionPitch,
                               const float* pfVolMask, unsigned int iVolMaskPitch,
                               const float* pfSinoMask, unsigned int iSinoMaskPitch)
@@ -258,11 +264,6 @@ bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPit
 	if (!ok)
 		return false;
 
-	// rescale sinogram
-	if (fSinogramScale != 1.0f)
-		processSino<opMul>(D_sinoData, fSinogramScale,
-		                   sinoPitch, dims);
-
 	ok = copyVolumeToDevice(pfReconstruction, iReconstructionPitch,
 	                        dims,
 	                        D_volumeData, volumePitch);
diff --git a/cuda/2d/cgls.cu b/cuda/2d/cgls.cu
index b6a9fae..905b960 100644
--- a/cuda/2d/cgls.cu
+++ b/cuda/2d/cgls.cu
@@ -102,14 +102,14 @@ bool CGLS::setBuffers(float* _D_volumeData, unsigned int _volumePitch,
 	return true;
 }
 
-bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
                          const float* pfReconstruction, unsigned int iReconstructionPitch,
                          const float* pfVolMask, unsigned int iVolMaskPitch,
                          const float* pfSinoMask, unsigned int iSinoMaskPitch)
 {
 	sliceInitialized = false;
 
-	return ReconAlgo::copyDataToGPU(pfSinogram, iSinogramPitch, fSinogramScale, pfReconstruction, iReconstructionPitch, pfVolMask, iVolMaskPitch, pfSinoMask, iSinoMaskPitch);
+	return ReconAlgo::copyDataToGPU(pfSinogram, iSinogramPitch, pfReconstruction, iReconstructionPitch, pfVolMask, iVolMaskPitch, pfSinoMask, iSinoMaskPitch);
 }
 
 bool CGLS::iterate(unsigned int iterations)
diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu
index a5b8a7a..65f1d68 100644
--- a/cuda/2d/fbp.cu
+++ b/cuda/2d/fbp.cu
@@ -321,15 +321,14 @@ bool FBP::iterate(unsigned int iterations)
 	if (fanProjs) {
 		float fOutputScale = 1.0 / (/*fPixelSize * fPixelSize * fPixelSize * */ fFanDetSize * fFanDetSize);
 
-		ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale);
+		ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale * this->fOutputScale);
 
 	} else {
 		// scale by number of angles. For the fan-beam case, this is already
 		// handled by FDK_PreWeight
 		float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles;
-		//fOutputScale /= fDetSize * fDetSize;
 
-		ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale);
+		ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale * this->fOutputScale);
 	}
 	if(!ok)
 	{
diff --git a/include/astra/cuda/2d/algo.h b/include/astra/cuda/2d/algo.h
index 3fccdef..b6ec231 100644
--- a/include/astra/cuda/2d/algo.h
+++ b/include/astra/cuda/2d/algo.h
@@ -56,6 +56,10 @@ public:
 
 	bool setSuperSampling(int raysPerDet, int raysPerPixelDim);
 
+	// Scale the final reconstruction.
+	// May be called at any time after setGeometry and before iterate(). Multiple calls stack.
+	bool setReconstructionScale(float fScale);
+
 	virtual bool enableVolumeMask();
 	virtual bool enableSinogramMask();
 
@@ -88,8 +92,7 @@ public:
 	// sinogram, reconstruction, volume mask and sinogram mask in system RAM,
 	// respectively. The corresponding pitch variables give the pitches
 	// of these buffers, measured in floats.
-	// The sinogram is multiplied by fSinogramScale after uploading it.
-	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
 	                           const float* pfReconstruction, unsigned int iReconstructionPitch,
 	                           const float* pfVolMask, unsigned int iVolMaskPitch,
 	                           const float* pfSinoMask, unsigned int iSinoMaskPitch);
diff --git a/include/astra/cuda/2d/cgls.h b/include/astra/cuda/2d/cgls.h
index 375a425..a854a74 100644
--- a/include/astra/cuda/2d/cgls.h
+++ b/include/astra/cuda/2d/cgls.h
@@ -47,7 +47,7 @@ public:
 	virtual bool setBuffers(float* D_volumeData, unsigned int volumePitch,
 	                        float* D_projData, unsigned int projPitch);
 
-	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
 	                           const float* pfReconstruction, unsigned int iReconstructionPitch,
 	                           const float* pfVolMask, unsigned int iVolMaskPitch,
 	                           const float* pfSinoMask, unsigned int iSinoMaskPitch);
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 88e235b..69140fc 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -151,6 +151,13 @@ void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm()
 	if (!ok) {
 		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set short-scan mode");
 	}
+
+	const CVolumeGeometry2D& volGeom = *m_pProjector->getVolumeGeometry();
+	float fPixelArea = volGeom.getPixelArea();
+	ok &= pFBP->setReconstructionScale(1.0f/fPixelArea);
+	if (!ok) {
+		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set reconstruction scale");
+	}
 }
 
 
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 939a026..6730cea 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -309,9 +309,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
 		m_bAlgoInit = true;
 	}
 
-	float fSinogramScale = 1.0f;
-
-	ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(), fSinogramScale,
+	ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(),
 	                            m_pReconstruction->getDataConst(), volgeom.getGridColCount(),
 	                            m_bUseReconstructionMask ? m_pReconstructionMask->getDataConst() : 0, volgeom.getGridColCount(),
 	                            m_bUseSinogramMask ? m_pSinogramMask->getDataConst() : 0, m_pSinogram->getGeometry()->getDetectorCount());
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index 423dc6c..95bef3c 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -264,8 +264,12 @@ void CFilteredBackProjectionAlgorithm::run(int _iNrIterations)
 	            DefaultBPPolicy(m_pReconstruction, &filteredSinogram));
 
 	// Scale data
-	int iAngleCount = m_pProjector->getProjectionGeometry()->getProjectionAngleCount();
-	(*m_pReconstruction) *= (PI/2)/iAngleCount;
+	const CVolumeGeometry2D& volGeom = *m_pProjector->getVolumeGeometry();
+	const CProjectionGeometry2D& projGeom = *m_pProjector->getProjectionGeometry();
+
+	int iAngleCount = projGeom.getProjectionAngleCount();
+	float fPixelArea = volGeom.getPixelArea();
+	(*m_pReconstruction) *= PI/(2*iAngleCount*fPixelArea);
 
 	m_pReconstruction->updateStatistics();
 }