X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=5ea532721ccbe6ca25f9409761dd513b180ea3f6;hp=2d0d60cc718263746866659b727ba5ed7b28feeb;hb=1a050c98763fbbc0662731b0b76953acede6f5d7;hpb=c8b19dfaffba9f06d8b6c40cb1bb83a8964867f7 diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 2d0d60c..5ea5327 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -1,7 +1,7 @@ /***************************************************************************** ** FILE IDENTIFICATION ** -** Name: backprojectors.cpp Classes for backprojection +** Name: backprojectors.cpp Classes for backprojection ** Programmer: Kevin Rosenberg ** Date Started: June 2000 ** @@ -33,7 +33,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 +41,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 +63,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 +79,7 @@ const char* const Backprojector::s_aszInterpName[] = #endif }; -const char* const Backprojector::s_aszInterpTitle[] = +const char* const Backprojector::s_aszInterpTitle[] = { "Nearest", "Linear", @@ -99,23 +99,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 +131,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 +153,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 +177,7 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; return false; } - + return true; } @@ -186,13 +186,13 @@ int Backprojector::convertBackprojectNameToID (const char* const backprojName) { int backprojID = BPROJ_INVALID; - + for (int i = 0; i < s_iBackprojectCount; i++) if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { backprojID = i; break; } - + return (backprojID); } @@ -200,10 +200,10 @@ const char* Backprojector::convertBackprojectIDToName (int bprojID) { static const char *bprojName = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) return (s_aszBackprojectName[bprojID]); - + return (bprojName); } @@ -211,10 +211,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,13 +223,13 @@ int Backprojector::convertInterpNameToID (const char* const interpName) { int interpID = INTERP_INVALID; - + for (int i = 0; i < s_iInterpCount; i++) if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { interpID = i; break; } - + return (interpID); } @@ -237,10 +237,10 @@ 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 +248,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,31 +263,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 - + iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection + if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) rotScale = PI / proj.nView(); // scale by number of PI rotations else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) 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 +310,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); @@ -352,7 +352,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 } @@ -368,33 +368,33 @@ 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 + + double x = xMin + xInc / 2; // Rectang coords of center of pixel for (int ix = 0; ix < nx; x += xInc, ix++) { 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); } @@ -403,7 +403,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; -} +} // CLASS IDENTICATION @@ -412,7 +412,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) { @@ -420,8 +420,8 @@ 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 + + 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++) { @@ -447,36 +447,36 @@ 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); - + 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 - + 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 + 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 + 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; @@ -490,16 +490,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(); } @@ -520,16 +520,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) + + // calculate detPosition for first point in image (ix=0, iy=0) double detPosColStart = 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]; @@ -539,31 +539,31 @@ 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; ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = nearest (curDetPos); // calc index in the filtered raysum vector - + int iDetPos = nearest (curDetPos); // 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 = static_cast(detPosFloor); - double frac = curDetPos - detPosFloor; // fraction distance from det + double frac = curDetPos - detPosFloor; // fraction distance from det if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]); + *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]); } else if (interpType == Backprojector::INTERP_CUBIC) { - double p = curDetPos; // position along detector + double p = curDetPos; // 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; @@ -591,14 +591,14 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou 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) + + // calculate L for first point in image (0, 0) long detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); - - double* deltaFilteredProj = NULL; + + double* deltaFilteredProj = NULL; CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_LINEAR) { // precalculate scaled difference for linear interpolation @@ -609,12 +609,12 @@ 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; 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; @@ -623,7 +623,7 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou else pImCol++; - } // end for iy + } // 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; @@ -631,23 +631,23 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou *pImCol++ += filteredProj[iDetPos]; else pImCol++; - } // end for iy + } // end for iy } else if (interpType == Backprojector::INTERP_LINEAR) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { const long iDetPos = curDetPos >> scaleShift; if (iDetPos >= 0 && iDetPos <= iLastDet) { - const long detRemainder = curDetPos & scaleBitmask; - *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); - } else + const long detRemainder = curDetPos & scaleBitmask; + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + } else pImCol++; - } // end for iy + } // 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 } // end for ix - + if (interpType == Backprojector::INTERP_LINEAR) delete deltaFilteredProj; else if (interpType == Backprojector::INTERP_CUBIC) @@ -659,15 +659,15 @@ 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); - + 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); @@ -675,24 +675,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; @@ -702,43 +702,43 @@ 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); - + 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 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) + int iDetPos = iDetCenter + nearest(dPos); // calc index in the filtered raysum vector + if (iDetPos >= 0 && iDetPos < nDet) pImCol[iy] += (filteredProj[iDetPos] / (dU * dU)); } 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); } 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); } - } // end for y - } // end for x + } // end for y + } // end for x if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp;