X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=8343314faca95ad774e4a6e3ba87e8cbd5f58a25;hp=ecdc3fdcf753bbe9c41dbc620b170500c375bf16;hb=23f5654dacb1952c15bda92c2606fae3a55e48ad;hpb=ee0105d74fec9d6bfd236e22e9e1d315e46c568e diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index ecdc3fd..8343314 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.17 2000/12/06 01:46:43 kevin Exp $ +** $Id: backprojectors.cpp,v 1.23 2001/01/04 21:28:41 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"}, @@ -101,7 +101,7 @@ Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char { m_fail = false; m_pBackprojectImplem = NULL; - + initBackprojector (proj, im, backprojName, interpName, interpFactor); } @@ -142,35 +142,35 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_failMessage = "Invalid interpolation name "; m_failMessage += interpName; } - + if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) { m_fail = true; return false; } - + if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) - m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor)); + 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)); + 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)); + 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; + m_fail = true; + m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; + return false; } - + return true; } @@ -179,24 +179,24 @@ int Backprojector::convertBackprojectNameToID (const char* const backprojName) { int backprojID = BPROJ_INVALID; - + for (int i = 0; i < s_iBackprojectCount; i++) - if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { - backprojID = i; - break; - } - - return (backprojID); + if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { + backprojID = i; + break; + } + + return (backprojID); } const char* Backprojector::convertBackprojectIDToName (int bprojID) { static const char *bprojName = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) - return (s_aszBackprojectName[bprojID]); - + return (s_aszBackprojectName[bprojID]); + return (bprojName); } @@ -204,10 +204,10 @@ const char* Backprojector::convertBackprojectIDToTitle (const int bprojID) { static const char *bprojTitle = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) - return (s_aszBackprojectTitle[bprojID]); - + return (s_aszBackprojectTitle[bprojID]); + return (bprojTitle); } @@ -216,24 +216,24 @@ int Backprojector::convertInterpNameToID (const char* const interpName) { int interpID = INTERP_INVALID; - + for (int i = 0; i < s_iInterpCount; i++) - if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { - interpID = i; - break; - } - - return (interpID); + if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { + interpID = i; + break; + } + + return (interpID); } const char* Backprojector::convertInterpIDToName (const int interpID) { static const char *interpName = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) - return (s_aszInterpName[interpID]); - + return (s_aszInterpName[interpID]); + return (interpName); } @@ -241,10 +241,10 @@ const char* Backprojector::convertInterpIDToTitle (const int interpID) { static const char *interpTitle = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) - return (s_aszInterpTitle[interpID]); - + return (s_aszInterpTitle[interpID]); + return (interpTitle); } @@ -257,27 +257,33 @@ Backprojector::convertInterpIDToTitle (const int interpID) // Pure virtual base class for all backprojectors. Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) - : proj(proj), im(im), interpType(interpType), m_interpFactor(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 rotScale = proj.rotInc(); - rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations - + + 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(); ny = im.ny(); im.arrayDataClear(); - + xMin = -proj.phmLen() / 2; // Retangular coords of phantom xMax = xMin + proj.phmLen(); yMin = -proj.phmLen() / 2; yMax = yMin + proj.phmLen(); - + xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; - + m_dFocalLength = proj.focalLength(); } @@ -294,20 +300,22 @@ Backproject::ScaleImageByRotIncrement () void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos) { - sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi); - errorIndexOutsideDetector (ix, iy, theta, L, iDetPos); + 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) { - ostringstream os; +#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 } @@ -321,7 +329,7 @@ void BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + 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; @@ -329,23 +337,23 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double 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 - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest (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 - v[ix][iy] += rotScale * filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest (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 + 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(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 - v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + double p = L / detInc; // position along detector + double pFloor = floor (p); + int iDetPos = iDetCenter + static_cast(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 + v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } } @@ -359,13 +367,13 @@ 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 (proj, im, interpType, interpFactor) +: Backproject (proj, im, interpType, interpFactor) { arrayR.initSetSize (im.nx(), im.ny()); arrayPhi.initSetSize (im.nx(), im.ny()); r = arrayR.getArray(); phi = arrayPhi.getArray(); - + double x, y; // Rectang coords of center of pixel int ix, iy; for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++) @@ -384,29 +392,29 @@ void BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++) { double L = r[ix][iy] * cos (theta - phi[ix][iy]); - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(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 - pImCol[iy] += filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest(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 + 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(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 - pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + double dPos = L / 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, theta, r[ix][iy], phi[ix][iy], L, iDetPos); + else + pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x @@ -421,14 +429,14 @@ 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 (proj, im, interpType, interpFactor) +: Backproject (proj, im, interpType, interpFactor) { // calculate center of first pixel v[0][0] double x = xMin + xInc / 2; double y = yMin + yInc / 2; start_r = sqrt (x * x + y * y); start_phi = atan2 (y, x); - + im.arrayDataClear(); } @@ -444,31 +452,31 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double 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) { double curDetPos = lColStart; ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { #ifdef DEBUG 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 - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - pImCol[iy] += filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest(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 + 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(detPosFloor); - double frac = detPos - detPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + double detPos = curDetPos / detInc; // position along detector + double detPosFloor = floor (detPos); + int iDetPos = iDetCenter + static_cast(detPosFloor); + double frac = detPos - detPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x @@ -485,40 +493,40 @@ void BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + // 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; - + // 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, rotScale=%8.5f\n", start_r, start_phi, rotScale); #endif 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) { #ifdef DEBUG 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 - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest (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 + *pImCol++ += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - double detPosFloor = floor (curDetPos); - int iDetPos = iDetCenter + static_cast(detPosFloor); - double frac = curDetPos - detPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); + double detPosFloor = floor (curDetPos); + int iDetPos = iDetCenter + static_cast(detPosFloor); + double frac = curDetPos - detPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); } } // end for y } // end for x @@ -534,43 +542,43 @@ void BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_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 * 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); - + for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { if (interpType == Backprojector::INTERP_NEAREST) { - 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 index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos]; + 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 index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - kint32 detPosFloor = curDetPos / scale; - kint32 detPosRemainder = curDetPos % scale; - if (detPosRemainder < 0) { - detPosFloor--; - detPosRemainder += scale; - } - int iDetPos = iDetCenter + detPosFloor; - double frac = detPosRemainder / dScale; - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + kint32 detPosFloor = curDetPos / scale; + kint32 detPosRemainder = curDetPos % scale; + if (detPosRemainder < 0) { + detPosFloor--; + detPosRemainder += scale; + } + int iDetPos = iDetCenter + detPosFloor; + double frac = detPosRemainder / dScale; + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x @@ -591,13 +599,13 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do static const kint32 scaleBitmask = scale - 1; static const kint32 halfScale = scale / 2; static const double dInvScale = 1. / 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 = new double [nDet]; if (interpType == Backprojector::INTERP_LINEAR) { @@ -605,34 +613,34 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do 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]; - + if (interpType == Backprojector::INTERP_NEAREST) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const int iDetPos = (curDetPos + halfScale) >> 16; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos]; + const int iDetPos = (curDetPos + halfScale) >> 16; + 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; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos]; + const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor; + 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; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + const kint32 iDetPos = curDetPos >> scaleShift; + const kint32 detRemainder = curDetPos & scaleBitmask; + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy } //end linear - } // end for ix - + } // end for ix + delete deltaFilteredProj; } @@ -641,10 +649,10 @@ 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); @@ -652,23 +660,23 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const 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; + 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; + 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 @@ -678,37 +686,37 @@ 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)); + 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); + 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