Merge pull request #41 from wjp/aniso

Add support for non-cube voxels

1 files changed, 64 insertions, 77 deletions

diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 5670873..dd6d551 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -67,34 +67,41 @@ template<typename ProjectionT>
 static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom,
                           unsigned int iProjectionAngleCount,
                           ProjectionT*& pProjs,
-                          float& fOutputScale)
+                          SProjectorParams3D& params)
 {
 	assert(pVolGeom);
 	assert(pProjs);
 
+#if 0
 	// TODO: Relative instead of absolute
 	const float EPS = 0.00001f;
 	if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthY()) > EPS)
 		return false;
 	if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthZ()) > EPS)
 		return false;
-
+#endif
 
 	// Translate
 	float dx = -(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2;
 	float dy = -(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2;
 	float dz = -(pVolGeom->getWindowMinZ() + pVolGeom->getWindowMaxZ()) / 2;
 
-	float factor = 1.0f / pVolGeom->getPixelLengthX();
+	float fx = 1.0f / pVolGeom->getPixelLengthX();
+	float fy = 1.0f / pVolGeom->getPixelLengthY();
+	float fz = 1.0f / pVolGeom->getPixelLengthZ();
 
 	for (int i = 0; i < iProjectionAngleCount; ++i) {
 		// CHECKME: Order of scaling and translation
 		pProjs[i].translate(dx, dy, dz);
-		pProjs[i].scale(factor);
+		pProjs[i].scale(fx, fy, fz);
 	}
 
+	params.fVolScaleX = pVolGeom->getPixelLengthX();
+	params.fVolScaleY = pVolGeom->getPixelLengthY();
+	params.fVolScaleZ = pVolGeom->getPixelLengthZ();
+
 	// CHECKME: Check factor
-	fOutputScale *= pVolGeom->getPixelLengthX();
+	//params.fOutputScale *= pVolGeom->getPixelLengthX();
 
 	return true;
 }
@@ -108,10 +115,8 @@ bool convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom,
 	dims.iVolY = pVolGeom->getGridRowCount();
 	dims.iVolZ = pVolGeom->getGridSliceCount();
 	dims.iProjAngles = pProjGeom->getProjectionCount();
-	dims.iProjU = pProjGeom->getDetectorColCount(),
-	dims.iProjV = pProjGeom->getDetectorRowCount(),
-	dims.iRaysPerDetDim = 1;
-	dims.iRaysPerVoxelDim = 1;
+	dims.iProjU = pProjGeom->getDetectorColCount();
+	dims.iProjV = pProjGeom->getDetectorRowCount();
 
 	if (dims.iVolX <= 0 || dims.iVolX <= 0 || dims.iVolX <= 0)
 		return false;
@@ -124,7 +129,7 @@ bool convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom,
 
 bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
                           const CParallelProjectionGeometry3D* pProjGeom,
-                          SPar3DProjection*& pProjs, float& fOutputScale)
+                          SPar3DProjection*& pProjs, SProjectorParams3D& params)
 {
 	assert(pVolGeom);
 	assert(pProjGeom);
@@ -141,16 +146,14 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 
 	bool ok;
 
-	fOutputScale = 1.0f;
-
-	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params);
 
 	return ok;
 }
 
 bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
                           const CParallelVecProjectionGeometry3D* pProjGeom,
-                          SPar3DProjection*& pProjs, float& fOutputScale)
+                          SPar3DProjection*& pProjs, SProjectorParams3D& params)
 {
 	assert(pVolGeom);
 	assert(pProjGeom);
@@ -164,16 +167,14 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 
 	bool ok;
 
-	fOutputScale = 1.0f;
-
-	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params);
 
 	return ok;
 }
 
 bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
                           const CConeProjectionGeometry3D* pProjGeom,
-                          SConeProjection*& pProjs, float& fOutputScale)
+                          SConeProjection*& pProjs, SProjectorParams3D& params)
 {
 	assert(pVolGeom);
 	assert(pProjGeom);
@@ -192,16 +193,14 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 
 	bool ok;
 
-	fOutputScale = 1.0f;
-
-	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params);
 
 	return ok;
 }
 
 bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
                           const CConeVecProjectionGeometry3D* pProjGeom,
