r508: no message

diff --git a/ChangeLog b/ChangeLog
index 71f6deee69af8fba2f6ac4a9ff483977144b7fa3..bfecccd0733b78aecb51b887d2b82b55b22d371f 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -22,6 +22,8 @@
        * ctsim: Added accelerator (hotkeys) to frames
 
        * ctsim: Online help added as well as HTML help
+
+       * backprojector: Added cubic interpolation
        
        * ctsim: Added icons to Frames on X-Window and MS Windows versions
 

diff --git a/NEWS b/NEWS
index 7a1b9959ab2b6e727df1b2d91331811e5b824537..78e00777aaee08f27763cad33907dcc5e6d81f27 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -2,6 +2,8 @@ Version 3.0alpha Features
 
 * Creation of arbitrary filter images
 
+* Now with cubic interpolation for reconstructions
+
 * All features of command line tools are now in graphical ctsim program!
 
 * Complex-valued image files now supported, lots of image math

diff --git a/include/backprojectors.h b/include/backprojectors.h
index 73dabe579cfe799e4c77e48ac6f339dba5c0796e..d9144deb6144504b91d5cb14f8241cc8602c2347 100644 (file)
--- a/include/backprojectors.h
+++ b/include/backprojectors.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: backprojectors.h,v 1.19 2001/01/28 19:10:18 kevin Exp $
+**  $Id: backprojectors.h,v 1.20 2001/02/09 01:54:20 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
@@ -52,6 +52,7 @@ class Backprojector
   static const int INTERP_INVALID;
   static const int INTERP_NEAREST;
   static const int INTERP_LINEAR;
+  static const int INTERP_CUBIC;
   static const int INTERP_FREQ_PREINTERPOLATION;
 #if HAVE_BSPLINE_INTERP
   static const int INTERP_BSPLINE;

diff --git a/include/cubicinterp.h b/include/cubicinterp.h
new file mode 100644 (file)
index 0000000..9e677cd
--- /dev/null
+++ b/include/cubicinterp.h
@@ -0,0 +1,35 @@
+/*****************************************************************************
+**  This is part of the CTSim program
+**  Copyright (c) 1983-2001 Kevin Rosenberg
+**
+**  $Id: cubicinterp.h,v 1.1 2001/02/09 01:54:20 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
+**  published by the Free Software Foundation.
+**
+**  This program is distributed in the hope that it will be useful,
+**  but WITHOUT ANY WARRANTY; without even the implied warranty of
+**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+**  GNU General Public License for more details.
+**
+**  You should have received a copy of the GNU General Public License
+**  along with this program; if not, write to the Free Software
+**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+******************************************************************************/
+
+class CubicInterpolator {
+private:
+  double *m_pdY2;  // second differential of y data
+
+  const double* const m_pdY;
+  const int m_n;
+
+public:
+  CubicInterpolator (const double* const y, int n);
+
+  ~CubicInterpolator ();
+
+  double interpolate (double x);
+};
+

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 0a386fb84eba92577e9b76937b518553d2df68cf..8afe1892054fb7969fdc878b290b1e587e1b55ec 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.24 2001/01/27 21:02:20 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.25 2001/02/09 01:54:20 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
@@ -59,18 +59,20 @@ const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / si
 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;
+const int Backprojector::INTERP_CUBIC = 2;
+const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
 #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;
