X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=2d0d60cc718263746866659b727ba5ed7b28feeb;hb=e99ef5538ce4756a0b3e77d347076ef5c041c629;hp=164b42eb78dca58d20ec61aebb3919c23cfce7c0;hpb=87f312d59cabca5080b481a20314601ea476c4be;p=ctsim.git diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 164b42e..2d0d60c 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.34 2003/07/04 21:39:40 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 @@ -270,12 +270,11 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType detInc = proj.detInc(); nDet = proj.nDet(); iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection - rotScale = proj.rotInc(); if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) - rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations + rotScale = PI / proj.nView(); // 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 + rotScale = (2 * PI) / proj.nView(); // scale by number of 2PI rotations else sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry()); @@ -527,37 +526,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])); + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]); } else if (interpType == Backprojector::INTERP_CUBIC) { - double p = iDetCenter + 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 - if (interpType == Backprojector::INTERP_CUBIC) + if (interpType == Backprojector::INTERP_LINEAR) + delete deltaFilteredProj; + else if (interpType == Backprojector::INTERP_CUBIC) delete pCubicInterp; } @@ -572,17 +582,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 = (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) - 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; @@ -598,31 +612,38 @@ 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]; + else + pImCol++; + } // 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]; + else + pImCol++; } // 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]); + const long iDetPos = curDetPos >> scaleShift; + if (iDetPos >= 0 && iDetPos <= iLastDet) { + const long detRemainder = curDetPos & scaleBitmask; + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + } else + pImCol++; } // 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) / 65536); + *pImCol++ += pCubicInterp->interpolate (static_cast(curDetPos) / scale); } } // end Cubic } // end for ix