Add flexible volgeom3d support to astraCudaBP

3 files changed, 24 insertions, 253 deletions

diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 2f7ea99..97bebf4 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -1331,173 +1331,30 @@ bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections,
 
 }
 
-bool astraCudaConeBP(float* pfVolume, const float* pfProjections,
-                     unsigned int iVolX,
-                     unsigned int iVolY,
-                     unsigned int iVolZ,
-                     unsigned int iProjAngles,
-                     unsigned int iProjU,
-                     unsigned int iProjV,
-                     float fOriginSourceDistance,
-                     float fOriginDetectorDistance,
-                     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;
-
-	SConeProjection* p = genConeProjections(iProjAngles,
-                                            iProjU, iProjV,
-                                            fOriginSourceDistance,
-                                            fOriginDetectorDistance,
-                                            fDetUSize, fDetVSize,
-                                            pfAngles);
-
-	bool ok;
-	ok = astraCudaConeBP(pfVolume, pfProjections, iVolX, iVolY, iVolZ,
-	                     iProjAngles, iProjU, iProjV, p, iGPUIndex, iVoxelSuperSampling);
-
-	delete[] p;
-
-	return ok;
-}
 
-bool astraCudaConeBP(float* pfVolume, const float* pfProjections,
-                     unsigned int iVolX,
-                     unsigned int iVolY,
-                     unsigned int iVolZ,
-                     unsigned int iProjAngles,
-                     unsigned int iProjU,
-                     unsigned int iProjV,
-                     const SConeProjection *pfAngles,
-                     int iGPUIndex, int iVoxelSuperSampling)
+bool astraCudaBP(float* pfVolume, const float* pfProjections,
+                 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)
-		return false;
-
-	dims.iProjAngles = iProjAngles;
-	dims.iProjU = iProjU;
-	dims.iProjV = iProjV;
-
-	if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0)
-		return false;
-
-	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
-
-	if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0)
-		return false;
-
-	if (iGPUIndex != -1) {
-		cudaSetDevice(iGPUIndex);
-		cudaError_t err = cudaGetLastError();
-
-		// Ignore errors caused by calling cudaSetDevice multiple times
-		if (err != cudaSuccess && err != cudaErrorSetOnActiveProcess)
-			return false;
-	}
-
-	cudaPitchedPtr D_volumeData = allocateVolumeData(dims);
-	bool ok = D_volumeData.ptr;
+	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
 	if (!ok)
 		return false;
 