+const int Backprojector::INTERP_BSPLINE = 4;
+const int Backprojector::INTERP_1BSPLINE = 5;
+const int Backprojector::INTERP_2BSPLINE = 6;
+const int Backprojector::INTERP_3BSPLINE = 7;
 #endif
 
 const char* Backprojector::s_aszInterpName[] = 
 {
   {"nearest"},
   {"linear"},
+  {"cubic"},
   {"freq_preinterpolationj"},
 #if HAVE_BSPLINE_INTERP
   {"bspline"},
@@ -84,6 +86,7 @@ const char* Backprojector::s_aszInterpTitle[] =
 {
   {"Nearest"},
   {"Linear"},
+  {"Cubic"},
   {"Frequency Preinterpolation"},
 #if HAVE_BSPLINE_INTERP
   {"B-Spline"},
@@ -330,6 +333,10 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 {
   double theta = view_angle;
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   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;
@@ -341,22 +348,25 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
       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
+        if (iDetPos >= 0 && iDetPos < nDet)
           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<int>(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
+        if (iDetPos >= 0 && iDetPos < nDet - 1)
           v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+      } else if (interpType = Backprojector::INTERP_CUBIC) {
+        double p = iDetCenter + (L / detInc);  // position along detector
+        if (p >= 0 && p < nDet)
+          v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
       }
     }
   }
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }  
 
 
@@ -393,6 +403,10 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
 {
   double theta = view_angle;
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
     
@@ -402,22 +416,25 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
       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
+        if (iDetPos >= 0 && iDetPos < nDet)
           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<int>(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
+        if (iDetPos >= 0 && iDetPos < nDet - 1)
           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+      } else if (interpType = Backprojector::INTERP_CUBIC) {
+        double p = iDetCenter + (L / detInc);  // position along detector
+        if (p >= 0 && p < nDet)
+          pImCol[iy] += pCubicInterp->interpolate (p);
       }
     }  // end for y 
   }    // end for x 
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 
 
@@ -453,6 +470,10 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
   double det_dy = yInc * sin (theta);
   double lColStart = start_r * cos (theta - start_phi);  // calculate L for first point in image
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
     double curDetPos = lColStart;
     ImageFileColumn pImCol = v[ix];
@@ -464,22 +485,25 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
       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
+        if (iDetPos >= 0 && iDetPos < nDet)
           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<int>(detPosFloor);
         double frac = detPos - detPosFloor;    // fraction distance from det 
-        if (iDetPos < 0 || iDetPos >= nDet - 1)
-          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-        else
+        if (iDetPos >= 0 && iDetPos < nDet - 1)
           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+      } else if (interpType = Backprojector::INTERP_CUBIC) {
+        double p = iDetCenter + (curDetPos / detInc);  // position along detector
+        if (p >= 0 && p < nDet)
+          pImCol[iy] += pCubicInterp->interpolate (p);
       }
     }  // end for y 
   }    // end for x 
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 
 
@@ -494,6 +518,10 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl
 {
   double theta = view_angle;
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   // 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;
@@ -515,21 +543,24 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl
       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
+        if (iDetPos >= 0 && iDetPos < nDet)
           *pImCol++ += filteredProj[iDetPos];
       } 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 
-        if (iDetPos < 0 || iDetPos >= nDet - 1)
-          errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
-        else
+        if (iDetPos > 0 && iDetPos < nDet - 1)
           *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
+      } else if (interpType = Backprojector::INTERP_CUBIC) {
+        double p = iDetCenter + curDetPos;     // position along detector
+        if (p >= 0 && p < nDet)
+          *pImCol++  += pCubicInterp->interpolate (p);
       }
     }  // end for y
   }    // end for x
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 
 // CLASS IDENTICATION
@@ -543,6 +574,10 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
 {
   double theta = view_angle;
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   static const kint32 scale = 1 << 16;
   static const double dScale = scale;
   static const kint32 halfScale = scale / 2;
@@ -579,9 +614,16 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
         else
           *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+      } else if (interpType = Backprojector::INTERP_CUBIC) {
+        double p = iDetCenter + (static_cast<double>(curDetPos) / scale);      // position along detector
+        if (p >= 0 && p < nDet)
+          *pImCol++  += pCubicInterp->interpolate (p);
       }
     }  // end for y
   }    // end for x
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 
 // CLASS IDENTICATION
@@ -606,13 +648,17 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
   