-                          SConeProjection*& pProjs, float& fOutputScale)
+                          SConeProjection*& pProjs, SProjectorParams3D& params)
 {
 	assert(pVolGeom);
 	assert(pProjGeom);
@@ -215,9 +214,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 
 	bool ok;
 
-	fOutputScale = 1.0f;
-
-	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+	ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, params);
 
 	return ok;
 }
@@ -227,7 +224,7 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
                           const CProjectionGeometry3D* pProjGeom,
                           SPar3DProjection*& pParProjs,
                           SConeProjection*& pConeProjs,
-                          float& fOutputScale)
+                          SProjectorParams3D& params)
 {
 	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(pProjGeom);
 	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(pProjGeom);
@@ -240,13 +237,13 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 	bool ok;
 
 	if (conegeom) {
-		ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, fOutputScale);
+		ok = convertAstraGeometry(pVolGeom, conegeom, pConeProjs, params);
 	} else if (conevec3dgeom) {
-		ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, fOutputScale);
+		ok = convertAstraGeometry(pVolGeom, conevec3dgeom, pConeProjs, params);
 	} else if (par3dgeom) {
-		ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, fOutputScale);
+		ok = convertAstraGeometry(pVolGeom, par3dgeom, pParProjs, params);
 	} else if (parvec3dgeom) {
-		ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, fOutputScale);
+		ok = convertAstraGeometry(pVolGeom, parvec3dgeom, pParProjs, params);
 	} else {
 		ok = false;
 	}
