r136: *** empty log message ***

diff --git a/ChangeLog b/ChangeLog
index c32573bb176ad6e0f2e1c0d8883498aedf5112a7..108313613f470aa8fd57f26230f6313376760556 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2.0.0-b2 - 7/07/00
+   Added zeropadding option to pjrec
+   Cleaned up SignalFilter class
+   Added zeropad options to html and cgi files
+       
 2.0.0-b1 - 7/05/00
    Updated trace level processing
    Added trace level to SignalFilter

diff --git a/cgi-bin/ctsim.cgi.in b/cgi-bin/ctsim.cgi.in
index 54a38a62d8bba21f3a6dac98fef1b8ba5aec81c2..8c9190b003ed0bd8eed36b227c7973f40df98f2a 100755 (executable)
--- a/cgi-bin/ctsim.cgi.in
+++ b/cgi-bin/ctsim.cgi.in
@@ -48,6 +48,7 @@ $error .= "Projection RotAngle must be between 0.1 and 2<br>" if ($PJ_RotAngle <
 my $IR_Nx = $Phantom_Nx;
 my $IR_Ny = $Phantom_Ny;
 my $IR_FilterMethod = FilterMetaChars($in{'IR_FilterMethod'});
+my $IR_Zeropad = FilterToNumber($in{'IR_Zeropad'});
 my $IR_Filter = FilterMetaChars($in{'IR_Filter'});
 my $IR_Filter_Param = FilterToNumber($in{'IR_Filter_Param'});
 my $IR_Interp = FilterMetaChars($in{'IR_Interp'});
@@ -108,7 +109,7 @@ $phm2if_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2if-lam" if $MPI;
 my $gp_cmd = "$phm2if_ver $phantom_fname $Phantom_Nx $Phantom_Ny --phantom $Phantom_Name --nsample $Phantom_NSample";
 my $pj_cmd = "$phm2pj_ver $pj_fname $PJ_NDet $PJ_NRot --phantom $Phantom_Name --nray $PJ_NRay --rotangle $PJ_RotAngle";
 my $pj_if_cmd = "$::bindir/pj2if $pj_fname $pj_if_fname";
-my $pjrec_cmd = "$pjrec_ver $pj_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj --filter-method $IR_FilterMethod";
+my $pjrec_cmd = "$pjrec_ver $pj_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj --filter-method $IR_FilterMethod --zeropad $IR_Zeropad";
 my $diff_cmd = "$diff_ver $phantom_fname $ir_fname $diff_fname --comp";
 my $compare_cmd = "$ifinfo_ver $phantom_fname $ir_fname";
 

diff --git a/configure b/configure
index 4ec38341a2c31ce34abb726528004c0699023e3f..960133d5b29287ab308d64a255e20bc953c77082 100755 (executable)
--- a/configure
+++ b/configure
@@ -710,7 +710,7 @@ fi
 
 PACKAGE=ctsim
 
-VERSION=2.0.0
+VERSION=2.0.0-b1
 
 if test "`cd $srcdir && pwd`" != "`pwd`" && test -f $srcdir/config.status; then
   { echo "configure: error: source directory already configured; run "make distclean" there first" 1>&2; exit 1; }

diff --git a/html/simulate.html.in b/html/simulate.html.in
index 632b4ef93427a505c01271fcf2ee7db1a3245c42..1af20f036c8c0270140f41b1ead3f4616cf6e30b 100644 (file)
--- a/html/simulate.html.in
+++ b/html/simulate.html.in
@@ -37,7 +37,11 @@ Rotation Angle<br>as a multiple of PI: <input type="text" name="PJ_RotAngle" siz
 Filter Method:<br>
 <input type="radio" name="IR_FilterMethod" value="convolution" checked>Convolution<br>
 <input type="radio" name="IR_FilterMethod" value="fourier">Fourier<br>
-<input type="radio" name="IR_FilterMethod" value="fft">FFT (Fast Fourier<br>
+<input type="radio" name="IR_FilterMethod" value="fft">FFT (Fast Fourier)<br>
+<input type="radio" name="IR_FilterMethod" value="fftw">FFTW (Fastest Fourier Library)<br>
+<p>
+Zeropad: (frequency-base filtering only)<br>
+<input type="text" name="IR_Zeropad" size="1" value="2"><br>
 <p>
 Interpolation:<br>
 <input type="radio" name="IR_Interp" value="linear" checked>Linear<br>

diff --git a/include/filter.h b/include/filter.h
index 8e0740b982e1bb34dc4592d090fb309c9b8e806c..7040a3bcabd2742ca0e94f0cfea2ed18cb6800dd 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.9 2000/07/04 22:21:01 kevin Exp $
+**  $Id: filter.h,v 1.10 2000/07/06 08:30:30 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
@@ -50,9 +50,7 @@ class SignalFilter {
        FILTER_METHOD_CONVOLUTION,
        FILTER_METHOD_FOURIER,
        FILTER_METHOD_FFT,
-       FILTER_METHOD_FFT_ZEROPAD_2,
-       FILTER_METHOD_FFT_ZEROPAD_4,
-       FILTER_METHOD_FFT_ZEROPAD_6
+       FILTER_METHOD_FFTW,
     } FilterMethodID;
 
     typedef enum {
@@ -75,17 +73,15 @@ class SignalFilter {
     static const char FILTER_METHOD_CONVOLUTION_STR[]=  "convolution";
     static const char FILTER_METHOD_FOURIER_STR[]=      "fourier";
     static const char FILTER_METHOD_FFT_STR[]=          "fft";
-    static const char FILTER_METHOD_FFT_ZEROPAD_2_STR[]="fft_zeropad2";
-    static const char FILTER_METHOD_FFT_ZEROPAD_4_STR[]="fft_zeropad4";
-    static const char FILTER_METHOD_FFT_ZEROPAD_6_STR[]="fft_zeropad6";
+    static const char FILTER_METHOD_FFTW_STR[]=         "fftw";
 
-    static const char DOMAIN_FREQUENCY_STR[]=    "frequency";
-    static const char DOMAIN_SPATIAL_STR[]= "spatial";
+    static const char DOMAIN_FREQUENCY_STR[]="frequency";
+    static const char DOMAIN_SPATIAL_STR[]="spatial";
 
 
-    SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalIncrement, int n, double param, const char* domainName, const int numIntegral = 0);
+    SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalIncrement, int n, double param, const char* domainName, const int zeropad = 0, const int numIntegral = 0);
 
-    SignalFilter (const FilterID filt_type, FilterMethodID filterMethodID, double bw, double signalIncrement, int n, double param, const DomainID domain, const int numIntegral = 0);
+    SignalFilter (const FilterID filt_type, FilterMethodID filterMethodID, double bw, double signalIncrement, int n, double param, const DomainID domain, const int zeropad = 0, const int numIntegral = 0);
 
     SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0);
 
@@ -104,11 +100,11 @@ class SignalFilter {
     void filterSignal (const double input[], double output[]) const;
     void filterSignal (const float input[], double output[]) const;
 
-    static void finiteFourierTransform (const float input[], complex<double> output[], const int n, const int direction);
+    static void finiteFourierTransform (const double input[], complex<double> output[], const int n, const int direction);
 
     static void finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, const int direction);
 
-    void finiteFourierTransform (const float input[], complex<double> output[], const int direction) const;
+    void finiteFourierTransform (const double input[], complex<double> output[], const int direction) const;
 
     void finiteFourierTransform (const complex<double> input[], complex<double> output[], const int direction) const;
 
@@ -151,6 +147,9 @@ class SignalFilter {
     complex<double>* m_complexVecFilter;
 #ifdef HAVE_FFTW
     fftw_plan m_planForward, m_planBackward;
+    fftw_complex* m_vecFftInput;
+#else
+    complex<double>* m_vecFftInput;
 #endif
 
     bool m_fail;
@@ -164,6 +163,7 @@ class SignalFilter {
     double m_filterParam;
     int m_numIntegral;
     int m_traceLevel;
+    int m_zeropad;
 
     static const FilterID convertFilterNameToID (const char* filterName);
     static const char* convertFilterIDToName (const FilterID filterID);
@@ -172,7 +172,7 @@ class SignalFilter {
     static const DomainID convertDomainNameToID (const char* domainName);
     static const char* convertDomainIDToName (const DomainID domainID);
 
-    void init (const FilterID filt_type, const FilterMethodID filterMethod, double bw, double signalIncrement, int n, double param, const DomainID domain, const int numInt);
+    void init (const FilterID filt_type, const FilterMethodID filterMethod, double bw, double signalIncrement, int n, double param, const DomainID domain, const int zeropad, const int numInt);
 
 double spatialResponseCalc (double x, double param, int n) const;
 

diff --git a/include/projections.h b/include/projections.h
index 1bbeafc0c8adebf234c223e587da96e2275d5354..f235e320ee578b198b387cb6c3ec68e7c59ab00b 100644 (file)
--- a/include/projections.h
+++ b/include/projections.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.h,v 1.7 2000/06/29 12:39:45 kevin Exp $
+**  $Id: projections.h,v 1.8 2000/07/06 08:30:30 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
@@ -53,7 +53,7 @@ class Projections
   bool detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int view_num);
   bool detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int view_num);
 
-  bool reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const char* const interpName, int interp_param, const char* const backprojName, const int trace);
+  bool reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* const interpName, int interp_param, const char* const backprojName, const int trace);
 
   void setNView (int nView);  // used in MPI to restrict # of views
   void setRotInc (double rotInc) { m_rotInc = rotInc;}

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 978a5025777a4ed4166e223242a01b853a31b0ef..0a24996805bba3a24f48bdb83c1bd1900c9b5e04 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.12 2000/07/05 17:59:26 kevin Exp $
+**  $Id: filter.cpp,v 1.13 2000/07/06 08:30:30 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
@@ -44,11 +44,12 @@
  *                     for spatial domain filters.  For ANALYTIC solutions, use numint = 0
  */
 
-SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, int numIntegral = 0)
+SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, int zeropad = 0, int numIntegral = 0)
 {
   m_vecFilter = NULL;
   m_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;
+  m_vecFftInput = NULL;
   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
     m_fail = true;
@@ -70,12 +71,12 @@ SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName
     m_failMessage += domainName;
     return;
   }
-  init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, numIntegral);
+  init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, zeropad, numIntegral);
 }
 
-SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int numIntegral = 0)
+SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int zeropad = 0, int numIntegral = 0)
 {
-  init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, numIntegral);
+  init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, zeropad, numIntegral);
 }
 
 SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0)
@@ -86,6 +87,7 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub
   m_vecFilter = NULL;
   m_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;
+  m_vecFftInput = NULL;
   m_filterParam = param;  
   m_numIntegral = numIntegral;
   m_idFilter = convertFilterNameToID (filterName);
@@ -105,7 +107,7 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub
 }
 
 void
-SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int numint)
+SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int zeropad, int numint)
 {
   m_bw = bw;
   m_idFilter = filterID;
@@ -123,37 +125,24 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   m_nSignalPoints = nSignalPoints;
   m_signalInc = signalIncrement;
   m_filterParam = param;  
-  
-  if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
-    int nFourier = m_nSignalPoints * m_nSignalPoints + 1;
-    double angleIncrement = (2. * PI) / m_nSignalPoints;
-    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);
-    }
-    m_nFilterPoints = m_nSignalPoints;
-    m_filterMin = -1. / (2 * m_signalInc);
-    m_filterMax = 1. / (2 * m_signalInc);
-    m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints;
-    m_vecFilter = new double [m_nFilterPoints];
-    int halfFilter = m_nFilterPoints / 2;
-    for (int i = 0; i < halfFilter; i++) 
-       m_vecFilter[i] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
-    for (int i = 0; i < halfFilter; i++)
-       m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i+1) / (halfFilter - 1) / (2 * m_signalInc);
-    if (halfFilter % 2) // odd
-      m_vecFilter[halfFilter] = 1 / (2 * m_signalInc);
-  } else if (m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
+  m_zeropad = zeropad;
+
+  m_vecFourierCosTable = NULL;
+  m_vecFourierSinTable = NULL;
+  m_vecFilter = NULL;
+  m_vecFftInput = NULL;
+
+  if (m_idFilterMethod == FILTER_METHOD_FFT)
+    m_idFilterMethod = FILTER_METHOD_FFTW;
+
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER || m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFTW) {
     m_nFilterPoints = m_nSignalPoints;
-    if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
+    if (m_zeropad > 0) {
       double logBase2 = log(m_nSignalPoints) / log(2);
-      int nextPowerOf2 = static_cast<int>(floor(logBase2)) + 1;
-      if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4)
-       nextPowerOf2++;
+      int nextPowerOf2 = static_cast<int>(floor(logBase2));
       if (logBase2 != floor(logBase2))
        nextPowerOf2++;
+      nextPowerOf2 += m_zeropad;
       m_nFilterPoints = 1 << nextPowerOf2;
       cout << "nFilterPoints = " << m_nFilterPoints << endl;
     }
