** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.12 2000/08/25 15:59:13 kevin Exp $
+** $Id: backprojectors.cpp,v 1.13 2000/08/31 08:38:58 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
return false;
}
- if (m_idBackproject == BPROJ_TRIG)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_TABLE)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_DIFF)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_DIFF2)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_IDIFF2)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_IDIFF3)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
- else {
- m_fail = true;
- m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
- return false;
+ if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
+ else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
+ else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
+ if (m_idBackproject == BPROJ_TRIG)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_TABLE)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_DIFF)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_DIFF2)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_IDIFF2)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_IDIFF3)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
+ } else {
+ m_fail = true;
+ m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
+ return false;
}
return true;
{
detInc = proj.detInc();
nDet = proj.nDet();
- iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
+ iDetCenter = nDet / 2; // index refering to L=0 projection
rotInc = proj.rotInc();
v = im.getArray();
xInc = (xMax - xMin) / nx; // size of cells
yInc = (yMax - yMin) / ny;
+
+ m_dFocalLength = proj.focalLength();
}
Backproject::~Backproject ()
void
BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
- int ix, iy;
- double x, y; // Rectang coords of center of pixel
+ double theta = view_angle;
- for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
- for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
+ 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;
+ for (int iy = 0; iy < ny; y += yInc, iy++) {
double r = sqrt (x * x + y * y); // distance of cell from center
double phi = atan2 (y, x); // angle of cell from center
double L = r * cos (theta - phi); // position on detector
v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
}
}
+ }
}
void
BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
+ double theta = view_angle;
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
void
BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
- double det_dx = xInc * sin (theta);
- double det_dy = yInc * cos (theta);
+ double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
+ double det_dx = xInc * cos (theta);
+ double det_dy = yInc * sin (theta);
double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
void
BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+ double theta = view_angle;
// Distance betw. detectors for an angle given in units of detectors
- double det_dx = xInc * sin (theta) / detInc;
- double det_dy = yInc * cos (theta) / detInc;
+ double det_dx = xInc * cos (theta) / detInc;
+ 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;
void
BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+ double theta = view_angle;
static const kint32 scale = 1 << 16;
static const double dScale = scale;
static const kint32 halfScale = scale / 2;
- const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+ const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (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 * scale);
void
BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+ double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
static const int scaleShift = 16;
static const kint32 scale = (1 << scaleShift);
static const kint32 scaleBitmask = scale - 1;
static const kint32 halfScale = scale / 2;
static const double dInvScale = 1. / scale;
- const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+ const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (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);
} //end linear
} // end for ix
}
+
+
+void
+BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+ double dFLPlusSin = m_dFocalLength + rsin_t;
+ double gamma = atan (rcos_t / dFLPlusSin);
+ 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
+ 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
+ pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / dL2;
+ }
+ } // end for y
+ } // end for x
+}
+
+void
+BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+
+ 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
+ dDetPos *= 2;
+
+ 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
+ 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] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / (dU * dU);
+ }
+ } // end for y
+ } // end for x
+}
+
projects / ctsim.git / blobdiff