X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=2d6600d37e6b60767bd416e05c216a494416ec24;hp=c40e2972ad606290afef9a33ec7dd0f9a8f96355;hb=d16eb37cbc73f67fc29a60645e0b1ac7fe32767e;hpb=1f05afcd917bbf8b56a9446e4f937864f09b33f4 diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index c40e297..2d6600d 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -7,9 +7,9 @@ ** Date Started: Aug 1984 ** ** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg +** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: procsignal.cpp,v 1.15 2001/01/12 03:49:07 kevin Exp $ +** $Id: procsignal.cpp,v 1.29 2001/03/18 18:08:25 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 @@ -27,6 +27,10 @@ #include "ct.h" +#ifdef HAVE_WXWINDOWS +#include "dlgezplot.h" +#endif + // FilterMethod ID/Names const int ProcessSignal::FILTER_METHOD_INVALID = -1; const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0; @@ -37,17 +41,17 @@ 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[] = { +const char* const ProcessSignal::s_aszFilterMethodName[] = { {"convolution"}, {"fourier"}, - {"fouier_table"}, + {"fouier-table"}, {"fft"}, #if HAVE_FFTW {"fftw"}, {"rfftw"}, #endif }; -const char* ProcessSignal::s_aszFilterMethodTitle[] = { +const char* const ProcessSignal::s_aszFilterMethodTitle[] = { {"Convolution"}, {"Fourier"}, {"Fouier Trigometric Table"}, @@ -63,11 +67,11 @@ 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[] = { +const char* const ProcessSignal::s_aszFilterGenerationName[] = { {"direct"}, - {"inverse_fourier"}, + {"inverse-fourier"}, }; -const char* ProcessSignal::s_aszFilterGenerationTitle[] = { +const char* const ProcessSignal::s_aszFilterGenerationTitle[] = { {"Direct"}, {"Inverse Fourier"}, }; @@ -78,10 +82,10 @@ 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) + 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); if (m_idFilterMethod == FILTER_METHOD_INVALID) { @@ -113,15 +117,16 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth } init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, - m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP); + 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) + 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; @@ -130,6 +135,7 @@ 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; @@ -145,12 +151,12 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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) { @@ -187,46 +193,40 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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 (); - pEZPlot->ezset ("title Filter Response: Natural Order"); - pEZPlot->ezset ("ylength 0.25"); - pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); - pEZPlot->plot (pSGP); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + 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 Fourier::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 (pSGP); + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } #endif ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD); 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 (pSGP); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } #endif Fourier::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 (pSGP); - delete pEZPlot; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } #endif for (i = 0; i < m_nFilterPoints; i++) { @@ -239,14 +239,18 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw } 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 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 Fourier Frequency: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif } // if (geometry) } // if (spatial filtering) @@ -254,22 +258,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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; -#ifdef DEBUG - if (m_traceLevel >= Trace::TRACE_CONSOLE) - std::cout << "nFilterPoints = " << m_nFilterPoints << endl; -#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); @@ -280,52 +272,54 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterMax -= m_dFilterInc; } - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + 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 - // Jan 2001: Direct seems to work for equilinear and equiangular - // however, inverse_fourier doesn't work for equiangular on all versions of CTSim tested - +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + 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 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); - double sinScale = sin (iDetFromZero * m_dSignalInc); - if (fabs(sinScale) < 1E-7) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc); 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; - 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); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } #endif Fourier::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 (pSGP); - delete pEZPlot; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + 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; @@ -342,32 +336,29 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_nFilterPoints = 1 << nextPowerOf2; } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; -#ifdef DEBUG +#if defined(DEBUG) || defined(_DEBUG) if (m_traceLevel >= Trace::TRACE_CONSOLE) - std::cout << "nFilterPoints = " << m_nFilterPoints << endl; + 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, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, + m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints); -#ifdef HAVE_SGP - EZPlot* pEZPlot = NULL; - if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot = new EZPlot; - 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; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot;; + dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints); + dlgEZPlot.ShowModal(); } #endif - + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - for (i = 0; i < m_nFilterPoints; i++) + for (i = 0; i < nSpatialPoints; 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); + for (i = 0; i < nSpatialPoints; i++) { + int iDetFromZero = i - ((nSpatialPoints - 1) / 2); double sinScale = sin (iDetFromZero * m_dSignalInc); if (fabs(sinScale) < 1E-7) sinScale = 1; @@ -377,6 +368,14 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw adSpatialFilter[i] *= dScale; } } +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot;; + dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints); + dlgEZPlot.ShowModal(); + } +#endif for (i = nSpatialPoints; i < m_nFilterPoints; i++) adSpatialFilter[i] = 0; @@ -387,15 +386,12 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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 (pSGP); - delete pEZPlot; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } #endif } @@ -868,3 +864,17 @@ 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)); + if (dLogBase2 != floor(dLogBase2)) + iLogBase2++; // raise up to next power of 2 + n = 1 << (iLogBase2 + (iZeropad - 1)); + } + + return n; +}