Adjust Parker weights to angular range

This fixes a global scaling in shortscan reconstruction.
Also add some debugging logs for the Parker weighting and a warning if
the angular range is too small.

Issue #229

1 files changed, 16 insertions, 5 deletions

diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 0b8d2ab..e549670 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -103,7 +103,7 @@ __global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsig
 	}
 }
 
-__global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fDetUSize, float fCentralFanAngle, const SDimensions3D dims)
+__global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fDetUSize, float fCentralFanAngle, float fScale, const SDimensions3D dims)
 {
 	float* projData = (float*)D_projData;
 	int angle = startAngle + blockIdx.y * g_anglesPerWeightBlock + threadIdx.y;
@@ -145,7 +145,7 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
 		fWeight = 0.0f;
 	}
 
-	fWeight *= 2; // adjust to effectively halved angular range
+	fWeight *= fScale;
 
 	for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV)
 	{
@@ -197,9 +197,11 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData,
 		if (fdA >= 0.0f) {
 			// going up from angles[0]
 			fAngleBase = angles[0];
+			ASTRA_DEBUG("Second angle >= first angle, so assuming angles are incrementing");
 		} else {
 			// going down from angles[0]
 			fAngleBase = angles[dims.iProjAngles - 1];
+			ASTRA_DEBUG("Second angle < first angle, so assuming angles are decrementing");
 		}
 
 		// We pick the lowest end of the range, and then
@@ -215,16 +217,25 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData,
 			fRelAngles[i] = f;
 		}
 
+		float fRange = fabs(fRelAngles[dims.iProjAngles-1] - fRelAngles[0]);
+		ASTRA_DEBUG("Assuming angles are linearly ordered and equally spaced for Parker weighting. Angular range %f radians", fRange);
+		float fScale = fRange / M_PI;
+
 		cudaError_t e1 = cudaMemcpyToSymbol(gC_angle, fRelAngles,
 		                                    dims.iProjAngles*sizeof(float), 0,
 		                                    cudaMemcpyHostToDevice);
 		assert(!e1);
 		delete[] fRelAngles;
 
-		float fCentralFanAngle = atanf(fDetUSize * (dims.iProjU*0.5f) /
-		                               (fSrcOrigin + fDetOrigin));
+		float fCentralFanAngle = fabs(atanf(fDetUSize * (dims.iProjU*0.5f) /
+		                               (fSrcOrigin + fDetOrigin)));
+
+		// Check range, but take possible discretisation and rounding into affect
+		if (fRange + 1e-3 < (M_PI + 2*fCentralFanAngle) * (dims.iProjAngles - 1)/ dims.iProjAngles) {
+			ASTRA_WARN("Angular range (%f rad) smaller than Parker weighting range (%f rad)", fRange, M_PI + 2*fCentralFanAngle);
+		}
 
-		devFDK_ParkerWeight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fDetUSize, fCentralFanAngle, dims);
+		devFDK_ParkerWeight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fDetUSize, fCentralFanAngle, fScale, dims);
 
 		if (!checkCuda(cudaThreadSynchronize(), "FDK_PreWeight ParkerWeight"))
 			return false;