r181: *** empty log message ***

diff --git a/ChangeLog b/ChangeLog
index 72fcec0bf7158fd16e3344a1b860f9c86f066384..50bb4f2cc59ec556f2e01d3db1a8d4ad2479abf5 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,10 +1,12 @@
-2.0.0-b9 - 8/15/00
+2.0.0-b9 - 8/22/00
    Added RCS Id strings to executable files
    Added RPM Spec file for RPM package creation
    Added loading of ASCII phanthom definitions from files
-   Added Filter-Generation option to reconstruction
-   Fixed compilation for non-SGP architectures
-   Decomposed SignalFilter class into ProcessSignal and SignalFilter classes   
+   Fixed compilation for non-SGP architectures 
+   Decomposed SignalFilter class into ProcessSignal and SignalFilter classes
+   Added Filter-Generation option to reconstruction to allow direct or 
+       inverse_fourier construction of filters
+
 2.0.0-b8 - 8/1/00
    Added line color support to SGP
    Fixed lineAbs bug

diff --git a/include/procsignal.h b/include/procsignal.h
index 1b08aea8fd5ce5d76877662acb92902df3c19dbf..8995dd675c519fd8be57001fc5da64c7bfb136ad 100644 (file)
--- a/include/procsignal.h
+++ b/include/procsignal.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.h,v 1.1 2000/08/19 23:00:05 kevin Exp $
+**  $Id: procsignal.h,v 1.2 2000/08/22 07:02:48 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
@@ -58,9 +58,9 @@ class ProcessSignal {
     static const int FILTER_GENERATION_DIRECT;
     static const int FILTER_GENERATION_INVERSE_FOURIER;
 
-    ProcessSignal (const char* szFilterName, const char* szFilterMethodName,double bw, double signalIncrement, int n, double param, const char* szDomainName, const char* szFilterGenerationName, const int zeropad = 0, const int preinterpolationFactor = 1);
+    ProcessSignal (const char* szFilterName, const char* szFilterMethodName,double bw, double signalIncrement, int n, double param, const char* szDomainName, const char* szFilterGenerationName, const int zeropad = 0, const int preinterpolationFactor = 1, const int iTraceLevel = TRACE_NONE);
 
-    ProcessSignal (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad = 0, const int iPreinterpolationFactor = 1);
+    ProcessSignal (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad = 0, const int iPreinterpolationFactor = 1, const int iTraceLevel = TRACE_NONE);
 
     ~ProcessSignal();
 
@@ -150,7 +150,7 @@ class ProcessSignal {
     fftw_plan m_complexPlanForward, m_complexPlanBackward;
 #endif
 
-    void init (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor);
+    void init (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, const int iTraceLevel);
 
     // transforms that use precalculated trig tables, therefore don't 
     // require number of data points (n) as an argument

diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp
index 348dd3f170d4fd06892c50d73b4e9bea8422f9cd..7038132c128db2a6ed0f4ef68457676594aadbda 100644 (file)
--- a/libctgraphics/ezplot.cpp
+++ b/libctgraphics/ezplot.cpp
@@ -6,7 +6,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ezplot.cpp,v 1.10 2000/08/19 22:59:06 kevin Exp $
+**  $Id: ezplot.cpp,v 1.11 2000/08/22 07:02:48 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
@@ -827,7 +827,7 @@ int
 EZPlot::axis_scale (double min, double max, int nint, double *minp, double *maxp, int *nintp)
 {
   if (min >= max || nint < 1) {
-    sys_error (ERR_WARNING, "Invalid params: min=%lf, min=%lf, num intervals=%d [axis_scale]", min, max, nint);
+    sys_error (ERR_WARNING, "Invalid params: min=%lf, max=%lf, num intervals=%d [axis_scale]", min, max, nint);
     return (FALSE);
   }
   

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index a88fad5a51c3cb1a663908a81b22dc95f8f7d07b..941d905ea76512f78814763a11ce62bd5de2b02a 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.25 2000/08/19 22:59:06 kevin Exp $
+**  $Id: filter.cpp,v 1.26 2000/08/22 07:02:48 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
@@ -401,31 +401,31 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param
 
   switch (filterID) {
   case FILTER_BANDLIMIT:
-    if (au >= bw / 2)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0.;
     else
       q = 1;
     break;
   case FILTER_ABS_BANDLIMIT:
-    if (au >= bw / 2)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0.;
     else
       q = au;
     break;
   case FILTER_TRIANGLE:
-    if (au >= bw)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0;
     else
       q = 1 - au / bw;
     break;
   case FILTER_COSINE:
-    if (au >= bw / 2)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0;
     else
       q = cos(PI * u / bw);
     break;
   case FILTER_ABS_COSINE:
-    if (au >= bw / 2)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0;
     else
       q = au * cos(PI * u / bw);
@@ -437,13 +437,13 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param
     q = au * bw * sinc (PI * bw * u, 1.);
     break;
   case FILTER_G_HAMMING:
-    if (au >= bw / 2)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0;
     else
       q = param + (1 - param) * cos (TWOPI * u / bw);
     break;
   case FILTER_ABS_G_HAMMING:
-    if (au >= bw / 2)
+    if (fabs(au) > fabs(bw / 2) + F_EPSILON)
       q = 0;
     else
       q = au * (param + (1 - param) * cos(TWOPI * u / bw));

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index 9fdbb682a1fbb07f9aa30022178bd25e5fafd1e7..70d759cbacc471f96ab39b3d8fb3ea3dd990e378 100644 (file)
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.cpp,v 1.1 2000/08/19 23:00:05 kevin Exp $
+**  $Id: procsignal.cpp,v 1.2 2000/08/22 07:02:48 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,8 +49,8 @@ const char* ProcessSignal::s_aszFilterMethodName[] = {
 };
 const char* ProcessSignal::s_aszFilterMethodTitle[] = {
   {"Convolution"},
-  {"Direct Fourier"},
-  {"Fouier Trigometric Table Lookout"},
+  {"Fourier"},
+  {"Fouier Trigometric Table"},
   {"FFT"},
 #if HAVE_FFTW
   {"FFTW"},
@@ -77,7 +77,7 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration
 // CLASS IDENTIFICATION
 //   ProcessSignal
 //
-ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad = 0, int iPreinterpolationFactor = 1)
+ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad = 0, int iPreinterpolationFactor = 1, int iTraceLevel = TRACE_NONE)
     : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
 {
   m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
@@ -109,12 +109,12 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth
     return;
   }
 
-  init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor);
+  init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel);
 }
 
 
 void
-ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor)
+ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel)
 {
   m_idFilter = idFilter;
   m_idDomain = idDomain;
@@ -124,7 +124,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
     m_fail = true;
     return;
   }
-  m_traceLevel = TRACE_NONE;
+  m_traceLevel = iTraceLevel;
   m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
   m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
   m_dBandwidth = dBandwidth;
@@ -152,7 +152,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
   if (! m_bFrequencyFiltering) {
 
     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
-       m_nFilterPoints = 2 * m_nSignalPoints - 1;
+       m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
        m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
        m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
        m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
@@ -160,63 +160,177 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        m_adFilter = new double[ m_nFilterPoints ];
        filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
-       m_nFilterPoints = m_nSignalPoints;
+       m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
        m_dFilterMin = -1. / (2 * m_dSignalInc);
        m_dFilterMax = 1. / (2 * m_dSignalInc);
        m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
        SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
        m_adFilter = new double[ m_nFilterPoints ];
        double adFrequencyFilter [m_nFilterPoints];
-       double adInverseFilter [m_nFilterPoints];
        filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
+       if (m_traceLevel >= TRACE_PLOT) {
+         SGPDriver sgpDriver ("Frequency Filter: Natural Order");
+         SGP sgp (sgpDriver);
+         EZPlot ezplot (sgp);
+
+         ezplot.ezset ("title Filter Response: Natural Order");
+         ezplot.addCurve (adFrequencyFilter, m_nFilterPoints);
+         ezplot.plot();
+         cio_put_str ("Press any key to continue");
+         cio_kb_getc ();
+       }
+           
        shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
-       ProcessSignal::finiteFourierTransform (adFrequencyFilter, adInverseFilter, m_nFilterPoints, 1);
-       for (int i = 0; i < m_nFilterPoints; i++)
-           m_adFilter [i] = adInverseFilter[i];
+       if (m_traceLevel >= TRACE_PLOT) {
+         SGPDriver sgpDriver ("Frequency Filter: Fourier Order");
+         SGP sgp (sgpDriver);
+         EZPlot ezplot (sgp);
+
+         ezplot.ezset ("title Filter Response: Fourier Order");
+         ezplot.addCurve (adFrequencyFilter, m_nFilterPoints);
+         ezplot.plot();
+         cio_put_str ("Press any key to continue");
+         cio_kb_getc ();
+       }
+       ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
+       if (m_traceLevel >= TRACE_PLOT) {
+         SGPDriver sgpDriver ("Inverse Fourier Frequency: Fourier Order");
+         SGP sgp (sgpDriver);
+         EZPlot ezplot (sgp);
+
+         ezplot.ezset ("title Inverse Fourier Frequency: Fourier Order");
+         ezplot.addCurve (m_adFilter, m_nFilterPoints);
+         ezplot.plot();
+         cio_put_str ("Press any key to continue");
+         cio_kb_getc ();
+       }
+       shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
+       if (m_traceLevel >= TRACE_PLOT) {
+         SGPDriver sgpDriver ("Inverse Fourier Frequency: Natural Order");
+         SGP sgp (sgpDriver);
+         EZPlot ezplot (sgp);
+
+         ezplot.ezset ("title Inverse Fourier Frequency: Natural Order");
+         ezplot.addCurve (m_adFilter, m_nFilterPoints);
+         ezplot.plot();
+         cio_put_str ("Press any key to continue");
+         cio_kb_getc ();
+       }
+       for (int i = 0; i < m_nFilterPoints; i++) {
+           m_adFilter[i] /= m_dSignalInc;
+       }
     }
   } 
 
   // Frequency-based filtering
   else if (m_bFrequencyFiltering) {
 
-    // calculate number of filter points with zeropadding
-    m_nFilterPoints = m_nSignalPoints;
-    if (m_iZeropad > 0) {
-      double logBase2 = log(m_nSignalPoints) / log(2);
-      int nextPowerOf2 = static_cast<int>(floor(logBase2));
-      if (logBase2 != floor(logBase2))
-       nextPowerOf2++;
-      nextPowerOf2 += (m_iZeropad - 1);
-      m_nFilterPoints = 1 << nextPowerOf2;
-      if (m_traceLevel >= TRACE_TEXT)
+    if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
+      // calculate number of filter points with zeropadding
+      m_nFilterPoints = m_nSignalPoints;
+      if (m_iZeropad > 0) {
+       double logBase2 = log(m_nFilterPoints) / log(2);
+       int nextPowerOf2 = static_cast<int>(floor(logBase2));
+       if (logBase2 != floor(logBase2))
+         nextPowerOf2++;
+       nextPowerOf2 += (m_iZeropad - 1);
+       m_nFilterPoints = 1 << nextPowerOf2;
+       if (m_traceLevel >= TRACE_TEXT)
        cout << "nFilterPoints = " << m_nFilterPoints << endl;
-    }
-    m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
+      }
+      m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 
-    if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
+      if (m_nFilterPoints % 2) { // Odd
        m_dFilterMin = -1. / (2 * m_dSignalInc);
        m_dFilterMax = 1. / (2 * m_dSignalInc);
        m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
-       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
-       m_adFilter = new double [m_nFilterPoints];
-       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
-       shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
+      } else { // Even
+       m_dFilterMin = -1. / (2 * m_dSignalInc);
+       m_dFilterMax = 1. / (2 * m_dSignalInc);
+       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
+       m_dFilterMax -= m_dFilterInc;
+      }
+
+      SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
+      m_adFilter = new double [m_nFilterPoints];
+      filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
+      if (m_traceLevel >= TRACE_PLOT) {
+       SGPDriver sgpDriver ("Frequency Filter: Natural Order");
+       SGP sgp (sgpDriver);
+       EZPlot ezplot (sgp);
+       
+       ezplot.ezset ("title Filter Filter: Natural Order");
+       ezplot.addCurve (m_adFilter, m_nFilterPoints);
+         ezplot.plot();
+         cio_put_str ("Press any key to continue");
+         cio_kb_getc ();
+      }
+      shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
+       if (m_traceLevel >= TRACE_PLOT) {
+         SGPDriver sgpDriver ("Frequency Filter: Fourier Order");
+         SGP sgp (sgpDriver);
+         EZPlot ezplot (sgp);
+
+         ezplot.ezset ("title Filter Filter: Fourier Order");
+         ezplot.addCurve (m_adFilter, m_nFilterPoints);
+         ezplot.plot();
+         cio_put_str ("Press any key to continue");
+         cio_kb_getc ();
+       }
     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
-       m_nFilterPoints = 2 * m_nSignalPoints - 1;
-       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
-       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
-       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
-       double adSpatialFilter [m_nFilterPoints];
-       double adInverseFilter [m_nFilterPoints];
-       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
-       filter.copyFilterData (adSpatialFilter, 0, m_nFilterPoints);
-       m_adFilter = new double [m_nFilterPoints];
-       finiteFourierTransform (adSpatialFilter, adInverseFilter, m_nFilterPoints, -1);
-       for (int i = 0; i < m_nFilterPoints; i++)
-           m_adFilter [i] = adInverseFilter[i];
+      // calculate number of filter points with zeropadding
+      int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
+      m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
+      m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
+      m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
+      m_nFilterPoints = nSpatialPoints;
+      if (m_iZeropad > 0) {
+       double logBase2 = log(nSpatialPoints) / log(2);
+       int nextPowerOf2 = static_cast<int>(floor(logBase2));
+       if (logBase2 != floor(logBase2))
+         nextPowerOf2++;
+       nextPowerOf2 += (m_iZeropad - 1);
+       m_nFilterPoints = 1 << nextPowerOf2;
       }
-    }
+      m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
+      if (m_traceLevel >= TRACE_TEXT)
+       cout << "nFilterPoints = " << m_nFilterPoints << endl;
+      double adSpatialFilter [m_nFilterPoints];
+      SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
+      filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
+      if (m_traceLevel >= TRACE_PLOT) {
+       SGPDriver sgpDriver ("Spatial Filter: Natural Order");
+       SGP sgp (sgpDriver);
+       EZPlot ezplot (sgp);
+       
+       ezplot.ezset ("title Spatial Filter: Natural Order");
+       ezplot.addCurve (adSpatialFilter, nSpatialPoints);
+       ezplot.plot();
+       cio_put_str ("Press any key to continue");
+       cio_kb_getc ();
+      }
+      for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
+       adSpatialFilter[i] = 0;
 
