r127: *** empty log message ***

diff --git a/ChangeLog b/ChangeLog
index 93df882c2a0a7052f1d2a1f52757ba0c3e159931..023c1c7d00bdfd6c3c3fb10faee8bbf145a2859c 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -3,6 +3,7 @@
    Fixed Array2dFile::labelsCopy()
    Added copy constructor and assignment for Array2dFileLabel class
    Added Timer to if-2.cpp and ifinfo.cpp
+   Added beginning of frequency (FFT) based filter to SignalFilter
        
 1.9.8 - 6/27/2000
    Rewrote Array2dFile class to be non-templated

diff --git a/include/ct.h b/include/ct.h
index 3788f9df5955e94b600228458f68acaab0defe2d..54dd3cfb52689f7c89c0bdce332489f6d19b7599 100644 (file)
--- a/include/ct.h
+++ b/include/ct.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ct.h,v 1.23 2000/06/26 21:15:24 kevin Exp $
+**  $Id: ct.h,v 1.24 2000/07/02 18:21:39 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
@@ -118,6 +118,7 @@ extern "C" {
 #include <exception>
 #include <stdexcept>
 #include <memory>
+#include <complex>
 
 using namespace std;
 

diff --git a/include/filter.h b/include/filter.h
index 19266a0edbce9e909d39bb127436f7de8a86fa68..5f5789b6fa2e857d1271ff042f4f8414b6adbb44 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.5 2000/06/30 02:03:27 kevin Exp $
+**  $Id: filter.h,v 1.6 2000/07/02 18:21:39 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
@@ -83,9 +83,9 @@ class SignalFilter {
     static const char DOMAIN_SPATIAL_STR[]= "spatial";
 
 
-    SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, const int numIntegral = 0);
+    SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalLength, int n, double param, const char* domainName, const int numIntegral = 0);
 
-    SignalFilter (const FilterID filt_type, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numIntegral = 0);
+    SignalFilter (const FilterID filt_type, FilterMethodID filterMethodID, double bw, double signalLength, int n, double param, const DomainID domain, const int numIntegral = 0);
 
     SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0);
 
@@ -94,10 +94,20 @@ class SignalFilter {
     double* getFilter (void) const
       { return m_vecFilter; }
 
+    int getNFilterPoints (void) const
+       { return m_nFilterPoints; }
+
     double convolve (const double f[], const double dx, const int n, const int np) const;
 
     double convolve (const float f[], const double dx, const int n, const int np) const;
 
+    void filterSignal (const double input[], double output[], double dx, const int n) const;
+    void filterSignal (const float input[], double output[], double dx, const int n) const;
+
+    static void finiteFourierTransform (const double input[], complex<double> output[], const int n, const int direction);
+
+    void finiteFourierTransform (const double input[], complex<double> output[], const int direction) const;
+
     bool fail(void) const      {return m_fail;}
     const string& failMessage(void) const {return m_failMessage;}
 
@@ -118,10 +128,15 @@ class SignalFilter {
 
  private:
     double m_bw;
-    int m_nPoints;
+    int m_nFilterPoints;
+    int m_nSignalPoints;
+    double m_signalLength;
     double m_xmin;
     double m_xmax;
     double* m_vecFilter;
+    double* m_vecFourierCosTable;
+    double* m_vecFourierSinTable;
+
     bool m_fail;
     string m_failMessage;
     string m_nameFilter;
@@ -133,12 +148,14 @@ class SignalFilter {
     double m_filterParam;
     int m_numIntegral;
 
-    static FilterID convertFilterNameToID (const char* filterName);
+    static const FilterID convertFilterNameToID (const char* filterName);
     static const char* convertFilterIDToName (const FilterID filterID);
+    static const FilterMethodID convertFilterMethodNameToID (const char* filterMethodName);
+    static const char* convertFilterMethodIDToName (const FilterMethodID filterMethodID);
     static const DomainID convertDomainNameToID (const char* domainName);
     static const char* convertDomainIDToName (const DomainID domainID);
 
-    void init (const FilterID filt_type, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numInt);
+    void init (const FilterID filt_type, const FilterMethodID filterMethod, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numInt);
 
 double spatialResponseCalc (double x, double param, int n) const;
 

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index d97ca2567f989d3018e42607241dc8c70e903854..bfb55e8dc4695ee5fb451c2847c3185d5e63990a 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.5 2000/06/30 02:03:27 kevin Exp $
+**  $Id: filter.cpp,v 1.6 2000/07/02 18:21:39 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
@@ -37,16 +37,18 @@
  *   int filt_type     Type of filter wanted
  *   double bw         Bandwidth of filter
  *   double xmin, xmax Filter limits
- *   int n             Number of points in filter
+ *   int n             Number of points in signal
  *   double param      General input parameter to filters
  *   int domain                FREQUENCY or SPATIAL domain wanted
  *   int numint                Number if intervals for calculating discrete inverse fourier xform
  *                     for spatial domain filters.  For ANALYTIC solutions, use numint = 0
  */
 
-SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numIntegral = 0)
+SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalLength, int n, double param, const char* domainName, int numIntegral = 0)
 {
   m_vecFilter = NULL;
+  m_vecFourierCosTable = NULL;
+  m_vecFourierSinTable = NULL;
   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
     m_fail = true;
@@ -54,6 +56,13 @@ SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, doub
     m_failMessage += filterName;
     return;
   }
+  m_idFilterMethod = convertFilterMethodNameToID (filterMethodName);
+  if (m_idFilterMethod == FILTER_METHOD_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid filter method name ";
+    m_failMessage += filterMethodName;
+    return;
+  }
   m_idDomain = convertDomainNameToID (domainName);
   if (m_idDomain == DOMAIN_INVALID) {
     m_fail = true;
@@ -61,12 +70,12 @@ SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, doub
     m_failMessage += domainName;
     return;
   }
-  init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numIntegral);
+  init (m_idFilter, m_idFilterMethod, bw, signalLength, n, param, m_idDomain, numIntegral);
 }
 
-SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numIntegral = 0)
+SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalLength, int n, double param, const DomainID domainID, int numIntegral = 0)
 {
-  init (filterID, bw, xmin, xmax, n, param, domainID, numIntegral);
+  init (filterID, filterMethodID, bw, signalLength, n, param, domainID, numIntegral);
 }
 
 SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0)
@@ -74,6 +83,8 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub
   m_bw = bw;
   m_nPoints = 0;
   m_vecFilter = NULL;
+  m_vecFourierCosTable = NULL;
+  m_vecFourierSinTable = NULL;
   m_filterParam = param;  
   m_numIntegral = numIntegral;
   m_idFilter = convertFilterNameToID (filterName);
@@ -93,24 +104,39 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub
 }
 
 void
-SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
+SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalLength, int n, double param, const DomainID domainID, int numint)
 {
   m_bw = bw;
   m_idFilter = filterID;
   m_idDomain = domainID;
-  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
+  m_idFilterMethod = filterMethodID;
+  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID) {
     m_fail = true;
     return;
   }
   m_nameFilter = convertFilterIDToName (m_idFilter);
   m_nameDomain = convertDomainIDToName (m_idDomain);
+  m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
   m_fail = false;
-  m_nPoints = n;
-  m_xmin = xmin;
-  m_xmax = xmax;
+  m_nSignalPoints = n;
+  m_nFilterPoints = 2 * m_nSignalPoints - 1;
+
+  m_signalLength = signalLength;
+  m_xmin = -signalLength;
+  m_xmax = signalLength;
   m_numIntegral = numint;
   m_filterParam = param;  
   m_vecFilter = new double[n];
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
+    int nFourier = n * n + 1;
+    double angleIncrement = (2. * PI) / n;
+    m_vecFourierCosTable = new double[ nFourier ];
+    m_vecFourierSinTable = new double[ nFourier ];
+    for (int i = 0; i < nFourier; i++) {
+      m_vecFourierCosTable[i] = cos (angleIncrement * i);
+      m_vecFourierSinTable[i] = sin (angleIncrement * i);
+    }
+  }
 
   double xinc = (m_xmax - m_xmin) / (m_nPoints - 1);
 
@@ -146,10 +172,12 @@ SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax
 SignalFilter::~SignalFilter (void)
 {
     delete m_vecFilter;
+    delete m_vecFourierSinTable;
+    delete m_vecFourierCosTable;
 }
 
 
-SignalFilter::FilterID
+const SignalFilter::FilterID
 SignalFilter::convertFilterNameToID (const char *filterName)
 {
   FilterID filterID = FILTER_INVALID;
@@ -207,6 +235,48 @@ SignalFilter::convertFilterIDToName (const FilterID filterID)
   return (name);
 }
       