-  // precalculate scaled difference for linear interpolation
-  double* deltaFilteredProj = new double [nDet];
+  double* deltaFilteredProj = NULL;  
+  CubicInterpolator* pCubicInterp = NULL;
   if (interpType == Backprojector::INTERP_LINEAR) {
+    // precalculate scaled difference for linear interpolation
+    deltaFilteredProj = new double [nDet];
     for (int i = 0; i < nDet - 1; i++)
       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
+    deltaFilteredProj[nDet - 1] = 0;  // last detector
+  } else if (interpType == Backprojector::INTERP_CUBIC) {
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
   }
-  deltaFilteredProj[nDet - 1] = 0;  // last detector
   
   int iLastDet = nDet - 1;
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
@@ -638,10 +684,17 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
         if (iDetPos >= 0 && iDetPos <= iLastDet)
           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
       }        // end for iy
-    } //end linear
+    } else if (interpType = Backprojector::INTERP_CUBIC) {
+      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+        *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
+      }
+    } // end cubic
   } // end for ix
   
-  delete deltaFilteredProj;
+  if (interpType == Backprojector::INTERP_LINEAR)
+    delete deltaFilteredProj;
+  else if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 
 
@@ -650,6 +703,10 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const
 {
   double beta = view_angle;
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
     
@@ -659,27 +716,29 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const
       double rsin_t = r[ix][iy] * sin (dAngleDiff);
       double dFLPlusSin = m_dFocalLength + rsin_t;
       double gamma =  atan (rcos_t / dFLPlusSin);
+      double dPos = gamma / detInc;  // position along detector
       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
       
       if (interpType == Backprojector::INTERP_NEAREST) {
-        int iDetPos =iDetCenter + nearest<int>(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
+        int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector      
+        if (iDetPos >= 0 && iDetPos < nDet)
           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<int>(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
+        if (iDetPos >= 0 && iDetPos < nDet - 1)
           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
+      } else if (interpType == Backprojector::INTERP_CUBIC) {
+        double d = iDetCenter + dPos;          // position along detector 
+        if (d >= 0 && d < nDet)
+          pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
       }
     }  // end for y 
   }    // end for x 
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 
 void
@@ -687,6 +746,10 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const
 {
   double beta = view_angle;
   
+  CubicInterpolator* pCubicInterp = NULL;
+  if (interpType == Backprojector::INTERP_CUBIC)
+    pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+  
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
     
@@ -698,27 +761,31 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const
       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
+      // of phantom, see Kak-Slaney Figure 3.22. This assumes that the detector is also
+      // located focal-length away from the origin.
       dDetPos *= 2; 
-      
+      double dPos = dDetPos / detInc;  // position along detector array 
+
       if (interpType == Backprojector::INTERP_NEAREST) {
-        int iDetPos = iDetCenter + nearest<int>(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
+        int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector 
+        if (iDetPos >= 0 && iDetPos < nDet)    
           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<int>(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);
+        if (iDetPos >= 0 && iDetPos < nDet - 1)
+          pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
+                           / (dU * dU);
+      } else if (interpType == Backprojector::INTERP_CUBIC) {
+        double d = iDetCenter + dPos;          // position along detector 
+        if (d >= 0 && d < nDet)
+          pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
       }
     }  // end for y 
   }    // end for x 
+
+  if (interpType == Backprojector::INTERP_CUBIC)
+    delete pCubicInterp;
 }
 

