X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=2b9eba45cdc4ec7ce4b300030c19a4fafaf5af1e;hp=06ad159b9d9577d73f099aee7a9586a3e19b6052;hb=352a0691b9bd67f6c93ea822b353be3c101f4adb;hpb=de411914da8b157958e9caae917bf1edeafbb713 diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 06ad159..2b9eba4 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -1,14 +1,12 @@ /***************************************************************************** ** FILE IDENTIFICATION ** -** Name: backprojectors.cpp Classes for backprojection +** Name: backprojectors.cpp Classes for backprojection ** Programmer: Kevin Rosenberg ** Date Started: June 2000 ** ** This is part of the CTSim program -** Copyright (c) 1983-2001 Kevin Rosenberg -** -** $Id: backprojectors.cpp,v 1.31 2001/03/11 15:27:30 kevin Exp $ +** Copyright (c) 1983-2009 Kevin Rosenberg ** ** 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 @@ -25,6 +23,7 @@ ******************************************************************************/ #include "ct.h" +#include "interpolator.h" const int Backprojector::BPROJ_INVALID = -1; const int Backprojector::BPROJ_TRIG = 0; @@ -32,20 +31,20 @@ const int Backprojector::BPROJ_TABLE = 1; const int Backprojector::BPROJ_DIFF = 2; const int Backprojector::BPROJ_IDIFF = 3; -const char* const Backprojector::s_aszBackprojectName[] = +const char* const Backprojector::s_aszBackprojectName[] = { - {"trig"}, - {"table"}, - {"diff"}, - {"idiff"}, + "trig", + "table", + "diff", + "idiff", }; -const char* const Backprojector::s_aszBackprojectTitle[] = +const char* const Backprojector::s_aszBackprojectTitle[] = { - {"Direct Trigometric"}, - {"Trigometric Table"}, - {"Difference Iteration"}, - {"Integer Difference Iteration"}, + "Direct Trigometric", + "Trigometric Table", + "Difference Iteration", + "Integer Difference Iteration", }; const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*); @@ -62,35 +61,35 @@ const int Backprojector::INTERP_2BSPLINE = 6; const int Backprojector::INTERP_3BSPLINE = 7; #endif -const char* const Backprojector::s_aszInterpName[] = +const char* const Backprojector::s_aszInterpName[] = { - {"nearest"}, - {"linear"}, - {"cubic"}, + "nearest", + "linear", + "cubic", #if HAVE_FREQ_PREINTERP - {"freq_preinterpolationj"}, + "freq_preinterpolationj", #endif #if HAVE_BSPLINE_INTERP - {"bspline"}, - {"1bspline"}, - {"2bspline"}, - {"3bspline"}, + "bspline", + "1bspline", + "2bspline", + "3bspline", #endif }; -const char* const Backprojector::s_aszInterpTitle[] = +const char* const Backprojector::s_aszInterpTitle[] = { - {"Nearest"}, - {"Linear"}, - {"Cubic"}, + "Nearest", + "Linear", + "Cubic", #if HAVE_FREQ_PREINTERP - {"Frequency Preinterpolation"}, + "Frequency Preinterpolation", #endif #if HAVE_BSPLINE_INTERP - {"B-Spline"}, - {"B-Spline 1st Order"}, - {"B-Spline 2nd Order"}, - {"B-Spline 3rd Order"}, + "B-Spline", + "B-Spline 1st Order", + "B-Spline 2nd Order", + "B-Spline 3rd Order", #endif }; @@ -98,23 +97,23 @@ const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const -Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, +Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor, const ReconstructionROI* pROI) { m_fail = false; m_pBackprojectImplem = NULL; - + initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI); } -void +void Backprojector::BackprojectView (const double* const viewData, const double viewAngle) { if (m_pBackprojectImplem != NULL) m_pBackprojectImplem->BackprojectView (viewData, viewAngle); } -void +void Backprojector::PostProcessing() { if (m_pBackprojectImplem != NULL) @@ -130,11 +129,11 @@ Backprojector::~Backprojector () // Backproject* projector = selectBackprojector (...) // // PURPOSE -// Selects a backprojector based on BackprojType +// Selects a backprojector based on BackprojType // and initializes the backprojector bool -Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, +Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor, const ReconstructionROI* pROI) { m_nameBackproject = backprojName; @@ -152,15 +151,15 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_failMessage = "Invalid interpolation name "; m_failMessage += interpName; } - + if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) { m_fail = true; return false; } - + if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI)); - else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) + else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI)); else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) { if (m_idBackproject == BPROJ_TRIG) @@ -176,7 +175,7 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; return false; } - + return true; } @@ -185,24 +184,24 @@ int Backprojector::convertBackprojectNameToID (const char* const backprojName) { int backprojID = BPROJ_INVALID; - - for (int i = 0; i < s_iBackprojectCount; i++) + + for (int i = 0; i < s_iBackprojectCount; i++) { if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { backprojID = i; break; } - - return (backprojID); + } + return (backprojID); } const char* Backprojector::convertBackprojectIDToName (int bprojID) { static const char *bprojName = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) return (s_aszBackprojectName[bprojID]); - + return (bprojName); } @@ -210,10 +209,10 @@ const char* Backprojector::convertBackprojectIDToTitle (const int bprojID) { static const char *bprojTitle = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) return (s_aszBackprojectTitle[bprojID]); - + return (bprojTitle); } @@ -222,24 +221,24 @@ int Backprojector::convertInterpNameToID (const char* const interpName) { int interpID = INTERP_INVALID; - - for (int i = 0; i < s_iInterpCount; i++) + + for (int i = 0; i < s_iInterpCount; i++) { if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { interpID = i; break; } - - return (interpID); + } + return (interpID); } const char* Backprojector::convertInterpIDToName (const int interpID) { static const char *interpName = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) return (s_aszInterpName[interpID]); - + return (interpName); } @@ -247,10 +246,10 @@ const char* Backprojector::convertInterpIDToTitle (const int interpID) { static const char *interpTitle = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) return (s_aszInterpTitle[interpID]); - + return (interpTitle); } @@ -262,32 +261,31 @@ Backprojector::convertInterpIDToTitle (const int interpID) // PURPOSE // Pure virtual base class for all backprojectors. -Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor, +Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor, const ReconstructionROI* pROI) : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false) { detInc = proj.detInc(); nDet = proj.nDet(); - iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection - rotScale = proj.rotInc(); - + iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection + if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) - rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations + rotScale = PI / proj.nView(); // scale by number of PI rotations else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) - rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations + rotScale = (2 * PI) / proj.nView(); // scale by number of 2PI rotations else sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry()); - + v = im.getArray(); nx = im.nx(); ny = im.ny(); im.arrayDataClear(); - + xMin = -proj.phmLen() / 2; // Retangular coords of phantom xMax = xMin + proj.phmLen(); yMin = -proj.phmLen() / 2; yMax = yMin + proj.phmLen(); - + if (pROI) { if (pROI->m_dXMin > xMin) xMin = pROI->m_dXMin; @@ -310,9 +308,12 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType } } - xInc = (xMax - xMin) / nx; // size of cells + xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; - + + im.setAxisIncrement (xInc, yInc); + im.setAxisExtent (xMin, xMax, yMin, yMax); + m_dFocalLength = proj.focalLength(); m_dSourceDetectorLength = proj.sourceDetectorLength(); } @@ -349,7 +350,7 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n"; os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n"; os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";; - + sys_error (ERR_WARNING, os.str().c_str()); #endif } @@ -365,33 +366,38 @@ void BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); - - double x = xMin + xInc / 2; // Rectang coords of center of pixel - for (int ix = 0; ix < nx; x += xInc, ix++) { + + double xstart = xMin + xInc / 2; // Rectang coords of center of pixel +#if HAVE_OPENMP + #pragma omp parallel for +#endif + for (int ix = 0; ix < nx; ix++) { + double x = xstart + (ix * xInc); + double y = yMin + yInc / 2; for (int iy = 0; iy < ny; y += yInc, iy++) { double r = sqrt (x * x + y * y); // distance of cell from center double phi = atan2 (y, x); // angle of cell from center double L = r * cos (theta - phi); // position on detector - + if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest (L / detInc); // calc'd index in the filter raysum array - + if (iDetPos >= 0 && iDetPos < nDet) v[ix][iy] += rotScale * filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - double p = L / detInc; // position along detector + double p = L / detInc; // position along detector double pFloor = floor (p); int iDetPos = iDetCenter + static_cast(pFloor); - double frac = p - pFloor; // fraction distance from det + double frac = p - pFloor; // fraction distance from det if (iDetPos >= 0 && iDetPos < nDet - 1) v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } else if (interpType == Backprojector::INTERP_CUBIC) { - double p = iDetCenter + (L / detInc); // position along detector + double p = iDetCenter + (L / detInc); // position along detector if (p >= 0 && p < nDet) v[ix][iy] += rotScale * pCubicInterp->interpolate (p); } @@ -400,7 +406,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; -} +} // CLASS IDENTICATION @@ -409,7 +415,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // PURPOSE // Precalculates trigometric function value for each point in image for backprojection. -BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor, const ReconstructionROI* pROI) : Backproject (proj, im, interpType, interpFactor, pROI) { @@ -417,14 +423,20 @@ BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int arrayPhi.initSetSize (im.nx(), im.ny()); r = arrayR.getArray(); phi = arrayPhi.getArray(); - - double x, y; // Rectang coords of center of pixel - int ix, iy; - for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++) - for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) { + + double xstart = xMin + xInc / 2; + +#if HAVE_OPENMP + #pragma omp parallel for +#endif + for (int ix = 0; ix < nx; ix++) { + double x = xstart + (ix * xInc); + double y = yMin + yInc / 2; + for (int iy = 0; iy < ny; iy++, y += yInc) { r[ix][iy] = sqrt (x * x + y * y); phi[ix][iy] = atan2 (y, x); } + } } BackprojectTable::~BackprojectTable () @@ -444,36 +456,42 @@ void BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); - + +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++) { double L = r[ix][iy] * cos (theta - phi[ix][iy]); - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(L / detInc); // calc index in the filtered raysum vector - - if (iDetPos >= 0 && iDetPos < nDet) + int iDetPos = iDetCenter + nearest(L / detInc); // calc index in the filtered raysum vector + + if (iDetPos >= 0 && iDetPos < nDet) { pImCol[iy] += filteredProj[iDetPos]; + } } else if (interpType == Backprojector::INTERP_LINEAR) { - double dPos = L / detInc; // position along detector + double dPos = L / detInc; // position along detector double dPosFloor = floor (dPos); int iDetPos = iDetCenter + static_cast(dPosFloor); - double frac = dPos - dPosFloor; // fraction distance from det - if (iDetPos >= 0 && iDetPos < nDet - 1) + double frac = dPos - dPosFloor; // fraction distance from det + if (iDetPos >= 0 && iDetPos < nDet - 1) { pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + } } else if (interpType == Backprojector::INTERP_CUBIC) { - double p = iDetCenter + (L / detInc); // position along detector - if (p >= 0 && p < nDet) + double p = iDetCenter + (L / detInc); // position along detector + if (p >= 0 && p < nDet) { pImCol[iy] += pCubicInterp->interpolate (p); + } } - } // end for y - } // end for x + } // end for y + } // end for x if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; @@ -487,16 +505,16 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // Backprojects by precalculating the change in L position for each x & y step in the image. // Iterates in x & y direction by adding difference in L position -BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, +BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor, const ReconstructionROI* pROI) : Backproject (proj, im, interpType, interpFactor, pROI) { - // calculate center of first pixel v[0][0] + // calculate center of first pixel v[0][0] double x = xMin + xInc / 2; double y = yMin + yInc / 2; start_r = sqrt (x * x + y * y); start_phi = atan2 (y, x); - + im.arrayDataClear(); } @@ -517,43 +535,60 @@ void BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - - // Distance between detectors for an angle given in units of detectors + + // Distance between detectors for an angle given in units of detectors double det_dx = xInc * cos (theta) / detInc; double det_dy = yInc * sin (theta) / detInc; - - // calculate detPosition for first point in image (ix=0, iy=0) - double detPosColStart = start_r * cos (theta - start_phi) / detInc; - + + // calculate detPosition for first point in image (ix=0, iy=0) + double detPosColBase = iDetCenter + start_r * cos (theta - start_phi) / detInc; + CubicPolyInterpolator* pCubicInterp = NULL; - if (interpType == Backprojector::INTERP_CUBIC) + double* deltaFilteredProj = NULL; + if (interpType == Backprojector::INTERP_LINEAR) { + // precalculate scaled difference for linear interpolation + deltaFilteredProj = new double [nDet]; + for (int i = 0; i < nDet - 1; i++) + deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i]; + deltaFilteredProj[nDet - 1] = 0; // last detector + } else if (interpType == Backprojector::INTERP_CUBIC) { pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); - - for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - double curDetPos = detPosColStart; + } + + int iLastDet = nDet - 1; + +#if HAVE_OPENMP + #pragma omp parallel for +#endif + for (int ix = 0; ix < nx; ix++) { + double detPos = detPosColBase + (ix * det_dx); ImageFileColumn pImCol = v[ix]; - - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { + + for (int iy = 0; iy < ny; iy++, detPos += det_dy) { if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest (curDetPos); // calc index in the filtered raysum vector - - if (iDetPos >= 0 && iDetPos < nDet) + int iDetPos = nearest (detPos); // calc index in the filtered raysum vector + if (iDetPos >= 0 && iDetPos < nDet) { *pImCol++ += filteredProj[iDetPos]; + } } else if (interpType == Backprojector::INTERP_LINEAR) { - double detPosFloor = floor (curDetPos); - int iDetPos = iDetCenter + static_cast(detPosFloor); - double frac = curDetPos - detPosFloor; // fraction distance from det - if (iDetPos > 0 && iDetPos < nDet - 1) - *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); + double detPosFloor = floor (detPos); + int iDetPos = static_cast(detPosFloor); + double frac = detPos - detPosFloor; // fraction distance from det + if (iDetPos >= 0 && iDetPos <= iLastDet) { + *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]); + } } else if (interpType == Backprojector::INTERP_CUBIC) { - double p = iDetCenter + curDetPos; // position along detector - if (p >= 0 && p < nDet) + double p = detPos; // position along detector + if (p >= 0 && p < nDet) { *pImCol++ += pCubicInterp->interpolate (p); + } } - } // end for y - } // end for x + } // end for y + } // end for x - if (interpType == Backprojector::INTERP_CUBIC) + if (interpType == Backprojector::INTERP_LINEAR) + delete deltaFilteredProj; + else if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; } @@ -568,19 +603,23 @@ void BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; // add half PI to view angle to get perpendicular theta angle +#if SIZEOF_LONG == 4 static const int scaleShift = 16; - static const kint32 scale = (1 << scaleShift); - static const kint32 scaleBitmask = scale - 1; - static const kint32 halfScale = scale / 2; +#elif SIZEOF_LONG == 8 + static const int scaleShift = 32; +#endif + static const long scale = (1L << scaleShift); + static const long scaleBitmask = scale - 1; + static const long halfScale = scale / 2; static const double dInvScale = 1. / scale; - - const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); - const kint32 det_dy = nearest (yInc * sin (theta) / detInc * scale); - - // calculate L for first point in image (0, 0) - kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); - - double* deltaFilteredProj = NULL; + + const long det_dx = nearest (xInc * cos (theta) / detInc * scale); + const long det_dy = nearest (yInc * sin (theta) / detInc * scale); + + // calculate L for first point in image (0, 0) + long detPosColBase = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + + double* deltaFilteredProj = NULL; CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_LINEAR) { // precalculate scaled difference for linear interpolation @@ -591,38 +630,41 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou } else if (interpType == Backprojector::INTERP_CUBIC) { pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); } - + int iLastDet = nDet - 1; - for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - kint32 curDetPos = detPosColStart; +#if HAVE_OPENMP + #pragma omp parallel for +#endif + for (int ix = 0; ix < nx; ix++) { + long detPos = detPosColBase + (ix * det_dx); ImageFileColumn pImCol = v[ix]; - - if (interpType == Backprojector::INTERP_NEAREST) { - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const int iDetPos = (curDetPos + halfScale) >> 16; - if (iDetPos >= 0 && iDetPos <= iLastDet) + + for (int iy = 0; iy < ny; iy++, detPos += det_dy) { + if (interpType == Backprojector::INTERP_NEAREST) { + const int iDetPos = (detPos + halfScale) >> scaleShift; + if (iDetPos >= 0 && iDetPos <= iLastDet) { *pImCol++ += filteredProj[iDetPos]; - } // end for iy - } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) { - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor; - if (iDetPos >= 0 && iDetPos <= iLastDet) + } else + pImCol++; + } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) { + const int iDetPos = ((detPos + halfScale) >> scaleShift) * m_interpFactor; + if (iDetPos >= 0 && iDetPos <= iLastDet) { *pImCol++ += filteredProj[iDetPos]; - } // end for iy - } else if (interpType == Backprojector::INTERP_LINEAR) { - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const kint32 iDetPos = curDetPos >> scaleShift; - const kint32 detRemainder = curDetPos & scaleBitmask; - if (iDetPos >= 0 && iDetPos <= iLastDet) + } else + pImCol++; + } else if (interpType == Backprojector::INTERP_LINEAR) { + const long iDetPos = detPos >> scaleShift; + if (iDetPos >= 0 && iDetPos <= iLastDet) { + const long detRemainder = detPos & scaleBitmask; *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); - } // end for iy - } else if (interpType == Backprojector::INTERP_CUBIC) { - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - *pImCol++ += pCubicInterp->interpolate (static_cast(curDetPos) / 65536); - } - } // end Cubic + } else + pImCol++; + } else if (interpType == Backprojector::INTERP_CUBIC) { + *pImCol++ += pCubicInterp->interpolate (static_cast(detPos) / scale); + } // end Cubic + } // end for iy } // end for ix - + if (interpType == Backprojector::INTERP_LINEAR) delete deltaFilteredProj; else if (interpType == Backprojector::INTERP_CUBIC) @@ -634,15 +676,18 @@ void BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle) { double beta = view_angle; - + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); - + +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - - for (int iy = 0; iy < ny; iy++) { + + for (int iy = 0; iy < ny; iy++) { double dAngleDiff = beta - phi[ix][iy]; double rcos_t = r[ix][iy] * cos (dAngleDiff); double rsin_t = r[ix][iy] * sin (dAngleDiff); @@ -650,24 +695,24 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double gamma = atan (rcos_t / dFLPlusSin); double dPos = gamma / detInc; // position along detector double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t); - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(dPos); // calc index in the filtered raysum vector + int iDetPos = iDetCenter + nearest(dPos); // calc index in the filtered raysum vector if (iDetPos >= 0 && iDetPos < nDet) pImCol[iy] += filteredProj[iDetPos] / dL2; } else if (interpType == Backprojector::INTERP_LINEAR) { double dPosFloor = floor (dPos); int iDetPos = iDetCenter + static_cast(dPosFloor); - double frac = dPos - dPosFloor; // fraction distance from det + double frac = dPos - dPosFloor; // fraction distance from det if (iDetPos >= 0 && iDetPos < nDet - 1) pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2; } else if (interpType == Backprojector::INTERP_CUBIC) { - double d = iDetCenter + dPos; // position along detector + double d = iDetCenter + dPos; // position along detector if (d >= 0 && d < nDet) pImCol[iy] += pCubicInterp->interpolate (d) / dL2; } - } // end for y - } // end for x + } // end for y + } // end for x if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; @@ -677,43 +722,47 @@ void BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle) { double beta = view_angle; - + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); - + +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++) { double dAngleDiff = beta - phi[ix][iy]; double rcos_t = r[ix][iy] * cos (dAngleDiff); double rsin_t = r[ix][iy] * sin (dAngleDiff); - + double dU = (m_dFocalLength + rsin_t) / m_dFocalLength; + double dU2 = dU * dU; + double dDetPos = rcos_t / dU; - // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22. - dDetPos *= m_dSourceDetectorLength / m_dFocalLength; - double dPos = dDetPos / detInc; // position along detector array + // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22. + dDetPos *= m_dSourceDetectorLength / m_dFocalLength; + double dPos = dDetPos / detInc; // position along detector array if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(dPos); // calc index in the filtered raysum vector - if (iDetPos >= 0 && iDetPos < nDet) - pImCol[iy] += (filteredProj[iDetPos] / (dU * dU)); + int iDetPos = iDetCenter + nearest(dPos); // calc index in the filtered raysum vector + if (iDetPos >= 0 && iDetPos < nDet) + pImCol[iy] += filteredProj[iDetPos] / dU2; } else if (interpType == Backprojector::INTERP_LINEAR) { double dPosFloor = floor (dPos); int iDetPos = iDetCenter + static_cast(dPosFloor); - double frac = dPos - dPosFloor; // fraction distance from det + double frac = dPos - dPosFloor; // fraction distance from det if (iDetPos >= 0 && iDetPos < nDet - 1) - pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) - / (dU * dU); + pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dU2; } else if (interpType == Backprojector::INTERP_CUBIC) { - double d = iDetCenter + dPos; // position along detector + double d = iDetCenter + dPos; // position along detector if (d >= 0 && d < nDet) - pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU); + pImCol[iy] += pCubicInterp->interpolate (d) / dU2; } - } // end for y - } // end for x + } // end for y + } // end for x if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp;