** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: procsignal.cpp,v 1.10 2000/12/16 06:12:47 kevin Exp $
+** $Id: procsignal.cpp,v 1.11 2000/12/29 15:45:06 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
// 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)
+: m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
{
m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
if (m_idFilterMethod == FILTER_METHOD_INVALID) {
m_failMessage += szDomainName;
return;
}
-
+
init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
}
m_idFilterGeneration = idFilterGeneration;
m_idGeometry = iGeometry;
m_dFocalLength = dFocalLength;
-
+
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;
m_dFilterParam = dFilterParam;
m_iZeropad = iZeropad;
m_iPreinterpolationFactor = iPreinterpolationFactor;
-
+
// scale signalInc/BW to signalInc/2 to adjust for imaginary detector
// through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
m_dSignalInc /= 2;
m_dBandwidth *= 2;
}
-
+
if (m_idFilterMethod == FILTER_METHOD_FFT) {
#if HAVE_FFTW
m_idFilterMethod = FILTER_METHOD_RFFTW;
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);
- m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
- 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_SPATIAL);
- m_adFilter = new double[ m_nFilterPoints ];
- filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
+ m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
+ m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
+ m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
+ 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_SPATIAL);
+ m_adFilter = new double[ m_nFilterPoints ];
+ filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
- m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
- 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_adFilter = new double[ m_nFilterPoints ];
- double* adFrequencyFilter = new double [m_nFilterPoints];
- filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
+ m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
+ 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_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 (*pSGP);
- pEZPlot->ezset ("title Filter Response: Natural Order");
- pEZPlot->ezset ("ylength 0.25");
- pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
- pEZPlot->plot();
- }
+ 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);
+ }
#endif
- shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
+ 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();
- }
+ 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);
+ }
#endif
- ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
- delete adFrequencyFilter;\r
+ ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
+ delete adFrequencyFilter;\r
#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();
- }
+ 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);
+ }
#endif
- shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
+ 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();
- delete pEZPlot;
- }
+ 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;
+ }
#endif
- for (i = 0; i < m_nFilterPoints; i++) {
- m_adFilter[i] /= m_dSignalInc;
- }
+ for (i = 0; i < m_nFilterPoints; i++) {
+ m_adFilter[i] /= m_dSignalInc;
+ }
}
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
- for (i = 0; i < m_nFilterPoints; i++)
- m_adFilter[i] *= 0.5;
+ 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;
- }
+ 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 (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 = 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;
+ 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;
+ if (m_traceLevel >= Trace::TRACE_CONSOLE)
+ std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
#endif
}
m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
-
+
if (m_nFilterPoints % 2) { // Odd
- m_dFilterMin = -1. / (2 * m_dSignalInc);
- m_dFilterMax = 1. / (2 * m_dSignalInc);
- m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
+ m_dFilterMin = -1. / (2 * m_dSignalInc);
+ m_dFilterMax = 1. / (2 * m_dSignalInc);
+ m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
} else { // Even
- m_dFilterMin = -1. / (2 * m_dSignalInc);
- m_dFilterMax = 1. / (2 * m_dSignalInc);
- m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
- m_dFilterMax -= m_dFilterInc;
+ m_dFilterMin = -1. / (2 * m_dSignalInc);
+ m_dFilterMax = 1. / (2 * m_dSignalInc);
+ m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
+ m_dFilterMax -= m_dFilterInc;
}
-
+
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
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
- for (i = 0; i < m_nFilterPoints; i++)
- m_adFilter[i] *= 0.5;
+ 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;
- }
+ 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;
+ }
}
#ifdef HAVE_SGP
EZPlot* pEZPlot = NULL;
if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
- pEZPlot = new EZPlot (*pSGP);
- 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();
+ 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);
}
#endif
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();
- delete pEZPlot;
+ 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;
}
#endif
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
m_nFilterPoints = nSpatialPoints;
if (m_iZeropad > 0) {
- double logBase2 = log(nSpatialPoints) / log(2);
- int nextPowerOf2 = static_cast<int>(floor(logBase2));
- if (logBase2 != floor(logBase2))
- nextPowerOf2++;
- nextPowerOf2 += (m_iZeropad - 1);
- m_nFilterPoints = 1 << nextPowerOf2;
+ double logBase2 = log(nSpatialPoints) / log(2);
+ int nextPowerOf2 = static_cast<int>(floor(logBase2));
+ if (logBase2 != floor(logBase2))
+ nextPowerOf2++;
+ nextPowerOf2 += (m_iZeropad - 1);
+ m_nFilterPoints = 1 << nextPowerOf2;
}
m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
#ifdef DEBUG
if (m_traceLevel >= Trace::TRACE_CONSOLE)
- std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
+ std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
#endif
double* adSpatialFilter = new double [m_nFilterPoints];\r
SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
#ifdef HAVE_SGP
EZPlot* pEZPlot = NULL;
if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
- pEZPlot = new EZPlot (*pSGP);
- pEZPlot->ezset ("title Spatial Filter: Natural Order");
- pEZPlot->ezset ("ylength 0.50");
- pEZPlot->ezset ("yporigin 0.00");
- pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
- pEZPlot->plot();
- delete pEZPlot;
+ 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;
}
#endif
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
- for (i = 0; i < m_nFilterPoints; i++)
- adSpatialFilter[i] *= 0.5;
+ for (i = 0; i < m_nFilterPoints; 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);
- 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;
- adSpatialFilter[i] *= dScale;
- }
+ 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;
+ adSpatialFilter[i] *= dScale;
+ }
}\r
for (i = nSpatialPoints; i < m_nFilterPoints; i++)
- adSpatialFilter[i] = 0;
-
+ adSpatialFilter[i] = 0;
+
m_adFilter = new double [m_nFilterPoints];
std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];\r
finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
- delete adSpatialFilter;\r
- for (i = 0; i < m_nFilterPoints; i++)
- m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc;
- delete acInverseFilter;\r
+ delete adSpatialFilter;\r
+ for (i = 0; i < m_nFilterPoints; i++)
+ m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc;
+ delete acInverseFilter;\r
#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();
- delete pEZPlot;\r
+ 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;\r
}
#endif
}
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);
m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
ProcessSignal::~ProcessSignal (void)
{
- delete [] m_adFourierSinTable;
- delete [] m_adFourierCosTable;
- delete [] m_adFilter;
-
+ 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;
- }
- if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
- rfftw_destroy_plan(m_realPlanForward);
- rfftw_destroy_plan(m_realPlanBackward);
- delete [] m_adRealFftInput;
- delete [] m_adRealFftSignal;
- }
+ if (m_idFilterMethod == FILTER_METHOD_FFTW) {
+ fftw_destroy_plan(m_complexPlanForward);
+ fftw_destroy_plan(m_complexPlanBackward);
+ delete [] m_adComplexFftInput;
+ delete [] m_adComplexFftSignal;
+ }
+ if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
+ rfftw_destroy_plan(m_realPlanForward);
+ rfftw_destroy_plan(m_realPlanBackward);
+ delete [] m_adRealFftInput;
+ delete [] m_adRealFftSignal;
+ }
#endif
}
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) {
+ 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 (s_aszFilterMethodName [fmID]);
+
return (name);
}
ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
{
static const char *title = "";
-
+
if (fmID >= 0 && fmID < s_iFilterMethodCount)
- return (s_aszFilterMethodTitle [fmID]);
-
+ return (s_aszFilterMethodTitle [fmID]);
+
return (title);
}
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) {
+ 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 (s_aszFilterGenerationName [fgID]);
+
return (name);
}
ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
{
static const char *name = "";
-
+
if (fgID >= 0 && fgID < s_iFilterGenerationCount)
- return (s_aszFilterGenerationTitle [fgID]);
-
+ return (s_aszFilterGenerationTitle [fgID]);
+
return (name);
}
int i;\r
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);
}
}\r
if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
- for (i = 0; i < m_nSignalPoints; i++)
- output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
+ for (i = 0; i < m_nSignalPoints; i++)
+ output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
} else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
double* inputSignal = new double [m_nFilterPoints];
for (i = 0; i < m_nSignalPoints; i++)
inputSignal[i] = 0; // zeropad
std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);\r
- delete inputSignal;
+ delete inputSignal;
for (i = 0; i < m_nFilterPoints; i++)
fftSignal[i] *= m_adFilter[i];
double* inverseFourier = new double [m_nFilterPoints];
finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);\r
- delete fftSignal;
+ delete fftSignal;
for (i = 0; i < m_nSignalPoints; i++)
output[i] = inverseFourier[i];\r
- delete inverseFourier;
+ delete inverseFourier;
} else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
double* inputSignal = new double [m_nFilterPoints];
for (i = 0; i < m_nSignalPoints; i++)
inputSignal[i] = 0; // zeropad
std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
finiteFourierTransform (inputSignal, fftSignal, -1);\r
- delete inputSignal;
+ delete inputSignal;
for (i = 0; i < m_nFilterPoints; i++)
fftSignal[i] *= m_adFilter[i];
double* inverseFourier = new double [m_nFilterPoints];
finiteFourierTransform (fftSignal, inverseFourier, 1);\r
- delete fftSignal;
+ delete fftSignal;
for (i = 0; i < m_nSignalPoints; i++)
output[i] = inverseFourier[i];\r
- delete inverseFourier;
+ delete inverseFourier;
}
#if HAVE_FFTW
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);
- for (i = 0; i < m_nFilterPoints; i++)
- m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];\r
- delete [] fftOutput;
- for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
+ 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);
+ for (i = 0; i < m_nFilterPoints; i++)
+ m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];\r
+ delete [] fftOutput;
+ 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);
- for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
- output[i] = ifftOutput[i];\r
- delete [] ifftOutput;
+
+ fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
+ rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
+ for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
+ output[i] = ifftOutput[i];\r
+ delete [] ifftOutput;
} 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);
- 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;
- }\r
- 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;\r
- delete [] ifftOutput;
+ 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);
+ 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;
+ }\r
+ 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;\r
+ delete [] ifftOutput;
}
#endif\r
delete input;
/* NAME
- * 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
- *
- * 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]
- * 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].
- */
+* 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
+*
+* 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]
+* 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].
+*/
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)];
for (int i = 0; i < np; i++)
sum += *func++ * *f2--;
#endif
-
+
return (sum * dx);
}
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)];
+ for (int i = 0; i < np; i++)
+ sum += func[i] * m_adFilter[n - i + (np - 1)];
#else
-double* f2 = m_adFilter + n + (np - 1);
-for (int i = 0; i < np; i++)
- sum += *func++ * *f2--;
+ double* f2 = m_adFilter + n + (np - 1);
+ for (int i = 0; i < np; i++)
+ sum += *func++ * *f2--;
#endif
-
+
return (sum * dx);
}
void
ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
{
- std::complex<double>* complexOutput = new std::complex<double> [n];
-
- finiteFourierTransform (input, complexOutput, n, direction);
- for (int i = 0; i < n; i++)
- output[i] = complexOutput[i].real();\r
- delete [] complexOutput;
+ std::complex<double>* complexOutput = new std::complex<double> [n];
+
+ finiteFourierTransform (input, complexOutput, n, direction);
+ for (int i = 0; i < n; i++)
+ output[i] = complexOutput[i].real();\r
+ delete [] complexOutput;
}
void
direction = -1;
else
direction = 1;
-
+
double angleIncrement = direction * 2 * PI / n;
for (int i = 0; i < n; i++) {
double sumReal = 0;
direction = -1;
else
direction = 1;
-
+
double angleIncrement = direction * 2 * PI / n;
for (int i = 0; i < n; i++) {
std::complex<double> sum (0,0);
direction = -1;
else
direction = 1;
-
+
double angleIncrement = direction * 2 * PI / n;
for (int i = 0; i < n; i++) {
- double sumReal = 0;
+ double sumReal = 0;
for (int j = 0; j < n; j++) {
double angle = i * j * angleIncrement;
sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
direction = -1;
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] * m_adFourierCosTable[tableIndex];
- sumImag += input[j] * m_adFourierSinTable[tableIndex];
+ sumReal += input[j] * m_adFourierCosTable[tableIndex];
+ sumImag += input[j] * m_adFourierSinTable[tableIndex];
} else {
- sumReal += input[j] * m_adFourierCosTable[tableIndex];
- sumImag -= input[j] * m_adFourierSinTable[tableIndex];
+ sumReal += input[j] * m_adFourierCosTable[tableIndex];
+ sumImag -= input[j] * m_adFourierSinTable[tableIndex];
}
}
if (direction < 0) {
direction = -1;
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]
- - input[j].imag() * m_adFourierSinTable[tableIndex];
- sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
- + input[j].imag() * 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]
- - input[j].imag() * -m_adFourierSinTable[tableIndex];
- sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
- + input[j].imag() * 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];
}
}
if (direction < 0) {
direction = -1;
else
direction = 1;
-
+
for (int i = 0; i < m_nFilterPoints; i++) {
- double sumReal = 0;
+ 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]
- - input[j].imag() * m_adFourierSinTable[tableIndex];
+ sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
+ - input[j].imag() * m_adFourierSinTable[tableIndex];
} else {
- sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
- - input[j].imag() * -m_adFourierSinTable[tableIndex];
+ sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
+ - input[j].imag() * -m_adFourierSinTable[tableIndex];
}
}
if (direction < 0) {
// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
-void
-ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
-{
+void\r
+ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)\r
+{\r
double* pdTemp = new double [n];\r
- int i;
- if (n % 2) { // Odd
- int iHalfN = (n - 1) / 2;
-
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i];
- } else { // Even
- int iHalfN = n / 2;
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
-}
-
-
-void
-ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
-{
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete pdTemp;\r
+}\r
+\r
+void\r
+ProcessSignal::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)\r
+{\r
+ std::complex<double>* pdTemp = new std::complex<double> [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + 1 + iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[0] = pdVector[iHalfN];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1] = pdVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete [] pdTemp;\r
+}\r
+\r
+
+void\r
+ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)\r
+{\r
double* pdTemp = new double [n];\r
- int i;
- if (n % 2) { // Odd
- int iHalfN = (n - 1) / 2;
-
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
pdTemp[iHalfN] = pdVector[0];\r
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i] = pdVector[i + iHalfN + 1];
- } else { // Even
- int iHalfN = n / 2;
- pdTemp[iHalfN] = pdVector[0];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i] = pdVector[i + iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i+1];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
-}
-
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN + 1];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete pdTemp;\r
+}\r
+\r
+void\r
+ProcessSignal::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)\r
+{\r
+ std::complex<double>* pdTemp = new std::complex<double> [n];\r
+ int i;\r
+ if (n % 2) { // Odd\r
+ int iHalfN = (n - 1) / 2;\r
+ \r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i + 1 + iHalfN] = pdVector[i + 1];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN + 1];\r
+ } else { // Even\r
+ int iHalfN = n / 2;\r
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)\r
+ pdTemp[i] = pdVector[i + iHalfN];\r
+ for (i = 0; i < iHalfN - 1; i++)\r
+ pdTemp[i + iHalfN + 1] = pdVector[i+1];\r
+ }\r
+ \r
+ for (i = 0; i < n; i++)\r
+ pdVector[i] = pdTemp[i];\r
+ delete [] pdTemp;\r
+}\r
+\r
projects / ctsim.git / blobdiff