** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.13 2000/08/31 08:38:58 kevin Exp $
+** $Id: backprojectors.cpp,v 1.18 2000/12/06 15:17:51 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
// PURPOSE
// Pure virtual base class for all backprojectors.
-Backproject::Backproject (const Projections& proj, ImageFile& im, const int interpType, const int interpFactor)
+Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
: proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
{
detInc = proj.detInc();
nDet = proj.nDet();
- iDetCenter = nDet / 2; // index refering to L=0 projection
- rotInc = proj.rotInc();
+ iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
+ rotScale = proj.rotInc();
+ rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
v = im.getArray();
nx = im.nx();
{
for (int ix = 0; ix < nx; ix++)
for (int iy = 0; iy < ny; iy++)
- v[ix][iy] *= rotInc;
+ v[ix][iy] *= rotScale;
}
void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
}
void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
-{
+{\r
+#if 1
ostringstream os;
os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
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());\r
+#endif
}
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] += rotInc * filteredProj[iDetPos];
+ v[ix][iy] += rotScale * filteredProj[iDetPos];
} else if (interpType == Backprojector::INTERP_LINEAR) {
double p = L / detInc; // position along detector
double pFloor = floor (p);
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] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
}
}
}
// Precalculates trigometric function value for each point in image for backprojection.
BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
- : Backproject::Backproject (proj, im, interpType, interpFactor)
+ : Backproject (proj, im, interpType, interpFactor)
{
- arrayR.initSetSize (nx, ny);
- arrayPhi.initSetSize (nx, ny);
+ arrayR.initSetSize (im.nx(), im.ny());
+ arrayPhi.initSetSize (im.nx(), im.ny());
r = arrayR.getArray();
phi = arrayPhi.getArray();
// Iterates in x & y direction by adding difference in L position
BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
- : Backproject::Backproject (proj, im, interpType, interpFactor)
+ : Backproject (proj, im, interpType, interpFactor)
{
// calculate center of first pixel v[0][0]
double x = xMin + xInc / 2;
double detPosColStart = start_r * cos (theta - start_phi) / detInc;
#ifdef DEBUG
- printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
+ 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;
kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
// precalculate scaled difference for linear interpolation
- double deltaFilteredProj [nDet];
+ double* deltaFilteredProj = new double [nDet];
if (interpType == Backprojector::INTERP_LINEAR) {
for (int i = 0; i < nDet - 1; i++)
deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
*pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
} // end for iy
} //end linear
- } // end for ix
+ } // end for ix\r
+\r
+ delete deltaFilteredProj;
}
if (iDetPos < 0 || iDetPos >= nDet - 1) {
; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
} else
- pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / dL2;
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
}
} // end for y
} // end for x
if (iDetPos < 0 || iDetPos >= nDet - 1)
; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
else
- pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / (dU * dU);
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU);
}
} // end for y
} // end for x