@@ -162,18 +151,38 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
     m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints;
     m_vecFilter = new double [m_nFilterPoints];
     int halfFilter = m_nFilterPoints / 2;
-    for (int i = 0; i < halfFilter; i++) 
-       m_vecFilter[i] = static_cast<double>(i) / (halfFilter - 1) /  (2 * m_signalInc) / m_nSignalPoints;
-    for (int i = 0; i < halfFilter; i++)
-       m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i+1) / (halfFilter - 1) /  (2 * m_signalInc) / m_nSignalPoints;
-    if (halfFilter % 2) // odd
-      m_vecFilter[halfFilter] = 1 / (2 * m_signalInc) / m_nSignalPoints;
+    for (int i = 0; i <= halfFilter; i++) 
+       m_vecFilter[i] = static_cast<double>(i) / halfFilter/ (2. * m_signalInc);
+    for (int i = 1; i <= halfFilter; i++)
+       m_vecFilter[m_nFilterPoints - i] = static_cast<double>(i) / halfFilter / (2. * m_signalInc);
+  }
+
+  // precalculate sin and cosine tables for fourier transform
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
+    int nFourier = m_nFilterPoints * m_nFilterPoints + 1;
+    double angleIncrement = (2. * PI) / m_nFilterPoints;
+    m_vecFourierCosTable = new double[ nFourier ];
+    m_vecFourierSinTable = new double[ nFourier ];
+    double angle = 0;
+    for (int i = 0; i < nFourier; i++) {
+      m_vecFourierCosTable[i] = cos (angle);
+      m_vecFourierSinTable[i] = sin (angle);
+      angle += angleIncrement;
+    }
+  }
 
 #if HAVE_FFTW
