projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r211: *** empty log message ***
[ctsim.git] / libctsim / procsignal.cpp
diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index 5c10c9e9ff2a44b6b9ba602c5352c52bf6c1094e..5a52c03c300ea3108787bc8ad17dc4356a04c557 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.5 2000/08/31 08:38:58 kevin Exp $
+**  $Id: procsignal.cpp,v 1.7 2000/09/07 14:29:05 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
@@ -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, int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength)
+ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength, SGP* pSGP)
     : 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, iTraceLevel, iGeometry, dFocalLength);
+  init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
 }
 
 
 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, int iTraceLevel, int iGeometry, double dFocalLength)
+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, int iGeometry, double dFocalLength, SGP* pSGP)
 {
   m_idFilter = idFilter;
   m_idDomain = idDomain;
@@ -178,54 +178,47 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        m_adFilter = new double[ m_nFilterPoints ];
        double adFrequencyFilter [m_nFilterPoints];
        filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
-       if (m_traceLevel >= Trace::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 ();
+#ifdef HAVE_SGP
+       EZPlot* pEZPlot = NULL;
+       if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
+         pEZPlot = new EZPlot (*pSGP);
+         pEZPlot->ezset ("title Filter Response: Natural Order");
+         pEZPlot->ezset ("ylength 0.25");
+         pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
+         pEZPlot->plot();
        }
-           
+#endif     
        shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
-       if (m_traceLevel >= Trace::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 ();
+#ifdef HAVE_SGP
+       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
+         pEZPlot->ezset ("title Filter Response: Fourier Order");
+         pEZPlot->ezset ("ylength 0.25");
+         pEZPlot->ezset ("yporigin 0.25");
+         pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
+         pEZPlot->plot();
        }
+#endif
        ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
-       if (m_traceLevel >= Trace::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 ();
+#ifdef HAVE_SGP
+       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
+         pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
+         pEZPlot->ezset ("ylength 0.25");
+         pEZPlot->ezset ("yporigin 0.50");
+         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
+         pEZPlot->plot();
        }
+#endif
        shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
-       if (m_traceLevel >= Trace::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 ();
+#ifdef HAVE_SGP
+       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
+         pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
+         pEZPlot->ezset ("ylength 0.25");
+         pEZPlot->ezset ("yporigin 0.75");
+         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
+         pEZPlot->plot();
+         delete pEZPlot;
        }
+#endif
        for (int i = 0; i < m_nFilterPoints; i++) {
            m_adFilter[i] /= m_dSignalInc;
        }
@@ -259,8 +252,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
          nextPowerOf2++;
        nextPowerOf2 += (m_iZeropad - 1);
        m_nFilterPoints = 1 << nextPowerOf2;
+#ifdef DEBUG
        if (m_traceLevel >= Trace::TRACE_CONSOLE)
-       cout << "nFilterPoints = " << m_nFilterPoints << endl;
+         cout << "nFilterPoints = " << m_nFilterPoints << endl;
+#endif
       }
       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 
@@ -279,6 +274,8 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       m_adFilter = new double [m_nFilterPoints];
       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 
+      // This doesn't work!
+      // Need to add filtering for divergent geometries & Frequency/Direct filtering
       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
        for (int i = 0; i < m_nFilterPoints; i++)
          m_adFilter[i] *= 0.5;
@@ -291,32 +288,31 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
          else
            sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
          double dScale = 0.5 * sinScale * sinScale;
-         //      m_adFilter[i] *= dScale;
+         m_adFilter[i] *= dScale;
        }
       }
-      if (m_traceLevel >= Trace::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 ();
+#ifdef HAVE_SGP
+      EZPlot* pEZPlot = NULL;
+      if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
+       pEZPlot = new EZPlot (*pSGP);
+       pEZPlot->ezset ("title Filter Filter: Natural Order");
+       pEZPlot->ezset ("ylength 0.50");
+       pEZPlot->ezset ("yporigin 0.00");
+       pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
+       pEZPlot->plot();
       }
+#endif
       shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
-       if (m_traceLevel >= Trace::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 ();
-       }
+#ifdef HAVE_SGP
+      if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
+       pEZPlot->ezset ("title Filter Filter: Fourier Order");
+       pEZPlot->ezset ("ylength 0.50");
+       pEZPlot->ezset ("yporigin 0.50");
+       pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
+       pEZPlot->plot();
+       delete pEZPlot;
+      }
+#endif
     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
       // calculate number of filter points with zeropadding
       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
@@ -333,22 +329,25 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        m_nFilterPoints = 1 << nextPowerOf2;
       }
       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
+#ifdef DEBUG
       if (m_traceLevel >= Trace::TRACE_CONSOLE)
        cout << "nFilterPoints = " << m_nFilterPoints << endl;
+#endif
       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::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 ();
+#ifdef HAVE_SGP
+      EZPlot* pEZPlot = NULL;
+      if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
+       pEZPlot = new EZPlot (*pSGP);
+       pEZPlot->ezset ("title Spatial Filter: Natural Order");
+       pEZPlot->ezset ("ylength 0.50");
+       pEZPlot->ezset ("yporigin 0.00");
+       pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
+       pEZPlot->plot();
+       delete pEZPlot;
       }
+#endif
       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
        for (int i = 0; i < m_nFilterPoints; i++)
          adSpatialFilter[i] *= 0.5;
@@ -372,17 +371,16 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       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::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 ();
+#ifdef HAVE_SGP
+      if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
+       pEZPlot->ezset ("title Spatial Filter: Inverse");
+       pEZPlot->ezset ("ylength 0.50");
+       pEZPlot->ezset ("yporigin 0.50");
+       pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
+       pEZPlot->plot();
+       delete pEZPlot;
       }
+#endif
     }
   }
   
Computed Tomography simulator and reconstructor