X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=71aa8d4824722d2940d874b0b87666007121954f;hp=ae31cf866967d5f5b1ebab1910e722100e5eb8c9;hb=1e88cf0f7fa4f690ea9f110e8ed3f2b5338d0a10;hpb=27a474e0622ebb7229fd5705552021f63d8f932d diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index ae31cf8..71aa8d4 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.9 2000/07/20 11:17:31 kevin Exp $ +** $Id: backprojectors.cpp,v 1.12 2000/08/25 15:59:13 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 @@ -26,29 +26,75 @@ #include "ct.h" -const char Backprojector::BPROJ_TRIG_STR[]= "trig"; -const char Backprojector::BPROJ_TABLE_STR[]= "table"; -const char Backprojector::BPROJ_DIFF_STR[]= "diff"; -const char Backprojector::BPROJ_DIFF2_STR[]= "diff2"; -const char Backprojector::BPROJ_IDIFF2_STR[]= "idiff2"; -const char Backprojector::BPROJ_IDIFF3_STR[]= "idiff3"; - -const char Backprojector::BPROJ_TRIG_TITLE_STR[]= "Direc Trigometric"; -const char Backprojector::BPROJ_TABLE_TITLE_STR[]= "Trig Table"; -const char Backprojector::BPROJ_DIFF_TITLE_STR[]= "Diff"; -const char Backprojector::BPROJ_DIFF2_TITLE_STR[]= "Diff2"; -const char Backprojector::BPROJ_IDIFF2_TITLE_STR[]= "Integer Diff2"; -const char Backprojector::BPROJ_IDIFF3_TITLE_STR[]= "Integer Diff3"; - -const char Backprojector::INTERP_NEAREST_STR[]= "nearest"; -const char Backprojector::INTERP_LINEAR_STR[]= "linear"; -const char Backprojector::INTERP_BSPLINE_STR[]= "bspline"; -const char Backprojector::INTERP_FREQ_PREINTERPOLATION_STR[]= "freq_preinterpolation"; - -const char Backprojector::INTERP_NEAREST_TITLE_STR[]= "Nearest"; -const char Backprojector::INTERP_LINEAR_TITLE_STR[]= "Linear"; -const char Backprojector::INTERP_BSPLINE_TITLE_STR[]= "B-Spline"; -const char Backprojector::INTERP_FREQ_PREINTERPOLATION_TITLE_STR[]= "Frequency Preinterpolation"; +const int Backprojector::BPROJ_INVALID = -1; +const int Backprojector::BPROJ_TRIG = 0; +const int Backprojector::BPROJ_TABLE = 1; +const int Backprojector::BPROJ_DIFF = 2; +const int Backprojector::BPROJ_DIFF2 = 3; +const int Backprojector::BPROJ_IDIFF2 = 4; +const int Backprojector::BPROJ_IDIFF3 = 5; + +const char* Backprojector::s_aszBackprojectName[] = +{ + {"trig"}, + {"table"}, + {"diff"}, + {"diff2"}, + {"idiff2"}, + {"idiff3"}, +}; + +const char* Backprojector::s_aszBackprojectTitle[] = +{ + {"Direct Trigometric"}, + {"Trigometric Table"}, + {"Difference Iteration"}, + {"Difference Iteration Optimized"}, + {"Integer Difference Iteration Optimized"}, + {"Integer Difference Iteration Highly-Optimized"}, +}; + +const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*); + +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; +#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; +#endif + +const char* Backprojector::s_aszInterpName[] = +{ + {"nearest"}, + {"linear"}, + {"freq_preinterpolationj"}, +#if HAVE_BSPLINE_INTERP + {"bspline"}, + {"1bspline"}, + {"2bspline"}, + {"3bspline"}, +#endif +}; + +const char* Backprojector::s_aszInterpTitle[] = +{ + {"Nearest"}, + {"Linear"}, + {"Frequency Preinterpolationj"}, +#if HAVE_BSPLINE_INTERP + {"B-Spline"}, + {"B-Spline 1st Order"}, + {"B-Spline 2nd Order"}, + {"B-Spline 3rd Order"}, +#endif +}; + +const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*); + Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) @@ -66,7 +112,7 @@ Backprojector::BackprojectView (const double* const viewData, const double viewA m_pBackprojectImplem->BackprojectView (viewData, viewAngle); } -Backprojector::~Backprojector (void) +Backprojector::~Backprojector () { delete m_pBackprojectImplem; } @@ -90,7 +136,7 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_failMessage = "Invalid backprojection name "; m_failMessage += backprojName; } - m_idInterpolation = convertInterpolationNameToID (interpName); + m_idInterpolation = convertInterpNameToID (interpName); if (m_idInterpolation == INTERP_INVALID) { m_fail = true; m_failMessage = "Invalid interpolation name "; @@ -124,107 +170,88 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const } -const Backprojector::BackprojectID +int Backprojector::convertBackprojectNameToID (const char* const backprojName) { - BackprojectID backprojID = BPROJ_INVALID; - - if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0) - backprojID = BPROJ_TRIG; - else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0) - backprojID = BPROJ_TABLE; - else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0) - backprojID = BPROJ_DIFF; - else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0) - backprojID = BPROJ_DIFF2; - else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0) - backprojID = BPROJ_IDIFF2; - else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0) - backprojID = BPROJ_IDIFF3; + int backprojID = BPROJ_INVALID; + + for (int i = 0; i < s_iBackprojectCount; i++) + if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { + backprojID = i; + break; + } return (backprojID); } const char* -Backprojector::convertBackprojectIDToName (const BackprojectID bprojID) +Backprojector::convertBackprojectIDToName (int bprojID) { - const char *bprojName = ""; - - if (bprojID == BPROJ_TRIG) - bprojName = BPROJ_TRIG_STR; - else if (bprojID == BPROJ_TABLE) - bprojName = BPROJ_TABLE_STR; - else if (bprojID == BPROJ_DIFF) - bprojName = BPROJ_DIFF_STR; - else if (bprojID == BPROJ_DIFF2) - bprojName = BPROJ_DIFF2_STR; - else if (bprojID == BPROJ_IDIFF2) - bprojName = BPROJ_IDIFF2_STR; - else if (bprojID == BPROJ_IDIFF3) - bprojName = BPROJ_IDIFF3_STR; + static const char *bprojName = ""; + + if (bprojID >= 0 && bprojID < s_iBackprojectCount) + return (s_aszBackprojectName[bprojID]); return (bprojName); } +const char* +Backprojector::convertBackprojectIDToTitle (const int bprojID) +{ + static const char *bprojTitle = ""; + + if (bprojID >= 0 && bprojID < s_iBackprojectCount) + return (s_aszBackprojectTitle[bprojID]); + return (bprojTitle); +} -const Backprojector::InterpolationID -Backprojector::convertInterpolationNameToID (const char* const interpName) + +int +Backprojector::convertInterpNameToID (const char* const interpName) { - InterpolationID interpID = INTERP_INVALID; - - if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0) - interpID = INTERP_NEAREST; - else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0) - interpID = INTERP_LINEAR; - else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATION_STR) == 0) - interpID = INTERP_FREQ_PREINTERPOLATION; -#if HAVE_BSPLINE_INTERP - else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0) - interpID = INTERP_BSPLINE; -#endif + int interpID = INTERP_INVALID; + + for (int i = 0; i < s_iInterpCount; i++) + if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { + interpID = i; + break; + } return (interpID); } +const char* +Backprojector::convertInterpIDToName (const int interpID) +{ + static const char *interpName = ""; -/* NAME - * name_of_interp Return name of interpolation method - * - * SYNOPSIS - * name = name_of_interp (interp_type) - * char *name Name of interpolation method - * int interp_type Method of interpolation - * - * NOTES - * Returns NULL if interp_type is invalid - */ + if (interpID >= 0 && interpID < s_iInterpCount) + return (s_aszInterpName[interpID]); + + return (interpName); +} const char* -Backprojector::convertInterpolationIDToName (const InterpolationID interpID) +Backprojector::convertInterpIDToTitle (const int interpID) { - if (interpID == INTERP_NEAREST) - return (INTERP_NEAREST_STR); - else if (interpID == INTERP_LINEAR) - return (INTERP_LINEAR_STR); - else if (interpID == INTERP_FREQ_PREINTERPOLATION) - return (INTERP_FREQ_PREINTERPOLATION_STR); -#if HAVE_BSPLINE_INTERP - else if (interpID == INTERP_BSPLINE) - return (INTERP_BSPLINE_STR); -#endif - else - return (""); + static const char *interpTitle = ""; + + if (interpID >= 0 && interpID < s_iInterpCount) + return (s_aszInterpTitle[interpID]); + + return (interpTitle); } + // CLASS IDENTICATION // Backproject // // PURPOSE // Pure virtual base class for all backprojectors. -Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType, const int interpFactor) +Backproject::Backproject (const Projections& proj, ImageFile& im, const int interpType, const int interpFactor) : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor) { detInc = proj.detInc(); @@ -246,11 +273,11 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Backproj yInc = (yMax - yMin) / ny; } -Backproject::~Backproject (void) +Backproject::~Backproject () {} void -Backproject::ScaleImageByRotIncrement (void) +Backproject::ScaleImageByRotIncrement () { for (int ix = 0; ix < nx; ix++) for (int iy = 0; iy < ny; iy++) @@ -259,18 +286,20 @@ Backproject::ScaleImageByRotIncrement (void) 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); + 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()); } @@ -320,7 +349,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // PURPOSE // Precalculates trigometric function value for each point in image for backprojection. -BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) : Backproject::Backproject (proj, im, interpType, interpFactor) { arrayR.initSetSize (nx, ny); @@ -337,7 +366,7 @@ BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Back } } -BackprojectTable::~BackprojectTable (void) +BackprojectTable::~BackprojectTable () { ScaleImageByRotIncrement(); } @@ -382,7 +411,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // Backprojects by precalculating the change in L position for each x & y step in the image. // Iterates in x & y direction by adding difference in L position -BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) +BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) : Backproject::Backproject (proj, im, interpType, interpFactor) { // calculate center of first pixel v[0][0] @@ -413,7 +442,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 @@ -464,7 +493,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 @@ -516,7 +545,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]; @@ -561,12 +590,14 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); // precalculate scaled difference for linear interpolation - double deltaFilteredProj [nDet - 1]; + double deltaFilteredProj [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]; @@ -574,21 +605,21 @@ 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