X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;ds=sidebyside;f=libctsim%2Fbackprojectors.cpp;h=ed388addbfb292718db205f1e5e13b2561a25c5f;hb=30e455abcd8cac05ce7afe43216ec9e26342e1cf;hp=d3cb2035dc125af3415006a2d204cc61b5cc7907;hpb=29ceda2bc7505f0496e1cb5d1e45b489cafa7891;p=ctsim.git diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index d3cb203..ed388ad 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.6 2000/07/09 08:16:17 kevin Exp $ +** $Id: backprojectors.cpp,v 1.7 2000/07/11 10:32:44 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 @@ -26,18 +26,18 @@ #include "ct.h" -Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName) +Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) { m_fail = false; m_pBackprojectImplem = NULL; - initBackprojector (proj, im, backprojName, interpName); + initBackprojector (proj, im, backprojName, interpName, interpFactor); } void Backprojector::BackprojectView (const double* const viewData, const double viewAngle) { - if (m_pBackprojectImplem) + if (m_pBackprojectImplem != NULL) m_pBackprojectImplem->BackprojectView (viewData, viewAngle); } @@ -54,7 +54,7 @@ Backprojector::~Backprojector (void) // and initializes the backprojector bool -Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName) +Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) { m_nameBackproject = backprojName; m_nameInterpolation = interpName; @@ -78,17 +78,17 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const } if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation)); + 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)); + 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)); + 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)); + 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)); + 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)); + m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor)); else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -152,8 +152,8 @@ Backprojector::convertInterpolationNameToID (const char* const interpName) interpID = INTERP_NEAREST; else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0) interpID = INTERP_LINEAR; - else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATE_STR) == 0) - interpID = INTERP_FREQ_PREINTERPOLATE; + else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATION_STR) == 0) + interpID = INTERP_FREQ_PREINTERPOLATION; #if HAVE_BSPLINE_INTERP else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0) interpID = INTERP_BSPLINE; @@ -182,8 +182,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID) return (INTERP_NEAREST_STR); else if (interpID == INTERP_LINEAR) return (INTERP_LINEAR_STR); - else if (interpID == INTERP_FREQ_PREINTERPOLATE) - return (INTERP_FREQ_PREINTERPOLATE_STR); + else if (interpID == INTERP_FREQ_PREINTERPOLATION) + return (INTERP_FREQ_PREINTERPOLATION_STR); #if HAVE_BSPLINE_INTERP else if (interpID == INTERP_BSPLINE) return (INTERP_BSPLINE_STR); @@ -199,8 +199,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID) // PURPOSE // Pure virtual base class for all backprojectors. -Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType) - : proj(proj), im(im), interpType(interpType) +Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType, const int interpFactor) + : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor) { detInc = proj.detInc(); nDet = proj.nDet(); @@ -219,9 +219,6 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Backproj xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; - - if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR) - sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType); } Backproject::~Backproject (void) @@ -298,8 +295,8 @@ 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, Backprojector::InterpolationID interpType) - : Backproject::Backproject (proj, im, interpType) +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : Backproject::Backproject (proj, im, interpType, interpFactor) { arrayR.initSetSize (nx, ny); arrayPhi.initSetSize (nx, ny); @@ -360,8 +357,8 @@ 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, Backprojector::InterpolationID interpType) - : Backproject::Backproject (proj, im, interpType) +BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : Backproject::Backproject (proj, im, interpType, interpFactor) { // calculate center of first pixel v[0][0] double x = xMin + xInc / 2; @@ -458,7 +455,7 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl if (iDetPos < 0 || iDetPos >= nDet - 1) errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else - *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); } } // end for y } // end for x @@ -539,29 +536,35 @@ 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 [nDet - 1]; + if (interpType == Backprojector::INTERP_LINEAR) { + for (int i = 0; i < nDet - 1; i++) + deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale; + } + for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; if (interpType == Backprojector::INTERP_NEAREST) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - int iDetPos = (curDetPos + halfScale) >> 16; + const int iDetPos = (curDetPos + halfScale) >> 16; + assert(iDetPos >= 0 && iDetPos < nDet); + *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; assert(iDetPos >= 0 && iDetPos < nDet); *pImCol++ += filteredProj[iDetPos]; } // end for iy } else if (interpType == Backprojector::INTERP_LINEAR) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - kint32 iDetPos = curDetPos >> scaleShift; - kint32 detRemainder = curDetPos & scaleBitmask; -#if 0 - double frac = detRemainder * (1. / 65536.); + const kint32 iDetPos = curDetPos >> scaleShift; + const kint32 detRemainder = curDetPos & scaleBitmask; assert(iDetPos >= 0 && iDetPos < nDet - 1); - *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]); -#else - assert(iDetPos >= 0 && iDetPos < nDet - 1); - const double* const detPointer = &filteredProj[iDetPos]; - *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale; -#endif + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy } //end linear } // end for ix