** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.24 2001/01/27 21:02:20 kevin Exp $
+** $Id: backprojectors.cpp,v 1.25 2001/02/09 01:54:20 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
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[] =
{
{"nearest"},
{"linear"},
+ {"cubic"},
{"freq_preinterpolationj"},
#if HAVE_BSPLINE_INTERP
{"bspline"},
{
{"Nearest"},
{"Linear"},
+ {"Cubic"},
{"Frequency Preinterpolation"},
#if HAVE_BSPLINE_INTERP
{"B-Spline"},
{
double theta = view_angle;
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
double x = xMin + xInc / 2; // Rectang coords of center of pixel
for (int ix = 0; ix < nx; x += xInc, ix++) {
double y = yMin + yInc / 2;
if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (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
+ 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<int>(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
+ 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;
}
{
double theta = view_angle;
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(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
+ 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<int>(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
+ 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
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
double det_dy = yInc * sin (theta);
double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
double curDetPos = lColStart;
ImageFileColumn pImCol = v[ix];
if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(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
+ if (iDetPos >= 0 && iDetPos < nDet)
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<int>(detPosFloor);
double frac = detPos - detPosFloor; // fraction distance from det
- if (iDetPos < 0 || iDetPos >= nDet - 1)
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
+ 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 + (curDetPos / detInc); // position along detector
+ if (p >= 0 && p < nDet)
+ pImCol[iy] += pCubicInterp->interpolate (p);
}
} // end for y
} // end for x
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
{
double theta = view_angle;
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
// Distance betw. detectors for an angle given in units of detectors
double det_dx = xInc * cos (theta) / detInc;
double det_dy = yInc * sin (theta) / detInc;
if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (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
+ 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);
double frac = curDetPos - detPosFloor; // fraction distance from det
- if (iDetPos < 0 || iDetPos >= nDet - 1)
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
+ 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 (p >= 0 && p < nDet)
+ *pImCol++ += pCubicInterp->interpolate (p);
}
} // end for y
} // end for x
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
// CLASS IDENTICATION
{
double theta = view_angle;
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
static const kint32 scale = 1 << 16;
static const double dScale = scale;
static const kint32 halfScale = scale / 2;
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
*pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ } else if (interpType = Backprojector::INTERP_CUBIC) {
+ double p = iDetCenter + (static_cast<double>(curDetPos) / scale); // position along detector
+ if (p >= 0 && p < nDet)
+ *pImCol++ += pCubicInterp->interpolate (p);
}
} // end for y
} // end for x
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
// CLASS IDENTICATION
// calculate L for first point in image (0, 0)
kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
- // precalculate scaled difference for linear interpolation
- double* deltaFilteredProj = new double [nDet];
+ double* deltaFilteredProj = NULL;
+ CubicInterpolator* 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 CubicInterpolator (filteredProj, nDet);
}
- deltaFilteredProj[nDet - 1] = 0; // last detector
int iLastDet = nDet - 1;
for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
if (iDetPos >= 0 && iDetPos <= iLastDet)
*pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
} // end for iy
- } //end linear
+ } else if (interpType = Backprojector::INTERP_CUBIC) {
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
+ }
+ } // end cubic
} // end for ix
- delete deltaFilteredProj;
+ if (interpType == Backprojector::INTERP_LINEAR)
+ delete deltaFilteredProj;
+ else if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
{
double beta = view_angle;
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
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<int>(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
+ int iDetPos = iDetCenter + nearest<int>(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<int>(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
+ 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
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
void
{
double beta = view_angle;
+ CubicInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
double dDetPos = rcos_t / dU;
// double to scale for imaginary detector that passes through origin
- // of phantom, see Kak-Slaney Figure 3.22
+ // of phantom, see Kak-Slaney Figure 3.22. This assumes that the detector is also
+ // located focal-length away from the origin.
dDetPos *= 2;
-
+ double dPos = dDetPos / detInc; // position along detector array
+
if (interpType == Backprojector::INTERP_NEAREST) {
- int iDetPos = iDetCenter + nearest<int>(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
+ int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
+ if (iDetPos >= 0 && iDetPos < nDet)
pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
} else if (interpType == Backprojector::INTERP_LINEAR) {
- double dPos = dDetPos / detInc; // position along detector
double dPosFloor = floor (dPos);
int iDetPos = iDetCenter + static_cast<int>(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);
+ if (iDetPos >= 0 && iDetPos < nDet - 1)
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
+ / (dU * dU);
+ } else if (interpType == Backprojector::INTERP_CUBIC) {
+ double d = iDetCenter + dPos; // position along detector
+ if (d >= 0 && d < nDet)
+ pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
}
} // end for y
} // end for x
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}