+      m_adFilter = new double [m_nFilterPoints];
+      complex<double> acInverseFilter [m_nFilterPoints];
+      finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
+      for (int i = 0; i < m_nFilterPoints; i++)
+       m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
+      if (m_traceLevel >= TRACE_PLOT) {
+       SGPDriver sgpDriver ("Spatial Filter: Inverse");
+       SGP sgp (sgpDriver);
+       EZPlot ezplot (sgp);
+       
+       ezplot.ezset ("title Spatial Filter: Inverse");
+       ezplot.addCurve (m_adFilter, m_nFilterPoints);
+       ezplot.plot();
+       cio_put_str ("Press any key to continue");
+       cio_kb_getc ();
+      }
+    }
+  }
+  
   // precalculate sin and cosine tables for fourier transform
   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
     int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
@@ -262,6 +376,7 @@ ProcessSignal::~ProcessSignal (void)
 {
     delete [] m_adFourierSinTable;
     delete [] m_adFourierCosTable;
+    delete [] m_adFilter;
 
 #if HAVE_FFTW
     if (m_idFilterMethod == FILTER_METHOD_FFTW) {
@@ -491,7 +606,7 @@ ProcessSignal::finiteFourierTransform (const double input[], double output[], co
 
     finiteFourierTransform (input, complexOutput, n, direction);
     for (int i = 0; i < n; i++)
-       output[i] = abs(complexOutput[n]);
+       output[i] = complexOutput[i].real();
 }
 
 void
@@ -670,13 +785,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
     int iHalfN = (n - 1) / 2;
     
     pdTemp[0] = pdVector[iHalfN];
-    for (int i = 1; i <= iHalfN; i++)
-      pdTemp[i] = pdVector[i+iHalfN];
-    for (int i = iHalfN+1; i < n; i++)
-      pdTemp[i] = pdVector[i-iHalfN];
+    for (int i = 0; i < iHalfN; i++)
+      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
+    for (int i = 0; i < iHalfN; i++)
+      pdTemp[i + iHalfN + 1] = pdVector[i];
   } else {     // Even
       int iHalfN = n / 2;
       pdTemp[0] = pdVector[iHalfN];
+      for (int i = 0; i < iHalfN; i++)
+       pdTemp[i + 1] = pdVector[i + iHalfN];
+      for (int i = 0; i < iHalfN - 1; i++)
+       pdTemp[i + iHalfN + 1] = pdVector[i];
   }
 
   for (int i = 0; i < n; i++)
@@ -693,13 +812,17 @@ ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
     int iHalfN = (n - 1) / 2;
     
     pdTemp[iHalfN] = pdVector[0];
-    for (int i = 1; i <= iHalfN; i++)
-      pdTemp[i] = pdVector[i+iHalfN];
-    for (int i = iHalfN+1; i < n; i++)
-      pdTemp[i] = pdVector[i-iHalfN];
+    for (int i = 0; i < iHalfN; i++)
+      pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
+    for (int i = 0; i < iHalfN; i++)
+      pdTemp[i] = pdVector[i + iHalfN + 1];
   } else {     // Even
       int iHalfN = n / 2;
       pdTemp[iHalfN] = pdVector[0];
+      for (int i = 0; i < iHalfN; i++)
+       pdTemp[i] = pdVector[i + iHalfN];
+      for (int i = 0; i < iHalfN - 1; i++)
+       pdTemp[i + iHalfN + 1] = pdVector[i+1];
   }
 
   for (int i = 0; i < n; i++)

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 813d7bc59f22fdc0973589eaffe354a655c270fa..d6b496d23cc27eeb6961955495fe66d56208903b 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.20 2000/08/19 22:59:06 kevin Exp $
+**  $Id: projections.cpp,v 1.21 2000/08/22 07:02:48 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
@@ -503,8 +503,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 #endif
 
   double filterBW = 1. / detInc;
