X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=cd44cf834d8bd0d2b776c1a0a48f958229114329;hp=db00df703287c9e9a045e3bcc8a2091b70d7d2da;hb=6afa21de8aa00b405de47584efe108c71df33e1b;hpb=ee0105d74fec9d6bfd236e22e9e1d315e46c568e

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index db00df7..cd44cf8 100644
--- 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)
-{
+{
+  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;
 	}
       }
-	  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];
+	  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];
+	  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 ];
+      }
+	  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;
+	  delete [] ifftOutput;
   }
 #endif
   delete input;