summaryrefslogtreecommitdiffstats
path: root/cuda/3d/fdk.cu
diff options
context:
space:
mode:
Diffstat (limited to 'cuda/3d/fdk.cu')
-rw-r--r--cuda/3d/fdk.cu21
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;