r137: *** empty log message ***

diff --git a/ChangeLog b/ChangeLog
index 108313613f470aa8fd57f26230f6313376760556..6cd4d1af0a6fd67c9b9ae3a33635143a66b20dc1 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -2,6 +2,8 @@
    Added zeropadding option to pjrec
    Cleaned up SignalFilter class
    Added zeropad options to html and cgi files
+   Added fourier_table filter method
+   Added FFTW routines to use real/half-complex transformations
        
 2.0.0-b1 - 7/05/00
    Updated trace level processing

diff --git a/configure b/configure
index 960133d5b29287ab308d64a255e20bc953c77082..889edc6163ba8fd8c12947d0bc3d656083ec479e 100755 (executable)
--- a/configure
+++ b/configure
@@ -710,7 +710,7 @@ fi
 
 PACKAGE=ctsim
 
-VERSION=2.0.0-b1
+VERSION=2.0.0-b2
 
 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/configure.in b/configure.in
index 63efb05ff93c00e6ef401841dc641377ebf8324d..121694d051f38d92e26fb4c24f88efc391d3218e 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-b1)
+AM_INIT_AUTOMAKE(ctsim,2.0.0-b2)
 AM_CONFIG_HEADER(config.h)
 
 dnl Checks for programs.

diff --git a/html/simulate.html.in b/html/simulate.html.in
index 1af20f036c8c0270140f41b1ead3f4616cf6e30b..c05887daee491cbd08442e5a23f9575144692490 100644 (file)
--- a/html/simulate.html.in
+++ b/html/simulate.html.in
@@ -37,8 +37,10 @@ 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="fourier_table">Fourier Table<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>
+<input type="radio" name="IR_FilterMethod" value="rfftw">RFFTW (Real-mode Fastest Fourier Library)<br>
 <p>
 Zeropad: (frequency-base filtering only)<br>
 <input type="text" name="IR_Zeropad" size="1" value="2"><br>