-  ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor);
-  processSignal.setTraceLevel(trace);
+  ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace);
 
   if (processSignal.fail()) {
       sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str());

diff --git a/src/dialogs.cpp b/src/dialogs.cpp
index 79913b4d28ac76de322a2164b8b78cf42afddf9b..4b2738ff2111cac5f5ef58fff26540812a38dff5 100644 (file)
--- a/src/dialogs.cpp
+++ b/src/dialogs.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.cpp,v 1.7 2000/07/23 01:49:03 kevin Exp $
+**  $Id: dialogs.cpp,v 1.8 2000/08/22 07:02:48 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
@@ -398,7 +398,7 @@ DialogGetProjectionParameters::getGeometry (void)
 /////////////////////////////////////////////////////////////////////
 
 
-DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1.,  int iDefaultFilterMethodID = SignalFilter::FILTER_METHOD_CONVOLUTION, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3)
+DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1.,  int iDefaultFilterMethodID = ProcessSignal::FILTER_METHOD_CONVOLUTION, int iDefaultFilterGenerationID = ProcessSignal::FILTER_GENERATION_INVALID, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3)
     : wxDialog (pParent, -1, "Set Reconstruction Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)
 {
   wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL);
@@ -411,9 +411,12 @@ DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* p
   m_pListBoxFilter->SetSelection (iDefaultFilterID);
   pTopSizer->Add (m_pListBoxFilter, 0, wxALL | wxALIGN_CENTER | wxEXPAND);
 
-  m_pListBoxFilterMethod = new StringValueAndTitleListBox (this, SignalFilter::getFilterMethodCount(), SignalFilter::getFilterMethodTitleArray(), SignalFilter::getFilterMethodNameArray());
+  m_pListBoxFilterMethod = new StringValueAndTitleListBox (this, ProcessSignal::getFilterMethodCount(), ProcessSignal::getFilterMethodTitleArray(), ProcessSignal::getFilterMethodNameArray());
   m_pListBoxFilterMethod->SetSelection (iDefaultFilterMethodID);
   pTopSizer->Add (m_pListBoxFilterMethod, 0, wxALL | wxALIGN_CENTER | wxEXPAND);
+  m_pListBoxFilterGeneration = new StringValueAndTitleListBox (this, ProcessSignal::getFilterGenerationCount(), ProcessSignal::getFilterGenerationTitleArray(), ProcessSignal::getFilterGenerationNameArray());
+  m_pListBoxFilterGeneration->SetSelection (iDefaultFilterGenerationID);
+  pTopSizer->Add (m_pListBoxFilterGeneration, 0, wxALL | wxALIGN_CENTER | wxEXPAND);
 
   m_pListBoxBackproject = new StringValueAndTitleListBox (this, Backprojector::getBackprojectCount(), Backprojector::getBackprojectTitleArray(), Backprojector::getBackprojectNameArray());
   m_pListBoxBackproject->SetSelection (iDefaultBackprojectID);
@@ -552,3 +555,9 @@ DialogGetReconstructionParameters::getBackprojectName (void)
 {
   return m_pListBoxBackproject->getSelectionStringValue();
 }
+
+const char*
+DialogGetReconstructionParameters::getFilterGenerationName (void)
+{
+  return m_pListBoxFilterGeneration->getSelectionStringValue();
+}

diff --git a/src/dialogs.h b/src/dialogs.h
index 6dcc2988cd5f5379114f045b8a741b67d64708e4..e4f890cbd44e8b491812db9d053416badcf2d7e1 100644 (file)
--- a/src/dialogs.h
+++ b/src/dialogs.h
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.h,v 1.8 2000/07/28 08:28:08 kevin Exp $
+**  $Id: dialogs.h,v 1.9 2000/08/22 07:02:48 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
@@ -33,6 +33,7 @@
 #include <string>
 #include "scanner.h"
 #include "phantom.h"
+#include "procsignal.h"
 #include "filter.h"
 
 // CLASS StringValueAndTitleListBox
@@ -136,7 +137,7 @@ class DialogGetProjectionParameters : public wxDialog
 class DialogGetReconstructionParameters : public wxDialog
 {
  public:
-    DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = SignalFilter::FILTER_METHOD_CONVOLUTION, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3);
+    DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = ProcessSignal::FILTER_METHOD_CONVOLUTION, int iDefaultFilterGeneration = ProcessSignal::FILTER_GENERATION_INVALID, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3);
     virtual ~DialogGetReconstructionParameters (void);
 
     unsigned int getXSize(void);