@@ -260,20 +257,17 @@ bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
 class AstraSIRT3d_internal {
 public:
 	SDimensions3D dims;
+	SProjectorParams3D params;
 	CUDAProjectionType3d projType;
 
 	float* angles;
 	float fOriginSourceDistance;
 	float fOriginDetectorDistance;
-	float fSourceZ;
-	float fDetSize;
 	float fRelaxation;
 
 	SConeProjection* projs;
 	SPar3DProjection* parprojs;
 
-	float fOutputScale;
-
 	bool initialized;
 	bool setStartReconstruction;
 
@@ -305,12 +299,9 @@ AstraSIRT3d::AstraSIRT3d()
 	pData->dims.iProjAngles = 0;
 	pData->dims.iProjU = 0;
 	pData->dims.iProjV = 0;
-	pData->dims.iRaysPerDetDim = 1;
-	pData->dims.iRaysPerVoxelDim = 1;
 
 	pData->projs = 0;
 	pData->parprojs = 0;
-	pData->fOutputScale = 1.0f;
 
 	pData->fRelaxation = 1.0f;
 
@@ -361,7 +352,7 @@ bool AstraSIRT3d::setGeometry(const CVolumeGeometry3D* pVolGeom,
 
 	ok = convertAstraGeometry(pVolGeom, pProjGeom,
 	                          pData->parprojs, pData->projs,
-	                          pData->fOutputScale);
+	                          pData->params);
 	if (!ok)
 		return false;
 
@@ -386,8 +377,8 @@ bool AstraSIRT3d::enableSuperSampling(unsigned int iVoxelSuperSampling,
 	if (iVoxelSuperSampling == 0 || iDetectorSuperSampling == 0)
 		return false;
 
-	pData->dims.iRaysPerVoxelDim = iVoxelSuperSampling;
-	pData->dims.iRaysPerDetDim = iDetectorSuperSampling;
+	pData->params.iRaysPerVoxelDim = iVoxelSuperSampling;
+	pData->params.iRaysPerDetDim = iDetectorSuperSampling;
 
 	return true;
 }
@@ -447,9 +438,9 @@ bool AstraSIRT3d::init()
 	bool ok;
 
 	if (pData->projType == PROJ_PARALLEL) {
-		ok = pData->sirt.setPar3DGeometry(pData->dims, pData->parprojs, pData->fOutputScale);
+		ok = pData->sirt.setPar3DGeometry(pData->dims, pData->parprojs, pData->params);
 	} else {
-		ok = pData->sirt.setConeGeometry(pData->dims, pData->projs, pData->fOutputScale);
+		ok = pData->sirt.setConeGeometry(pData->dims, pData->projs, pData->params);
 	}
 
 	if (!ok)
@@ -652,19 +643,16 @@ float AstraSIRT3d::computeDiffNorm()
 class AstraCGLS3d_internal {
 public:
 	SDimensions3D dims;
+	SProjectorParams3D params;
 	CUDAProjectionType3d projType;
 
 	float* angles;
 	float fOriginSourceDistance;
 	float fOriginDetectorDistance;
-	float fSourceZ;
-	float fDetSize;
 
 	SConeProjection* projs;
 	SPar3DProjection* parprojs;
 
-	float fOutputScale;
-
 	bool initialized;
 	bool setStartReconstruction;
 
@@ -696,12 +684,9 @@ AstraCGLS3d::AstraCGLS3d()
 	pData->dims.iProjAngles = 0;
 	pData->dims.iProjU = 0;
 	pData->dims.iProjV = 0;
-	pData->dims.iRaysPerDetDim = 1;
-	pData->dims.iRaysPerVoxelDim = 1;
 
 	pData->projs = 0;
 	pData->parprojs = 0;
-	pData->fOutputScale = 1.0f;
 
 	pData->initialized = false;
 	pData->setStartReconstruction = false;
@@ -750,7 +735,7 @@ bool AstraCGLS3d::setGeometry(const CVolumeGeometry3D* pVolGeom,
 
 	ok = convertAstraGeometry(pVolGeom, pProjGeom,
 	                          pData->parprojs, pData->projs,
-	                          pData->fOutputScale);
+	                          pData->params);
 	if (!ok)
 		return false;
 
@@ -774,8 +759,8 @@ bool AstraCGLS3d::enableSuperSampling(unsigned int iVoxelSuperSampling,
 	if (iVoxelSuperSampling == 0 || iDetectorSuperSampling == 0)
 		return false;
 
-	pData->dims.iRaysPerVoxelDim = iVoxelSuperSampling;
-	pData->dims.iRaysPerDetDim = iDetectorSuperSampling;
+	pData->params.iRaysPerVoxelDim = iVoxelSuperSampling;
+	pData->params.iRaysPerDetDim = iDetectorSuperSampling;
 
 	return true;
 }
@@ -829,9 +814,9 @@ bool AstraCGLS3d::init()
 	bool ok;
 
 	if (pData->projType == PROJ_PARALLEL) {
-		ok = pData->cgls.setPar3DGeometry(pData->dims, pData->parprojs, pData->fOutputScale);
+		ok = pData->cgls.setPar3DGeometry(pData->dims, pData->parprojs, pData->params);
 	} else {
-		ok = pData->cgls.setConeGeometry(pData->dims, pData->projs, pData->fOutputScale);
+		ok = pData->cgls.setConeGeometry(pData->dims, pData->projs, pData->params);
 	}
 
 	if (!ok)
@@ -1039,23 +1024,23 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections,
                  Cuda3DProjectionKernel projKernel)
 {
 	SDimensions3D dims;
+	SProjectorParams3D params;
+
+	params.iRaysPerDetDim = iDetectorSuperSampling;
 
 	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
 	if (!ok)
 		return false;
 
-	dims.iRaysPerDetDim = iDetectorSuperSampling;
 	if (iDetectorSuperSampling == 0)
 		return false;
 
 	SPar3DProjection* pParProjs;
 	SConeProjection* pConeProjs;
 
-	float outputScale;
-
 	ok = convertAstraGeometry(pVolGeom, pProjGeom,
 	                          pParProjs, pConeProjs,
-	                          outputScale);
+	                          params);
 
 
 	if (iGPUIndex != -1) {
@@ -1093,10 +1078,10 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections,
 	if (pParProjs) {
 		switch (projKernel) {
 		case ker3d_default:
-			ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, outputScale);
+			ok &= Par3DFP(D_volumeData, D_projData, dims, pParProjs, params);
 			break;
 		case ker3d_sum_square_weights:
-			ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, outputScale*outputScale);
+			ok &= Par3DFP_SumSqW(D_volumeData, D_projData, dims, pParProjs, params);
 			break;
 		default:
 			assert(false);
@@ -1104,7 +1089,7 @@ bool astraCudaFP(const float* pfVolume, float* pfProjections,
 	} else {
 		switch (projKernel) {
 		case ker3d_default:
-			ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, outputScale);
+			ok &= ConeFP(D_volumeData, D_projData, dims, pConeProjs, params);
 			break;
 		default:
 			assert(false);
@@ -1129,21 +1114,20 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections,
                  int iGPUIndex, int iVoxelSuperSampling)
 {
 	SDimensions3D dims;
+	SProjectorParams3D params;
+
+	params.iRaysPerVoxelDim = iVoxelSuperSampling;
 
 	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
 	if (!ok)
 		return false;
 
-	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
-
 	SPar3DProjection* pParProjs;
 	SConeProjection* pConeProjs;
 
-	float outputScale;
-
 	ok = convertAstraGeometry(pVolGeom, pProjGeom,
 	                          pParProjs, pConeProjs,
-	                          outputScale);
+	                          params);
 
 	if (iGPUIndex != -1) {
 		cudaSetDevice(iGPUIndex);
@@ -1179,9 +1163,9 @@ bool astraCudaBP(float* pfVolume, const float* pfProjections,
 	}
 
 	if (pParProjs)
-		ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, outputScale);
+		ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, params);
 	else
-		ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, outputScale);
+		ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, params);
 
 	ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX);
 