diff --git a/include/ct.h b/include/ct.h
index 5cbb36134ae1b98c76d41833a588fe5d0b1e4a1b..5a0dc9b5a7b8a6ed69d4d58eee6de644d4ecb3d6 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.25 2000/07/04 18:33:35 kevin Exp $
+**  $Id: ct.h,v 1.26 2000/07/06 18:37:24 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
@@ -110,6 +110,7 @@ extern "C" {
 #endif
 
 #ifdef HAVE_FFTW
+#include <rfftw.h>
 #include <fftw.h>
 #endif
 

diff --git a/include/filter.h b/include/filter.h
index 7040a3bcabd2742ca0e94f0cfea2ed18cb6800dd..e4c4f3f9db30df46d6da0077fe6b2b381a4333cf 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.10 2000/07/06 08:30:30 kevin Exp $
+**  $Id: filter.h,v 1.11 2000/07/06 18:37:24 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
@@ -28,6 +28,7 @@
 #ifndef FILTER_H
 #define FILTER_H
 
+
 class SignalFilter {
  public:
 
@@ -49,8 +50,10 @@ class SignalFilter {
        FILTER_METHOD_INVALID,
        FILTER_METHOD_CONVOLUTION,
        FILTER_METHOD_FOURIER,
+       FILTER_METHOD_FOURIER_TABLE,
        FILTER_METHOD_FFT,
        FILTER_METHOD_FFTW,
+       FILTER_METHOD_RFFTW
     } FilterMethodID;
 
     typedef enum {
@@ -72,8 +75,10 @@ class SignalFilter {
     
     static const char FILTER_METHOD_CONVOLUTION_STR[]=  "convolution";
     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";
     static const char FILTER_METHOD_FFTW_STR[]=         "fftw";
+    static const char FILTER_METHOD_RFFTW_STR[]=        "rfftw";
 
     static const char DOMAIN_FREQUENCY_STR[]="frequency";
     static const char DOMAIN_SPATIAL_STR[]="spatial";
@@ -101,12 +106,12 @@ class SignalFilter {
     void filterSignal (const float input[], double output[]) const;
 
     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);
+    static void finiteFourierTransform (const complex<double> input[], double output[], const int n, const int direction);
 
     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;
+    void finiteFourierTransform (const complex<double> input[], double output[], const int direction) const;
 
     void setTraceLevel (int traceLevel) {m_traceLevel = traceLevel; }
 
@@ -146,8 +151,10 @@ class SignalFilter {
     double* m_vecFourierSinTable;
     complex<double>* m_complexVecFilter;
 #ifdef HAVE_FFTW
-    fftw_plan m_planForward, m_planBackward;
-    fftw_complex* m_vecFftInput;
+    fftw_real* m_vecRealFftInput;
+    rfftw_plan m_realPlanForward, m_realPlanBackward;
+    fftw_complex* m_vecComplexFftInput;
+    fftw_plan m_complexPlanForward, m_complexPlanBackward;
 #else
     complex<double>* m_vecFftInput;
 #endif

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 0a24996805bba3a24f48bdb83c1bd1900c9b5e04..6a344d42a6273957a8d847c264bcdebd438c0357 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.13 2000/07/06 08:30:30 kevin Exp $
+**  $Id: filter.cpp,v 1.14 2000/07/06 18:37:24 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
@@ -49,7 +49,6 @@ SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName
   m_vecFilter = NULL;
   m_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;
-  m_vecFftInput = NULL;
   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
     m_fail = true;
@@ -87,7 +86,6 @@ 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);
@@ -130,19 +128,18 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   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) {
+  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) {
     m_nFilterPoints = m_nSignalPoints;
     if (m_zeropad > 0) {
       double logBase2 = log(m_nSignalPoints) / log(2);
       int nextPowerOf2 = static_cast<int>(floor(logBase2));
       if (logBase2 != floor(logBase2))
        nextPowerOf2++;
-      nextPowerOf2 += m_zeropad;
+      nextPowerOf2 += (m_zeropad - 1);
       m_nFilterPoints = 1 << nextPowerOf2;
       cout << "nFilterPoints = " << m_nFilterPoints << endl;
     }
@@ -158,7 +155,7 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   }
 
   // precalculate sin and cosine tables for fourier transform
-  if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
     int nFourier = m_nFilterPoints * m_nFilterPoints + 1;
     double angleIncrement = (2. * PI) / m_nFilterPoints;
     m_vecFourierCosTable = new double[ nFourier ];
@@ -172,15 +169,25 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   }
 
 #if HAVE_FFTW
-  if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+  if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
       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);
-    m_vecFftInput = new fftw_complex [ m_nFilterPoints ];
-    for (int i = 0; i < m_nFilterPoints; i++) 
-       m_vecFftInput[i].re = m_vecFftInput[i].im = 0;
+  if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
+      m_complexPlanForward = m_complexPlanBackward = NULL;
+      m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
+      m_realPlanBackward = rfftw_create_plan (m_nFilterPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
+      m_vecRealFftInput = new fftw_real [ m_nFilterPoints ];
+      for (int i = 0; i < m_nFilterPoints; i++) 
+         m_vecRealFftInput[i] = 0;
+  } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+      m_realPlanForward = m_realPlanBackward = NULL;
+      m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
+      m_complexPlanBackward = fftw_create_plan (m_nFilterPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
+      m_vecComplexFftInput = new fftw_complex [ m_nFilterPoints ];
+      for (int i = 0; i < m_nFilterPoints; i++) 
+         m_vecComplexFftInput[i].re = m_vecComplexFftInput[i].im = 0;
   }
 #endif
 
@@ -227,11 +234,17 @@ SignalFilter::~SignalFilter (void)
     delete [] m_vecFilter;
     delete [] m_vecFourierSinTable;
     delete [] m_vecFourierCosTable;
-    delete [] m_vecFftInput;
+
 #if HAVE_FFTW
     if (m_idFilterMethod == FILTER_METHOD_FFTW) {
-       fftw_destroy_plan(m_planForward);
-       fftw_destroy_plan(m_planBackward);
+       fftw_destroy_plan(m_complexPlanForward);
+       fftw_destroy_plan(m_complexPlanBackward);
+       delete [] m_vecComplexFftInput;
+    }
+    if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
+       rfftw_destroy_plan(m_realPlanForward);
+       rfftw_destroy_plan(m_realPlanBackward);
+       delete [] m_vecRealFftInput;
     }
 #endif
 }
@@ -304,10 +317,14 @@ SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
     fmID = FILTER_METHOD_CONVOLUTION;
   else if (strcasecmp (filterMethodName, FILTER_METHOD_FOURIER_STR) == 0)
     fmID = FILTER_METHOD_FOURIER;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FOURIER_TABLE_STR) == 0)
+    fmID = FILTER_METHOD_FOURIER_TABLE;
   else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_STR) == 0)
     fmID = FILTER_METHOD_FFT;
   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;
 
   return (fmID);
 }
