projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r246: More modifications for MSVC
[ctsim.git] / libctsim / procsignal.cpp
diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index 5a52c03c300ea3108787bc8ad17dc4356a04c557..db00df703287c9e9a045e3bcc8a2091b70d7d2da 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.7 2000/09/07 14:29:05 kevin Exp $
+**  $Id: procsignal.cpp,v 1.8 2000/12/06 01:46:43 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
@@ -176,7 +176,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        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* adFrequencyFilter = new double [m_nFilterPoints];
        filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
 #ifdef HAVE_SGP
        EZPlot* pEZPlot = NULL;
@@ -199,6 +199,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        }
 #endif
        ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
+       delete adFrequencyFilter;\r
 #ifdef HAVE_SGP
        if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
          pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
@@ -333,7 +334,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       if (m_traceLevel >= Trace::TRACE_CONSOLE)
        cout << "nFilterPoints = " << m_nFilterPoints << endl;
 #endif
-      double adSpatialFilter [m_nFilterPoints];
+      double* adSpatialFilter = new double [m_nFilterPoints];\r
       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
       filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
 #ifdef HAVE_SGP
@@ -362,15 +363,18 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
          double dScale = 0.5 * sinScale * sinScale;
          adSpatialFilter[i] *= dScale;
        }
-      }
-      for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
+      }\r
+         int i;
+      for (i = nSpatialPoints; i < m_nFilterPoints; i++)
        adSpatialFilter[i] = 0;
 
       m_adFilter = new double [m_nFilterPoints];
-      complex<double> acInverseFilter [m_nFilterPoints];
+      complex<double>* acInverseFilter = new complex<double> [m_nFilterPoints];\r
       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
-      for (int i = 0; i < m_nFilterPoints; i++)
-       m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
+       delete adSpatialFilter;\r
+       for (i = 0; i < m_nFilterPoints; i++)
+         m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
+       delete acInverseFilter;\r
 #ifdef HAVE_SGP
       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
        pEZPlot->ezset ("title Spatial Filter: Inverse");
@@ -378,7 +382,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        pEZPlot->ezset ("yporigin 0.50");
        pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
        pEZPlot->plot();
-       delete pEZPlot;
+       delete pEZPlot;\r
       }
 #endif
     }
@@ -386,7 +390,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
   
   // 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;
+    int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
     double angleIncrement = (2. * PI) / m_nFilterPoints;
     m_adFourierCosTable = new double[ nFourier ];
     m_adFourierSinTable = new double[ nFourier ];
@@ -523,8 +527,9 @@ ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 void
 ProcessSignal::filterSignal (const float constInput[], double output[]) const
 {
-  double input [m_nSignalPoints];
-  for (int i = 0; i < m_nSignalPoints; i++)
+  double* input = new double [m_nSignalPoints];
+  int i;\r
+  for (i = 0; i < m_nSignalPoints; i++)
     input[i] = constInput[i];
 
   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
@@ -537,38 +542,44 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
       int iDetFromCenter = i - (m_nSignalPoints / 2);
       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
     }
-  }
+  }\r
   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
-      for (int i = 0; i < m_nSignalPoints; i++)
+      for (i = 0; i < m_nSignalPoints; i++)
        output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
-    double inputSignal[m_nFilterPoints];
-    for (int i = 0; i < m_nSignalPoints; i++)
+    double* inputSignal = new double [m_nFilterPoints];
+    for (i = 0; i < m_nSignalPoints; i++)
       inputSignal[i] = input[i];
-    for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
+    for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
       inputSignal[i] = 0;  // zeropad
-    complex<double> fftSignal[m_nFilterPoints];
-    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
-    for (int i = 0; i < m_nFilterPoints; i++)
+    complex<double>* fftSignal = new complex<double> [m_nFilterPoints];
+    finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);\r
+       delete inputSignal;
+    for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
-    double inverseFourier[m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
-    for (int i = 0; i < m_nSignalPoints; i++) 
-      output[i] = inverseFourier[i];
+    double* inverseFourier = new double [m_nFilterPoints];
+    finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);\r
+       delete fftSignal;
+    for (i = 0; i < m_nSignalPoints; i++) 
+      output[i] = inverseFourier[i];\r
+       delete inverseFourier;
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
-    double inputSignal[m_nFilterPoints];
-    for (int i = 0; i < m_nSignalPoints; i++)
+    double* inputSignal = new double [m_nFilterPoints];
+    for (i = 0; i < m_nSignalPoints; i++)
       inputSignal[i] = input[i];