@@ -145,6 +146,7 @@ class DialogGetReconstructionParameters : public wxDialog
     double getFilterParam(void);
     const char* getFilterMethodName(void);
     unsigned int getZeropad(void);
+    const char* getFilterGenerationName(void);
     const char* getInterpName(void);
     unsigned int getInterpParam(void);
     const char* getBackprojectName(void);
@@ -158,6 +160,7 @@ class DialogGetReconstructionParameters : public wxDialog
 
     StringValueAndTitleListBox* m_pListBoxFilter;
     StringValueAndTitleListBox* m_pListBoxFilterMethod;
+    StringValueAndTitleListBox* m_pListBoxFilterGeneration;
     StringValueAndTitleListBox* m_pListBoxInterp;
     StringValueAndTitleListBox* m_pListBoxBackproject;
 

diff --git a/src/views.cpp b/src/views.cpp
index 1ca02ebf624e834c812079968c83138b65907c2d..0ccbc09f47e7782559a52c8ae22cec859c78d2c2 100644 (file)
--- a/src/views.cpp
+++ b/src/views.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.13 2000/08/02 18:06:00 kevin Exp $
+**  $Id: views.cpp,v 1.14 2000/08/22 07:02:48 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
@@ -625,7 +625,7 @@ ProjectionFileView::OnProperties (wxCommandEvent& event)
 void
 ProjectionFileView::OnReconstruct (wxCommandEvent& event)
 {
-  DialogGetReconstructionParameters dialogReconstruction (m_frame, 256, 256, SignalFilter::FILTER_ABS_BANDLIMIT, 1., SignalFilter::FILTER_METHOD_CONVOLUTION, 3, Backprojector::INTERP_LINEAR, 1, Backprojector::BPROJ_IDIFF3);
+  DialogGetReconstructionParameters dialogReconstruction (m_frame, 256, 256, SignalFilter::FILTER_ABS_BANDLIMIT, 1., ProcessSignal::FILTER_METHOD_CONVOLUTION, ProcessSignal::FILTER_GENERATION_INVALID, 3, Backprojector::INTERP_LINEAR, 1, Backprojector::BPROJ_IDIFF3);
   int retVal = dialogReconstruction.ShowModal();
   if (retVal == wxID_OK) {
     int xSize = dialogReconstruction.getXSize();
@@ -634,6 +634,7 @@ ProjectionFileView::OnReconstruct (wxCommandEvent& event)
     double optFilterParam = dialogReconstruction.getFilterParam();
     wxString optFilterMethodName = dialogReconstruction.getFilterMethodName();
     int optZeropad = dialogReconstruction.getZeropad();
+    wxString optFilterGenerationName = dialogReconstruction.getFilterGenerationName();
     wxString optInterpName = dialogReconstruction.getInterpName();
     int optInterpParam = dialogReconstruction.getInterpParam();
     wxString optBackprojectName = dialogReconstruction.getBackprojectName();
@@ -642,7 +643,7 @@ ProjectionFileView::OnReconstruct (wxCommandEvent& event)
       ImageFile& imageFile = pReconDoc->getImageFile();
       const Projections& rProj = GetDocument()->getProjections();
       imageFile.setArraySize (xSize, ySize);
-      rProj.reconstruct (imageFile, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeropad, optInterpName.c_str(), optInterpParam, optBackprojectName.c_str(), TRACE_NONE);
+      rProj.reconstruct (imageFile, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeropad, optFilterGenerationName.c_str(), optInterpName.c_str(), optInterpParam, optBackprojectName.c_str(), TRACE_NONE);
       pReconDoc->Modify(true);
       pReconDoc->UpdateAllViews(this);
       ostringstream os;

diff --git a/tools/pjrec.cpp b/tools/pjrec.cpp
index 1638dbea45023bd5856b75c1598e4781b115a17e..10adbff4a2d929eb60cfbe85cc5c39487f298386 100644 (file)
--- a/tools/pjrec.cpp
+++ b/tools/pjrec.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pjrec.cpp,v 1.12 2000/08/19 22:59:06 kevin Exp $
+**  $Id: pjrec.cpp,v 1.13 2000/08/22 07:02:48 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,7 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.12 2000/08/19 22:59:06 kevin Exp $";
+static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.13 2000/08/22 07:02:48 kevin Exp $";
 
 void 
 pjrec_usage (const char *program)
@@ -72,12 +72,12 @@ pjrec_usage (const char *program)
   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_cosine    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 << "    cosine        Cosine" << endl;
   cout << "    triangle      Triangle" << endl;
   cout << "    hamming       Hamming" << endl;
   cout << "  --filter-method  Filter method before backprojections\n";;