r256: *** empty log message ***

[ctsim.git] / libctsim / procsignal.cpp

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index db00df703287c9e9a045e3bcc8a2091b70d7d2da..cd44cf834d8bd0d2b776c1a0a48f958229114329 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.8 2000/12/06 01:46:43 kevin Exp $
+**  $Id: procsignal.cpp,v 1.9 2000/12/16 02:44:26 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
@@ -115,7 +115,8 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth
 
 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, SGP* pSGP)
-{
+{\r
+  int i;
   m_idFilter = idFilter;
   m_idDomain = idDomain;
   m_idFilterMethod = idFilterMethod;
@@ -220,15 +221,15 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
          delete pEZPlot;
        }
 #endif
-       for (int i = 0; i < m_nFilterPoints; i++) {
+       for (i = 0; i < m_nFilterPoints; i++) {
            m_adFilter[i] /= m_dSignalInc;
        }
     }
     if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
-       for (int i = 0; i < m_nFilterPoints; i++)
+       for (i = 0; i < m_nFilterPoints; i++)
            m_adFilter[i] *= 0.5;
     } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
-       for (int i = 0; i < m_nFilterPoints; i++) {
+       for (i = 0; i < m_nFilterPoints; i++) {
          int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
          double sinScale = sin (iDetFromZero * m_dSignalInc);
          if (fabs(sinScale) < 1E-7)
@@ -278,10 +279,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       // 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++)
+       for (i = 0; i < m_nFilterPoints; i++)
          m_adFilter[i] *= 0.5;
       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
-       for (int i = 0; i < m_nFilterPoints; i++) {
+       for (i = 0; i < m_nFilterPoints; i++) {
          int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
          double sinScale = sin (iDetFromZero * m_dSignalInc);
          if (fabs(sinScale) < 1E-7)
@@ -350,10 +351,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       }
 #endif
       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
-       for (int i = 0; i < m_nFilterPoints; i++)
+       for (i = 0; i < m_nFilterPoints; i++)
          adSpatialFilter[i] *= 0.5;
       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
-       for (int i = 0; i < m_nFilterPoints; i++) {
+       for (i = 0; i < m_nFilterPoints; i++) {
          int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
          double sinScale = sin (iDetFromZero * m_dSignalInc);
          if (fabs(sinScale) < 1E-7)
@@ -364,7 +365,6 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
          adSpatialFilter[i] *= dScale;
        }
       }\r
-         int i;
       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
        adSpatialFilter[i] = 0;
 
@@ -395,7 +395,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
     m_adFourierCosTable = new double[ nFourier ];
     m_adFourierSinTable = new double[ nFourier ];
     double angle = 0;
-    for (int i = 0; i < nFourier; i++) {
+    for (i = 0; i < nFourier; i++) {
       m_adFourierCosTable[i] = cos (angle);
       m_adFourierSinTable[i] = sin (angle);
       angle += angleIncrement;
@@ -404,7 +404,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
 
 #if HAVE_FFTW
   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
-    for (int i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
+    for (i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
       m_adFilter[i] /= m_nFilterPoints;
   }
 
@@ -413,16 +413,16 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
-    for (int i = 0; i < m_nFilterPoints; i++) 
+    for (i = 0; i < m_nFilterPoints; i++) 
       m_adRealFftInput[i] = 0;
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
-    for (int i = 0; i < m_nFilterPoints; i++) 
+    for (i = 0; i < m_nFilterPoints; i++) 
       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
-    for (int i = 0; i < m_nOutputPoints; i++) 
+    for (i = 0; i < m_nOutputPoints; i++) 
       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
   }
 #endif
@@ -583,34 +583,38 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
   }
 #if HAVE_FFTW
   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
-      for (int i = 0; i < m_nSignalPoints; i++)
+      for (i = 0; i < m_nSignalPoints; i++)
          m_adRealFftInput[i] = input[i];
 
-      fftw_real fftOutput [ m_nFilterPoints ];
+      fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
       rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
-      for (int i = 0; i < m_nFilterPoints; i++)
-         m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
-      for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
-       m_adRealFftSignal[i] = 0;
+      for (i = 0; i < m_nFilterPoints; i++)
+           m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];\r
+         delete [] fftOutput;
+      for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
+            m_adRealFftSignal[i] = 0;
 
-      fftw_real ifftOutput [ m_nOutputPoints ];
+      fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
       rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
-      for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
-         output[i] = ifftOutput[i];
+      for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
+           output[i] = ifftOutput[i];\r
+         delete [] ifftOutput;
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
-      for (int i = 0; i < m_nSignalPoints; i++)
+      for (i = 0; i < m_nSignalPoints; i++)
          m_adComplexFftInput[i].re = input[i];
 
-      fftw_complex fftOutput [ m_nFilterPoints ];
+      fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
       fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
-      for (int i = 0; i < m_nFilterPoints; i++) {
+      for (i = 0; i < m_nFilterPoints; i++) {
          m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
          m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
-      }
-      fftw_complex ifftOutput [ m_nOutputPoints ];
+      }\r
+         delete [] fftOutput;
+      fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
       fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
-      for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) 
-         output[i] = ifftOutput[i].re;
+      for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) 
+         output[i] = ifftOutput[i].re;\r
+         delete [] ifftOutput;
   }
 #endif\r
   delete input;