X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=ebfe2a4358326b16aafafb2b30ef8fe84b3240f9;hp=8343314faca95ad774e4a6e3ba87e8cbd5f58a25;hb=9776c9a12ba53419d34563a5ec57c90e3d6798f4;hpb=23f5654dacb1952c15bda92c2606fae3a55e48ad diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 8343314..ebfe2a4 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -6,9 +6,9 @@ ** Date Started: June 2000 ** ** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg +** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.23 2001/01/04 21:28:41 kevin Exp $ +** $Id: backprojectors.cpp,v 1.27 2001/02/22 18:22:40 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 @@ -30,28 +30,22 @@ const int Backprojector::BPROJ_INVALID = -1; const int Backprojector::BPROJ_TRIG = 0; const int Backprojector::BPROJ_TABLE = 1; const int Backprojector::BPROJ_DIFF = 2; -const int Backprojector::BPROJ_DIFF2 = 3; -const int Backprojector::BPROJ_IDIFF2 = 4; -const int Backprojector::BPROJ_IDIFF3 = 5; +const int Backprojector::BPROJ_IDIFF = 3; -const char* Backprojector::s_aszBackprojectName[] = +const char* const Backprojector::s_aszBackprojectName[] = { {"trig"}, {"table"}, {"diff"}, - {"diff2"}, - {"idiff2"}, - {"idiff3"}, + {"idiff"}, }; -const char* Backprojector::s_aszBackprojectTitle[] = +const char* const Backprojector::s_aszBackprojectTitle[] = { {"Direct Trigometric"}, {"Trigometric Table"}, {"Difference Iteration"}, - {"Difference Iteration Optimized"}, - {"Integer Difference Iteration Optimized"}, - {"Integer Difference Iteration Highly-Optimized"}, + {"Integer Difference Iteration"}, }; const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*); @@ -59,19 +53,23 @@ const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / si const int Backprojector::INTERP_INVALID = -1; const int Backprojector::INTERP_NEAREST = 0; const int Backprojector::INTERP_LINEAR = 1; -const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 2; +const int Backprojector::INTERP_CUBIC = 2; +const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3; #if HAVE_BSPLINE_INTERP -const int Backprojector::INTERP_BSPLINE = 3; -const int Backprojector::INTERP_1BSPLINE = 4; -const int Backprojector::INTERP_2BSPLINE = 5; -const int Backprojector::INTERP_3BSPLINE = 6; +const int Backprojector::INTERP_BSPLINE = 4; +const int Backprojector::INTERP_1BSPLINE = 5; +const int Backprojector::INTERP_2BSPLINE = 6; +const int Backprojector::INTERP_3BSPLINE = 7; #endif -const char* Backprojector::s_aszInterpName[] = +const char* const Backprojector::s_aszInterpName[] = { {"nearest"}, {"linear"}, + {"cubic"}, +#if HAVE_FREQ_PREINTERP {"freq_preinterpolationj"}, +#endif #if HAVE_BSPLINE_INTERP {"bspline"}, {"1bspline"}, @@ -80,11 +78,14 @@ const char* Backprojector::s_aszInterpName[] = #endif }; -const char* Backprojector::s_aszInterpTitle[] = +const char* const Backprojector::s_aszInterpTitle[] = { {"Nearest"}, {"Linear"}, + {"Cubic"}, +#if HAVE_FREQ_PREINTERP {"Frequency Preinterpolation"}, +#endif #if HAVE_BSPLINE_INTERP {"B-Spline"}, {"B-Spline 1st Order"}, @@ -159,12 +160,8 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); else if (m_idBackproject == BPROJ_DIFF) m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_DIFF2) - m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF2) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF3) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_IDIFF) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor)); } else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -330,6 +327,10 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double { 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 y = yMin + yInc / 2; @@ -341,22 +342,25 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest (L / detInc); // calc'd index in the filter raysum array - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos); - else + 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 pFloor = floor (p); int iDetPos = iDetCenter + static_cast(pFloor); double frac = p - pFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos); - else + 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 + if (p >= 0 && p < nDet) + v[ix][iy] += rotScale * pCubicInterp->interpolate (p); } } } + + if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; } @@ -393,6 +397,10 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl { 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]; @@ -402,22 +410,25 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest(L / detInc); // calc index in the filtered raysum vector - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); - else + if (iDetPos >= 0 && iDetPos < nDet) pImCol[iy] += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { 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) - errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); - else + 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) + pImCol[iy] += pCubicInterp->interpolate (p); } } // end for y } // end for x + + if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; } @@ -445,153 +456,60 @@ BackprojectDiff::~BackprojectDiff() ScaleImageByRotIncrement(); } -void -BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle) -{ - double theta = view_angle; // add half PI to view angle to get perpendicular theta angle - double det_dx = xInc * cos (theta); - double det_dy = yInc * sin (theta); - double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image - - for (int ix = 0; ix < nx; ix++, lColStart += det_dx) { - double curDetPos = lColStart; - ImageFileColumn pImCol = v[ix]; - - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { -#ifdef DEBUG - printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos); -#endif - if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(curDetPos / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - pImCol[iy] += filteredProj[iDetPos]; - } else if (interpType == Backprojector::INTERP_LINEAR) { - double detPos = curDetPos / detInc; // position along detector - double detPosFloor = floor (detPos); - int iDetPos = iDetCenter + static_cast(detPosFloor); - double frac = detPos - detPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); - } - } // end for y - } // end for x -} - - -// CLASS IDENTICATION -// BackprojectDiff2 -// -// PURPOSE -// Optimized version of BackprojectDiff void -BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle) +BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - // Distance betw. 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; -#ifdef DEBUG - printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale); -#endif + CubicPolyInterpolator* pCubicInterp = NULL; + if (interpType == Backprojector::INTERP_CUBIC) + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); + 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) { -#ifdef DEBUG - printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(curDetPos)]); -#endif if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest (curDetPos); // calc index in the filtered raysum vector - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else + 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) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else + if (iDetPos > 0 && iDetPos < nDet - 1) *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); + } else if (interpType = Backprojector::INTERP_CUBIC) { + double p = iDetCenter + curDetPos; // position along detector + if (p >= 0 && p < nDet) + *pImCol++ += pCubicInterp->interpolate (p); } } // end for y } // end for x -} -// CLASS IDENTICATION -// BackprojectIntDiff2 -// -// PURPOSE -// Integer version of BackprojectDiff2 - -void -BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle) -{ - double theta = view_angle; - - static const kint32 scale = 1 << 16; - static const double dScale = scale; - static const kint32 halfScale = scale / 2; - - 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 * scale); - - for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - kint32 curDetPos = detPosColStart; - ImageFileColumn pImCol = v[ix]; - - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - if (interpType == Backprojector::INTERP_NEAREST) { - int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale)); - int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos]; - } else if (interpType == Backprojector::INTERP_LINEAR) { - kint32 detPosFloor = curDetPos / scale; - kint32 detPosRemainder = curDetPos % scale; - if (detPosRemainder < 0) { - detPosFloor--; - detPosRemainder += scale; - } - int iDetPos = iDetCenter + detPosFloor; - double frac = detPosRemainder / dScale; - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); - } - } // end for y - } // end for x + if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; } + // CLASS IDENTICATION -// BackprojectIntDiff3 +// BackprojectIntDiff // // PURPOSE -// Highly optimized version of BackprojectIntDiff2 +// Highly optimized and integer version of BackprojectDiff void -BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle) +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 static const int scaleShift = 16; @@ -606,13 +524,17 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do // calculate L for first point in image (0, 0) kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); - // precalculate scaled difference for linear interpolation - double* deltaFilteredProj = new double [nDet]; + double* deltaFilteredProj = NULL; + CubicPolyInterpolator* pCubicInterp = 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]) * dInvScale; + deltaFilteredProj[nDet - 1] = 0; // last detector + } else if (interpType == Backprojector::INTERP_CUBIC) { + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); } - deltaFilteredProj[nDet - 1] = 0; // last detector int iLastDet = nDet - 1; for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { @@ -638,10 +560,17 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do if (iDetPos >= 0 && iDetPos <= iLastDet) *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy - } //end linear + } 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 } // end for ix - delete deltaFilteredProj; + if (interpType == Backprojector::INTERP_LINEAR) + delete deltaFilteredProj; + else if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; } @@ -650,6 +579,10 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const { 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]; @@ -659,27 +592,29 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double rsin_t = r[ix][iy] * sin (dAngleDiff); double dFLPlusSin = m_dFocalLength + rsin_t; 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(gamma / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos - ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos); - } else + 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 dPos = gamma / 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) { - ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos); - } else + 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 + if (d >= 0 && d < nDet) + pImCol[iy] += pCubicInterp->interpolate (d) / dL2; } } // end for y } // end for x + + if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; } void @@ -687,6 +622,10 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const { 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]; @@ -698,27 +637,31 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double dU = (m_dFocalLength + rsin_t) / m_dFocalLength; double dDetPos = rcos_t / dU; // double to scale for imaginary detector that passes through origin - // of phantom, see Kak-Slaney Figure 3.22 + // of phantom, see Kak-Slaney Figure 3.22. This assumes that the detector is also + // located focal-length away from the origin. dDetPos *= 2; - + double dPos = dDetPos / detInc; // position along detector array + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(dDetPos / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - ; /// errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos); - else + 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 dPos = dDetPos / 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) - ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos); - else - pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU); + 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 + if (d >= 0 && d < nDet) + pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU); } } // end for y } // end for x + + if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; }