X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;fp=libctsim%2Fprocsignal.cpp;h=be88ecc31e34cee618746eddd3537a4f6fe85d98;hp=2206ae1838cf41ebd3f2e13eb34353fdb060670e;hb=1a050c98763fbbc0662731b0b76953acede6f5d7;hpb=c8b19dfaffba9f06d8b6c40cb1bb83a8964867f7 diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 2206ae1..be88ecc 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -1,9 +1,9 @@ /***************************************************************************** ** 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 @@ -81,9 +81,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 +115,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 +136,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 +147,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 +168,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 +200,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 +253,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 +271,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 +285,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 +317,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 +341,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 +352,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 +378,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 +396,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,13 +410,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_adRealFftInput = static_cast(fftw_malloc (sizeof(double) * m_nFilterPoints)); m_adRealFftOutput = static_cast(fftw_malloc (sizeof(double) * m_nFilterPoints)); @@ -424,7 +424,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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++) + for (i = 0; i < m_nFilterPoints; i++) m_adRealFftInput[i] = 0; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { m_adComplexFftInput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints)); @@ -434,13 +434,13 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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++) + for (i = 0; i < m_nFilterPoints; i++) m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0; - for (i = 0; i < m_nOutputPoints; i++) + for (i = 0; i < m_nOutputPoints; i++) m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0; } #endif - + } ProcessSignal::~ProcessSignal (void) @@ -448,7 +448,7 @@ 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); @@ -473,13 +473,13 @@ 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) { fmID = i; break; } - + return (fmID); } @@ -487,10 +487,10 @@ const char * ProcessSignal::convertFilterMethodIDToName (const int fmID) { static const char *name = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) return (s_aszFilterMethodName [fmID]); - + return (name); } @@ -498,10 +498,10 @@ const char * ProcessSignal::convertFilterMethodIDToTitle (const int fmID) { static const char *title = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) return (s_aszFilterMethodTitle [fmID]); - + return (title); } @@ -510,13 +510,13 @@ 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) { fgID = i; break; } - + return (fgID); } @@ -524,10 +524,10 @@ const char * ProcessSignal::convertFilterGenerationIDToName (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) return (s_aszFilterGenerationName [fgID]); - + return (name); } @@ -535,10 +535,10 @@ const char * ProcessSignal::convertFilterGenerationIDToTitle (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) return (s_aszFilterGenerationTitle [fgID]); - + return (name); } @@ -549,7 +549,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); @@ -578,7 +578,7 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const 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) { @@ -595,7 +595,7 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const 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; } @@ -603,27 +603,27 @@ 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_execute (m_realPlanForward); for (i = 0; i < m_nFilterPoints; i++) m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i]; for (i = m_nFilterPoints; i < m_nOutputPoints; i++) - m_adRealFftSignal[i] = 0; - + m_adRealFftSignal[i] = 0; + fftw_execute (m_realPlanBackward); for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) output[i] = m_adRealFftBackwardOutput[i]; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { for (i = 0; i < m_nSignalPoints; i++) m_adComplexFftInput[i][0] = input[i]; - + fftw_execute (m_complexPlanForward); for (i = 0; i < m_nFilterPoints; i++) { m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0]; m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1]; } fftw_execute (m_complexPlanBackward); - for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) output[i] = m_adComplexFftBackwardOutput[i][0]; } #endif @@ -632,36 +632,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)]; @@ -670,16 +670,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)]; @@ -688,7 +688,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); } @@ -697,7 +697,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(); @@ -709,9 +709,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); @@ -758,9 +758,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; @@ -782,9 +782,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]; @@ -843,18 +843,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]; } }