@@ -321,10 +338,14 @@ SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID)
     return (FILTER_METHOD_CONVOLUTION_STR);
   else if (fmID == FILTER_METHOD_FOURIER)
     return (FILTER_METHOD_FOURIER_STR);
+  else if (fmID == FILTER_METHOD_FOURIER_TABLE)
+    return (FILTER_METHOD_FOURIER_TABLE_STR);
   else if (fmID == FILTER_METHOD_FFT)
     return (FILTER_METHOD_FFT_STR);
   else if (fmID == FILTER_METHOD_FFTW)
     return (FILTER_METHOD_FFTW_STR);
+  else if (fmID == FILTER_METHOD_RFFTW)
+    return (FILTER_METHOD_RFFTW_STR);
 
   return (name);
 }
@@ -363,33 +384,60 @@ 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) {
+    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
     complex<double> fftSignal[m_nFilterPoints];
-    complex<double> complexOutput[m_nFilterPoints];
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
     complex<double> filteredSignal[m_nFilterPoints];
+    dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
+    double inverseFourier[m_nFilterPoints];
+    finiteFourierTransform (filteredSignal, inverseFourier, m_nFilterPoints, 1);
+    for (int i = 0; i < m_nSignalPoints; i++) 
+      output[i] = inverseFourier[i];
+  } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
     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);
+    complex<double> fftSignal[m_nFilterPoints];
+    finiteFourierTransform (inputSignal, fftSignal, -1);
+    complex<double> filteredSignal[m_nFilterPoints];
     dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints);
-    finiteFourierTransform (filteredSignal, complexOutput, m_nFilterPoints, 1);
+    double inverseFourier[m_nFilterPoints];
+    finiteFourierTransform (filteredSignal, inverseFourier, 1);
     for (int i = 0; i < m_nSignalPoints; i++) 
-      output[i] = complexOutput[i].real();
+      output[i] = inverseFourier[i];
   }
 #if HAVE_FFTW
-  else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+  else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
       for (int i = 0; i < m_nSignalPoints; i++)
-         m_vecFftInput[i].re = input[i];
+         m_vecRealFftInput[i] = input[i];
+
+      fftw_real out[m_nFilterPoints];
+      rfftw_one (m_realPlanForward, m_vecRealFftInput, out);
+      for (int i = 0; i < m_nFilterPoints; i++) {
+         out[i] *= m_vecFilter[i];
+      }
+      fftw_real outFiltered[m_nFilterPoints];
+      rfftw_one(m_realPlanBackward, out, outFiltered);
+      for (int i = 0; i < m_nSignalPoints; i++) 
+         output[i] = outFiltered[i];
+  } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+      for (int i = 0; i < m_nSignalPoints; i++)
+         m_vecComplexFftInput[i].re = input[i];
 
       fftw_complex out[m_nFilterPoints];
-      fftw_one(m_planForward, m_vecFftInput, out);
+      fftw_one(m_complexPlanForward, m_vecComplexFftInput, 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;
+         out[i].re *= m_vecFilter[i];
+         out[i].im *= m_vecFilter[i];
       }
       fftw_complex outFiltered[m_nFilterPoints];
-      fftw_one(m_planBackward, out, outFiltered);
+      fftw_one(m_complexPlanBackward, out, outFiltered);
       for (int i = 0; i < m_nSignalPoints; i++) 
          output[i] = outFiltered[i].re;
   }
