X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;fp=libctsim%2Fprocsignal.cpp;h=b717f8037e104bebc5a453e4571249e5f81173e3;hp=67641fd52f5b306fde9469c25845c107bb749c07;hb=ed576a069b0ba9de2e9153707012eba620ac606e;hpb=7f8f356151b0c8db0dbbf1c1896cc22630d6c774 diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 67641fd..b717f80 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: procsignal.cpp,v 1.19 2001/01/12 16:41:56 kevin Exp $ +** $Id: procsignal.cpp,v 1.20 2001/01/12 21:53:27 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 @@ -82,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, SGP* pSGP) + : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) { m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); if (m_idFilterMethod == FILTER_METHOD_INVALID) { @@ -117,15 +117,15 @@ 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, 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, SGP* pSGP) { int i; m_idFilter = idFilter; @@ -191,46 +191,41 @@ 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)) + 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 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++) { @@ -243,14 +238,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) @@ -266,9 +265,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw nextPowerOf2++; nextPowerOf2 += (m_iZeropad - 1); m_nFilterPoints = 1 << nextPowerOf2; -#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 } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; @@ -285,52 +284,54 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw } SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, - m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + 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 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 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; @@ -347,33 +348,31 @@ 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); + m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints); -#if defined(HAVE_WXWINDOWS) && defined(DEBUG) - EZPlotDialog pEZPlotDlg = NULL; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) if (g_bRunningWXWindows && m_traceLevel > 0) { - pEZPlotDlg = new EZPlotDialog; - pEZPlot->getEZPlot()->ezset ("title Spatial Filter: Natural Order"); - pEZPlot->getEZPlot()->ezset ("ylength 0.50"); - pEZPlot->getEZPlot()->ezset ("yporigin 0.00"); - pEZPlot->getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints); + EZPlotDialog dlgEZPlot;; + dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints); + dlgEZPlot.ShowModal(); } #endif - -// #define PRE_JAN_2001 1 + +#define PRE_JAN_2001 1 #ifdef PRE_JAN_2001 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; @@ -383,6 +382,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; @@ -393,46 +400,74 @@ 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; +#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 + #else for (i = nSpatialPoints; i < m_nFilterPoints; i++) adSpatialFilter[i] = 0; - -// for (i = 0; i < m_nFilterPoints; i++) -// adSpatialFilter[i] /= m_dSignalInc; - + 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]); + 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); - double sinScale = sin (iDetFromZero * m_dSignalInc); - if (fabs(sinScale) < 1E-7) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - m_adFilter[i] *= dScale; + if (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) - if (g_bRunningWXWindows && pEZPlotDlg && m_traceLevel > 0) { - pEZPlotDlg->getEZPlot()->ezset ("title Spatial Filter: Inverse"); - pEZPlotDlg->getEZPlot()->ezset ("ylength 0.50"); - pEZPlotDlg->getEZPlot()->ezset ("yporigin 0.50"); - pEZPlotDlg->getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); - pEZPlotDlg->ShowModal(); - delete pEZPlotDlg; +#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 }