@@ -1204,21 +1188,21 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume,
                       int iGPUIndex, int iVoxelSuperSampling)
 {
 	SDimensions3D dims;
+	SProjectorParams3D params;
+
+	params.iRaysPerVoxelDim = iVoxelSuperSampling;
 
 	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
 	if (!ok)
 		return false;
 
-	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
 
 	SPar3DProjection* pParProjs;
 	SConeProjection* pConeProjs;
 
-	float outputScale;
-
 	ok = convertAstraGeometry(pVolGeom, pProjGeom,
 	                          pParProjs, pConeProjs,
-	                          outputScale);
+	                          params);
 
 	if (iGPUIndex != -1) {
 		cudaSetDevice(iGPUIndex);
@@ -1255,9 +1239,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume,
 	processSino3D<opSet>(D_projData, 1.0f, dims);
 
 	if (pParProjs)
-		ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, outputScale);
+		ok &= Par3DBP(D_pixelWeight, D_projData, dims, pParProjs, params);
 	else
-		ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, outputScale);
+		ok &= ConeBP(D_pixelWeight, D_projData, dims, pConeProjs, params);
 
 	processVol3D<opInvert>(D_pixelWeight, dims);
 	if (!ok) {
@@ -1272,9 +1256,9 @@ bool astraCudaBP_SIRTWeighted(float* pfVolume,
 	ok &= zeroVolumeData(D_volumeData, dims);
 	// Do BP into D_volumeData
 	if (pParProjs)
-		ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, outputScale);
+		ok &= Par3DBP(D_volumeData, D_projData, dims, pParProjs, params);
 	else
-		ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, outputScale);
+		ok &= ConeBP(D_volumeData, D_projData, dims, pConeProjs, params);
 
 	// Multiply with weights
 	processVol3D<opMul>(D_volumeData, D_pixelWeight, dims);
@@ -1314,6 +1298,9 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections,
                   int iGPUIndex, int iVoxelSuperSampling)
 {
 	SDimensions3D dims;
+	SProjectorParams3D params;
+
+	params.iRaysPerVoxelDim = iVoxelSuperSampling;
 
 	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
 
@@ -1323,7 +1310,6 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections,
 	if (!ok)
 		return false;
 
-	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
 
 	if (iVoxelSuperSampling == 0)
 		return false;
@@ -1367,6 +1353,7 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections,
 
 
 	// TODO: Offer interface for SrcZ, DetZ
+	// TODO: Voxel scaling
 	ok &= FDK(D_volumeData, D_projData, fOriginSourceDistance,
 	          fOriginDetectorDistance, 0, 0, fDetUSize, fDetVSize,
 	          dims, pfAngles, bShortScan);