X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=2b9eba45cdc4ec7ce4b300030c19a4fafaf5af1e;hp=7926b0d27301b9946ea575369bd41fca88a80f6f;hb=352a0691b9bd67f6c93ea822b353be3c101f4adb;hpb=ebb615bb35cf344a765aa65414982904e59f961b diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 7926b0d..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-2000 Kevin Rosenberg -** -** $Id: backprojectors.cpp,v 1.20 2000/12/18 09:31:26 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,33 +23,28 @@ ******************************************************************************/ #include "ct.h" +#include "interpolator.h" 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 char* Backprojector::s_aszBackprojectName[] = -{ - {"trig"}, - {"table"}, - {"diff"}, - {"diff2"}, - {"idiff2"}, - {"idiff3"}, +const int Backprojector::BPROJ_IDIFF = 3; + +const char* const Backprojector::s_aszBackprojectName[] = +{ + "trig", + "table", + "diff", + "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"}, + "Direct Trigometric", + "Trigometric Table", + "Difference Iteration", + "Integer Difference Iteration", }; const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*); @@ -59,37 +52,44 @@ 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"}, - {"freq_preinterpolationj"}, + "nearest", + "linear", + "cubic", +#if HAVE_FREQ_PREINTERP + "freq_preinterpolationj", +#endif #if HAVE_BSPLINE_INTERP - {"bspline"}, - {"1bspline"}, - {"2bspline"}, - {"3bspline"}, + "bspline", + "1bspline", + "2bspline", + "3bspline", #endif }; -const char* Backprojector::s_aszInterpTitle[] = +const char* const Backprojector::s_aszInterpTitle[] = { - {"Nearest"}, - {"Linear"}, - {"Frequency Preinterpolationj"}, + "Nearest", + "Linear", + "Cubic", +#if HAVE_FREQ_PREINTERP + "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 }; @@ -97,21 +97,29 @@ 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 +void Backprojector::BackprojectView (const double* const viewData, const double viewAngle) { if (m_pBackprojectImplem != NULL) m_pBackprojectImplem->BackprojectView (viewData, viewAngle); } +void +Backprojector::PostProcessing() +{ + if (m_pBackprojectImplem != NULL) + m_pBackprojectImplem->PostProcessing(); +} + Backprojector::~Backprojector () { delete m_pBackprojectImplem; @@ -121,11 +129,12 @@ 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, 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; @@ -149,26 +158,22 @@ 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)); - else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) - m_pBackprojectImplem = static_cast(new BackprojectEquiangular(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, pROI)); else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) { - if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_TABLE) - 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)); + if (m_idBackproject == BPROJ_TRIG) + 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, pROI)); + else if (m_idBackproject == BPROJ_DIFF) + 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, pROI)); } else { - m_fail = true; - m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; - return false; + m_fail = true; + m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; + return false; } return true; @@ -180,12 +185,12 @@ 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; - } - + for (int i = 0; i < s_iBackprojectCount; i++) { + if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { + backprojID = i; + break; + } + } return (backprojID); } @@ -195,7 +200,7 @@ Backprojector::convertBackprojectIDToName (int bprojID) static const char *bprojName = ""; if (bprojID >= 0 && bprojID < s_iBackprojectCount) - return (s_aszBackprojectName[bprojID]); + return (s_aszBackprojectName[bprojID]); return (bprojName); } @@ -206,7 +211,7 @@ Backprojector::convertBackprojectIDToTitle (const int bprojID) static const char *bprojTitle = ""; if (bprojID >= 0 && bprojID < s_iBackprojectCount) - return (s_aszBackprojectTitle[bprojID]); + return (s_aszBackprojectTitle[bprojID]); return (bprojTitle); } @@ -217,12 +222,12 @@ 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; - } - + for (int i = 0; i < s_iInterpCount; i++) { + if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { + interpID = i; + break; + } + } return (interpID); } @@ -232,7 +237,7 @@ Backprojector::convertInterpIDToName (const int interpID) static const char *interpName = ""; if (interpID >= 0 && interpID < s_iInterpCount) - return (s_aszInterpName[interpID]); + return (s_aszInterpName[interpID]); return (interpName); } @@ -243,7 +248,7 @@ Backprojector::convertInterpIDToTitle (const int interpID) static const char *interpTitle = ""; if (interpID >= 0 && interpID < s_iInterpCount) - return (s_aszInterpTitle[interpID]); + return (s_aszInterpTitle[interpID]); return (interpTitle); } @@ -256,20 +261,20 @@ 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) - : proj(proj), im(im), interpType(interpType), m_interpFactor(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 - 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 - else - sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry()); + 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(); @@ -281,15 +286,47 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType yMin = -proj.phmLen() / 2; yMax = yMin + proj.phmLen(); - xInc = (xMax - xMin) / nx; // size of cells + 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(); } Backproject::~Backproject () {} +void +Backproject::PostProcessing() +{ + m_bPostProcessingDone = true; +} + void Backproject::ScaleImageByRotIncrement () { @@ -300,12 +337,12 @@ Backproject::ScaleImageByRotIncrement () void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos) { - sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi); - errorIndexOutsideDetector (ix, iy, theta, L, iDetPos); + sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi); + errorIndexOutsideDetector (ix, iy, theta, L, iDetPos); } void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos) -{ +{ #if 1 std::ostringstream os; os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n"; @@ -314,7 +351,7 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n"; os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";; - sys_error (ERR_WARNING, os.str().c_str()); + sys_error (ERR_WARNING, os.str().c_str()); #endif } @@ -330,8 +367,17 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double { double theta = view_angle; - double x = xMin + xInc / 2; // Rectang coords of center of pixel - for (int ix = 0; ix < nx; x += xInc, ix++) { + CubicPolyInterpolator* pCubicInterp = NULL; + if (interpType == Backprojector::INTERP_CUBIC) + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); + + 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 @@ -339,25 +385,28 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double 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 + 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 - v[ix][iy] += rotScale * filteredProj[iDetPos]; + 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 - v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + 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) + 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; +} // CLASS IDENTICATION @@ -366,26 +415,41 @@ 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()); 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 () { - ScaleImageByRotIncrement(); +} + +void +BackprojectTable::PostProcessing() +{ + if (! m_bPostProcessingDone) { + ScaleImageByRotIncrement(); + m_bPostProcessingDone = true; + } } void @@ -393,6 +457,13 @@ 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); + +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; @@ -400,24 +471,30 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl 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) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); - else - pImCol[iy] += filteredProj[iDetPos]; + 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 - pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + 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) { + 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 + } // end for y + } // end for x + + if (interpType == Backprojector::INTERP_CUBIC) + delete pCubicInterp; } @@ -428,10 +505,11 @@ 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] + // 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); @@ -440,208 +518,157 @@ BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int in im.arrayDataClear(); } -BackprojectDiff::~BackprojectDiff() +BackprojectDiff::~BackprojectDiff () { - ScaleImageByRotIncrement(); } void -BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle) +BackprojectDiff::PostProcessing() { - 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 + if (! m_bPostProcessingDone) { + ScaleImageByRotIncrement(); + m_bPostProcessingDone = true; + } } - -// 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 - for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - double curDetPos = detPosColStart; - ImageFileColumn pImCol = v[ix]; + // calculate detPosition for first point in image (ix=0, iy=0) + double detPosColBase = iDetCenter + start_r * cos (theta - start_phi) / detInc; - 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 - *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 - *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); - } - } // 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; + CubicPolyInterpolator* pCubicInterp = NULL; + 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); + } - const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); - const kint32 det_dy = nearest (yInc * sin (theta) / detInc * scale); + int iLastDet = nDet - 1; - // 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; +#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 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]; + 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) { - 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]); + 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 = 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; + else 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 +#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); + 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); - // 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) { - 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) - *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) - *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) - *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); - } // end for iy - } //end linear - } // end for ix - - delete deltaFilteredProj; + 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]; + } 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]; + } 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) + delete pCubicInterp; } @@ -650,36 +677,45 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const { 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); 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 - pImCol[iy] += filteredProj[iDetPos] / dL2; + 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 - pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2; + double dPosFloor = floor (dPos); + int iDetPos = iDetCenter + static_cast(dPosFloor); + 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 + 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; } void @@ -687,6 +723,13 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const { 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]; @@ -696,29 +739,32 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const 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; - // double to scale for imaginary detector that passes through origin - // of phantom, see Kak-Slaney Figure 3.22 - dDetPos *= 2; + // 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(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 - 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 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); + double dPosFloor = floor (dPos); + int iDetPos = iDetCenter + static_cast(dPosFloor); + double frac = dPos - dPosFloor; // fraction distance from det + if (iDetPos >= 0 && iDetPos < nDet - 1) + pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dU2; + } else if (interpType == Backprojector::INTERP_CUBIC) { + double d = iDetCenter + dPos; // position along detector + if (d >= 0 && d < nDet) + 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; }