diff --git a/libctsupport/cubicinterp.cpp b/libctsupport/cubicinterp.cpp
new file mode 100644 (file)
index 0000000..8392dec
--- /dev/null
+++ b/libctsupport/cubicinterp.cpp
@@ -0,0 +1,75 @@
+/*****************************************************************************
+**  This is part of the CTSim program
+**  Copyright (c) 1983-2001 Kevin Rosenberg
+**
+**  $Id: cubicinterp.cpp,v 1.1 2001/02/09 01:54:21 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
+**  published by the Free Software Foundation.
+**
+**  This program is distributed in the hope that it will be useful,
+**  but WITHOUT ANY WARRANTY; without even the implied warranty of
+**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+**  GNU General Public License for more details.
+**
+**  You should have received a copy of the GNU General Public License
+**  along with this program; if not, write to the Free Software
+**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+******************************************************************************/
+
+
+#include "ctsupport.h"
+#include "cubicinterp.h"
+
+
+CubicInterpolator::CubicInterpolator (const double* const y, const int n)
+  : m_pdY(y), m_n(n)
+{
+  m_pdY2 = new double [n];
+  m_pdY2[0] = 0.0;   // second deriviative = 0 at beginning and end
+  m_pdY2[n-1] = 0;
+
+  double* temp = new double [n - 1];
+  temp[0] = 0.0
+  for (int i = 1; i < n - 1; i++) { 
+    double t = 2 + (0.5 * m_pdY2[i-1]);
+    temp[i] = y[i+1] + y[i-1] - y[i] - y[i];
+    temp[i] = (3 * temp[i] - 0.5 * temp[i-1]) / t;
+    m_pdY2[i] = -0.5 / t;
+  }
+
+  for (i = n - 2; i >= 0; i--) 
+    m_pdY2[i] = temp[i] + m_pdY2[i] * m_pdY2[i + 1];
+
+  delete temp;
+}
+
+CubicInterpolator::~CubicInterpolator ()
+{
+  delete m_pdY2;
+}
+
+
+double
+CubicInterpolator::interpolate (double x)
+{
+  const static double oneSixth = (1. / 6.);
+  int lo = static_cast<int>(floor(x));
+  int hi = lo + 1;
+
+  if (lo < 0 || hi >= m_n) {
+    sys_error (ERR_SEVERE, "X range out of bounds [CubicInterpolator::interpolate]");
+    return (0);
+  }
+
+  double loFr = hi - x;
+  double hiFr = 1 - a;
+  double y = loFr * m_pdY[lo] + hiFr * m_pdY[hi];
+  y += oneSixth * ((loFr*loFr*loFr - loFr) * m_pdY2[lo] + (hiFr*hiFr*hiFr - hiFr) * m_pdY2[hi]);
+  
+  return y;
+}
+
+
+

diff --git a/msvc/ctsim/ctsim.plg b/msvc/ctsim/ctsim.plg
index c131702410dc6cf2a26f8bc18c37f3a1d474c2f2..72e9dd791ffdcacbad8d540835be8e82555b0a44 100644 (file)
--- a/msvc/ctsim/ctsim.plg
+++ b/msvc/ctsim/ctsim.plg
@@ -3,16 +3,60 @@
 <pre>
 <h1>Build Log</h1>
 <h3>
---------------------Configuration: ctsim - Win32 Debug--------------------
+--------------------Configuration: libctsim - Win32 Debug--------------------
 </h3>
 <h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP10B.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP108.tmp" with contents
