From 52ebf0670ae39abea04344c32a1dc054c6816a03 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Wed, 24 May 2017 11:52:05 +0200
Subject: Document FDK weighting slightly better

---
 cuda/3d/cone_bp.cu | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

(limited to 'cuda')

diff --git a/cuda/3d/cone_bp.cu b/cuda/3d/cone_bp.cu
index c3405b8..02242b0 100644
--- a/cuda/3d/cone_bp.cu
+++ b/cuda/3d/cone_bp.cu
@@ -126,6 +126,9 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng
 			float fVNum = fCv.w + fX * fCv.x + fY * fCv.y + fZ * fCv.z;
 			float fDen  = fCd.w + fX * fCd.x + fY * fCd.y + fZ * fCd.z;
 
+			// fCd.w = -|| u v s || (determinant of 3x3 matrix with cols u,v,s)
+			// fDen =  || u v (x-s) ||
+
 			float fU,fV, fr;
 
 			for (int idx = 0; idx < ZSIZE; idx++)
@@ -134,9 +137,17 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng
 				fU = fUNum * fr;
 				fV = fVNum * fr;
 				float fVal = tex3D(gT_coneProjTexture, fU, fAngle, fV);
-				if (FDKWEIGHT)
+				if (FDKWEIGHT) {
+					// The correct factor here is this one:
+					// Z[idx] += (fr*fCd.w)*(fr*fCd.w)*fVal;
+					// This is the square of the inverse magnification factor
+					// from fX,fY,fZ to the detector.
+
+					// Since we are assuming we have a circular cone
+					// beam trajectory, fCd.w is constant, and we instead
+					// multiply by fCd.w*fCd.w in the FDK preweighting step.
 					Z[idx] += fr*fr*fVal;
-				else
+				} else
 					Z[idx] += fVal;
 
 				fUNum += fCu.z;
-- 
cgit v1.2.3


From 9fa2ffdc5348f8f19de48d06a72b82bdc1ba8f22 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 30 Nov 2017 16:59:18 +0100
Subject: Start on fixing FDK and BP voxel-size weighting factors

---
 cuda/2d/astra.cu   |  2 +-
 cuda/3d/cone_bp.cu |  6 +++++-
 cuda/3d/fdk.cu     | 23 ++++++++++++++++-------
 cuda/3d/fdk.h      |  3 ++-
 4 files changed, 24 insertions(+), 10 deletions(-)

(limited to 'cuda')

diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index c0132b2..0069bc3 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -343,7 +343,7 @@ bool AstraFBP::run()
 
 		astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance,
 		              pData->fOriginDetectorDistance, 0.0f,
-		              pData->dims.fDetScale, 1.0f, // TODO: Are these correct?
+		              pData->dims.fDetScale, 1.0f, 1.0f, // TODO: Are these correct?
 		              pData->bShortScan, dims3d, pData->angles);
 	}
 
