From 9fa2ffdc5348f8f19de48d06a72b82bdc1ba8f22 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 30 Nov 2017 16:59:18 +0100 Subject: Start on fixing FDK and BP voxel-size weighting factors --- cuda/3d/fdk.cu | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) (limited to 'cuda/3d/fdk.cu') 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<<>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, dims); + devFDK_preweight<<>>(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; -- cgit v1.2.3 From d61ce8cb50cee2145d66a209cbeb2b07ae645355 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 30 Nov 2017 17:41:13 +0100 Subject: Adapt FBP_CUDA voxel-size weighting factors --- cuda/3d/fdk.cu | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'cuda/3d/fdk.cu') 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 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/3d/fdk.cu') 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