X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=2b9eba45cdc4ec7ce4b300030c19a4fafaf5af1e;hp=d38b7189f6b6c878b09ec0549dca5daf242ce24e;hb=352a0691b9bd67f6c93ea822b353be3c101f4adb;hpb=e46615378520c3e801074560ebdc36be9558c6dd diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index d38b718..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$ +** 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 @@ -33,7 +31,7 @@ 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", @@ -41,7 +39,7 @@ const char* const Backprojector::s_aszBackprojectName[] = "idiff", }; -const char* const Backprojector::s_aszBackprojectTitle[] = +const char* const Backprojector::s_aszBackprojectTitle[] = { "Direct Trigometric", "Trigometric Table", @@ -63,7 +61,7 @@ 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", @@ -79,7 +77,7 @@ const char* const Backprojector::s_aszInterpName[] = #endif }; -const char* const Backprojector::s_aszInterpTitle[] = +const char* const Backprojector::s_aszInterpTitle[] = { "Nearest", "Linear", @@ -99,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) @@ -131,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; @@ -153,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) @@ -177,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; } @@ -186,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); } @@ -211,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); } @@ -223,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); } @@ -248,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); } @@ -263,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; @@ -311,9 +308,9 @@ 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); @@ -353,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 } @@ -369,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); } @@ -404,7 +406,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; -} +} // CLASS IDENTICATION @@ -413,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) { @@ -421,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 () @@ -448,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; @@ -491,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(); } @@ -521,16 +535,16 @@ 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 = iDetCenter + 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; - double* deltaFilteredProj = NULL; + double* deltaFilteredProj = NULL; if (interpType == Backprojector::INTERP_LINEAR) { // precalculate scaled difference for linear interpolation deltaFilteredProj = new double [nDet]; @@ -540,31 +554,37 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double } 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) { - double curDetPos = detPosColStart; + +#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 = 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); + double detPosFloor = floor (detPos); int iDetPos = static_cast(detPosFloor); - double frac = curDetPos - detPosFloor; // fraction distance from det - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]); + 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 = 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_LINEAR) delete deltaFilteredProj; @@ -588,18 +608,18 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou #elif SIZEOF_LONG == 8 static const int scaleShift = 32; #endif - static const long scale = (1 << scaleShift); + 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 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 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); - - double* deltaFilteredProj = NULL; + + // 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 @@ -610,39 +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) { - long 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) >> scaleShift; - 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) >> scaleShift) * 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 long iDetPos = curDetPos >> scaleShift; - const long detRemainder = curDetPos & scaleBitmask; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *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) / scale); - } - } // end Cubic + } 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]); + } 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) @@ -654,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); @@ -670,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; @@ -697,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;