Merge branch 'consistent_scaling'

7 files changed, 32 insertions, 20 deletions

diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index 1319a87..822f746 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -1462,12 +1462,10 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 		Cuda3DProjectionKernel projKernel = ker3d_default;
 		int detectorSuperSampling = 1;
 		int voxelSuperSampling = 1;
-		bool densityWeighting = false;
 		if (projector) {
 			projKernel = projector->getProjectionKernel();
 			detectorSuperSampling = projector->getDetectorSuperSampling();
 			voxelSuperSampling = projector->getVoxelSuperSampling();
-			densityWeighting = projector->getDensityWeighting();
 		}
 
 		size_t inx, iny, inz;
@@ -1513,7 +1511,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP");
 
-			ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, voxelSuperSampling, densityWeighting);
+			ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, voxelSuperSampling);
 			if (!ok) ASTRA_ERROR("Error performing sub-BP");
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done");
 		}
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 88e235b..c1d3dc8 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_pReconstruction->getGeometry();
+	float fPixelArea = volGeom.getPixelArea();
+	ok &= pFBP->setReconstructionScale(1.0f/fPixelArea);
+	if (!ok) {
+		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set reconstruction scale");
+	}
 }
 
 
diff --git a/src/CudaProjector3D.cpp b/src/CudaProjector3D.cpp
index 3ea7043..e5c55cc 100644
--- a/src/CudaProjector3D.cpp
+++ b/src/CudaProjector3D.cpp
@@ -67,7 +67,6 @@ void CCudaProjector3D::_clear()
 	m_iVoxelSuperSampling = 1;
 	m_iDetectorSuperSampling = 1;
 	m_iGPUIndex = -1;
-	m_bDensityWeighting = false;
 }
 
 //----------------------------------------------------------------------------------------
@@ -132,13 +131,6 @@ bool CCudaProjector3D::initialize(const Config& _cfg)
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1);
 	CC.markOptionParsed("DetectorSuperSampling");
 
-	if (dynamic_cast<CConeProjectionGeometry3D*>(m_pProjectionGeometry) ||
-	    dynamic_cast<CConeVecProjectionGeometry3D*>(m_pProjectionGeometry))
-	{
-		m_bDensityWeighting = _cfg.self.getOptionBool("DensityWeighting", false);
-		CC.markOptionParsed("DensityWeighting");
-	}
-
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
 	CC.markOptionParsed("GPUIndex");
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 1e81390..6730cea 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -309,10 +309,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
 		m_bAlgoInit = true;
 	}
 
-	float fPixelSize = volgeom.getPixelLengthX();
-	float fSinogramScale = 1.0f/(fPixelSize*fPixelSize);
-
-	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/Features.cpp b/src/Features.cpp
index 9114131..09a3499 100644
--- a/src/Features.cpp
+++ b/src/Features.cpp
@@ -34,6 +34,12 @@ _AstraExport bool hasFeature(const std::string &flag) {
 	if (flag == "cuda") {
 		return cudaEnabled();
 	}
+	if (flag == "projectors_scaled_as_line_integrals") {
+		return true;
+	}
+	if (flag == "fan_cone_BP_density_weighting_by_default") {
+		return true;
+	}
 
 	return false;
 }
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index 423dc6c..6b4093d 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -167,6 +167,11 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 
 	m_filterConfig = getFilterConfigForAlgorithm(_cfg, this);
 
+	const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
+	if (!parprojgeom) {
+		ASTRA_ERROR("FBP currently only supports parallel projection geometries.");
+		return false;
+	}
 
 	// TODO: check that the angles are linearly spaced between 0 and pi
 
@@ -264,8 +269,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();
 }
diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp
index e09a3bc..806572f 100644
--- a/src/GeometryUtil2D.cpp
+++ b/src/GeometryUtil2D.cpp
@@ -28,6 +28,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/GeometryUtil2D.h"
 
 #include <cmath>
+#include <cstdio>
 
 namespace astra {
 
@@ -158,14 +159,16 @@ bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float
 	// project origin on detector line ( == project source on detector line)
 
 	double t = (- proj.fDetSX) * proj.fDetUX + (- proj.fDetSY) * proj.fDetUY;
+	t /= (proj.fDetUX * proj.fDetUX + proj.fDetUY * proj.fDetUY);
 
 	fOffset = (float)t - 0.5*iProjDets;
 
-	// TODO: CHECKME
 	fOriginDetector = sqrt((proj.fDetSX + t * proj.fDetUX)*(proj.fDetSX + t * proj.fDetUX) + (proj.fDetSY + t * proj.fDetUY)*(proj.fDetSY + t * proj.fDetUY));
 
-	//float fAngle = atan2(proj.fDetSX + t * proj.fDetUX - proj.fSrcX, proj.fDetSY + t * proj.fDetUY); // TODO: Fix order + sign
-	fAngle = atan2(proj.fDetUY, proj.fDetUX); // TODO: Check order + sign
+	fAngle = atan2(proj.fDetUY, proj.fDetUX);
+
+	//fprintf(stderr, "getFanParams: s = (%f,%f) d = (%f,%f) u = (%f,%f)\n", proj.fSrcX, proj.fSrcY, proj.fDetSX, proj.fDetSY, proj.fDetUX, proj.fDetUY);
+	//fprintf(stderr, "getFanParams: fOS = %f, fOD = %f, detsize = %f, offset = %f (t = %f), angle = %f\n", fOriginSource, fOriginDetector, fDetSize, fOffset, t, fAngle);
 
 	return true;
 }