X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=d3cb2035dc125af3415006a2d204cc61b5cc7907;hb=6850134e4711a842cc7c0a306a9e4243df06b952;hp=1f6e92c1bece4c526d80abdad76e47d1d70a7f75;hpb=08f34bf3ba14d4f436f4d2ef0ee5af1d6eb266ac;p=ctsim.git diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 1f6e92c..d3cb203 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.4 2000/06/25 17:32:24 kevin Exp $ +** $Id: backprojectors.cpp,v 1.6 2000/07/09 08:16:17 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -87,6 +87,8 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation)); else if (m_idBackproject == BPROJ_IDIFF2) m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_IDIFF3) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation)); else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -112,6 +114,8 @@ Backprojector::convertBackprojectNameToID (const char* const backprojName) backprojID = BPROJ_DIFF2; else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0) backprojID = BPROJ_IDIFF2; + else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0) + backprojID = BPROJ_IDIFF3; return (backprojID); } @@ -131,6 +135,8 @@ Backprojector::convertBackprojectIDToName (const BackprojectID bprojID) bprojName = BPROJ_DIFF2_STR; else if (bprojID == BPROJ_IDIFF2) bprojName = BPROJ_IDIFF2_STR; + else if (bprojID == BPROJ_IDIFF3) + bprojName = BPROJ_IDIFF3_STR; return (bprojName); } @@ -146,6 +152,8 @@ Backprojector::convertInterpolationNameToID (const char* const interpName) interpID = INTERP_NEAREST; else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0) interpID = INTERP_LINEAR; + else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATE_STR) == 0) + interpID = INTERP_FREQ_PREINTERPOLATE; #if HAVE_BSPLINE_INTERP else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0) interpID = INTERP_BSPLINE; @@ -174,6 +182,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID) return (INTERP_NEAREST_STR); else if (interpID == INTERP_LINEAR) return (INTERP_LINEAR_STR); + else if (interpID == INTERP_FREQ_PREINTERPOLATE) + return (INTERP_FREQ_PREINTERPOLATE_STR); #if HAVE_BSPLINE_INTERP else if (interpID == INTERP_BSPLINE) return (INTERP_BSPLINE_STR); @@ -465,12 +475,12 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do { double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle - kint32 scale = 1 << 16; - double dScale = scale; - kint32 halfScale = scale / 2; + static const kint32 scale = 1 << 16; + static const double dScale = scale; + static const kint32 halfScale = scale / 2; - kint32 det_dx = nearest (xInc * sin (theta) / detInc * scale); - kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); + const kint32 det_dx = nearest (xInc * sin (theta) / detInc * scale); + const kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); // calculate L for first point in image (0, 0) kint32 detPosColStart = nearest (start_r * cos (theta - start_phi) / detInc * scale); @@ -500,8 +510,59 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do if (iDetPos < 0 || iDetPos >= nDet - 1) errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else - *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x } + +// CLASS IDENTICATION +// BackprojectIntDiff3 +// +// PURPOSE +// Highly optimized version of BackprojectIntDiff2 + +void +BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle) +{ + double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle + static const int scaleShift = 16; + static const kint32 scale = (1 << scaleShift); + static const kint32 scaleBitmask = scale - 1; + static const double dScale = scale; + static const kint32 halfScale = scale / 2; + static const double dInvScale = 1. / scale; + + const kint32 det_dx = nearest (xInc * sin (theta) / detInc * scale); + const kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); + + // calculate L for first point in image (0, 0) + kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + + for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { + kint32 curDetPos = detPosColStart; + ImageFileColumn pImCol = v[ix]; + + if (interpType == Backprojector::INTERP_NEAREST) { + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { + int iDetPos = (curDetPos + halfScale) >> 16; + assert(iDetPos >= 0 && iDetPos < nDet); + *pImCol++ += filteredProj[iDetPos]; + } // end for iy + } else if (interpType == Backprojector::INTERP_LINEAR) { + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { + kint32 iDetPos = curDetPos >> scaleShift; + kint32 detRemainder = curDetPos & scaleBitmask; +#if 0 + double frac = detRemainder * (1. / 65536.); + assert(iDetPos >= 0 && iDetPos < nDet - 1); + *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]); +#else + assert(iDetPos >= 0 && iDetPos < nDet - 1); + const double* const detPointer = &filteredProj[iDetPos]; + *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale; +#endif + } // end for iy + } //end linear + } // end for ix +}