Add flexible volgeom3d support to astraCudaBP_SIRTWeighted

3 files changed, 66 insertions, 139 deletions

diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index eff928d..2f7ea99 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -306,7 +306,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 
 
 bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
-                          const CConeVecProjectionGeometry3D* pProjGeom,
+                          const CProjectionGeometry3D* pProjGeom,
                           SPar3DProjection*& pParProjs,
                           SConeProjection*& pConeProjs,
                           float& fOutputScale)
@@ -322,13 +322,13 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 	bool ok;
 
 	if (conegeom) {
-		ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, outputScale);
+		ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, fOutputScale);
 	} else if (conevec3dgeom) {
-		ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, outputScale);
+		ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, fOutputScale);
 	} else if (par3dgeom) {
-		ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, outputScale);
+		ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, fOutputScale);
 	} else if (parvec3dgeom) {
-		ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, outputScale);
+		ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, fOutputScale);
 	} else {
 		ok = false;
 	}
@@ -1471,40 +1471,6 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections,
 	return ok;
 }
 
-// This computes the column weights, divides by them, and adds the
-// result to the current volume. This is both more expensive and more
-// GPU memory intensive than the regular BP, but allows saving system RAM.
-bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, const float* pfProjections,
-                      unsigned int iVolX,
-                      unsigned int iVolY,
-                      unsigned int iVolZ,
-                      unsigned int iProjAngles,
-                      unsigned int iProjU,
-                      unsigned int iProjV,
-                      float fDetUSize,
-                      float fDetVSize,
-                      const float *pfAngles,
-                      int iGPUIndex, int iVoxelSuperSampling)
-{
-	if (iVolX == 0 || iVolY == 0 || iVolZ == 0)
-		return false;
-	if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0)
-		return false;
-
-	SPar3DProjection* p = genPar3DProjections(iProjAngles,
-                                             iProjU, iProjV,
-                                             fDetUSize, fDetVSize,
-                                             pfAngles);
-
-	bool ok;
-	ok = astraCudaPar3DBP_SIRTWeighted(pfVolume, pfProjections, iVolX, iVolY, iVolZ,
-	                      iProjAngles, iProjU, iProjV, p, iGPUIndex, iVoxelSuperSampling);
-
-	delete[] p;
-
-	return ok;
-}
-
 
 bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections,
                       unsigned int iVolX,
@@ -1582,33 +1548,30 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections,
 // This computes the column weights, divides by them, and adds the
 // result to the current volume. This is both more expensive and more
 // GPU memory intensive than the regular BP, but allows saving system RAM.
-bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume,
+bool astraCudaBP_SIRTWeighted(float* pfVolume,
                       const float* pfProjections,
-                      unsigned int iVolX,
-                      unsigned int iVolY,
-                      unsigned int iVolZ,
-                      unsigned int iProjAngles,
-                      unsigned int iProjU,
-                      unsigned int iProjV,
-                      const SPar3DProjection *pfAngles,
+                      const CVolumeGeometry3D* pVolGeom,
+                      const CProjectionGeometry3D* pProjGeom,
                       int iGPUIndex, int iVoxelSuperSampling)
 {
 	SDimensions3D dims;
 
-	dims.iVolX = iVolX;
-	dims.iVolY = iVolY;
-	dims.iVolZ = iVolZ;
-	if (iVolX == 0 || iVolY == 0 || iVolZ == 0)
+	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
+	if (!ok)
 		return false;
 
-	dims.iProjAngles = iProjAngles;
-	dims.iProjU = iProjU;
-	dims.iProjV = iProjV;
+	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
 
-	if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0)
-		return false;
+	SPar3DProjection* pParProjs;
+	SConeProjection* pConeProjs;
 
