diff options
| -rw-r--r-- | cuda/2d/fbp.cu | 16 | ||||
| -rw-r--r-- | cuda/3d/cone_bp.cu | 21 | ||||
| -rw-r--r-- | cuda/3d/fdk.cu | 26 | ||||
| -rw-r--r-- | cuda/3d/fdk.h | 3 |
4 files changed, 50 insertions, 16 deletions
diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index e7ed277..ecaf544 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -260,6 +260,7 @@ bool FBP::iterate(unsigned int iterations) bool ok = false; + float fFanDetSize = 0.0f; if (fanProjs) { // Call FDK_PreWeight to handle fan beam geometry. We treat // this as a cone beam setup of a single slice: @@ -271,12 +272,12 @@ bool FBP::iterate(unsigned int iterations) float *pfAngles = new float[dims.iProjAngles]; - float fOriginSource, fOriginDetector, fDetSize, fOffset; + float fOriginSource, fOriginDetector, fOffset; for (unsigned int i = 0; i < dims.iProjAngles; ++i) { bool ok = astra::getFanParameters(fanProjs[i], dims.iProjDets, pfAngles[i], fOriginSource, fOriginDetector, - fDetSize, fOffset); + fFanDetSize, fOffset); if (!ok) { ASTRA_ERROR("FBP_CUDA: Failed to extract circular fan beam parameters from fan beam geometry"); return false; @@ -300,7 +301,7 @@ bool FBP::iterate(unsigned int iterations) astraCUDA3d::FDK_PreWeight(tmp, fOriginSource, fOriginDetector, 0.0f, - fDetSize, 1.0f, + fFanDetSize, 1.0f, /* fPixelSize */ 1.0f, m_bShortScan, dims3d, pfAngles); } else { // TODO: How should different detector pixel size in different @@ -326,12 +327,17 @@ bool FBP::iterate(unsigned int iterations) } - float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles; - if (fanProjs) { + float fOutputScale = 1.0 / (/*fPixelSize * fPixelSize * fPixelSize * */ fFanDetSize * fFanDetSize); + ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, 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)dims.iProjAngles; + //fOutputScale /= fDetSize * fDetSize; + ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale); } if(!ok) diff --git a/cuda/3d/cone_bp.cu b/cuda/3d/cone_bp.cu index 3dd0c97..2d12d00 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; @@ -255,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 4e926f2..48194c4 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,18 @@ __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; + // 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 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) { @@ -142,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) { @@ -156,7 +166,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 +179,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 +321,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 f4049e6..6f6e73b 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, |