+  if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+      for (int i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
+       m_vecFilter[i] /= m_nFilterPoints;
+
     m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
     m_planBackward = fftw_create_plan (m_nFilterPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
-#endif
+    m_vecFftInput = new fftw_complex [ m_nFilterPoints ];
+    for (int i = 0; i < m_nFilterPoints; i++) 
+       m_vecFftInput[i].re = m_vecFftInput[i].im = 0;
   }
+#endif
 
  if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
     m_nFilterPoints = 2 * m_nSignalPoints - 1;
@@ -215,11 +224,12 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
 
 SignalFilter::~SignalFilter (void)
 {
-    delete m_vecFilter;
-    delete m_vecFourierSinTable;
-    delete m_vecFourierCosTable;
+    delete [] m_vecFilter;
+    delete [] m_vecFourierSinTable;
+    delete [] m_vecFourierCosTable;
+    delete [] m_vecFftInput;
 #if HAVE_FFTW
-    if (m_idFilterMethod == FILTER_METHOD_FFT) {
+    if (m_idFilterMethod == FILTER_METHOD_FFTW) {
        fftw_destroy_plan(m_planForward);
        fftw_destroy_plan(m_planBackward);
     }
@@ -296,12 +306,8 @@ SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
     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;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFTW_STR) == 0)
+    fmID = FILTER_METHOD_FFTW;
 
   return (fmID);
 }
@@ -317,12 +323,8 @@ SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID)
     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);
+  else if (fmID == FILTER_METHOD_FFTW)
+    return (FILTER_METHOD_FFTW_STR);
 
   return (name);
 }
