/*****************************************************************************
** File IDENTIFICATION
**
-** Name: filter.cpp
-** Purpose: Routines for signal-procesing filters
-** Progammer: Kevin Rosenberg
-** Date Started: Aug 1984
+** Name: procsignal.cpp
+** Purpose: Routines for processing signals and projections
+** Progammer: Kevin Rosenberg
+** 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.16 2001/01/12 04:28:37 kevin Exp $
+** $Id: procsignal.cpp,v 1.32 2001/03/30 19:17:32 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
#include "ct.h"
+#ifdef HAVE_WXWINDOWS
+#include "nographics.h"
+#endif
+
// FilterMethod ID/Names
const int ProcessSignal::FILTER_METHOD_INVALID = -1;
const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
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"},
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"},
};
// 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) {
}
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;
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;
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) {
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++) {
} 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)
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<int>(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);
}
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 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;
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);
-#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;
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;
m_adFilter = new double [m_nFilterPoints];
std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
-#define PRE_JAN_2001 1
-#ifdef PRE_JAN_2001
finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
delete adSpatialFilter;
for (i = 0; i < m_nFilterPoints; i++)
m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
delete acInverseFilter;
-#else
- finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, FORWARD);
- delete adSpatialFilter;
- for (i = 0; i < m_nFilterPoints; i++)
- m_adFilter[i] = std::abs(acInverseFilter[i]);
- delete acInverseFilter;
-#endif
-
-#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
}
}
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_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[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_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++)
}
}
+
+int
+ProcessSignal::addZeropadFactor (int n, int iZeropad)
+{
+ if (iZeropad > 0) {
+ double dLogBase2 = log(n) / log(2);
+ int iLogBase2 = static_cast<int>(floor (dLogBase2));
+ int iPaddedN = 1 << (iLogBase2 + iZeropad);
+#ifdef DEBUG
+ sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);
+#endif
+ return iPaddedN;
+ }
+
+ return n;
+}
projects / ctsim.git / blobdiff