Add cone_fp kernel support for anisotropic voxels

2 files changed, 31 insertions, 10 deletions

diff --git a/cuda/3d/cone_fp.cu b/cuda/3d/cone_fp.cu
index 5a31b65..2feec06 100644
--- a/cuda/3d/cone_fp.cu
+++ b/cuda/3d/cone_fp.cu
@@ -139,7 +139,8 @@ template<class COORD>
 __global__ void cone_FP_t(float* D_projData, unsigned int projPitch,
                           unsigned int startSlice,
                           unsigned int startAngle, unsigned int endAngle,
-                          const SDimensions3D dims, float fOutputScale)
+                          const SDimensions3D dims,
+                          float fScale1, float fScale2, float fOutputScale)
 {
 	COORD c;
 
@@ -188,7 +189,7 @@ __global__ void cone_FP_t(float* D_projData, unsigned int projPitch,
 		const float b1 = c.c1(fSrcX,fSrcY,fSrcZ) - a1 * c.c0(fSrcX,fSrcY,fSrcZ);
 		const float b2 = c.c2(fSrcX,fSrcY,fSrcZ) - a2 * c.c0(fSrcX,fSrcY,fSrcZ);
 
-		const float fDistCorr = sqrt(a1*a1+a2*a2+1.0f) * fOutputScale;
+		const float fDistCorr = sqrt(a1*a1*fScale1+a2*a2*fScale2+1.0f) * fOutputScale;
 
 		float fVal = 0.0f;
 
@@ -214,7 +215,8 @@ template<class COORD>
 __global__ void cone_FP_SS_t(float* D_projData, unsigned int projPitch,
                              unsigned int startSlice,
                              unsigned int startAngle, unsigned int endAngle,
-                             const SDimensions3D dims, int iRaysPerDetDim, float fOutputScale)
+                             const SDimensions3D dims, int iRaysPerDetDim,
+                             float fScale1, float fScale2, float fOutputScale)
 {
 	COORD c;
 
@@ -272,7 +274,7 @@ __global__ void cone_FP_SS_t(float* D_projData, unsigned int projPitch,
 		const float b1 = c.c1(fSrcX,fSrcY,fSrcZ) - a1 * c.c0(fSrcX,fSrcY,fSrcZ);
 		const float b2 = c.c2(fSrcX,fSrcY,fSrcZ) - a2 * c.c0(fSrcX,fSrcY,fSrcZ);
 
-		const float fDistCorr = sqrt(a1*a1+a2*a2+1.0f) * fOutputScale;
+		const float fDistCorr = sqrt(a1*a1*fScale1+a2*a2*fScale2+1.0f) * fOutputScale;
 
 		float fVal = 0.0f;
 
@@ -372,23 +374,38 @@ bool ConeFP_Array_internal(cudaPitchedPtr D_projData,
 				// printf("angle block: %d to %d, %d (%dx%d, %dx%d)\n", blockStart, blockEnd, blockDirection, dimGrid.x, dimGrid.y, dimBlock.x, dimBlock.y);
 
 				if (blockDirection == 0) {
+					float fS1 = params.fVolScaleY / params.fVolScaleX;
+					fS1 *= fS1;
+					float fS2 = params.fVolScaleZ / params.fVolScaleX;
+					fS2 *= fS2;
+					float fS0 = params.fOutputScale * params.fVolScaleX;
 					for (unsigned int i = 0; i < dims.iVolX; i += g_blockSlices)
 						if (params.iRaysPerDetDim == 1)
-							cone_FP_t<DIR_X><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.fOutputScale);
+							cone_FP_t<DIR_X><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, fS1, fS2, fS0);
 						else
-							cone_FP_SS_t<DIR_X><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.iRaysPerDetDim, params.fOutputScale);
+							cone_FP_SS_t<DIR_X><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.iRaysPerDetDim, fS1, fS2, fS0);
 				} else if (blockDirection == 1) {
+					float fS1 = params.fVolScaleX / params.fVolScaleY;
+					fS1 *= fS1;
+					float fS2 = params.fVolScaleZ / params.fVolScaleY;
+					fS2 *= fS2;
+					float fS0 = params.fOutputScale * params.fVolScaleY;
 					for (unsigned int i = 0; i < dims.iVolY; i += g_blockSlices)
 						if (params.iRaysPerDetDim == 1)
-							cone_FP_t<DIR_Y><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.fOutputScale);
+							cone_FP_t<DIR_Y><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, fS1, fS2, fS0);
 						else
-							cone_FP_SS_t<DIR_Y><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.iRaysPerDetDim, params.fOutputScale);
+							cone_FP_SS_t<DIR_Y><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.iRaysPerDetDim, fS1, fS2, fS0);
 				} else if (blockDirection == 2) {
+					float fS1 = params.fVolScaleX / params.fVolScaleZ;
+					fS1 *= fS1;
+					float fS2 = params.fVolScaleY / params.fVolScaleZ;
+					fS2 *= fS2;
+					float fS0 = params.fOutputScale * params.fVolScaleZ;
 					for (unsigned int i = 0; i < dims.iVolZ; i += g_blockSlices)
 						if (params.iRaysPerDetDim == 1)
-							cone_FP_t<DIR_Z><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.fOutputScale);
+							cone_FP_t<DIR_Z><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, fS1, fS2, fS0);
 						else
-							cone_FP_SS_t<DIR_Z><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.iRaysPerDetDim, params.fOutputScale);
+							cone_FP_SS_t<DIR_Z><<<dimGrid, dimBlock, 0, stream>>>((float*)D_projData.ptr, D_projData.pitch/sizeof(float), i, blockStart, blockEnd, dims, params.iRaysPerDetDim, fS1, fS2, fS0);
 				}
 
 			}
diff --git a/cuda/3d/dims3d.h b/cuda/3d/dims3d.h
index 569b395..a15c67a 100644
--- a/cuda/3d/dims3d.h
+++ b/cuda/3d/dims3d.h
@@ -57,12 +57,16 @@ struct SProjectorParams3D {
 	SProjectorParams3D() :
 	    iRaysPerDetDim(1), iRaysPerVoxelDim(1),
 	    fOutputScale(1.0f),
+	    fVolScaleX(1.0f), fVolScaleY(1.0f), fVolScaleZ(1.0f),
 	    ker(ker3d_default)
 	{ }
 
 	unsigned int iRaysPerDetDim;
 	unsigned int iRaysPerVoxelDim;
 	float fOutputScale;
+	float fVolScaleX;
+	float fVolScaleY;
+	float fVolScaleZ;
 	Cuda3DProjectionKernel ker;
 };