@@ -361,32 +363,35 @@ SignalFilter::filterSignal (const float input[], double output[]) const
     for (int i = 0; i < m_nSignalPoints; i++)
       output[i] = convolve (input, m_signalInc, i, m_nSignalPoints);
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
-    complex<double> fftSignal[m_nSignalPoints];
-    complex<double> complexOutput[m_nSignalPoints];
-    complex<double> filteredSignal[m_nSignalPoints];
-    finiteFourierTransform (input, fftSignal, m_nSignalPoints, -1);
-    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nSignalPoints);
-    finiteFourierTransform (filteredSignal, complexOutput, m_nSignalPoints, 1);
+    complex<double> fftSignal[m_nFilterPoints];
+    complex<double> complexOutput[m_nFilterPoints];
+    complex<double> filteredSignal[m_nFilterPoints];
+    double inputSignal[m_nFilterPoints];
+    for (int i = 0; i < m_nSignalPoints; i++)
+      inputSignal[i] = input[i];
+    for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
+      inputSignal[i] = 0;  // zeropad
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
+    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
+    finiteFourierTransform (filteredSignal, complexOutput, m_nFilterPoints, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
       output[i] = complexOutput[i].real();
-  } else if (m_idFilterMethod == FILTER_METHOD_FFT || FILTER_METHOD_FFT_ZEROPAD_2 || FILTER_METHOD_FFT_ZEROPAD_4) {
+  }
 #if HAVE_FFTW
-    fftw_complex in[m_nFilterPoints], out[m_nFilterPoints];
-    for (int i = 0; i < m_nSignalPoints; i++) {
-      in[i].re = input[i];
-      in[i].im = 0;
-    }
-    for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) {
-      in[i].re = in[i].im = 0;      // ZeroPad
-    }
-    fftw_one(m_planForward, in, out);
-    for (int i = 0; i < m_nFilterPoints; i++) {
-      out[i].re = m_vecFilter[i] * out[i].re;
-      out[i].im = m_vecFilter[i] * out[i].im;
-    }
-    fftw_one(m_planBackward, out, in);
-    for (int i = 0; i < m_nSignalPoints; i++) 
-      output[i] = in[i].re;
+  else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+      for (int i = 0; i < m_nSignalPoints; i++)
+         m_vecFftInput[i].re = input[i];
+
+      fftw_complex out[m_nFilterPoints];
+      fftw_one(m_planForward, m_vecFftInput, out);
+      for (int i = 0; i < m_nFilterPoints; i++) {
+         out[i].re = m_vecFilter[i] * out[i].re;
+         out[i].im = m_vecFilter[i] * out[i].im;
+      }
+      fftw_complex outFiltered[m_nFilterPoints];
+      fftw_one(m_planBackward, out, outFiltered);
+      for (int i = 0; i < m_nSignalPoints; i++) 
+         output[i] = outFiltered[i].re;
   }
 #endif
 }
@@ -718,19 +723,19 @@ for (int i = 0; i < np; i++)
 
 
 void
-SignalFilter::finiteFourierTransform (const float input[], complex<double> output[], const int n, int direction)
+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;
+  double angleIncrement = direction * 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;
+      double angle = i * j * angleIncrement;
       sumReal += input[j] * cos(angle);
       sumImag += input[j] * sin(angle);
     }
@@ -767,16 +772,16 @@ SignalFilter::finiteFourierTransform (const complex<double> input[], complex<dou
 }
 
 void
