X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=3a28676ff606bbb94c5f4b9753f8b189d2c5268d;hp=d6701ca41ae9aabacd8a638f4752189a76287ec9;hb=24e04db129360d73d3a43f024108086f51a2e45d;hpb=728e9fcbe0b785e56e7dcd2bd305786c647050af diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index d6701ca..3a28676 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -1,15 +1,13 @@ /***************************************************************************** ** File IDENTIFICATION -** +** ** Name: procsignal.cpp ** Purpose: Routines for processing signals and projections -** Progammer: Kevin Rosenberg +** Progammer: Kevin Rosenberg ** Date Started: Aug 1984 ** ** This is part of the CTSim program -** Copyright (c) 1983-2001 Kevin Rosenberg -** -** $Id: procsignal.cpp,v 1.31 2001/03/29 21:25:49 kevin Exp $ +** 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 @@ -42,23 +40,23 @@ const int ProcessSignal::FILTER_METHOD_FFTW = 4; const int ProcessSignal::FILTER_METHOD_RFFTW =5 ; #endif const char* const ProcessSignal::s_aszFilterMethodName[] = { - {"convolution"}, - {"fourier"}, - {"fouier-table"}, - {"fft"}, + "convolution", + "fourier", + "fouier-table", + "fft", #if HAVE_FFTW - {"fftw"}, - {"rfftw"}, + "fftw", + "rfftw", #endif }; const char* const ProcessSignal::s_aszFilterMethodTitle[] = { - {"Convolution"}, - {"Fourier"}, - {"Fouier Trigometric Table"}, - {"FFT"}, + "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*); @@ -68,12 +66,12 @@ const int ProcessSignal::FILTER_GENERATION_INVALID = -1; const int ProcessSignal::FILTER_GENERATION_DIRECT = 0; const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1; const char* const ProcessSignal::s_aszFilterGenerationName[] = { - {"direct"}, - {"inverse-fourier"}, + "direct", + "inverse-fourier", }; const char* const ProcessSignal::s_aszFilterGenerationTitle[] = { - {"Direct"}, - {"Inverse Fourier"}, + "Direct", + "Inverse Fourier", }; const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*); @@ -81,9 +79,9 @@ 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, +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) { @@ -115,17 +113,17 @@ 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, + + 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, +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; @@ -136,7 +134,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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; @@ -147,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 adjust for imaginary detector through origin of phantom + + // 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) { 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; @@ -168,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); @@ -200,7 +198,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints); dlgEZPlot.ShowModal(); } -#endif +#endif Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); #ifdef HAVE_SGP if (g_bRunningWXWindows && m_traceLevel > 0) { @@ -253,14 +251,14 @@ 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 = addZeropadFactor (m_nSignalPoints, m_iZeropad); m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - + if (isOdd (m_nFilterPoints)) { // Odd m_dFilterMin = -1. / (2 * m_dSignalInc); m_dFilterMax = 1. / (2 * m_dSignalInc); @@ -271,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; @@ -285,7 +283,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.ShowModal(); } #endif - + // 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 @@ -317,9 +315,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.ShowModal(); } #endif - + // FILTERING: FREQUENCY - INVERSE FOURIER - + } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { // calculate number of filter points with zeropadding int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1; @@ -341,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)) @@ -352,7 +350,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.ShowModal(); } #endif - + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { for (i = 0; i < nSpatialPoints; i++) adSpatialFilter[i] *= 0.5; @@ -378,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); @@ -396,7 +394,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw #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; @@ -410,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 | FFTW_USE_WISDOM); - m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_USE_WISDOM); - 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 | FFTW_USE_WISDOM); - m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM); - 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) @@ -443,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 } @@ -464,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); } @@ -489,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); } @@ -501,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); } @@ -526,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); } @@ -538,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) { @@ -564,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) { @@ -581,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; } @@ -594,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; @@ -631,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)]; @@ -669,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)]; @@ -687,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); } @@ -696,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(); @@ -708,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); @@ -757,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; @@ -781,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]; @@ -842,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]; } } @@ -871,9 +892,11 @@ ProcessSignal::addZeropadFactor (int n, int iZeropad) if (iZeropad > 0) { double dLogBase2 = log(n) / log(2); int iLogBase2 = static_cast(floor (dLogBase2)); - if (dLogBase2 != floor(dLogBase2)) - iLogBase2++; // raise up to next power of 2 - n = 1 << (iLogBase2 + (iZeropad - 1)); + int iPaddedN = 1 << (iLogBase2 + iZeropad); +#ifdef DEBUG + sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN); +#endif + return iPaddedN; } return n;