X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=d38b7189f6b6c878b09ec0549dca5daf242ce24e;hb=e46615378520c3e801074560ebdc36be9558c6dd;hp=bab1dc7b18fc75e349cc961c43165f884545ed87;hpb=c953cbb6ffc2fd50e736230f4e6976a025983cff;p=ctsim.git diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index bab1dc7..d38b718 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.29 2001/03/01 07:30:49 kevin Exp $ +** $Id$ ** ** 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 +25,7 @@ ******************************************************************************/ #include "ct.h" +#include "interpolator.h" const int Backprojector::BPROJ_INVALID = -1; const int Backprojector::BPROJ_TRIG = 0; @@ -34,18 +35,18 @@ const int Backprojector::BPROJ_IDIFF = 3; const char* const Backprojector::s_aszBackprojectName[] = { - {"trig"}, - {"table"}, - {"diff"}, - {"idiff"}, + "trig", + "table", + "diff", + "idiff", }; 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*); @@ -64,33 +65,33 @@ const int Backprojector::INTERP_3BSPLINE = 7; 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[] = { - {"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,12 +99,13 @@ const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const -Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) +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); + initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI); } void @@ -133,7 +135,8 @@ Backprojector::~Backprojector () // and initializes the backprojector bool -Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) +Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, + const char* const interpName, const int interpFactor, const ReconstructionROI* pROI) { m_nameBackproject = backprojName; m_nameInterpolation = interpName; @@ -157,18 +160,18 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const } if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) - m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI)); else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) - m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI)); else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) { if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI)); else if (m_idBackproject == BPROJ_TABLE) - m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI)); else if (m_idBackproject == BPROJ_DIFF) - m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI)); else if (m_idBackproject == BPROJ_IDIFF) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI)); } else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -260,7 +263,8 @@ 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(); @@ -285,9 +289,34 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType yMin = -proj.phmLen() / 2; yMax = yMin + proj.phmLen(); + if (pROI) { + if (pROI->m_dXMin > xMin) + xMin = pROI->m_dXMin; + if (pROI->m_dXMax < xMax) + xMax = pROI->m_dXMax; + if (pROI->m_dYMin > yMin) + yMin = pROI->m_dYMin; + if (pROI->m_dYMax < yMax) + yMax = pROI->m_dYMax; + + if (xMin > xMax) { + double temp = xMin; + xMin = xMax; + xMax = temp; + } + if (yMin > yMax) { + double temp = yMin; + yMin = yMax; + yMax = temp; + } + } + 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(); } @@ -365,7 +394,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double 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) { + } 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); @@ -384,8 +413,9 @@ 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, const int interpFactor) -: Backproject (proj, im, interpType, interpFactor) +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, + const int interpFactor, const ReconstructionROI* pROI) +: Backproject (proj, im, interpType, interpFactor, pROI) { arrayR.initSetSize (im.nx(), im.ny()); arrayPhi.initSetSize (im.nx(), im.ny()); @@ -441,7 +471,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl 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) { + } else if (interpType == Backprojector::INTERP_CUBIC) { double p = iDetCenter + (L / detInc); // position along detector if (p >= 0 && p < nDet) pImCol[iy] += pCubicInterp->interpolate (p); @@ -461,8 +491,9 @@ 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, const int interpFactor) -: Backproject (proj, im, interpType, interpFactor) +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] double x = xMin + xInc / 2; @@ -496,37 +527,48 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double 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; + double detPosColStart = 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); + } + 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 = iDetCenter + 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 = iDetCenter + static_cast(detPosFloor); + int iDetPos = 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])); - } else if (interpType = Backprojector::INTERP_CUBIC) { - double p = iDetCenter + curDetPos; // position along detector + 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) *pImCol++ += pCubicInterp->interpolate (p); } } // 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; } @@ -541,17 +583,21 @@ 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 = (1 << 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); + 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) - kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + long detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); double* deltaFilteredProj = NULL; CubicPolyInterpolator* pCubicInterp = NULL; @@ -567,31 +613,32 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou int iLastDet = nDet - 1; for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - kint32 curDetPos = detPosColStart; + 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) >> 16; + const int iDetPos = (curDetPos + 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; + const int iDetPos = ((curDetPos + 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; + const long iDetPos = curDetPos >> scaleShift; + const long detRemainder = curDetPos & scaleBitmask; if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + } // end for iy - } else if (interpType = Backprojector::INTERP_CUBIC) { + } else if (interpType == Backprojector::INTERP_CUBIC) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - *pImCol++ += pCubicInterp->interpolate (static_cast(curDetPos) / 65536); + *pImCol++ += pCubicInterp->interpolate (static_cast(curDetPos) / scale); } } // end Cubic } // end for ix