r138: *** empty log message ***

diff --git a/ChangeLog b/ChangeLog
index 6cd4d1af0a6fd67c9b9ae3a33635143a66b20dc1..d48b1ae75cf99daf8f98aa4beecb4ed1c52d9299 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,8 +1,10 @@
+2.0.0-b3 - 7/09/00
+       
 2.0.0-b2 - 7/07/00
-   Added zeropadding option to pjrec
    Cleaned up SignalFilter class
+   Added zeropad option to pjrec
    Added zeropad options to html and cgi files
-   Added fourier_table filter method
+   Added fourier_table and rfttw filter methods
    Added FFTW routines to use real/half-complex transformations
        
 2.0.0-b1 - 7/05/00

diff --git a/configure.in b/configure.in
index 121694d051f38d92e26fb4c24f88efc391d3218e..776c723fe12af6a03f60b5bb7f6b8ecc4df1794c 100644 (file)
--- a/configure.in
+++ b/configure.in
@@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout
 dnl CDPATH=
 
 AC_INIT(src/pjrec.cpp)
-AM_INIT_AUTOMAKE(ctsim,2.0.0-b2)
+AM_INIT_AUTOMAKE(ctsim,2.0.0-b3)
 AM_CONFIG_HEADER(config.h)
 
 dnl Checks for programs.

diff --git a/include/backprojectors.h b/include/backprojectors.h
index 6b8b9c4635b76c183161167f14ec2f60e3f2b9b9..0a2a8e12c68b75ed0c91ffc17b42a97cb1e3a659 100644 (file)
--- a/include/backprojectors.h
+++ b/include/backprojectors.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.h,v 1.5 2000/06/25 17:32:24 kevin Exp $
+**  $Id: backprojectors.h,v 1.6 2000/07/07 15:30:59 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
@@ -39,7 +39,8 @@ class Backprojector
     BPROJ_TABLE,
     BPROJ_DIFF,
     BPROJ_DIFF2,
-    BPROJ_IDIFF2
+    BPROJ_IDIFF2,
+    BPROJ_IDIFF3
 } BackprojectID;
 
   typedef enum {
@@ -59,6 +60,7 @@ class Backprojector
   static const char BPROJ_DIFF_STR[]=     "diff";
   static const char BPROJ_DIFF2_STR[]=    "diff2";
   static const char BPROJ_IDIFF2_STR[]=   "idiff2";
+  static const char BPROJ_IDIFF3_STR[]=   "idiff3";
  
   static const char INTERP_NEAREST_STR[]=  "nearest";
   static const char INTERP_LINEAR_STR[]=   "linear";
@@ -187,3 +189,14 @@ class BackprojectIntDiff2 : public BackprojectDiff
 };
 
 
+class BackprojectIntDiff3 : public BackprojectDiff
+{
+ public:
+  BackprojectIntDiff3 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
+    :  BackprojectDiff::BackprojectDiff (proj, im, interpType)
+    {}
+  
+  void BackprojectView (const double* const t, double view_angle);
+};
+
+

diff --git a/include/filter.h b/include/filter.h
index e4c4f3f9db30df46d6da0077fe6b2b381a4333cf..9ee7390a4acd83c96bd9a72b7f9990fbfe7bede9 100644 (file)
--- a/include/filter.h
+++ b/include/filter.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: filter.h,v 1.11 2000/07/06 18:37:24 kevin Exp $
+**  $Id: filter.h,v 1.12 2000/07/07 15:30:59 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,8 +52,10 @@ class SignalFilter {
        FILTER_METHOD_FOURIER,
        FILTER_METHOD_FOURIER_TABLE,
        FILTER_METHOD_FFT,
+#if HAVE_FFTW
        FILTER_METHOD_FFTW,
-       FILTER_METHOD_RFFTW
+       FILTER_METHOD_RFFTW,
+#endif
     } FilterMethodID;
 
     typedef enum {
@@ -77,8 +79,10 @@ class SignalFilter {
     static const char FILTER_METHOD_FOURIER_STR[]=      "fourier";
     static const char FILTER_METHOD_FOURIER_TABLE_STR[]="fourier_table";
     static const char FILTER_METHOD_FFT_STR[]=          "fft";
+#if HAVE_FFTW
     static const char FILTER_METHOD_FFTW_STR[]=         "fftw";
     static const char FILTER_METHOD_RFFTW_STR[]=        "rfftw";
+#endif
 
     static const char DOMAIN_FREQUENCY_STR[]="frequency";
     static const char DOMAIN_SPATIAL_STR[]="spatial";
@@ -136,8 +140,6 @@ class SignalFilter {
 
     static double spatialResponseAnalytic (FilterID fType, double bw, double x, double param);
 
-    static void dotProduct (const double v1[], const complex<double> v2[], complex<double> output[], const int n);
-
  private:
     double m_bw;
     int m_nFilterPoints;
@@ -155,8 +157,6 @@ class SignalFilter {
     rfftw_plan m_realPlanForward, m_realPlanBackward;
     fftw_complex* m_vecComplexFftInput;
     fftw_plan m_complexPlanForward, m_complexPlanBackward;
-#else
-    complex<double>* m_vecFftInput;
 #endif
 
     bool m_fail;

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 1f6e92c1bece4c526d80abdad76e47d1d70a7f75..0684c16c5c638e95e01eec79a8d5988968b5ed04 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
 **
-**  $Id: backprojectors.cpp,v 1.4 2000/06/25 17:32:24 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.5 2000/07/07 15:30:59 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
@@ -87,6 +87,8 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
     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 if (m_idBackproject == BPROJ_IDIFF3)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation));
   else {
     m_fail = true;
     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
@@ -112,6 +114,8 @@ Backprojector::convertBackprojectNameToID (const char* const backprojName)
     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);
 }
@@ -131,6 +135,8 @@ Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
     bprojName = BPROJ_DIFF2_STR;
   else if (bprojID == BPROJ_IDIFF2)
     bprojName = BPROJ_IDIFF2_STR;
+  else if (bprojID == BPROJ_IDIFF3)
+    bprojName = BPROJ_IDIFF3_STR;
 
   return (bprojName);
 }