-SignalFilter::finiteFourierTransform (const float input[], complex<double> output[], int direction) const
+SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
 {
   if (direction < 0)
     direction = -1;
   else 
     direction = 1;
     
-  for (int i = 0; i < m_nSignalPoints; i++) {
+  for (int i = 0; i < m_nFilterPoints; i++) {
     double sumReal = 0, sumImag = 0;
-    for (int j = 0; j < m_nSignalPoints; j++) {
+    for (int j = 0; j < m_nFilterPoints; j++) {
       int tableIndex = i * j;
       if (direction > 0) {
        sumReal += input[i] * m_vecFourierCosTable[tableIndex];
@@ -787,8 +792,8 @@ SignalFilter::finiteFourierTransform (const float input[], complex<double> outpu
       }
     }
     if (direction < 0) {
-      sumReal /= m_nSignalPoints;
-      sumImag /= m_nSignalPoints;
+      sumReal /= m_nFilterPoints;
+      sumImag /= m_nFilterPoints;
     }
     output[i] = complex<double> (sumReal, sumImag);
   }

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 0182eec8d4af2ff6bae3ec189f2a3ad942974c6f..db6431411cff6a435d78a2dfb31a9407a7a1f95b 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.11 2000/07/04 22:21:01 kevin Exp $
+**  $Id: projections.cpp,v 1.12 2000/07/06 08:30:30 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
@@ -487,7 +487,7 @@ Projections::printScanInfo (void) const
  */
 
 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)
+Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* const interpName, int interp_param, const char* const backprojectName, const int trace)
 {
   int nview = m_nView;
   double detInc = m_detInc;
@@ -505,7 +505,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 #endif
 
   double filterBW = 1. / detInc;
-  SignalFilter filter (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", 0);
+  SignalFilter filter (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", zeropad);
   filter.setTraceLevel(trace);
 
   if (filter.fail()) {

diff --git a/src/pjrec.cpp b/src/pjrec.cpp
index 782563cb1dd76d40a48bf0003b6427ae9fc1e4b9..c22b472eb7212d69d9c57439d9d0ec050d478e76 100644 (file)
--- a/src/pjrec.cpp
+++ b/src/pjrec.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pjrec.cpp,v 1.6 2000/07/04 22:21:01 kevin Exp $
+**  $Id: pjrec.cpp,v 1.7 2000/07/06 08:30:30 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
@@ -29,13 +29,14 @@
 #include "timer.h"
 
 
-enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
 
 static struct option my_options[] =
 {
   {"interp", 1, 0, O_INTERP},
   {"filter", 1, 0, O_FILTER},
   {"filter-method", 1, 0, O_FILTER_METHOD},
+  {"zeropad", 1, 0, O_ZEROPAD},
   {"filter-param", 1, 0, O_FILTER_PARAM},
   {"backproj", 1, 0, O_BACKPROJ},
   {"trace", 1, 0, O_TRACE},
@@ -53,48 +54,54 @@ pjrec_usage (const char *program)
   cout << "usage: " << fileBasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl;
   cout << "Image reconstruction from raysum projections" << endl;
   cout << endl;
-  cout << "   raysum-file     Input raysum file" << endl;
-  cout << "   image-file      Output image file in SDF2D format" << endl;
-  cout << "   nx-image        Number of columns in output image" << endl;
-  cout << "   ny-image        Number of rows in output image" << endl;
-  cout << "   --interp        Interpolation method during backprojection" << endl;
-  cout << "       nearest     Nearest neighbor interpolation" << endl;
-  cout << "       linear      Linear interpolation" << endl;
+  cout << "  raysum-file     Input raysum file" << endl;
+  cout << "  image-file      Output image file in SDF2D format" << endl;
+  cout << "  nx-image        Number of columns in output image" << endl;
+  cout << "  ny-image        Number of rows in output image" << endl;
+  cout << "  --interp        Interpolation method during backprojection" << endl;
+  cout << "    nearest     Nearest neighbor interpolation" << endl;
+  cout << "    linear      Linear interpolation" << endl;
 #if HAVE_BSPLINE_INTERP
-  cout << "       bspline     B-spline interpolation" << endl;
+  cout << "    bspline     B-spline interpolation" << endl;
 #endif
-  cout << "    --filter       Filter name" << endl;
-  cout << "       abs_bandlimit Abs * Bandlimiting (default)" << endl;
-  cout << "       abs_sinc      Abs * Sinc" << endl;
-  cout << "       abs_cos       Abs * Cosine" << endl;
-  cout << "       abs_hamming   Abs * Hamming" << endl;
-  cout << "       shepp         Shepp-Logan" << endl;
-  cout << "       bandlimit     Bandlimiting" << endl;
-  cout << "       sinc          Sinc" << endl;
-  cout << "       cos           Cosine" << endl;
-  cout << "       triangle      Triangle" << endl;
-  cout << "       hamming       Hamming" << endl;
-  cout << "    --filter-method  Filter method before backprojections\n";;
-  cout << "       convolution      Spatial filtering (default)\n";
-  cout << "       fourier          Frequency filtering with discete fourier\n";
-  cout << "    --backproj     Backprojection Method" << endl;
-  cout << "       trig        Trigometric functions at every point" << endl;
-  cout << "       table       Trigometric functions with precalculated table" << endl;
-  cout << "       diff        Difference method" << endl;
-  cout << "       diff2       Optimized difference method (default)" << endl;
-  cout << "       idiff2      Optimized difference method with integer math" << endl;
-  cout << "    --filter-param Alpha level for Hamming filter" << endl;
-  cout << "    --trace        Set tracing to level" << endl;
-  cout << "         none      No tracing (default)" << endl;
-  cout << "         text      Text level tracing" << endl;
-  cout << "         phm       Trace phantom" << endl;
-  cout << "         rays      Trace allrays" << endl;
-  cout << "         plot      Trace plotting" << endl;
-  cout << "         clipping  Trace clipping" << endl;
-  cout << "    --verbose      Turn on verbose mode" << endl;
-  cout << "    --debug        Turn on debug mode" << endl;
-  cout << "    --version      Print version" << endl;
-  cout << "    --help         Print this help message" << endl;
+  cout << "  --filter       Filter name" << endl;
+  cout << "    abs_bandlimit Abs * Bandlimiting (default)" << endl;
+  cout << "    abs_sinc      Abs * Sinc" << endl;
+  cout << "    abs_cos       Abs * Cosine" << endl;
+  cout << "    abs_hamming   Abs * Hamming" << endl;
+  cout << "    shepp         Shepp-Logan" << endl;
+  cout << "    bandlimit     Bandlimiting" << endl;
+  cout << "    sinc          Sinc" << endl;
+  cout << "    cos           Cosine" << endl;
+  cout << "    triangle      Triangle" << endl;
+  cout << "    hamming       Hamming" << endl;
+  cout << "  --filter-method  Filter method before backprojections\n";;
+  cout << "    convolution      Spatial filtering (default)\n";
+  cout << "    fourier          Frequency filtering with discete fourier\n";
+  cout << "    fft              Fast Fourier Transform\n";
+#if HAVE_FFTW
+  cout << "    fftw             Fast Fourier Transform in the West library\n";
+#endif
+  cout << "  --zeropad n   Set zeropad level (default = 0)\n";
+  cout << "                set n to number of powers to two to pad\n";
+  cout << "  --backproj    Backprojection Method" << endl;
+  cout << "    trig        Trigometric functions at every point" << endl;
+  cout << "    table       Trigometric functions with precalculated table" << endl;
+  cout << "    diff        Difference method" << endl;
+  cout << "    diff2       Optimized difference method (default)" << endl;
+  cout << "    idiff2      Optimized difference method with integer math" << endl;
+  cout << "  --filter-param Alpha level for Hamming filter" << endl;
+  cout << "  --trace        Set tracing to level" << endl;
+  cout << "     none      No tracing (default)" << endl;
+  cout << "     text      Text level tracing" << endl;
+  cout << "     phm       Trace phantom" << endl;
+  cout << "     rays      Trace allrays" << endl;
+  cout << "     plot      Trace plotting" << endl;
+  cout << "     clipping  Trace clipping" << endl;
+  cout << "  --verbose      Turn on verbose mode" << endl;
+  cout << "  --debug        Turn on debug mode" << endl;
+  cout << "  --version      Print version" << endl;
+  cout << "  --help         Print this help message" << endl;
 }
 
 
@@ -115,6 +122,7 @@ pjrec_main (int argc, char * argv[])
   char *endptr;
   int optVerbose = 0;
   int optDebug = 0;
+  int optZeroPad = 0;
   int optTrace = TRACE_NONE;
   double optFilterParam = -1;
   string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR;
@@ -163,6 +171,12 @@ pjrec_main (int argc, char * argv[])
            pjrec_usage(argv[0]);
          }
          break;
+       case O_ZEROPAD:
+         optZeroPad = strtol(optarg, &endptr, 10);
+         if (endptr != optarg + strlen(optarg)) {
+           pjrec_usage(argv[0]);
+         }
+         break;
        case O_VERBOSE:
          optVerbose = 1;
          break;
@@ -236,11 +250,13 @@ pjrec_main (int argc, char * argv[])
   TimerCollectiveMPI timerBcast (mpiWorld.getComm());
   mpiWorld.BcastString (optBackprojName);
   mpiWorld.BcastString (optFilterName);
+  mpiWorld.BcastString (optFilterMethodName);
   mpiWorld.BcastString (optInterpName);
   mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
+  mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&optInterpParam, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0);
   mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
@@ -279,7 +295,7 @@ pjrec_main (int argc, char * argv[])
 
 #ifdef HAVE_MPI
   TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
-  projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
+  projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
   if (optVerbose)
       timerReconstruct.timerEndAndReport ("Time to reconstruct");
 
@@ -288,7 +304,7 @@ pjrec_main (int argc, char * argv[])
   if (optVerbose)
       timerReduce.timerEndAndReport ("Time to reduce image");
 #else
-  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
+  projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace);
 #endif
 
 #ifdef HAVE_MPI