r117: *** empty log message ***

[ctsim.git] / libctsim / backprojectors.cpp

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index d05bccae2cbcaeb2fac8322d3452a03fc04de44b..9e350923ab4594287c868e21ee4c8f4b371bf8a0 100644 (file)
--- 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
 **
 **  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.3 2000/06/22 10:17:28 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
 **
 **  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"
 
 
 #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 (...)
 
 // FUNCTION IDENTIFICATION
 //     Backproject* projector = selectBackprojector (...)


@@ -34,24 +53,122 @@
 //     Selects a backprojector based on BackprojType 
 //     and initializes the backprojector
 
 //     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)
+{
+  m_nameBackproject = backprojName;
+  m_nameInterpolation = interpName;
+  m_idBackproject = convertBackprojectNameToID (backprojName);
+  m_idInterpolation = convertInterpolationNameToID (interpName);
+  m_pBackprojectImplem = NULL;
+
+  if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
+    m_fail = true;
+    return false;
+  }
+
+  if (m_idBackproject == BPROJ_TRIG)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_TABLE)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_DIFF)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_DIFF2)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_IDIFF2)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
+  else {
+    m_fail = true;
+    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;
+
+  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;
+
+  return (bprojName);
+}
+
+
+
+const Backprojector::InterpolationID
+Backprojector::convertInterpolationNameToID (const char* const interpName)

 {
 {

-    Backproject* bj = NULL;
-
-    if (bjType == O_BPROJ_TRIG)
-       bj = static_cast<Backproject*>(new BackprojectTrig (proj, im, interpType));
-    else if (bjType == O_BPROJ_TABLE)
-       bj = static_cast<Backproject*>(new BackprojectTable (proj, im, interpType));
-    else if (bjType == O_BPROJ_DIFF)
-       bj = static_cast<Backproject*>(new BackprojectDiff (proj, im, interpType));
-    else if (bjType == O_BPROJ_DIFF2)
-       bj = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, interpType));
-    else if (bjType == O_BPROJ_IDIFF2)
-       bj = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, interpType));
-    else 
-      sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
-
-    return (bj);
+  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;
+#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);
+#if HAVE_BSPLINE_INTERP
+  else if (interpID == INTERP_BSPLINE)
+    return (INTERP_BSPLINE_STR);
+#endif
+  else
+    return ("");

 }
 
 
 }
 
 


@@ -61,7 +178,7 @@ Backproject* selectBackprojector (BackprojType bjType, const Projections& proj,
 // PURPOSE
 //   Pure virtual base class for all backprojectors.
 
 // 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();
   : proj(proj), im(im), interpType(interpType)
 {
   detInc = proj.detInc();


@@ -82,7 +199,7 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Interpol
   xInc = (xMax - xMin) / nx;   // size of cells
   yInc = (yMax - yMin) / ny;
 
   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);
 }
 
     sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
 }
 


@@ -133,14 +250,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
 
       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<int> (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];
        int iDetPos = iDetCenter + nearest<int> (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<int>(pFloor);
          double p = L / detInc;        // position along detector
          double pFloor = floor (p);
          int iDetPos = iDetCenter + static_cast<int>(pFloor);


@@ -160,7 +277,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 // PURPOSE
 //   Precalculates trigometric function value for each point in image for backprojection.
 
 // 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);
   : Backproject::Backproject (proj, im, interpType)
 {
   arrayR.initSetSize (nx, ny);


@@ -193,14 +310,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]);
 
     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<int>(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<int>(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<int>(dPosFloor);
        double dPos = L / detInc;               // position along detector 
        double dPosFloor = floor (dPos);
        int iDetPos = iDetCenter + static_cast<int>(dPosFloor);


@@ -222,7 +339,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
 
 //   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] 
   :  Backproject::Backproject (proj, im, interpType)
 {
   // calculate center of first pixel v[0][0] 


@@ -255,14 +372,14 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
 #ifdef DEBUG
       printf ("[%2d,%2d]:  %8.5lf  ", ix, iy, curDetPos);
 #endif
 #ifdef DEBUG
       printf ("[%2d,%2d]:  %8.5lf  ", ix, iy, curDetPos);
 #endif

-      if (interpType == I_NEAREST) {
+      if (interpType == Backprojector::INTERP_NEAREST) {

        int iDetPos = iDetCenter + nearest<int>(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<int>(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<int>(detPosFloor);
        double detPos = curDetPos / detInc;             // position along detector 
        double detPosFloor = floor (detPos);
        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);


@@ -306,14 +423,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<int>(L))]);
 #endif
 #ifdef DEBUG
       printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
 #endif

-      if (interpType == I_NEAREST) {
+      if (interpType == Backprojector::INTERP_NEAREST) {

        int iDetPos = iDetCenter + nearest<int> (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<int> (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<int>(detPosFloor);
        double frac = curDetPos - detPosFloor;  // fraction distance from det 
        double detPosFloor = floor (curDetPos);
        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
        double frac = curDetPos - detPosFloor;  // fraction distance from det 


@@ -337,26 +454,22 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
 {
   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
 

-#if SIZEOF_LONG == 8
-  long int scale = 1 << 32;
-#else
-  long int scale = 1 << 16;
-#endif
+  kint32 scale = 1 << 16;

   double dScale = scale;
   double dScale = scale;

-  long int halfScale = scale / 2;
+  kint32 halfScale = scale / 2;

 
 

-  long int det_dx = nearest<long int> (xInc * sin (theta) / detInc * scale);
-  long int det_dy = nearest<long int> (yInc * cos (theta) / detInc * scale);
+  kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
+  kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);

 
   // calculate L for first point in image (0, 0) 
 
   // calculate L for first point in image (0, 0) 

-  long int detPosColStart = nearest<long int> (start_r * cos (theta - start_phi) / detInc * scale);
+  kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);

        
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
        
   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) {
     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 
 
        int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
        int iDetPos = iDetCenter + detPosNearest;       // calc index in the filtered raysum vector 
 


@@ -364,9 +477,9 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
            errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
        else
          *pImCol++ += filteredProj[iDetPos];
            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;
        if (detPosRemainder < 0) {
          detPosFloor--;
          detPosRemainder += scale;