-    for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
+    for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
       inputSignal[i] = 0;  // zeropad
-    complex<double> fftSignal[m_nFilterPoints];
-    finiteFourierTransform (inputSignal, fftSignal, -1);
-    for (int i = 0; i < m_nFilterPoints; i++)
+    complex<double>* fftSignal = new complex<double> [m_nFilterPoints];
+    finiteFourierTransform (inputSignal, fftSignal, -1);\r
+       delete inputSignal;
+    for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
-    double inverseFourier[m_nFilterPoints];
-    finiteFourierTransform (fftSignal, inverseFourier, 1);
-    for (int i = 0; i < m_nSignalPoints; i++) 
-      output[i] = inverseFourier[i];
+    double* inverseFourier = new double [m_nFilterPoints];
+    finiteFourierTransform (fftSignal, inverseFourier, 1);\r
+       delete fftSignal;
+    for (i = 0; i < m_nSignalPoints; i++) 
+      output[i] = inverseFourier[i];\r
+       delete inverseFourier;
   }
 #if HAVE_FFTW
   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
@@ -601,7 +612,8 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
       for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) 
          output[i] = ifftOutput[i].re;
   }
-#endif
+#endif\r
+  delete input;
 }
 
 
@@ -670,11 +682,12 @@ for (int i = 0; i < np; i++)
 void
 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 {
-    complex<double> complexOutput[n];
+    complex<double>* complexOutput = new complex<double> [n];
 
     finiteFourierTransform (input, complexOutput, n, direction);
     for (int i = 0; i < n; i++)
-       output[i] = complexOutput[i].real();
+       output[i] = complexOutput[i].real();\r
+       delete [] complexOutput;
 }
 
 void
@@ -848,25 +861,26 @@ ProcessSignal::finiteFourierTransform (const complex<double> input[], double out
 void
 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
 {
-  double* pdTemp = new double [n];
+  double* pdTemp = new double [n];\r
+  int i;
   if (n % 2) { // Odd
     int iHalfN = (n - 1) / 2;
     
     pdTemp[0] = pdVector[iHalfN];
-    for (int i = 0; i < iHalfN; i++)
+    for (i = 0; i < iHalfN; i++)
       pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
-    for (int i = 0; i < iHalfN; i++)
+    for (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++)
+      for (i = 0; i < iHalfN; i++)
        pdTemp[i + 1] = pdVector[i + iHalfN];
-      for (int i = 0; i < iHalfN - 1; i++)
+      for (i = 0; i < iHalfN - 1; i++)
        pdTemp[i + iHalfN + 1] = pdVector[i];
   }
 
-  for (int i = 0; i < n; i++)
+  for (i = 0; i < n; i++)
     pdVector[i] = pdTemp[i];
   delete pdTemp;
 }
@@ -875,25 +889,26 @@ ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
 void
 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
 {
-  double* pdTemp = new double [n];
+  double* pdTemp = new double [n];\r
+  int i;
   if (n % 2) { // Odd
     int iHalfN = (n - 1) / 2;
     
-    pdTemp[iHalfN] = pdVector[0];
-    for (int i = 0; i < iHalfN; i++)
+    pdTemp[iHalfN] = pdVector[0];\r
+       for (i = 0; i < iHalfN; i++)
       pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
-    for (int i = 0; i < iHalfN; i++)
+    for (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++)
+      for (i = 0; i < iHalfN; i++)
        pdTemp[i] = pdVector[i + iHalfN];
-      for (int i = 0; i < iHalfN - 1; i++)
+      for (i = 0; i < iHalfN - 1; i++)
        pdTemp[i + iHalfN + 1] = pdVector[i+1];
   }
 
-  for (int i = 0; i < n; i++)
+  for (i = 0; i < n; i++)
     pdVector[i] = pdTemp[i];
   delete pdTemp;
 }
Computed Tomography simulator and reconstructor