X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=726e5b8d30189c14e472f5e4eee802930742964e;hp=4983256ffd0e57846690d9bf7bb5ef9c8a1f544a;hb=7ec2cd66921180a624813dff9f8bac76c6b268cc;hpb=286d655a25df2668bd65ad365676c6ecc94415a1 diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 4983256..726e5b8 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.10 2000/07/22 15:45:33 kevin Exp $ +** $Id: backprojectors.cpp,v 1.21 2001/01/01 10:14:34 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 @@ -84,7 +84,7 @@ const char* Backprojector::s_aszInterpTitle[] = { {"Nearest"}, {"Linear"}, - {"Frequency Preinterpolationj"}, + {"Frequency Preinterpolation"}, #if HAVE_BSPLINE_INTERP {"B-Spline"}, {"B-Spline 1st Order"}, @@ -148,22 +148,27 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const return false; } - if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_TABLE) - m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_DIFF) - m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_DIFF2) - m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF2) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF3) - m_pBackprojectImplem = static_cast(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(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor)); + else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) + m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor)); + else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) { + if (m_idBackproject == BPROJ_TRIG) + m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_TABLE) + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_DIFF) + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_DIFF2) + m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_IDIFF2) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_IDIFF3) + m_pBackprojectImplem = static_cast(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; @@ -251,13 +256,20 @@ Backprojector::convertInterpIDToTitle (const int interpID) // 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 - 1) / 2; // index refering to L=0 projection - rotInc = proj.rotInc(); + rotScale = proj.rotInc(); + + if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) + rotScale /= (proj.nView() * proj.rotInc() / PI); // 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 + else + sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry()); v = im.getArray(); nx = im.nx(); @@ -271,6 +283,8 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const int inte xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; + + m_dFocalLength = proj.focalLength(); } Backproject::~Backproject () @@ -281,23 +295,27 @@ Backproject::ScaleImageByRotIncrement () { 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) { - printf ("r=%f, phi=%f\n", r, phi); + sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi); errorIndexOutsideDetector (ix, iy, theta, L, iDetPos); } void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos) -{ - printf ("ix=%d, iy=%d\n", ix, iy); - printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc); - printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter); - printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc); - printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc); - sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos); +{ +#if 1 + std::ostringstream os; + os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n"; + os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n"; + os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n"; + os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n"; + os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";; + + sys_error (ERR_WARNING, os.str().c_str()); +#endif } @@ -310,12 +328,12 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl 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 @@ -326,7 +344,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double 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); @@ -335,9 +353,10 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double 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]); } } + } } @@ -348,10 +367,10 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // 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(); @@ -372,7 +391,7 @@ BackprojectTable::~BackprojectTable () 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]; @@ -410,7 +429,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // 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; @@ -429,9 +448,9 @@ BackprojectDiff::~BackprojectDiff() 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) { @@ -440,7 +459,7 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { #ifdef DEBUG - printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos); + printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos); #endif if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest(curDetPos / detInc); // calc index in the filtered raysum vector @@ -473,17 +492,17 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double 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; #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; @@ -491,7 +510,7 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { #ifdef DEBUG - printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(L))]); + printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(curDetPos)]); #endif if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest (curDetPos); // calc index in the filtered raysum vector @@ -522,14 +541,14 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl 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 (xInc * sin (theta) / detInc * scale); - const kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); + const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); + const kint32 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 * scale); @@ -543,7 +562,7 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale)); int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else *pImCol++ += filteredProj[iDetPos]; @@ -574,26 +593,28 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do 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 (xInc * sin (theta) / detInc * scale); - const kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); + const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); + const kint32 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); // precalculate scaled difference for linear interpolation - double deltaFilteredProj [nDet - 1]; + 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; } + deltaFilteredProj[nDet - 1] = 0; // last detector + int iLastDet = nDet - 1; for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; @@ -601,22 +622,103 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do if (interpType == Backprojector::INTERP_NEAREST) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { const int iDetPos = (curDetPos + halfScale) >> 16; - assert(iDetPos >= 0 && iDetPos < nDet); - *pImCol++ += filteredProj[iDetPos]; + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos]; } // 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; - assert(iDetPos >= 0 && iDetPos < nDet); + if (iDetPos >= 0 && iDetPos <= iLastDet) *pImCol++ += filteredProj[iDetPos]; } // 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; - assert(iDetPos >= 0 && iDetPos < nDet - 1); - *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy } //end linear - } // end for ix + } // end for ix + + delete deltaFilteredProj; +} + + +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(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(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] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / 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(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(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); + } + } // end for y + } // end for x } +