+const SignalFilter::FilterMethodID
+SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
+{
+  FilterMethodID fmID = FILTER_METHOD_INVALID;
+
+  if (strcasecmp (filterMethodName, FILTER_METHOD_CONVOLUTION_STR) == 0)
+    fmID = FILTER_METHOD_CONVOLUTION;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FOURIER_STR) == 0)
+    fmID = FILTER_METHOD_FOURIER;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_STR) == 0)
+    fmID = FILTER_METHOD_FFT;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_ZEROPAD_2_STR) == 0)
+    fmID = FILTER_METHOD_FFT_ZEROPAD_2;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_ZEROPAD_4_STR) == 0)
+    fmID = FILTER_METHOD_FFT_ZEROPAD_4;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_ZEROPAD_6_STR) == 0)
+    fmID = FILTER_METHOD_FFT_ZEROPAD_6;
+
+  return (fmID);
+}
+
+const char *
+SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID)
+{
+  const char *name = "";
+
+  if (fmID == FILTER_METHOD_CONVOLUTION)
+    return (FILTER_METHOD_CONVOLUTION_STR);
+  else if (fmID == FILTER_METHOD_FOURIER)
+    return (FILTER_METHOD_FOURIER_STR);
+  else if (fmID == FILTER_METHOD_FFT)
+    return (FILTER_METHOD_FFT_STR);
+  else if (fmID == FILTER_METHOD_FFT_ZEROPAD_2)
+    return (FILTER_METHOD_FFT_ZEROPAD_2_STR);
+  else if (fmID == FILTER_METHOD_FFT_ZEROPAD_4)
+    return (FILTER_METHOD_FFT_ZEROPAD_4_STR);
+  else if (fmID == FILTER_METHOD_FFT_ZEROPAD_6)
+    return (FILTER_METHOD_FFT_ZEROPAD_6_STR);
+
+  return (name);
+}
+
 const SignalFilter::DomainID
 SignalFilter::convertDomainNameToID (const char* const domainName)
 {
@@ -234,6 +304,32 @@ SignalFilter::convertDomainIDToName (const DomainID domain)
 }
 
 
+void
+SignalFilter::filterSignal (const double input[], double output[], double dx, const int n) const
+{
+  if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
+    for (int i = 0; i < n; i++)
+      output[i] = convolve (input, dx, i, n);
+  } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
+    complex<double> fftSignal[n];
+    complex<double> complexOutput;
+    finiteFourierTransform (input, fftSignal, 1);
+    finiteFourierTransform (fftSignal, complexOutput, -1);
+    for (int i = 0; i < n; i++)
+      output[i] = complexOutput[i].mag();
+  }
+}
+
+void
+SignalFilter::filterSignal (const float input[], double output[], double dx, const int n) const
+{
+  if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
+    for (int i = 0; i < n; i++)
+      output[i] = convolve (input, dx, i, n);
+  }
+}
+
+
 double
 SignalFilter::response (double x)
 {
@@ -492,10 +588,9 @@ SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, d
 double 
 SignalFilter::integral_abscos (double u, double w)
 {
-  if (fabs (u) > F_EPSILON)
-    return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
-  else
-    return (w * w / 2);
+  return (fabs (u) > F_EPSILON 
+     ? (cos(u * w) - 1) / (u * u) + w / u * sin (u * w) 
+     : (w * w / 2));
 }
 
 
@@ -560,3 +655,57 @@ for (int i = 0; i < np; i++)
   return (sum * dx);
 }
 