+[
+/nologo /G6 /MTd /W3 /Gm /Gi /GR /GX /Zi /Od /Gy /I "..\..\..\wx2.2.5\src\png" /I "..\..\..\wx2.2.5\src\zlib" /I "..\..\INCLUDE" /I "..\..\getopt" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2.2.5\include" /D "_DEBUG" /D "HAVE_WXWIN" /D "HAVE_STRING_H" /D "HAVE_GETOPT_H" /D "WIN32" /D "_MBCS" /D "_LIB" /D "MSVC" /D "HAVE_FFTW" /D "HAVE_PNG" /D "HAVE_SGP" /D "HAVE_WXWINDOWS" /D "__WXMSW__" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /FR"Debug/" /Fp"Debug/libctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c 
+"C:\ctsim\libctsim\backprojectors.cpp"
+]
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP108.tmp" 
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP109.tmp" with contents
 [
-/nologo /G6 /MTd /W3 /Gm /Gi /GR /GX /Zi /Od /Gy /I "\wx2.2.5\include" /I "..\..\..\fftw-2.1.3\fftw" /I "\wx2.2.5\src\png" /I "\wx2.2.5\src\zlib" /I "..\..\include" /I "..\..\getopt" /I "..\..\..\fftw-2.1.3\rfftw" /D VERSION=\"3.0.0beta1\" /D "_DEBUG" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.0.0alpha5\" /FR"Debug/" /Fp"Debug/ctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c 
-"C:\ctsim\src\views.cpp"
+/nologo /out:"Debug\libctsim.lib" 
+.\Debug\array2dfile.obj
+.\Debug\backprojectors.obj
+.\Debug\clip.obj
+.\Debug\consoleio.obj
+.\Debug\dlgezplot.obj
+.\Debug\ezplot.obj
+.\Debug\ezset.obj
+.\Debug\ezsupport.obj
+.\Debug\filter.obj
+.\Debug\fnetorderstream.obj
+.\Debug\fourier.obj
+.\Debug\getopt.obj
+.\Debug\getopt1.obj
+.\Debug\globalvars.obj
+.\Debug\hashtable.obj
+.\Debug\imagefile.obj
+.\Debug\mathfuncs.obj
+.\Debug\phantom.obj
+.\Debug\plotfile.obj
+.\Debug\pol.obj
+.\Debug\procsignal.obj
+.\Debug\projections.obj
+.\Debug\reconstruct.obj
+.\Debug\scanner.obj
+.\Debug\sgp.obj
+.\Debug\strfuncs.obj
+.\Debug\syserror.obj
+.\Debug\trace.obj
+.\Debug\transformmatrix.obj
+.\Debug\xform.obj
+.\Debug\cubicinterp.obj
 ]
-Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP10B.tmp" 
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP10C.tmp" with contents
+Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP109.tmp"
+<h3>Output Window</h3>
+Compiling...
+backprojectors.cpp
+Creating library...
+<h3>
+--------------------Configuration: ctsim - Win32 Debug--------------------
+</h3>
+<h3>Command Lines</h3>
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP10A.tmp" with contents
 [
 comctl32.lib winmm.lib rpcrt4.lib ws2_32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ../libctsim/Debug/libctsim.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib libcmtd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib wxd.lib xpmd.lib tiffd.lib zlibd.lib pngd.lib opengl32.lib glu32.lib /nologo /subsystem:windows /incremental:yes /pdb:"Debug/ctsim.pdb" /debug /machine:I386 /out:"Debug/ctsim.exe" /pdbtype:sept /libpath:"\wx2.2.5\lib" 
 .\Debug\ctsim.obj
@@ -33,10 +77,8 @@ comctl32.lib winmm.lib rpcrt4.lib ws2_32.lib kernel32.lib user32.lib gdi32.lib w
 \wx2.2.5\lib\zlibd.lib
 \wx2.2.5\lib\tiffd.lib
 ]
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP10C.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP10A.tmp"
 <h3>Output Window</h3>
-Compiling...
-views.cpp
 Linking...
 
 

diff --git a/src/views.cpp b/src/views.cpp
index c1debac9f66427f8f1718f5113e8791f75876840..2f4d443235751acc7fcebbddc8f4830f7c58cbfd 100644 (file)
--- a/src/views.cpp
+++ b/src/views.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.99 2001/02/08 06:25:07 kevin Exp $
+**  $Id: views.cpp,v 1.100 2001/02/09 01:54:21 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
@@ -268,7 +268,7 @@ ImageFileView::OnProperties (wxCommandEvent& event)
   if (rIF.nx() == 0 || rIF.ny() == 0)
     *theApp->getLog() << "Properties: empty imagefile\n";
   else {
-    const std::string& rFilename = rIF.getFilename();
+    const std::string rFilename = m_pFrame->GetTitle().c_str();
     std::ostringstream os;
     double min, max, mean, mode, median, stddev;
     rIF.statistics (rIF.getArray(), min, max, mean, mode, median, stddev);