diff options
| -rw-r--r-- | cuda/3d/fdk.cu | 21 |
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; |