diff options
| -rw-r--r-- | cuda/2d/fan_bp.cu | 11 | ||||
| -rw-r--r-- | cuda/2d/fbp.cu | 6 | ||||
| -rw-r--r-- | src/GeometryUtil2D.cpp | 9 |
3 files changed, 15 insertions, 11 deletions
diff --git a/cuda/2d/fan_bp.cu b/cuda/2d/fan_bp.cu index 0d21897..428485c 100644 --- a/cuda/2d/fan_bp.cu +++ b/cuda/2d/fan_bp.cu @@ -271,11 +271,14 @@ __global__ void devFanBP_FBPWeighted(float* D_volData, unsigned int volPitch, un const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX; const float fDen = fDetUX * fYD - fDetUY * fXD; + const float fr = __fdividef(1.0f, fDen); + + // fDen = || u (x-s) || + // Required scale factor is ( || u s || / || u (x-s) || ) ^ 2. + // The factor || u s || ^ 2 is handled by the preweighting - const float fWeight = fXD*fXD + fYD*fYD; - - const float fT = fNum / fDen; - fVal += tex2D(gT_FanProjTexture, fT, fA) / fWeight; + const float fT = fNum * fr; + fVal += tex2D(gT_FanProjTexture, fT, fA) * fr * fr; fA += 1.0f; } diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index f0edc19..28cdd92 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -301,7 +301,7 @@ bool FBP::iterate(unsigned int iterations) astraCUDA3d::FDK_PreWeight(tmp, fOriginSource, fOriginDetector, 0.0f, - fFanDetSize, 1.0f, /* fPixelSize */ 1.0f, + fFanDetSize, 1.0f, /* fPixelSize, but is normalized */ 1.0f, m_bShortScan, dims3d, pfAngles); } else { // TODO: How should different detector pixel size in different @@ -328,9 +328,7 @@ bool FBP::iterate(unsigned int iterations) } if (fanProjs) { - float fOutputScale = 1.0 / (/*fPixelSize * fPixelSize * fPixelSize * */ fFanDetSize * fFanDetSize); - - ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale * fProjectorScale * fReconstructionScale); + ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fProjectorScale * fReconstructionScale); } else { // scale by number of angles. For the fan-beam case, this is already diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp index e09a3bc..806572f 100644 --- a/src/GeometryUtil2D.cpp +++ b/src/GeometryUtil2D.cpp @@ -28,6 +28,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include "astra/GeometryUtil2D.h" #include <cmath> +#include <cstdio> namespace astra { @@ -158,14 +159,16 @@ bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float // project origin on detector line ( == project source on detector line) double t = (- proj.fDetSX) * proj.fDetUX + (- proj.fDetSY) * proj.fDetUY; + t /= (proj.fDetUX * proj.fDetUX + proj.fDetUY * proj.fDetUY); fOffset = (float)t - 0.5*iProjDets; - // TODO: CHECKME fOriginDetector = sqrt((proj.fDetSX + t * proj.fDetUX)*(proj.fDetSX + t * proj.fDetUX) + (proj.fDetSY + t * proj.fDetUY)*(proj.fDetSY + t * proj.fDetUY)); - //float fAngle = atan2(proj.fDetSX + t * proj.fDetUX - proj.fSrcX, proj.fDetSY + t * proj.fDetUY); // TODO: Fix order + sign - fAngle = atan2(proj.fDetUY, proj.fDetUX); // TODO: Check order + sign + fAngle = atan2(proj.fDetUY, proj.fDetUX); + + //fprintf(stderr, "getFanParams: s = (%f,%f) d = (%f,%f) u = (%f,%f)\n", proj.fSrcX, proj.fSrcY, proj.fDetSX, proj.fDetSY, proj.fDetUX, proj.fDetUY); + //fprintf(stderr, "getFanParams: fOS = %f, fOD = %f, detsize = %f, offset = %f (t = %f), angle = %f\n", fOriginSource, fOriginDetector, fDetSize, fOffset, t, fAngle); return true; } |