-	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
+	float outputScale;
+
+	ok = convertAstraGeometry(pVolGeom, pProjGeom,
+	                          pParProjs, pConeProjs,
+	                          outputScale);
+
+	// TODO: OutputScale
 
 	if (iGPUIndex != -1) {
 		cudaSetDevice(iGPUIndex);
@@ -1621,7 +1584,7 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume,
 
 
 	cudaPitchedPtr D_pixelWeight = allocateVolumeData(dims);
-	bool ok = D_pixelWeight.ptr;
+	ok = D_pixelWeight.ptr;
 	if (!ok)
 		return false;
 
@@ -1643,7 +1606,12 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume,
 	// Compute weights
 	ok &= zeroVolumeData(D_pixelWeight, dims);
 	processSino3D<opSet>(D_projData, 1.0f, dims);
-	ok &= Par3DBP(D_pixelWeight, D_projData, dims, pfAngles);
+
+	if (pParProjs)
+		ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs);
+	else
+		ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs);
+
 	processVol3D<opInvert>(D_pixelWeight, dims);
 	if (!ok) {
 		cudaFree(D_pixelWeight.ptr);
@@ -1656,7 +1624,11 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume,
 	                              dims, dims.iProjU);
 	ok &= zeroVolumeData(D_volumeData, dims);
 	// Do BP into D_volumeData
-	ok &= Par3DBP(D_volumeData, D_projData, dims, pfAngles);
+	if (pParProjs)
+		ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs);
+	else
+		ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs);
+
 	// Multiply with weights
 	processVol3D<opMul>(D_volumeData, D_pixelWeight, dims);
 
@@ -1679,6 +1651,9 @@ bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume,
 	cudaFree(D_volumeData.ptr);
 	cudaFree(D_projData.ptr);
 
+	delete[] pParProjs;
+	delete[] pConeProjs;
+
 	return ok;
 
 }
diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h
index 6bac8b2..b2e4e08 100644
--- a/cuda/3d/astra3d.h
+++ b/cuda/3d/astra3d.h
@@ -378,26 +378,9 @@ _AstraExport bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections,
                       const SPar3DProjection *pfAngles,
                       int iGPUIndex, int iVoxelSuperSampling);
 
-_AstraExport bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, const float* pfProjections,
-                      unsigned int iVolX,
-                      unsigned int iVolY,
-                      unsigned int iVolZ,
-                      unsigned int iProjAngles,
-                      unsigned int iProjU,
-                      unsigned int iProjV,
-                      float fDetUSize,
-                      float fDetVSize,
-                      const float *pfAngles,
-                      int iGPUIndex, int iVoxelSuperSampling);
-
-_AstraExport bool astraCudaPar3DBP_SIRTWeighted(float* pfVolume, const float* pfProjections,
-                      unsigned int iVolX,
-                      unsigned int iVolY,
-                      unsigned int iVolZ,
-                      unsigned int iProjAngles,
-                      unsigned int iProjU,
-                      unsigned int iProjV,
-                      const SPar3DProjection *pfAngles,
+_AstraExport bool astraCudaBP_SIRTWeighted(float* pfVolume, const float* pfProjections,
+                      const CVolumeGeometry3D* pVolGeom,
+                      const CProjectionGeometry3D* pProjGeom,
                       int iGPUIndex, int iVoxelSuperSampling);
 
 _AstraExport bool astraCudaFDK(float* pfVolume, const float* pfProjections,
diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp
index abcf096..7117cfc 100644
--- a/src/CudaBackProjectionAlgorithm3D.cpp
+++ b/src/CudaBackProjectionAlgorithm3D.cpp
@@ -107,16 +107,8 @@ 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");
-	}
+	m_bSIRTWeighting = _cfg.self->getOptionBool("SIRTWeighting", false);
+	CC.markOptionParsed("SIRTWeighting");
 
 	// success
 	m_bIsInitialized = _check();
@@ -178,7 +170,12 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
 	const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry();
 
-	if (conegeom) {
+	if (m_bSIRTWeighting) {
+		astraCudaBP_SIRTWeighted(pReconMem->getData(),
+		                         pSinoMem->getDataConst(),
+		                         &volgeom, projgeom,
+		                         m_iGPUIndex, m_iVoxelSuperSampling);
+	} else if (conegeom) {
 		astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(),
 		                volgeom.getGridColCount(),
 		                volgeom.getGridRowCount(),
@@ -193,55 +190,27 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 		                conegeom->getProjectionAngles(),
 		                m_iGPUIndex, m_iVoxelSuperSampling);
 	} else if (par3dgeom) {
-		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);
-		}
+		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 if (parvec3dgeom) {
-		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);
-		}
+		astraCudaPar3DBP(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(),