+
+void
+SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
+{
+  if (direction < 0)
+    direction = -1;
+  else 
+    direction = 1;
+    
+  double angleIncrement = 2 * PI / n;
+  for (int i = 0; i < n; i++) {
+    double sumReal = 0;
+    double sumImag = 0;
+    for (int j = 0; j < n; j++) {
+      double angle = i * j * angleIncrement * direction;
+      sumReal += input[i] * cos(angle);
+      sumImag += input[i] * sin(angle);
+    }
+    if (direction > 0) {
+      sumReal /= n;
+      sumImag /= n;
+    }
+    output[i] = complex<double> (sumReal, sumImag);
+  }
+}
+
+void
+SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
+{
+  if (direction < 0)
+    direction = -1;
+  else 
+    direction = 1;
+    
+  double angleIncrement = 2 * PI / m_nPoints;
+  for (int i = 0; i < m_nPoints; i++) {
+    double sumReal = 0, sumImag = 0;
+    for (int j = 0; j < m_nPoints; j++) {
+      int tableIndex = i * j;
+      if (direction > 0) {
+       sumReal += input[i] * m_vecFourierCosTable[tableIndex];
+       sumImag += input[i] * m_vecFourierSinTable[tableIndex];
+      } else {
+       sumReal += input[i] * m_vecFourierCosTable[tableIndex];
+       sumImag -= input[i] * m_vecFourierSinTable[tableIndex];
+      }
+    }
+    if (direction > 0) {
+      sumReal /= m_nPoints;
+      sumImag /= m_nPoints;
+    }
+    output[i] = complex<double> (sumReal, sumImag);
+  }
+}

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 3a071b478bceb2f25761068260791512d3f32ca7..6e21b3568285ded1049dccf2d2ca9f7052e26a8f 100644 (file)
--- a/libctsim/projections.cpp
+++ b/libctsim/projections.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.7 2000/06/29 12:39:46 kevin Exp $
+**  $Id: projections.cpp,v 1.8 2000/07/02 18:21:39 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
@@ -490,10 +490,10 @@ bool
 Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const char* const interpName, int interp_param, const char* const backprojectName, const int trace)
 {
   int nview = m_nView;
-  double det_inc = m_detInc;
-  double detlen = (m_nDet - 1) * det_inc;
-  int n_filtered_proj = m_nDet;
-  double filtered_proj [n_filtered_proj];   // convolved result
+  double detInc = m_detInc;
+  double detlen = (m_nDet - 1) * detInc;
+  int n_filteredProj = m_nDet;
+  double filteredProj [n_filteredProj];   // filtered projections
 
 #ifdef HAVE_BSPLINE_INTERP
   int spline_order = 0, zoom_factor = 0;
@@ -501,16 +501,13 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     zoom_factor = interp_param;
     spline_order = 3;
     zoom_factor = 3;
-    n_filtered_proj = (m_nDet - 1) * (zoom_factor + 1) + 1;
+    n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1;
   }
 #endif
 
-  int n_vec_filter = 2 * m_nDet - 1;
-  double filterBW = 1. / det_inc;
-  double filterMin = -detlen;
-  double filterMax = detlen;
+  double filterBW = 1. / detInc;
 
-  SignalFilter filter (filterName, filterBW, filterMin, filterMax, n_vec_filter, filt_param, "spatial", 0);
+  SignalFilter filter (filterName, filterMethodName, filterBW, detlen, m_nDet, filt_param, "spatial", 0);
   if (filter.fail()) {
       sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", filter.failMessage().c_str());
       return false;
@@ -521,13 +518,14 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 
 #if HAVE_SGP
   SGP_ID gid;
+  int n_vec_filter = filter.getNFilterPoints();
   double plot_xaxis [n_vec_filter];                    // array for plotting 
 
   if (trace > TRACE_TEXT)  {
     int i;
     double f;
-    double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
-    for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
+    double filterInc = (detlen * 2) / (n_vec_filter - 1);
+    for (i = 0, f = -detlen; i < n_vec_filter; i++, f += filterInc)
       plot_xaxis[i] = f;
       
     gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
@@ -536,7 +534,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     sgp2_close (gid);
   }
   if (trace >= TRACE_TEXT) {
-    printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
+    printf ("nview=%d, ndet=%d, det_start=%.4f, detInc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
   }
 #endif  //HAVE_SGP
 
@@ -553,8 +551,10 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
     DetectorArray& darray = getDetectorArray (iview);
     DetectorValue* detval = darray.detValues();
 
-    for (int j = 0; j < m_nDet; j++)
-      filtered_proj[j] = filter.convolve (detval, det_inc, j, m_nDet);
+    filter.filterSignal (detval, filteredProj, detInc, m_nDet);
+
+    //    for (int j = 0; j < m_nDet; j++)
+    //      filteredProj[j] = filter.convolve (detval, detInc, j, m_nDet);
 
 #ifdef HAVE_SGP
     if (trace >= TRACE_PLOT)  {
@@ -576,23 +576,23 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
       ezset  ("xlength .5.");
       ezset ("box");
       ezset ("grid");
-      gid = ezplot (filtered_proj, plot_xaxis, m_nDet);
+      gid = ezplot (filteredProj, plot_xaxis, m_nDet);
     }
 #endif  //HAVE_SGP
 
 #ifdef HAVE_BSPLINE_INTERP
     if (interp_type == I_BSPLINE) 
-       bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+       bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
     
 #ifdef HAVE_SGP
     if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
-       bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
-      ezplot_1d (filtered_proj, n_filtered_proj);
+       bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
+      ezplot_1d (filteredProj, n_filteredProj);
     }
 #endif
 #endif
 
-    bj.BackprojectView (filtered_proj, darray.viewAngle());
+    bj.BackprojectView (filteredProj, darray.viewAngle());
 
 #ifdef HAVE_SGP
     if (trace >= TRACE_PLOT) {