X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=a7933b21eb2a4ffb508c88d5caa35ffc834d3723;hp=7bf2711f0a626190dd01fa297d478b0d5c1a1431;hb=5c6b29ab4885308cc3381af6e0a68f4804956d2e;hpb=c24c1c0721df40e77822ad2b9ec01a944012ff42 diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 7bf2711..a7933b2 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.10 2000/12/16 06:12:47 kevin Exp $ +** $Id: procsignal.cpp,v 1.11 2000/12/29 15:45:06 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 @@ -78,7 +78,7 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration // ProcessSignal // ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength, SGP* pSGP) - : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) +: m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) { m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); if (m_idFilterMethod == FILTER_METHOD_INVALID) { @@ -108,7 +108,7 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth m_failMessage += szDomainName; return; } - + init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP); } @@ -123,7 +123,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_idFilterGeneration = idFilterGeneration; m_idGeometry = iGeometry; m_dFocalLength = dFocalLength; - + if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) { m_fail = true; return; @@ -137,14 +137,14 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterParam = dFilterParam; m_iZeropad = iZeropad; m_iPreinterpolationFactor = iPreinterpolationFactor; - + // scale signalInc/BW to signalInc/2 to adjust for imaginary detector // through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { m_dSignalInc /= 2; m_dBandwidth *= 2; } - + if (m_idFilterMethod == FILTER_METHOD_FFT) { #if HAVE_FFTW m_idFilterMethod = FILTER_METHOD_RFFTW; @@ -154,165 +154,165 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw return; #endif } - + bool m_bFrequencyFiltering = true; if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) m_bFrequencyFiltering = false; - + // Spatial-based filtering if (! m_bFrequencyFiltering) { - + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { - m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; - m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); - m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); - 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_SPATIAL); - m_adFilter = new double[ m_nFilterPoints ]; - filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); + 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_SPATIAL); + m_adFilter = new double[ m_nFilterPoints ]; + filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { - m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - 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 = new double [m_nFilterPoints]; - filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + 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 = 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); - pEZPlot->ezset ("title Filter Response: Natural Order"); - pEZPlot->ezset ("ylength 0.25"); - pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); - pEZPlot->plot(); - } + EZPlot* pEZPlot = NULL; + if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot = new EZPlot (); + pEZPlot->ezset ("title Filter Response: Natural Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + } #endif - shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); + 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"); - pEZPlot->ezset ("yporigin 0.25"); - pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); - pEZPlot->plot(); - } + if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot->ezset ("title Filter Response: Fourier Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->ezset ("yporigin 0.25"); + pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + } #endif - ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1); - delete adFrequencyFilter; + 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"); - pEZPlot->ezset ("yporigin 0.50"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - } + if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->ezset ("yporigin 0.50"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + } #endif - shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); + 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"); - pEZPlot->ezset ("yporigin 0.75"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - delete pEZPlot; - } + if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->ezset ("yporigin 0.75"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; + } #endif - for (i = 0; i < m_nFilterPoints; i++) { - m_adFilter[i] /= m_dSignalInc; - } + for (i = 0; i < m_nFilterPoints; i++) { + m_adFilter[i] /= m_dSignalInc; + } } if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - for (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] *= 0.5; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] *= 0.5; } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { - 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) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - m_adFilter[i] *= dScale; - } + 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) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + m_adFilter[i] *= dScale; + } } // if (geometry) } // if (spatial filtering) - + else if (m_bFrequencyFiltering) { // Frequency-based filtering - + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { // calculate number of filter points with zeropadding m_nFilterPoints = m_nSignalPoints; if (m_iZeropad > 0) { - double logBase2 = log(m_nFilterPoints) / log(2); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_iZeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; + double logBase2 = log(m_nFilterPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; #ifdef DEBUG - if (m_traceLevel >= Trace::TRACE_CONSOLE) - std::cout << "nFilterPoints = " << m_nFilterPoints << endl; + if (m_traceLevel >= Trace::TRACE_CONSOLE) + std::cout << "nFilterPoints = " << m_nFilterPoints << endl; #endif } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - + if (m_nFilterPoints % 2) { // Odd - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); } else { // Even - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; - m_dFilterMax -= m_dFilterInc; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; + m_dFilterMax -= m_dFilterInc; } - + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); 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 (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] *= 0.5; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] *= 0.5; } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { - 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) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - m_adFilter[i] *= dScale; - } + 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) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + m_adFilter[i] *= dScale; + } } #ifdef HAVE_SGP EZPlot* pEZPlot = NULL; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot = new EZPlot (*pSGP); - pEZPlot->ezset ("title Filter Filter: Natural Order"); - pEZPlot->ezset ("ylength 0.50"); - pEZPlot->ezset ("yporigin 0.00"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); + pEZPlot = new EZPlot; + pEZPlot->ezset ("title Filter Filter: Natural Order"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.00"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); } #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"); - pEZPlot->ezset ("yporigin 0.50"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - delete pEZPlot; + pEZPlot->ezset ("title Filter Filter: Fourier Order"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.50"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; } #endif } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { @@ -323,17 +323,17 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1); m_nFilterPoints = nSpatialPoints; if (m_iZeropad > 0) { - double logBase2 = log(nSpatialPoints) / log(2); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_iZeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; + double logBase2 = log(nSpatialPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; #ifdef DEBUG if (m_traceLevel >= Trace::TRACE_CONSOLE) - std::cout << "nFilterPoints = " << m_nFilterPoints << endl; + std::cout << "nFilterPoints = " << m_nFilterPoints << endl; #endif double* adSpatialFilter = new double [m_nFilterPoints]; SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); @@ -341,48 +341,48 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw #ifdef HAVE_SGP EZPlot* pEZPlot = NULL; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot = new EZPlot (*pSGP); - pEZPlot->ezset ("title Spatial Filter: Natural Order"); - pEZPlot->ezset ("ylength 0.50"); - pEZPlot->ezset ("yporigin 0.00"); - pEZPlot->addCurve (adSpatialFilter, nSpatialPoints); - pEZPlot->plot(); - delete pEZPlot; + pEZPlot = new EZPlot; + pEZPlot->ezset ("title Spatial Filter: Natural Order"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.00"); + pEZPlot->addCurve (adSpatialFilter, nSpatialPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; } #endif if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - for (i = 0; i < m_nFilterPoints; i++) - adSpatialFilter[i] *= 0.5; + for (i = 0; i < m_nFilterPoints; i++) + adSpatialFilter[i] *= 0.5; } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { - 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) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - adSpatialFilter[i] *= dScale; - } + 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) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + adSpatialFilter[i] *= dScale; + } } for (i = nSpatialPoints; i < m_nFilterPoints; i++) - adSpatialFilter[i] = 0; - + adSpatialFilter[i] = 0; + m_adFilter = new double [m_nFilterPoints]; std::complex* acInverseFilter = new std::complex [m_nFilterPoints]; finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1); - delete adSpatialFilter; - for (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc; - delete acInverseFilter; + delete adSpatialFilter; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] = std::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; + pEZPlot->ezset ("title Spatial Filter: Inverse"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.50"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; } #endif } @@ -401,13 +401,13 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw angle += angleIncrement; } } - + #if HAVE_FFTW if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) { for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft m_adFilter[i] /= m_nFilterPoints; } - + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE); m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); @@ -431,23 +431,23 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw ProcessSignal::~ProcessSignal (void) { - delete [] m_adFourierSinTable; - delete [] m_adFourierCosTable; - delete [] m_adFilter; - + delete [] m_adFourierSinTable; + delete [] m_adFourierCosTable; + delete [] m_adFilter; + #if HAVE_FFTW - if (m_idFilterMethod == FILTER_METHOD_FFTW) { - fftw_destroy_plan(m_complexPlanForward); - fftw_destroy_plan(m_complexPlanBackward); - delete [] m_adComplexFftInput; - delete [] m_adComplexFftSignal; - } - if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - rfftw_destroy_plan(m_realPlanForward); - rfftw_destroy_plan(m_realPlanBackward); - delete [] m_adRealFftInput; - delete [] m_adRealFftSignal; - } + if (m_idFilterMethod == FILTER_METHOD_FFTW) { + fftw_destroy_plan(m_complexPlanForward); + fftw_destroy_plan(m_complexPlanBackward); + delete [] m_adComplexFftInput; + delete [] m_adComplexFftSignal; + } + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { + rfftw_destroy_plan(m_realPlanForward); + rfftw_destroy_plan(m_realPlanBackward); + delete [] m_adRealFftInput; + delete [] m_adRealFftSignal; + } #endif } @@ -455,24 +455,24 @@ int ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName) { int fmID = FILTER_METHOD_INVALID; - + for (int i = 0; i < s_iFilterMethodCount; i++) - if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { + if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { fmID = i; break; } - - return (fmID); + + return (fmID); } const char * ProcessSignal::convertFilterMethodIDToName (const int fmID) { static const char *name = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodName [fmID]); - + return (s_aszFilterMethodName [fmID]); + return (name); } @@ -480,10 +480,10 @@ const char * ProcessSignal::convertFilterMethodIDToTitle (const int fmID) { static const char *title = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodTitle [fmID]); - + return (s_aszFilterMethodTitle [fmID]); + return (title); } @@ -492,24 +492,24 @@ int ProcessSignal::convertFilterGenerationNameToID (const char* const fgName) { int fgID = FILTER_GENERATION_INVALID; - + for (int i = 0; i < s_iFilterGenerationCount; i++) - if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { + if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { fgID = i; break; } - - return (fgID); + + return (fgID); } const char * ProcessSignal::convertFilterGenerationIDToName (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) - return (s_aszFilterGenerationName [fgID]); - + return (s_aszFilterGenerationName [fgID]); + return (name); } @@ -517,10 +517,10 @@ const char * ProcessSignal::convertFilterGenerationIDToTitle (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) - return (s_aszFilterGenerationTitle [fgID]); - + return (s_aszFilterGenerationTitle [fgID]); + return (name); } @@ -531,7 +531,7 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const int i; for (i = 0; i < m_nSignalPoints; i++) input[i] = constInput[i]; - + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { for (int i = 0; i < m_nSignalPoints; i++) { int iDetFromCenter = i - (m_nSignalPoints / 2); @@ -544,8 +544,8 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const } } if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { - for (i = 0; i < m_nSignalPoints; i++) - output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints); + 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 = new double [m_nFilterPoints]; for (i = 0; i < m_nSignalPoints; i++) @@ -554,15 +554,15 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const inputSignal[i] = 0; // zeropad std::complex* fftSignal = new std::complex [m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); - delete inputSignal; + delete inputSignal; for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; double* inverseFourier = new double [m_nFilterPoints]; finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1); - delete fftSignal; + delete fftSignal; for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; - delete inverseFourier; + delete inverseFourier; } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { double* inputSignal = new double [m_nFilterPoints]; for (i = 0; i < m_nSignalPoints; i++) @@ -571,50 +571,50 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const inputSignal[i] = 0; // zeropad std::complex* fftSignal = new std::complex [m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, -1); - delete inputSignal; + delete inputSignal; for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; double* inverseFourier = new double [m_nFilterPoints]; finiteFourierTransform (fftSignal, inverseFourier, 1); - delete fftSignal; + delete fftSignal; for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; - delete inverseFourier; + delete inverseFourier; } #if HAVE_FFTW else if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - for (i = 0; i < m_nSignalPoints; i++) - m_adRealFftInput[i] = input[i]; - - fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ]; - rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); - 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++) + for (i = 0; i < m_nSignalPoints; i++) + m_adRealFftInput[i] = input[i]; + + fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ]; + rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); + 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 = new fftw_real [ m_nOutputPoints ]; - rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); - for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i]; - delete [] ifftOutput; + + fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ]; + rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = ifftOutput[i]; + delete [] ifftOutput; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - for (i = 0; i < m_nSignalPoints; i++) - m_adComplexFftInput[i].re = input[i]; - - fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ]; - fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); - 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; - } - delete [] fftOutput; - fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ]; - fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); - for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i].re; - delete [] ifftOutput; + for (i = 0; i < m_nSignalPoints; i++) + m_adComplexFftInput[i].re = input[i]; + + fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ]; + fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); + 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; + } + delete [] fftOutput; + fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ]; + fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = ifftOutput[i].re; + delete [] ifftOutput; } #endif delete input; @@ -622,36 +622,36 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const /* NAME - * convolve Discrete convolution of two functions - * - * SYNOPSIS - * r = convolve (f1, f2, dx, n, np, func_type) - * double r Convolved result - * double f1[], f2[] Functions to be convolved - * double dx Difference between successive x values - * int n Array index to center convolution about - * int np Number of points in f1 array - * int func_type EVEN or ODD or EVEN_AND_ODD function f2 - * - * NOTES - * f1 is the projection data, its indices range from 0 to np - 1. - * The index for f2, the filter, ranges from -(np-1) to (np-1). - * There are 3 ways to handle the negative vertices of f2: - * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. - * All filters used in reconstruction are even. - * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] - * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 - * for negative indices. Since f2 must range from -(np-1) to (np-1), - * if we add (np - 1) to f2's array index, then f2's index will - * range from 0 to 2 * (np - 1), and the origin, x = 0, will be - * stored at f2[np-1]. - */ +* convolve Discrete convolution of two functions +* +* SYNOPSIS +* r = convolve (f1, f2, dx, n, np, func_type) +* double r Convolved result +* double f1[], f2[] Functions to be convolved +* double dx Difference between successive x values +* int n Array index to center convolution about +* int np Number of points in f1 array +* int func_type EVEN or ODD or EVEN_AND_ODD function f2 +* +* NOTES +* f1 is the projection data, its indices range from 0 to np - 1. +* The index for f2, the filter, ranges from -(np-1) to (np-1). +* There are 3 ways to handle the negative vertices of f2: +* 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. +* All filters used in reconstruction are even. +* 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] +* 3. If f2 is both ODD AND EVEN, then we must store the value of f2 +* for negative indices. Since f2 must range from -(np-1) to (np-1), +* if we add (np - 1) to f2's array index, then f2's index will +* range from 0 to 2 * (np - 1), and the origin, x = 0, will be +* stored at f2[np-1]. +*/ double ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const { double sum = 0.0; - + #if UNOPTIMIZED_CONVOLUTION for (int i = 0; i < np; i++) sum += func[i] * m_adFilter[n - i + (np - 1)]; @@ -660,7 +660,7 @@ ProcessSignal::convolve (const double func[], const double dx, const int n, cons for (int i = 0; i < np; i++) sum += *func++ * *f2--; #endif - + return (sum * dx); } @@ -669,16 +669,16 @@ double ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const { double sum = 0.0; - + #if UNOPTIMIZED_CONVOLUTION -for (int i = 0; i < np; i++) - sum += func[i] * m_adFilter[n - i + (np - 1)]; + for (int i = 0; i < np; i++) + sum += func[i] * m_adFilter[n - i + (np - 1)]; #else -double* f2 = m_adFilter + n + (np - 1); -for (int i = 0; i < np; i++) - sum += *func++ * *f2--; + double* f2 = m_adFilter + n + (np - 1); + for (int i = 0; i < np; i++) + sum += *func++ * *f2--; #endif - + return (sum * dx); } @@ -686,12 +686,12 @@ for (int i = 0; i < np; i++) void ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction) { - std::complex* complexOutput = new std::complex [n]; - - finiteFourierTransform (input, complexOutput, n, direction); - for (int i = 0; i < n; i++) - output[i] = complexOutput[i].real(); - delete [] complexOutput; + std::complex* complexOutput = new std::complex [n]; + + finiteFourierTransform (input, complexOutput, n, direction); + for (int i = 0; i < n; i++) + output[i] = complexOutput[i].real(); + delete [] complexOutput; } void @@ -701,7 +701,7 @@ ProcessSignal::finiteFourierTransform (const double input[], std::complex input[], std:: direction = -1; else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { std::complex sum (0,0); @@ -750,10 +750,10 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl direction = -1; else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { - double sumReal = 0; + double sumReal = 0; for (int j = 0; j < n; j++) { double angle = i * j * angleIncrement; sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle); @@ -774,17 +774,17 @@ ProcessSignal::finiteFourierTransform (const double input[], std::complex 0) { - sumReal += input[j] * m_adFourierCosTable[tableIndex]; - sumImag += input[j] * m_adFourierSinTable[tableIndex]; + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag += input[j] * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j] * m_adFourierCosTable[tableIndex]; - sumImag -= input[j] * m_adFourierSinTable[tableIndex]; + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag -= input[j] * m_adFourierSinTable[tableIndex]; } } if (direction < 0) { @@ -803,21 +803,21 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], std:: direction = -1; else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { double sumReal = 0, sumImag = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * m_adFourierSinTable[tableIndex]; - sumImag += input[j].real() * m_adFourierSinTable[tableIndex] - + input[j].imag() * m_adFourierCosTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * -m_adFourierSinTable[tableIndex]; - sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] - + input[j].imag() * m_adFourierCosTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; } } if (direction < 0) { @@ -835,17 +835,17 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl direction = -1; else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { - double sumReal = 0; + double sumReal = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; } } if (direction < 0) { @@ -862,58 +862,112 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1) // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1 -void -ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) -{ +void +ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) +{ double* pdTemp = new double [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } else { // Even - int iHalfN = n / 2; - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - - -void -ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) -{ + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + +void +ProcessSignal::shuffleNaturalToFourierOrder (std::complex* pdVector, const int n) +{ + std::complex* pdTemp = new std::complex [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete [] pdTemp; +} + + +void +ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) +{ double* pdTemp = new double [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN + 1]; - } else { // Even - int iHalfN = n / 2; - pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i+1]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + +void +ProcessSignal::shuffleFourierToNaturalOrder (std::complex* pdVector, const int n) +{ + std::complex* pdTemp = new std::complex [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete [] pdTemp; +} +