@@ -756,11 +804,11 @@ SignalFilter::finiteFourierTransform (const complex<double> input[], complex<dou
   else 
     direction = 1;
     
-  double angleIncrement = 2 * PI / n;
+  double angleIncrement = direction * 2 * PI / n;
   for (int i = 0; i < n; i++) {
     complex<double> sum (0,0);
     for (int j = 0; j < n; j++) {
-      double angle = i * j * angleIncrement * direction;
+      double angle = i * j * angleIncrement;
       complex<double> exponentTerm (cos(angle), sin(angle));
       sum += input[j] * exponentTerm;
     }
@@ -771,6 +819,28 @@ SignalFilter::finiteFourierTransform (const complex<double> input[], complex<dou
   }
 }
 
+void
+SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
+{
+  if (direction < 0)
+    direction = -1;
+  else 
+    direction = 1;
+    
+  double angleIncrement = direction * 2 * PI / n;
+  for (int i = 0; i < n; i++) {
+      double sumReal = 0;
+    for (int j = 0; j < n; j++) {
+      double angle = i * j * angleIncrement;
+      sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
+    }
+    if (direction < 0) {
+      sumReal /= n;
+    }
+    output[i] = sumReal;
+  }
+}
+
 void
 SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
 {
@@ -784,11 +854,11 @@ SignalFilter::finiteFourierTransform (const double input[], complex<double> outp
     for (int j = 0; j < m_nFilterPoints; j++) {
       int tableIndex = i * j;
       if (direction > 0) {
-       sumReal += input[i] * m_vecFourierCosTable[tableIndex];
-       sumImag += input[i] * m_vecFourierSinTable[tableIndex];
+       sumReal += input[j] * m_vecFourierCosTable[tableIndex];
+       sumImag += input[j] * m_vecFourierSinTable[tableIndex];
       } else {
-       sumReal += input[i] * m_vecFourierCosTable[tableIndex];
-       sumImag -= input[i] * m_vecFourierSinTable[tableIndex];
+       sumReal += input[j] * m_vecFourierCosTable[tableIndex];
+       sumImag -= input[j] * m_vecFourierSinTable[tableIndex];
       }
     }
     if (direction < 0) {
@@ -799,8 +869,7 @@ SignalFilter::finiteFourierTransform (const double input[], complex<double> outp
   }
 }
 
-// (a+bi) * (c + di) = (ac - db) + (bc + da)i
-#if 0
+// (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 void
 SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
 {
@@ -809,26 +878,57 @@ SignalFilter::finiteFourierTransform (const complex<double> input[], complex<dou
   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];
-       sumImag += input[i] * m_vecFourierSinTable[tableIndex];
+       sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] 
+         - input[j].imag() * m_vecFourierSinTable[tableIndex];
+       sumImag += input[j].real() * m_vecFourierSinTable[tableIndex]
+         + input[j].imag() * m_vecFourierCosTable[tableIndex];
       } else {
-       sumReal += input[i] * m_vecFourierCosTable[tableIndex];
-       sumImag -= input[i] * m_vecFourierSinTable[tableIndex];
+       sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] 
+         - input[j].imag() * -m_vecFourierSinTable[tableIndex];
+       sumImag += input[j].real() * -m_vecFourierSinTable[tableIndex]
+         + input[j].imag() * m_vecFourierCosTable[tableIndex];
       }
     }
-    if (direction > 0) {
-      sumReal /= m_nSignalPoints;
-      sumImag /= m_nSignalPoints;
+    if (direction < 0) {
+      sumReal /= m_nFilterPoints;
+      sumImag /= m_nFilterPoints;
     }
     output[i] = complex<double> (sumReal, sumImag);
   }
 }
-#endif
+
+void
+SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
+{
+  if (direction < 0)
+    direction = -1;
+  else 
+    direction = 1;
+    
+  for (int i = 0; i < m_nFilterPoints; i++) {
+      double sumReal = 0;
+    for (int j = 0; j < m_nFilterPoints; j++) {
+      int tableIndex = i * j;
+      if (direction > 0) {
+       sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] 
+         - input[j].imag() * m_vecFourierSinTable[tableIndex];
+      } else {
+       sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] 
+         - input[j].imag() * -m_vecFourierSinTable[tableIndex];
+      }
+    }
+    if (direction < 0) {
+      sumReal /= m_nFilterPoints;
+    }
+    output[i] = sumReal;
+  }
+}
+
 
 void 
 SignalFilter::dotProduct (const double v1[], const complex<double> v2[], complex<double> output[], const int n)

diff --git a/src/pjrec.cpp b/src/pjrec.cpp
index c22b472eb7212d69d9c57439d9d0ec050d478e76..c64f6ae849c958e6e228e0907550d80099f4607a 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.7 2000/07/06 08:30:30 kevin Exp $
+**  $Id: pjrec.cpp,v 1.8 2000/07/06 18:37:24 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
@@ -78,9 +78,11 @@ pjrec_usage (const char *program)
   cout << "  --filter-method  Filter method before backprojections\n";;
   cout << "    convolution      Spatial filtering (default)\n";
   cout << "    fourier          Frequency filtering with discete fourier\n";
+  cout << "    fourier_table    Frequency filtering with table lookup fourier\n";
   cout << "    fft              Fast Fourier Transform\n";
 #if HAVE_FFTW
-  cout << "    fftw             Fast Fourier Transform in the West library\n";
+  cout << "    fftw             Fast Fourier Transform West library\n";
+  cout << "    rfftw            Fast Fourier Transform West (real-mode) library\n";
 #endif
   cout << "  --zeropad n   Set zeropad level (default = 0)\n";
   cout << "                set n to number of powers to two to pad\n";