From 1dd79f23f783564719a52de7d9b54b17005c32d7 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 2 May 2014 09:20:54 +0000
Subject: Add SIRT-Weighted BP3D (par3d-only) for use in large BP

---
 src/CudaBackProjectionAlgorithm3D.cpp | 80 ++++++++++++++++++++++++++---------
 1 file changed, 60 insertions(+), 20 deletions(-)

(limited to 'src')

diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp
index b096756..891d7fe 100644
--- a/src/CudaBackProjectionAlgorithm3D.cpp
+++ b/src/CudaBackProjectionAlgorithm3D.cpp
@@ -53,6 +53,7 @@ CCudaBackProjectionAlgorithm3D::CCudaBackProjectionAlgorithm3D()
 	m_bIsInitialized = false;
 	m_iGPUIndex = -1;
 	m_iVoxelSuperSampling = 1;
+	m_bSIRTWeighting = false;
 }
 
 //----------------------------------------------------------------------------------------
@@ -106,6 +107,17 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg)
 	m_iVoxelSuperSampling = (int)_cfg.self->getOptionNumerical("VoxelSuperSampling", 1);
 	CC.markOptionParsed("VoxelSuperSampling");
 
+	CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
+	ASTRA_ASSERT(pSinoMem);
+	const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry();
+const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom);
+	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
+	if (parvec3dgeom || par3dgeom) {
+		// This option is only supported for Par3D currently
+		m_bSIRTWeighting = _cfg.self->getOptionBool("SIRTWeighting", false);
+		CC.markOptionParsed("SIRTWeighting");
+	}
+
 	// success
 	m_bIsInitialized = _check();
 	return m_bIsInitialized;
@@ -181,27 +193,55 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 		                conegeom->getProjectionAngles(),
 		                m_iGPUIndex, m_iVoxelSuperSampling);
 	} else if (par3dgeom) {
-		astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(),
-		                 volgeom.getGridColCount(),
-		                 volgeom.getGridRowCount(),
-		                 volgeom.getGridSliceCount(),
-		                 par3dgeom->getProjectionCount(),
-		                 par3dgeom->getDetectorColCount(),
-		                 par3dgeom->getDetectorRowCount(),
-		                 par3dgeom->getDetectorSpacingX(),
-		                 par3dgeom->getDetectorSpacingY(),
-		                 par3dgeom->getProjectionAngles(),
-		                 m_iGPUIndex, m_iVoxelSuperSampling);
+		if (!m_bSIRTWeighting) {
+			astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(),
+			                 volgeom.getGridColCount(),
+			                 volgeom.getGridRowCount(),
+			                 volgeom.getGridSliceCount(),
+			                 par3dgeom->getProjectionCount(),
+			                 par3dgeom->getDetectorColCount(),
+			                 par3dgeom->getDetectorRowCount(),
+			                 par3dgeom->getDetectorSpacingX(),
+			                 par3dgeom->getDetectorSpacingY(),
+			                 par3dgeom->getProjectionAngles(),
+			                 m_iGPUIndex, m_iVoxelSuperSampling);
+		} else {
+			astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(),
+			                 pSinoMem->getDataConst(),
+			                 volgeom.getGridColCount(),
+			                 volgeom.getGridRowCount(),
+			                 volgeom.getGridSliceCount(),
+			                 par3dgeom->getProjectionCount(),
+			                 par3dgeom->getDetectorColCount(),
+			                 par3dgeom->getDetectorRowCount(),
+			                 par3dgeom->getDetectorSpacingX(),
+			                 par3dgeom->getDetectorSpacingY(),
+			                 par3dgeom->getProjectionAngles(),
+			                 m_iGPUIndex, m_iVoxelSuperSampling);
+		}
 	} else if (parvec3dgeom) {
-		astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(),
-		                 volgeom.getGridColCount(),
-		                 volgeom.getGridRowCount(),
-		                 volgeom.getGridSliceCount(),
-		                 parvec3dgeom->getProjectionCount(),
-		                 parvec3dgeom->getDetectorColCount(),
-		                 parvec3dgeom->getDetectorRowCount(),
-		                 parvec3dgeom->getProjectionVectors(),
-		                 m_iGPUIndex, m_iVoxelSuperSampling);
+		if (!m_bSIRTWeighting) {
+			astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(),
+			                 volgeom.getGridColCount(),
+			                 volgeom.getGridRowCount(),
+			                 volgeom.getGridSliceCount(),
+			                 parvec3dgeom->getProjectionCount(),
+			                 parvec3dgeom->getDetectorColCount(),
+			                 parvec3dgeom->getDetectorRowCount(),
+			                 parvec3dgeom->getProjectionVectors(),
+			                 m_iGPUIndex, m_iVoxelSuperSampling);
+		} else {
+			astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(),
+			                 pSinoMem->getDataConst(),
+			                 volgeom.getGridColCount(),
+			                 volgeom.getGridRowCount(),
+			                 volgeom.getGridSliceCount(),
+			                 parvec3dgeom->getProjectionCount(),
+			                 parvec3dgeom->getDetectorColCount(),
+			                 parvec3dgeom->getDetectorRowCount(),
+			                 parvec3dgeom->getProjectionVectors(),
+			                 m_iGPUIndex, m_iVoxelSuperSampling);
+		}
 	} else if (conevecgeom) {
 		astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(),
 		                volgeom.getGridColCount(),
-- 
cgit v1.2.3