** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.32 2001/03/18 18:08:25 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
******************************************************************************/
#include "ct.h"
+#include "interpolator.h"
const int Backprojector::BPROJ_INVALID = -1;
const int Backprojector::BPROJ_TRIG = 0;
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*);
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
};
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());
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<int> (curDetPos); // calc index in the filtered raysum vector
+ int iDetPos = nearest<int> (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<int>(detPosFloor);
+ int iDetPos = static_cast<int>(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;
}
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<kint32> (xInc * cos (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
+ const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
+ const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
// calculate L for first point in image (0, 0)
- kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
+ long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
double* deltaFilteredProj = NULL;
CubicPolyInterpolator* pCubicInterp = NULL;
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<double>(curDetPos) / 65536);
+ *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
}
} // end Cubic
} // end for ix