X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=db00df703287c9e9a045e3bcc8a2091b70d7d2da;hb=7a50ad8860d015c3e06038475929bc62011b561c;hp=a9ab2c7956053165a89eff72eea2ad5ec813a5e5;hpb=bfcc769cf8019eabc8c65c07257c8dbee4b4c977;p=ctsim.git diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index a9ab2c7..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.6 2000/09/02 05:10:39 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,8 +176,9 @@ 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; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { pEZPlot = new EZPlot (*pSGP); @@ -186,8 +187,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); pEZPlot->plot(); } - +#endif shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); +#ifdef HAVE_SGP if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { pEZPlot->ezset ("title Filter Response: Fourier Order"); pEZPlot->ezset ("ylength 0.25"); @@ -195,7 +197,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); pEZPlot->plot(); } +#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"); pEZPlot->ezset ("ylength 0.25"); @@ -203,7 +208,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw pEZPlot->addCurve (m_adFilter, m_nFilterPoints); pEZPlot->plot(); } +#endif shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); +#ifdef HAVE_SGP if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order"); pEZPlot->ezset ("ylength 0.25"); @@ -212,6 +219,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw pEZPlot->plot(); delete pEZPlot; } +#endif for (int i = 0; i < m_nFilterPoints; i++) { m_adFilter[i] /= m_dSignalInc; } @@ -267,6 +275,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; @@ -279,9 +289,10 @@ 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; } } +#ifdef HAVE_SGP EZPlot* pEZPlot = NULL; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { pEZPlot = new EZPlot (*pSGP); @@ -291,7 +302,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw pEZPlot->addCurve (m_adFilter, m_nFilterPoints); pEZPlot->plot(); } +#endif shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints); +#ifdef HAVE_SGP if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { pEZPlot->ezset ("title Filter Filter: Fourier Order"); pEZPlot->ezset ("ylength 0.50"); @@ -300,6 +313,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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; @@ -320,9 +334,10 @@ 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 EZPlot* pEZPlot = NULL; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { pEZPlot = new EZPlot (*pSGP); @@ -333,6 +348,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw pEZPlot->plot(); delete pEZPlot; } +#endif if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { for (int i = 0; i < m_nFilterPoints; i++) adSpatialFilter[i] *= 0.5; @@ -347,29 +363,34 @@ 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"); pEZPlot->ezset ("ylength 0.50"); pEZPlot->ezset ("yporigin 0.50"); pEZPlot->addCurve (m_adFilter, m_nFilterPoints); pEZPlot->plot(); - delete pEZPlot; + delete pEZPlot; } +#endif } } // 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 ]; @@ -506,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) { @@ -520,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) { @@ -584,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; } @@ -653,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 @@ -831,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; } @@ -858,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; }