-	cudaPitchedPtr D_projData = allocateProjectionData(dims);
-	ok = D_projData.ptr;
-	if (!ok) {
-		cudaFree(D_volumeData.ptr);
-		return false;
-	}
-
-	ok &= copyProjectionsToDevice(pfProjections, D_projData,
-	                              dims, dims.iProjU);
-
-	ok &= zeroVolumeData(D_volumeData, dims);
-
-	if (!ok) {
-		cudaFree(D_volumeData.ptr);
-		cudaFree(D_projData.ptr);
-		return false;
-	}
-
-	ok &= ConeBP(D_volumeData, D_projData, dims, pfAngles);
-
-	ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX);
-
-
-	cudaFree(D_volumeData.ptr);
-	cudaFree(D_projData.ptr);
-
-	return ok;
-
-}
-
-bool astraCudaPar3DBP(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(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,
-                      unsigned int iVolY,
-                      unsigned int iVolZ,
-                      unsigned int iProjAngles,
-                      unsigned int iProjU,
-                      unsigned int iProjV,
-                      const SPar3DProjection *pfAngles,
-                      int iGPUIndex, int iVoxelSuperSampling)
-{
-	SDimensions3D dims;
+	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
 
-	dims.iVolX = iVolX;
-	dims.iVolY = iVolY;
-	dims.iVolZ = iVolZ;
-	if (iVolX == 0 || iVolY == 0 || iVolZ == 0)
-		return false;
+	SPar3DProjection* pParProjs;
+	SConeProjection* pConeProjs;
 
-	dims.iProjAngles = iProjAngles;
-	dims.iProjU = iProjU;
-	dims.iProjV = iProjV;
+	float outputScale;
 
-	if (iProjAngles == 0 || iProjU == 0 || iProjV == 0 || pfAngles == 0)
-		return false;
+	ok = convertAstraGeometry(pVolGeom, pProjGeom,
+	                          pParProjs, pConeProjs,
+	                          outputScale);
 
-	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
+	// TODO: OutputScale
 
 	if (iGPUIndex != -1) {
 		cudaSetDevice(iGPUIndex);
@@ -1510,7 +1367,7 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections,
 
 
 	cudaPitchedPtr D_volumeData = allocateVolumeData(dims);
-	bool ok = D_volumeData.ptr;
+	ok = D_volumeData.ptr;
 	if (!ok)
 		return false;
 
@@ -1532,7 +1389,10 @@ bool astraCudaPar3DBP(float* pfVolume, const float* pfProjections,
 		return false;
 	}
 
-	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);
 
 	ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX);
 
diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h
index b2e4e08..5464d2f 100644
--- a/cuda/3d/astra3d.h
+++ b/cuda/3d/astra3d.h
@@ -332,50 +332,9 @@ _AstraExport bool astraCudaPar3DFP(const float* pfVolume, float* pfProjections,
                       Cuda3DProjectionKernel projKernel);
 
 
-_AstraExport bool astraCudaConeBP(float* pfVolume, const float* pfProjections,
-                     unsigned int iVolX,
-                     unsigned int iVolY,
-                     unsigned int iVolZ,
-                     unsigned int iProjAngles,
-                     unsigned int iProjU,
-                     unsigned int iProjV,
-                     float fOriginSourceDistance,
-                     float fOriginDetectorDistance,
-                     float fDetUSize,
-                     float fDetVSize,
-                     const float *pfAngles,
-                     int iGPUIndex, int iVoxelSuperSampling);
-
-_AstraExport bool astraCudaConeBP(float* pfVolume, const float* pfProjections,
-                     unsigned int iVolX,
-                     unsigned int iVolY,
-                     unsigned int iVolZ,
-                     unsigned int iProjAngles,
-                     unsigned int iProjU,
-                     unsigned int iProjV,
-                     const SConeProjection *pfAngles,
-                     int iGPUIndex, int iVoxelSuperSampling);
-
-_AstraExport bool astraCudaPar3DBP(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(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(float* pfVolume, const float* pfProjections,
+                      const CVolumeGeometry3D* pVolGeom,
+                      const CProjectionGeometry3D* pProjGeom,
                       int iGPUIndex, int iVoxelSuperSampling);
 
 _AstraExport bool astraCudaBP_SIRTWeighted(float* pfVolume, const float* pfProjections,
diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp
index 7117cfc..a8a1b0a 100644
--- a/src/CudaBackProjectionAlgorithm3D.cpp
+++ b/src/CudaBackProjectionAlgorithm3D.cpp
@@ -164,10 +164,6 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 	ASTRA_ASSERT(pReconMem);
 
 	const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry();
-	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom);
-	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom);
-	const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom);
-	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
 	const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry();
 
 	if (m_bSIRTWeighting) {
@@ -175,54 +171,10 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 		                         pSinoMem->getDataConst(),
 		                         &volgeom, projgeom,
 		                         m_iGPUIndex, m_iVoxelSuperSampling);
-	} else if (conegeom) {
-		astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(),
-		                volgeom.getGridColCount(),
-		                volgeom.getGridRowCount(),
-		                volgeom.getGridSliceCount(),
-		                conegeom->getProjectionCount(),
-		                conegeom->getDetectorColCount(),
-		                conegeom->getDetectorRowCount(),
-		                conegeom->getOriginSourceDistance(),
-		                conegeom->getOriginDetectorDistance(),
-		                conegeom->getDetectorSpacingX(),
-		                conegeom->getDetectorSpacingY(),
-		                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);
-	} 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);
-	} else if (conevecgeom) {
-		astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(),
-		                volgeom.getGridColCount(),
-		                volgeom.getGridRowCount(),
-		                volgeom.getGridSliceCount(),
-		                conevecgeom->getProjectionCount(),
-		                conevecgeom->getDetectorColCount(),
-		                conevecgeom->getDetectorRowCount(),
-		                conevecgeom->getProjectionVectors(),
-		                m_iGPUIndex, m_iVoxelSuperSampling);
 	} else {
-		ASTRA_ASSERT(false);
+		astraCudaBP(pReconMem->getData(), pSinoMem->getDataConst(),
+		            &volgeom, projgeom,
+		            m_iGPUIndex, m_iVoxelSuperSampling);
 	}
 
 }