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 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
        
 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
 
 
 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; }
 
 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)
 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.
 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>
 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="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>
 <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
 **
 **  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
 **
 **  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
 #endif
 
 #ifdef HAVE_FFTW

+#include <rfftw.h>

 #include <fftw.h>
 #endif
 
 #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
 **
 **  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
 **
 **  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
 
 #ifndef FILTER_H
 #define FILTER_H
 

+

 class SignalFilter {
  public:
 
 class SignalFilter {
  public:
 


@@ -49,8 +50,10 @@ class SignalFilter {
        FILTER_METHOD_INVALID,
        FILTER_METHOD_CONVOLUTION,
        FILTER_METHOD_FOURIER,
        FILTER_METHOD_INVALID,
        FILTER_METHOD_CONVOLUTION,
        FILTER_METHOD_FOURIER,

+       FILTER_METHOD_FOURIER_TABLE,

        FILTER_METHOD_FFT,
        FILTER_METHOD_FFTW,
        FILTER_METHOD_FFT,
        FILTER_METHOD_FFTW,

+       FILTER_METHOD_RFFTW

     } FilterMethodID;
 
     typedef enum {
     } 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_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_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";
 
     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);
     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[], 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 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[], 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; }
 
 
     void setTraceLevel (int traceLevel) {m_traceLevel = traceLevel; }
 


@@ -146,8 +151,10 @@ class SignalFilter {
     double* m_vecFourierSinTable;
     complex<double>* m_complexVecFilter;
 #ifdef HAVE_FFTW
     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
 #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
 **
 **  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
 **
 **  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_vecFilter = NULL;
   m_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;

-  m_vecFftInput = NULL;

   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
     m_fail = true;
   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_vecFilter = NULL;
   m_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;

-  m_vecFftInput = NULL;

   m_filterParam = param;  
   m_numIntegral = numIntegral;
   m_idFilter = convertFilterNameToID (filterName);
   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_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_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++;
     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;
     }
       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
   }
 
   // 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 ];
     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 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;
       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
 
   }
 #endif
 


@@ -227,11 +234,17 @@ 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_FFTW) {
 #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
 }
     }
 #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;
     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_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);
 }
 
   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);
     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_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);
 }
 
   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) {
     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> fftSignal[m_nFilterPoints];

-    complex<double> complexOutput[m_nFilterPoints];
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);

     complex<double> filteredSignal[m_nFilterPoints];
     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
     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);
     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++) 
     for (int i = 0; i < m_nSignalPoints; i++) 

-      output[i] = complexOutput[i].real();
+      output[i] = inverseFourier[i];

   }
 #if HAVE_FFTW
   }
 #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++)
       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_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++) {
       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_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;
   }
       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;
     
   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++) {
   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;
     }
       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
 {
 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) {
     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 {
       } 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) {
       }
     }
     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
 {
 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;
     
   else 
     direction = 1;
     

-  for (int i = 0; i < m_nSignalPoints; i++) {
+  for (int i = 0; i < m_nFilterPoints; i++) {

     double sumReal = 0, sumImag = 0;
     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) {
       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 {
       } 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);
   }
 }
     }
     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)
 
 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
 **
 **  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
 **
 **  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 << "  --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 << "    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";
 #endif
   cout << "  --zeropad n   Set zeropad level (default = 0)\n";
   cout << "                set n to number of powers to two to pad\n";