@@ -500,7 +506,58 @@ 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 double dScale = scale;
+  static const kint32 halfScale = scale / 2;
+
+  const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
+  const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+
+  // calculate L for first point in image (0, 0) 
+  kint32 detPosColStart = nearest<kint32> ((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];
+
+    for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+      if (interpType == Backprojector::INTERP_NEAREST) {
+       int iDetPos = (curDetPos + halfScale) >> scaleShift;
+
+       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) {
+       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

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 6a344d42a6273957a8d847c264bcdebd438c0357..65163f3185b0be2de8beb7342ad643f0701ff875 100644 (file)
--- a/libctsim/filter.cpp
+++ b/libctsim/filter.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: filter.cpp,v 1.14 2000/07/06 18:37:24 kevin Exp $
+**  $Id: filter.cpp,v 1.15 2000/07/07 15:30:59 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
@@ -129,10 +129,21 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   m_vecFourierSinTable = NULL;
   m_vecFilter = NULL;
 
-  if (m_idFilterMethod == FILTER_METHOD_FFT)
-    m_idFilterMethod = FILTER_METHOD_FFTW;
+  if (m_idFilterMethod == FILTER_METHOD_FFT) {
+#if HAVE_FFTW
+    m_idFilterMethod = FILTER_METHOD_RFFTW;
+#else
+    m_fail = true;
+    m_failMessage = "FFT not yet implemented";
+    return;
+#endif
+  }
 
-  if (m_idFilterMethod == FILTER_METHOD_FOURIER || FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER || FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT
+#if HAVE_FFTW
+      || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW
+#endif
+      ) {
     m_nFilterPoints = m_nSignalPoints;
     if (m_zeropad > 0) {
       double logBase2 = log(m_nSignalPoints) / log(2);
@@ -321,10 +332,12 @@ SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
     fmID = FILTER_METHOD_FOURIER_TABLE;
   else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_STR) == 0)
     fmID = FILTER_METHOD_FFT;
+#if HAVE_FFTW
   else if (strcasecmp (filterMethodName, FILTER_METHOD_FFTW_STR) == 0)
     fmID = FILTER_METHOD_FFTW;
   else if (strcasecmp (filterMethodName, FILTER_METHOD_RFFTW_STR) == 0)
     fmID = FILTER_METHOD_RFFTW;
+#endif
 
   return (fmID);
 }
@@ -342,10 +355,12 @@ SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID)
     return (FILTER_METHOD_FOURIER_TABLE_STR);
   else if (fmID == FILTER_METHOD_FFT)
     return (FILTER_METHOD_FFT_STR);
+#if HAVE_FFTW
   else if (fmID == FILTER_METHOD_FFTW)
     return (FILTER_METHOD_FFTW_STR);
   else if (fmID == FILTER_METHOD_RFFTW)
     return (FILTER_METHOD_RFFTW_STR);
+#endif
 
   return (name);
 }
@@ -391,10 +406,10 @@ SignalFilter::filterSignal (const float input[], double output[]) const
       inputSignal[i] = 0;  // zeropad
     complex<double> fftSignal[m_nFilterPoints];
     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
-    complex<double> filteredSignal[m_nFilterPoints];
-    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
+    for (int i = 0; i < m_nFilterPoints; i++)
+      fftSignal[i] *= m_vecFilter[i];
     double inverseFourier[m_nFilterPoints];
-    finiteFourierTransform (filteredSignal, inverseFourier, m_nFilterPoints, 1);
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
@@ -405,10 +420,10 @@ SignalFilter::filterSignal (const float input[], double output[]) const
       inputSignal[i] = 0;  // zeropad
     complex<double> fftSignal[m_nFilterPoints];
     finiteFourierTransform (inputSignal, fftSignal, -1);
-    complex<double> filteredSignal[m_nFilterPoints];
-    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
+    for (int i = 0; i < m_nFilterPoints; i++)
+      fftSignal[i] *= m_vecFilter[i];
     double inverseFourier[m_nFilterPoints];
-    finiteFourierTransform (filteredSignal, inverseFourier, 1);
+    finiteFourierTransform (fftSignal, inverseFourier, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
       output[i] = inverseFourier[i];
   }
@@ -930,9 +945,3 @@ SignalFilter::finiteFourierTransform (const complex<double> input[], double outp
 }
 
 
-void 
-SignalFilter::dotProduct (const double v1[], const complex<double> v2[], complex<double> output[], const int n)
-{
-    for (int i = 0; i < n; i++)
-       output[i] = v1[i] * v2[i];
-}