diff --git a/cuda/3d/cone_bp.cu b/cuda/3d/cone_bp.cu
index 02242b0..aff0834 100644
--- a/cuda/3d/cone_bp.cu
+++ b/cuda/3d/cone_bp.cu
@@ -266,7 +266,11 @@ bool ConeBP_Array(cudaPitchedPtr D_volumeData,
 {
 	bindProjDataTexture(D_projArray);
 
-	float fOutputScale = params.fOutputScale * params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ;
+	float fOutputScale;
+	if (params.bFDKWeighting)
+		fOutputScale = params.fOutputScale / (params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ);
+	else
+		fOutputScale = params.fOutputScale * (params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ);
 
 	for (unsigned int th = 0; th < dims.iProjAngles; th += g_MaxAngles) {
 		unsigned int angleCount = g_MaxAngles;
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 27357ad..0630afe 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -59,7 +59,7 @@ __constant__ float gC_angle[g_MaxAngles];
 // per-detector u/v shifts?
 
 
-__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, const SDimensions3D dims)
+__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, float fVoxSize, const SDimensions3D dims)
 {
 	float* projData = (float*)D_projData;
 	int angle = startAngle + blockIdx.y * g_anglesPerWeightBlock + threadIdx.y;
@@ -83,10 +83,17 @@ __global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsig
 
 	float fV = (startDetectorV - 0.5f*dims.iProjV + 0.5f) * fDetVSize + fZShift;
 
-	//const float fW = fCentralRayLength;
-	//const float fW = fCentralRayLength * (M_PI / 2.0f) / (float)dims.iProjAngles;
+	// Four contributions to the weighting factors:
+	// fCentralRayLength / fRayLength   : the main FDK preweighting factor
+	// fSrcOrigin / (fDetUSize * fCentralRayLength)
+	//                                  : to adjust the filter to the det width
+	// || u v s || ^ 2                  : see cone_bp.cu, FDKWEIGHT
+	// pi / (2 * iProjAngles)           : scaling of the integral over angles
+
 	const float fW1 = fSrcOrigin * fDetUSize * fDetVSize;
-	const float fW = fCentralRayLength * fW1 * fW1 * (M_PI / 2.0f) / (float)dims.iProjAngles;
+	const float fW2 = fCentralRayLength / (fDetUSize * fSrcOrigin);
+	const float fW3 = fVoxSize * fVoxSize;
+	const float fW = fCentralRayLength * fW1 * fW1 * fW2 * fW3 * (M_PI / 2.0f) / (float)dims.iProjAngles;
 
 	for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV)
 	{
@@ -156,7 +163,8 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
 bool FDK_PreWeight(cudaPitchedPtr D_projData,
                 float fSrcOrigin, float fDetOrigin,
                 float fZShift,
-                float fDetUSize, float fDetVSize, bool bShortScan,
+                float fDetUSize, float fDetVSize, float fVoxSize,
+				bool bShortScan,
                 const SDimensions3D& dims, const float* angles)
 {
 	// The pre-weighting factor for a ray is the cosine of the angle between
@@ -168,7 +176,7 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData,
 
 	int projPitch = D_projData.pitch/sizeof(float);
 
-	devFDK_preweight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, dims);
+	devFDK_preweight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, fVoxSize, dims);
 
 	cudaTextForceKernelsCompletion();
 
@@ -310,8 +318,9 @@ bool FDK(cudaPitchedPtr D_volumeData,
 
 
 #if 1
+	// NB: assuming cube voxels (params.fVolScaleX)
 	ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin,
-	                fZShift, fDetUSize, fDetVSize,
+	                fZShift, fDetUSize, fDetVSize, params.fVolScaleX,
 	                bShortScan, dims, pfAngles);
 #else
 	ok = true;
diff --git a/cuda/3d/fdk.h b/cuda/3d/fdk.h
index 0165af9..0f8d871 100644
--- a/cuda/3d/fdk.h
+++ b/cuda/3d/fdk.h
@@ -35,7 +35,8 @@ namespace astraCUDA3d {
 bool FDK_PreWeight(cudaPitchedPtr D_projData,
                 float fSrcOrigin, float fDetOrigin,
                 float fZShift,
-                float fDetUSize, float fDetVSize, bool bShortScan,
+                float fDetUSize, float fDetVSize, float fVoxSize,
+                bool bShortScan,
                 const SDimensions3D& dims, const float* angles);
 
 bool FDK(cudaPitchedPtr D_volumeData,
-- 
cgit v1.2.3


From d61ce8cb50cee2145d66a209cbeb2b07ae645355 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 30 Nov 2017 17:41:13 +0100
Subject: Adapt FBP_CUDA voxel-size weighting factors

---
 cuda/2d/astra.cu | 11 ++++++++---
 cuda/3d/fdk.cu   |  3 ++-
 2 files changed, 10 insertions(+), 4 deletions(-)

(limited to 'cuda')

diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index 0069bc3..41ce0c8 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -343,7 +343,7 @@ bool AstraFBP::run()
 
 		astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance,
 		              pData->fOriginDetectorDistance, 0.0f,
-		              pData->dims.fDetScale, 1.0f, 1.0f, // TODO: Are these correct?
+		              pData->dims.fDetScale, 1.0f, pData->fPixelSize, // TODO: Are these correct?
 		              pData->bShortScan, dims3d, pData->angles);
 	}
 
@@ -366,12 +366,17 @@ bool AstraFBP::run()
 
 	}
 
-	float fOutputScale = (M_PI / 2.0f) / (float)pData->dims.iProjAngles;
 
 	if (pData->bFanBeam) {
-		ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->fanProjections, fOutputScale);
+		float fOutputScale = 1.0 / (pData->fPixelSize * pData->fPixelSize * pData->fPixelSize * pData->dims.fDetScale * pData->dims.fDetScale);
 
+		ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->fanProjections, fOutputScale);
 	} else {
+		// scale by number of angles. For the fan-beam case, this is already
+		// handled by FDK_PreWeight
+		float fOutputScale = (M_PI / 2.0f) / (float)pData->dims.iProjAngles;
+		fOutputScale /= pData->dims.fDetScale * pData->dims.fDetScale;
+
 		ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets, fOutputScale);
 	}
 	if(!ok)
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 0630afe..11c68e1 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -83,12 +83,13 @@ __global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsig
 
 	float fV = (startDetectorV - 0.5f*dims.iProjV + 0.5f) * fDetVSize + fZShift;
 
-	// Four contributions to the weighting factors:
+	// Contributions to the weighting factors:
 	// fCentralRayLength / fRayLength   : the main FDK preweighting factor
 	// fSrcOrigin / (fDetUSize * fCentralRayLength)
 	//                                  : to adjust the filter to the det width
 	// || u v s || ^ 2                  : see cone_bp.cu, FDKWEIGHT
 	// pi / (2 * iProjAngles)           : scaling of the integral over angles
+	// fVoxSize ^ 2                     : ...
 
 	const float fW1 = fSrcOrigin * fDetUSize * fDetVSize;
 	const float fW2 = fCentralRayLength / (fDetUSize * fSrcOrigin);
-- 
cgit v1.2.3


From a9ee92fa1f5d72cf5fa178fbf6a6a547c95f7cc8 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 30 Nov 2017 17:54:15 +0100
Subject: Fix shortscan (FBP/FDK) scaling factor

---
 cuda/3d/fdk.cu | 2 ++
 1 file changed, 2 insertions(+)

(limited to 'cuda')

diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 11c68e1..08e1ebf 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -150,6 +150,8 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
 		fWeight = 0.0f;
 	}
 
+	fWeight *= 2; // adjust to effectively halved angular range
+
 	for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV)
 	{
 
-- 
cgit v1.2.3