X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=db00df703287c9e9a045e3bcc8a2091b70d7d2da;hp=5a52c03c300ea3108787bc8ad17dc4356a04c557;hb=ee0105d74fec9d6bfd236e22e9e1d315e46c568e;hpb=806adf54f5b8d061662696b3b498bfab3cd8b2e6 diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 5a52c03..db00df7 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.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; #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]; 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++) + } + int i; + for (i = nSpatialPoints; i < m_nFilterPoints; i++) adSpatialFilter[i] = 0; m_adFilter = new double [m_nFilterPoints]; - complex acInverseFilter [m_nFilterPoints]; + complex* acInverseFilter = new complex [m_nFilterPoints]; finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1); - for (int i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc; + delete adSpatialFilter; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc; + delete acInverseFilter; #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; } #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; + 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); } - } + } 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 fftSignal[m_nFilterPoints]; - finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); - for (int i = 0; i < m_nFilterPoints; i++) + complex* fftSignal = new complex [m_nFilterPoints]; + finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); + 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); + delete fftSignal; + for (i = 0; i < m_nSignalPoints; i++) + output[i] = inverseFourier[i]; + 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 fftSignal[m_nFilterPoints]; - finiteFourierTransform (inputSignal, fftSignal, -1); - for (int i = 0; i < m_nFilterPoints; i++) + complex* fftSignal = new complex [m_nFilterPoints]; + finiteFourierTransform (inputSignal, fftSignal, -1); + 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); + delete fftSignal; + for (i = 0; i < m_nSignalPoints; i++) + output[i] = inverseFourier[i]; + 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 + 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 complexOutput[n]; + complex* complexOutput = new complex [n]; finiteFourierTransform (input, complexOutput, n, direction); for (int i = 0; i < n; i++) - output[i] = complexOutput[i].real(); + output[i] = complexOutput[i].real(); + delete [] complexOutput; } void @@ -848,25 +861,26 @@ ProcessSignal::finiteFourierTransform (const complex input[], double out void ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) { - double* pdTemp = new double [n]; + double* pdTemp = new double [n]; + 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]; + 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]; + 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; }