X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=3a28676ff606bbb94c5f4b9753f8b189d2c5268d;hp=b717f8037e104bebc5a453e4571249e5f81173e3;hb=24e04db129360d73d3a43f024108086f51a2e45d;hpb=ed576a069b0ba9de2e9153707012eba620ac606e diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index b717f80..3a28676 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -1,15 +1,13 @@ /***************************************************************************** ** File IDENTIFICATION -** -** Name: filter.cpp -** Purpose: Routines for signal-procesing filters -** Progammer: Kevin Rosenberg -** Date Started: Aug 1984 ** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg +** Name: procsignal.cpp +** Purpose: Routines for processing signals and projections +** Progammer: Kevin Rosenberg +** Date Started: Aug 1984 ** -** $Id: procsignal.cpp,v 1.20 2001/01/12 21:53:27 kevin Exp $ +** This is part of the CTSim program +** Copyright (c) 1983-2009 Kevin Rosenberg ** ** 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 @@ -28,7 +26,7 @@ #include "ct.h" #ifdef HAVE_WXWINDOWS -#include "../src/dlgezplot.h" +#include "nographics.h" #endif // FilterMethod ID/Names @@ -41,24 +39,24 @@ const int ProcessSignal::FILTER_METHOD_FFT = 3; const int ProcessSignal::FILTER_METHOD_FFTW = 4; const int ProcessSignal::FILTER_METHOD_RFFTW =5 ; #endif -const char* ProcessSignal::s_aszFilterMethodName[] = { - {"convolution"}, - {"fourier"}, - {"fouier_table"}, - {"fft"}, +const char* const ProcessSignal::s_aszFilterMethodName[] = { + "convolution", + "fourier", + "fouier-table", + "fft", #if HAVE_FFTW - {"fftw"}, - {"rfftw"}, + "fftw", + "rfftw", #endif }; -const char* ProcessSignal::s_aszFilterMethodTitle[] = { - {"Convolution"}, - {"Fourier"}, - {"Fouier Trigometric Table"}, - {"FFT"}, +const char* const ProcessSignal::s_aszFilterMethodTitle[] = { + "Convolution", + "Fourier", + "Fouier Trigometric Table", + "FFT", #if HAVE_FFTW - {"FFTW"}, - {"Real/Half-Complex FFTW"}, + "FFTW", + "Real/Half-Complex FFTW", #endif }; const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*); @@ -67,13 +65,13 @@ const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / const int ProcessSignal::FILTER_GENERATION_INVALID = -1; const int ProcessSignal::FILTER_GENERATION_DIRECT = 0; const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1; -const char* ProcessSignal::s_aszFilterGenerationName[] = { - {"direct"}, - {"inverse_fourier"}, +const char* const ProcessSignal::s_aszFilterGenerationName[] = { + "direct", + "inverse-fourier", }; -const char* ProcessSignal::s_aszFilterGenerationTitle[] = { - {"Direct"}, - {"Inverse Fourier"}, +const char* const ProcessSignal::s_aszFilterGenerationTitle[] = { + "Direct", + "Inverse Fourier", }; const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*); @@ -81,10 +79,10 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration // CLASS IDENTIFICATION // 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) +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, double dSourceDetectorLength, SGP* pSGP) : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) { m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); @@ -115,17 +113,18 @@ 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); + + init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, + m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, + dSourceDetectorLength, pSGP); } 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) +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, double dSourceDetectorLength, SGP* pSGP) { int i; m_idFilter = idFilter; @@ -134,7 +133,8 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_idFilterGeneration = idFilterGeneration; m_idGeometry = iGeometry; m_dFocalLength = dFocalLength; - + m_dSourceDetectorLength = dSourceDetectorLength; + 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; @@ -145,18 +145,18 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dBandwidth = dBandwidth; m_nSignalPoints = nSignalPoints; m_dSignalInc = dSignalIncrement; - m_dFilterParam = dFilterParam; + m_dFilterParam = dFilterParam; m_iZeropad = iZeropad; m_iPreinterpolationFactor = iPreinterpolationFactor; - - // scale signalInc/BW to signalInc/2 to adjust for imaginary detector - // through origin of phantom rather than 2 times distance to detector, + + // scale signalInc/BW to adjust for imaginary detector through origin of phantom // see Kak-Slaney Fig 3.22, for Collinear diagram if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - m_dSignalInc /= 2; - m_dBandwidth *= 2; + double dEquilinearScale = m_dSourceDetectorLength / m_dFocalLength; + m_dSignalInc /= dEquilinearScale; + m_dBandwidth *= dEquilinearScale; } - + if (m_idFilterMethod == FILTER_METHOD_FFT) { #if HAVE_FFTW m_idFilterMethod = FILTER_METHOD_RFFTW; @@ -166,14 +166,14 @@ 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); @@ -192,14 +192,13 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw double* adFrequencyFilter = new double [m_nFilterPoints]; filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) - EZPlotDialog* pEZPlotDlg = NULL; if (g_bRunningWXWindows && m_traceLevel > 0) { EZPlotDialog dlgEZPlot; dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Natural Order"); dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints); dlgEZPlot.ShowModal(); } -#endif +#endif Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); #ifdef HAVE_SGP if (g_bRunningWXWindows && m_traceLevel > 0) { @@ -252,27 +251,15 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw #endif } // 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; -#if defined(DEBUG) || defined(_DEBUG) - if (m_traceLevel >= Trace::TRACE_CONSOLE) - sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints); -#endif - } + m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad); m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - - if (m_nFilterPoints % 2) { // Odd + + if (isOdd (m_nFilterPoints)) { // Odd m_dFilterMin = -1. / (2 * m_dSignalInc); m_dFilterMax = 1. / (2 * m_dSignalInc); m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); @@ -282,12 +269,12 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; m_dFilterMax -= m_dFilterInc; } - - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, + + 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); - + #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) if (g_bRunningWXWindows && m_traceLevel > 0) { EZPlotDialog dlgEZPlot; @@ -297,10 +284,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw } #endif - // This doesn't work: Need to add filtering for divergent geometries & Frequency/Direct filtering - // Jan 2001: Direct seems to work for equilinear and equiangular - // however, inverse_fourier doesn't work for equiangular on all versions of CTSim tested - // Scaling is done with data in frequency space, natural order + // This works fairly well. I'm not sure why since scaling for geometries is done on + // frequency filter rather than spatial filter as it should be. + // It gives values slightly off than freq/inverse filtering if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { for (i = 0; i < m_nFilterPoints; i++) m_adFilter[i] *= 0.5; @@ -353,7 +339,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints); #endif double* adSpatialFilter = new double [m_nFilterPoints]; - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints); #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) @@ -364,9 +350,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.ShowModal(); } #endif - -#define PRE_JAN_2001 1 -#ifdef PRE_JAN_2001 + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { for (i = 0; i < nSpatialPoints; i++) adSpatialFilter[i] *= 0.5; @@ -392,7 +376,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw #endif for (i = nSpatialPoints; i < m_nFilterPoints; i++) adSpatialFilter[i] = 0; - + m_adFilter = new double [m_nFilterPoints]; std::complex* acInverseFilter = new std::complex [m_nFilterPoints]; finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD); @@ -408,71 +392,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.ShowModal(); } #endif - -#else - for (i = nSpatialPoints; i < m_nFilterPoints; i++) - adSpatialFilter[i] = 0; - - std::complex* acInverseFilter = new std::complex [m_nFilterPoints]; - finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, FORWARD); - delete adSpatialFilter; - m_adFilter = new double [m_nFilterPoints]; - for (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] = std::abs(acInverseFilter[i]); - delete acInverseFilter; - -#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) - if (g_bRunningWXWindows && m_traceLevel > 0) { - EZPlotDialog dlgEZPlot; - dlgEZPlot.getEZPlot()->ezset ("title Inverse Spatial Filter: Fourier order"); - dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); - dlgEZPlot.ShowModal(); - } -#endif - Fourier::shuffleFourierToNaturalOrder(m_adFilter, m_nFilterPoints); -#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) - if (g_bRunningWXWindows && m_traceLevel > 0) { - EZPlotDialog dlgEZPlot; - dlgEZPlot.getEZPlot()->ezset ("title Inverse Spatial Filter: Natural order"); - dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); - dlgEZPlot.ShowModal(); - } -#endif - if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - 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); - if (abs(iDetFromZero) < m_nSignalPoints) { - double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc); - double dScale = 0.5 * sinScale * sinScale; - m_adFilter[i] *= dScale; - } - } - } -#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) - if (g_bRunningWXWindows && m_traceLevel > 0) { - EZPlotDialog dlgEZPlot; - dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Spatial Filter: Natural order"); - dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); - dlgEZPlot.ShowModal(); - } -#endif - Fourier::shuffleNaturalToFourierOrder(m_adFilter, m_nFilterPoints); -#endif - -#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) - if (g_bRunningWXWindows && m_traceLevel > 0) { - EZPlotDialog dlgEZPlot; - dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter Inverse Post Geometry Filtering"); - dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); - dlgEZPlot.ShowModal(); - } -#endif } } - + // precalculate sin and cosine tables for fourier transform if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1; @@ -486,32 +408,37 @@ 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); - m_adRealFftInput = new fftw_real [ m_nFilterPoints ]; - m_adRealFftSignal = new fftw_real [ m_nOutputPoints ]; - for (i = 0; i < m_nFilterPoints; i++) + m_adRealFftInput = static_cast(fftw_malloc (sizeof(double) * m_nFilterPoints)); + m_adRealFftOutput = static_cast(fftw_malloc (sizeof(double) * m_nFilterPoints)); + m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE); + m_adRealFftSignal = static_cast(fftw_malloc (sizeof(double) * m_nOutputPoints)); + m_adRealFftBackwardOutput = static_cast(fftw_malloc (sizeof(double) * m_nOutputPoints)); + m_realPlanBackward = fftw_plan_r2r_1d (m_nOutputPoints, m_adRealFftSignal, m_adRealFftBackwardOutput, FFTW_HC2R, FFTW_ESTIMATE); + 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 (i = 0; i < m_nFilterPoints; i++) - m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0; - for (i = 0; i < m_nOutputPoints; i++) - m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0; + m_adComplexFftInput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints)); + m_adComplexFftOutput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints)); + m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD, FFTW_ESTIMATE); + m_adComplexFftSignal = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints)); + m_adComplexFftBackwardOutput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints)); + m_complexPlanBackward = fftw_plan_dft_1d (m_nOutputPoints, m_adComplexFftSignal, m_adComplexFftBackwardOutput, FFTW_BACKWARD, FFTW_ESTIMATE); + + for (i = 0; i < m_nFilterPoints; i++) + m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0; + for (i = 0; i < m_nOutputPoints; i++) + m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0; } #endif - + } ProcessSignal::~ProcessSignal (void) @@ -519,19 +446,23 @@ ProcessSignal::~ProcessSignal (void) 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; + fftw_free (m_adComplexFftInput); + fftw_free (m_adComplexFftOutput); + fftw_free (m_adComplexFftSignal); + fftw_free (m_adComplexFftBackwardOutput); } if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - rfftw_destroy_plan(m_realPlanForward); - rfftw_destroy_plan(m_realPlanBackward); - delete [] m_adRealFftInput; - delete [] m_adRealFftSignal; + fftw_destroy_plan(m_realPlanForward); + fftw_destroy_plan(m_realPlanBackward); + fftw_free (m_adRealFftInput); + fftw_free (m_adRealFftOutput); + fftw_free (m_adRealFftSignal); + fftw_free (m_adRealFftBackwardOutput); } #endif } @@ -540,24 +471,24 @@ int ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName) { int fmID = FILTER_METHOD_INVALID; - - for (int i = 0; i < s_iFilterMethodCount; i++) + + for (int i = 0; i < s_iFilterMethodCount; i++) { 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 (name); } @@ -565,10 +496,10 @@ const char * ProcessSignal::convertFilterMethodIDToTitle (const int fmID) { static const char *title = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) return (s_aszFilterMethodTitle [fmID]); - + return (title); } @@ -577,24 +508,24 @@ int ProcessSignal::convertFilterGenerationNameToID (const char* const fgName) { int fgID = FILTER_GENERATION_INVALID; - - for (int i = 0; i < s_iFilterGenerationCount; i++) + + for (int i = 0; i < s_iFilterGenerationCount; i++) { 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 (name); } @@ -602,10 +533,10 @@ const char * ProcessSignal::convertFilterGenerationIDToTitle (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) return (s_aszFilterGenerationTitle [fgID]); - + return (name); } @@ -614,21 +545,34 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const { double* input = new double [m_nSignalPoints]; int i; + +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (i = 0; i < m_nSignalPoints; i++) input[i] = constInput[i]; - + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (int i = 0; i < m_nSignalPoints; i++) { int iDetFromCenter = i - (m_nSignalPoints / 2); input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc); } } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (int i = 0; i < m_nSignalPoints; i++) { int iDetFromCenter = i - (m_nSignalPoints / 2); input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc); } } if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (i = 0; i < m_nSignalPoints; i++) output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints); } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) { @@ -640,12 +584,15 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const std::complex* fftSignal = new std::complex [m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD); delete inputSignal; +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; double* inverseFourier = new double [m_nFilterPoints]; finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD); delete fftSignal; - for (i = 0; i < m_nSignalPoints; i++) + for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; delete inverseFourier; } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { @@ -657,12 +604,15 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const std::complex* fftSignal = new std::complex [m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, FORWARD); delete inputSignal; +#if HAVE_OPENMP + #pragma omp parallel for +#endif for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; double* inverseFourier = new double [m_nFilterPoints]; finiteFourierTransform (fftSignal, inverseFourier, BACKWARD); delete fftSignal; - for (i = 0; i < m_nSignalPoints; i++) + for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; delete inverseFourier; } @@ -670,36 +620,31 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const 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); + + fftw_execute (m_realPlanForward); for (i = 0; i < m_nFilterPoints; i++) - m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i]; - delete [] fftOutput; + m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i]; 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); + m_adRealFftSignal[i] = 0; + + fftw_execute (m_realPlanBackward); for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i]; - delete [] ifftOutput; + output[i] = m_adRealFftBackwardOutput[i]; } 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); + m_adComplexFftInput[i][0] = input[i]; + + fftw_execute (m_complexPlanForward); +#if HAVE_OPENMP + #pragma omp parallel for +#endif 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; + m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0]; + m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1]; } - 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; + fftw_execute (m_complexPlanBackward); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = m_adComplexFftBackwardOutput[i][0]; } #endif delete input; @@ -707,36 +652,36 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const /* NAME -* convolve Discrete convolution of two functions +* 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 +* 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] +* 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]. +* 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 +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)]; @@ -745,16 +690,16 @@ 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); } -double +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)]; @@ -763,7 +708,7 @@ ProcessSignal::convolve (const float func[], const double dx, const int n, const for (int i = 0; i < np; i++) sum += *func++ * *f2--; #endif - + return (sum * dx); } @@ -772,7 +717,7 @@ 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(); @@ -784,9 +729,9 @@ ProcessSignal::finiteFourierTransform (const double input[], std::complex input[], std:: { if (direction < 0) direction = -1; - else + else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { std::complex sum (0,0); @@ -833,9 +778,9 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl { if (direction < 0) direction = -1; - else + else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { double sumReal = 0; @@ -857,9 +802,9 @@ ProcessSignal::finiteFourierTransform (const double input[], std::complex input[], std:: { if (direction < 0) direction = -1; - else + 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] + 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] + 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]; @@ -918,18 +863,18 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl { if (direction < 0) direction = -1; - else + else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { 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] + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - input[j].imag() * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - input[j].imag() * -m_adFourierSinTable[tableIndex]; } } @@ -940,3 +885,19 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl } } + +int +ProcessSignal::addZeropadFactor (int n, int iZeropad) +{ + if (iZeropad > 0) { + double dLogBase2 = log(n) / log(2); + int iLogBase2 = static_cast(floor (dLogBase2)); + int iPaddedN = 1 << (iLogBase2 + iZeropad); +#ifdef DEBUG + sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN); +#endif + return iPaddedN; + } + + return n; +}