X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=d3cb2035dc125af3415006a2d204cc61b5cc7907;hb=6850134e4711a842cc7c0a306a9e4243df06b952;hp=d05bccae2cbcaeb2fac8322d3452a03fc04de44b;hpb=99dd1d6ed10db1f669a5fe6af71225a50fc0ddfb;p=ctsim.git diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index d05bcca..d3cb203 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.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: backprojectors.cpp,v 1.6 2000/07/09 08:16:17 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,6 +26,25 @@ #include "ct.h" +Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName) +{ + m_fail = false; + m_pBackprojectImplem = NULL; + + initBackprojector (proj, im, backprojName, interpName); +} + +void +Backprojector::BackprojectView (const double* const viewData, const double viewAngle) +{ + if (m_pBackprojectImplem) + m_pBackprojectImplem->BackprojectView (viewData, viewAngle); +} + +Backprojector::~Backprojector (void) +{ + delete m_pBackprojectImplem; +} // FUNCTION IDENTIFICATION // Backproject* projector = selectBackprojector (...) @@ -34,24 +53,143 @@ // Selects a backprojector based on BackprojType // and initializes the backprojector -Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType) +bool +Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName) { - Backproject* bj = NULL; - - if (bjType == O_BPROJ_TRIG) - bj = static_cast(new BackprojectTrig (proj, im, interpType)); - else if (bjType == O_BPROJ_TABLE) - bj = static_cast(new BackprojectTable (proj, im, interpType)); - else if (bjType == O_BPROJ_DIFF) - bj = static_cast(new BackprojectDiff (proj, im, interpType)); - else if (bjType == O_BPROJ_DIFF2) - bj = static_cast(new BackprojectDiff2 (proj, im, interpType)); - else if (bjType == O_BPROJ_IDIFF2) - bj = static_cast(new BackprojectIntDiff2 (proj, im, interpType)); - else - sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]"); - - return (bj); + m_nameBackproject = backprojName; + m_nameInterpolation = interpName; + m_pBackprojectImplem = NULL; + m_idBackproject = convertBackprojectNameToID (backprojName); + if (m_idBackproject == BPROJ_INVALID) { + m_fail = true; + m_failMessage = "Invalid backprojection name "; + m_failMessage += backprojName; + } + m_idInterpolation = convertInterpolationNameToID (interpName); + if (m_idInterpolation == INTERP_INVALID) { + m_fail = true; + 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 (m_idBackproject == BPROJ_TRIG) + m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_TABLE) + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_DIFF) + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_DIFF2) + m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_IDIFF2) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_IDIFF3) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation)); + else { + m_fail = true; + m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; + return false; + } + + return true; +} + + +const Backprojector::BackprojectID +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; + + return (backprojID); +} + +const char* +Backprojector::convertBackprojectIDToName (const BackprojectID 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; + + return (bprojName); +} + + + +const Backprojector::InterpolationID +Backprojector::convertInterpolationNameToID (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_PREINTERPOLATE_STR) == 0) + interpID = INTERP_FREQ_PREINTERPOLATE; +#if HAVE_BSPLINE_INTERP + else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0) + interpID = INTERP_BSPLINE; +#endif + + return (interpID); +} + + +/* 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 + */ + +const char* +Backprojector::convertInterpolationIDToName (const InterpolationID interpID) +{ + if (interpID == INTERP_NEAREST) + return (INTERP_NEAREST_STR); + else if (interpID == INTERP_LINEAR) + return (INTERP_LINEAR_STR); + else if (interpID == INTERP_FREQ_PREINTERPOLATE) + return (INTERP_FREQ_PREINTERPOLATE_STR); +#if HAVE_BSPLINE_INTERP + else if (interpID == INTERP_BSPLINE) + return (INTERP_BSPLINE_STR); +#endif + else + return (""); } @@ -61,7 +199,7 @@ Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, // PURPOSE // Pure virtual base class for all backprojectors. -Backproject::Backproject (const Projections& proj, ImageFile& im, const InterpolationType interpType) +Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType) : proj(proj), im(im), interpType(interpType) { detInc = proj.detInc(); @@ -82,7 +220,7 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Interpol xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; - if (interpType != I_NEAREST && interpType != I_LINEAR) + if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR) sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType); } @@ -133,14 +271,14 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double double phi = atan2 (y, x); // angle of cell from center double L = r * cos (theta - phi); // position on detector - if (interpType == I_NEAREST) { + 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] += rotInc * filteredProj[iDetPos]; - } else if (interpType == I_LINEAR) { + } else if (interpType == Backprojector::INTERP_LINEAR) { double p = L / detInc; // position along detector double pFloor = floor (p); int iDetPos = iDetCenter + static_cast(pFloor); @@ -160,7 +298,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, InterpolationType interpType) +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) : Backproject::Backproject (proj, im, interpType) { arrayR.initSetSize (nx, ny); @@ -193,14 +331,14 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl for (int iy = 0; iy < ny; iy++) { double L = r[ix][iy] * cos (theta - phi[ix][iy]); - if (interpType == I_NEAREST) { + 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]; - } else if (interpType == I_LINEAR) { + } else if (interpType == Backprojector::INTERP_LINEAR) { double dPos = L / detInc; // position along detector double dPosFloor = floor (dPos); int iDetPos = iDetCenter + static_cast(dPosFloor); @@ -222,7 +360,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, InterpolationType interpType) +BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) : Backproject::Backproject (proj, im, interpType) { // calculate center of first pixel v[0][0] @@ -255,14 +393,14 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double #ifdef DEBUG printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos); #endif - if (interpType == I_NEAREST) { + 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]; - } else if (interpType == I_LINEAR) { + } else if (interpType == Backprojector::INTERP_LINEAR) { double detPos = curDetPos / detInc; // position along detector double detPosFloor = floor (detPos); int iDetPos = iDetCenter + static_cast(detPosFloor); @@ -306,14 +444,14 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl #ifdef DEBUG printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(L))]); #endif - if (interpType == I_NEAREST) { + 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]; - } else if (interpType == I_LINEAR) { + } else if (interpType == Backprojector::INTERP_LINEAR) { double detPosFloor = floor (curDetPos); int iDetPos = iDetCenter + static_cast(detPosFloor); double frac = curDetPos - detPosFloor; // fraction distance from det @@ -337,26 +475,22 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do { double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle -#if SIZEOF_LONG == 8 - long int scale = 1 << 32; -#else - long int scale = 1 << 16; -#endif - double dScale = scale; - long int halfScale = scale / 2; + static const kint32 scale = 1 << 16; + static const double dScale = scale; + static const kint32 halfScale = scale / 2; - long int det_dx = nearest (xInc * sin (theta) / detInc * scale); - long int det_dy = nearest (yInc * cos (theta) / detInc * scale); + const kint32 det_dx = nearest (xInc * sin (theta) / detInc * scale); + const kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); // calculate L for first point in image (0, 0) - long int detPosColStart = nearest (start_r * cos (theta - start_phi) / detInc * scale); + kint32 detPosColStart = nearest (start_r * cos (theta - start_phi) / detInc * scale); for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - long int curDetPos = detPosColStart; + kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - if (interpType == I_NEAREST) { + 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 @@ -364,9 +498,9 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else *pImCol++ += filteredProj[iDetPos]; - } else if (interpType == I_LINEAR) { - long int detPosFloor = curDetPos / scale; - long int detPosRemainder = curDetPos % scale; + } else if (interpType == Backprojector::INTERP_LINEAR) { + kint32 detPosFloor = curDetPos / scale; + kint32 detPosRemainder = curDetPos % scale; if (detPosRemainder < 0) { detPosFloor--; detPosRemainder += scale; @@ -376,8 +510,59 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do if (iDetPos < 0 || iDetPos >= nDet - 1) errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else - *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x } + +// CLASS IDENTICATION +// BackprojectIntDiff3 +// +// PURPOSE +// Highly optimized version of BackprojectIntDiff2 + +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 + static const int scaleShift = 16; + static const kint32 scale = (1 << scaleShift); + static const kint32 scaleBitmask = scale - 1; + static const double dScale = scale; + 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); + + // calculate L for first point in image (0, 0) + kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + + 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) { + int iDetPos = (curDetPos + halfScale) >> 16; + assert(iDetPos >= 0 && iDetPos < nDet); + *pImCol++ += filteredProj[iDetPos]; + } // end for iy + } else if (interpType == Backprojector::INTERP_LINEAR) { + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { + kint32 iDetPos = curDetPos >> scaleShift; + kint32 detRemainder = curDetPos & scaleBitmask; +#if 0 + double frac = detRemainder * (1. / 65536.); + assert(iDetPos >= 0 && iDetPos < nDet - 1); + *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]); +#else + assert(iDetPos >= 0 && iDetPos < nDet - 1); + const double* const detPointer = &filteredProj[iDetPos]; + *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale; +#endif + } // end for iy + } //end linear + } // end for ix +}