/*****************************************************************************
** File IDENTIFICATION
-**
-** Name: filter.cpp
-** Purpose: Routines for signal-procesing filters
-** Progammer: Kevin Rosenberg
-** Date Started: Aug 1984
**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
+** Name: procsignal.cpp
+** Purpose: Routines for processing signals and projections
+** Progammer: Kevin Rosenberg
+** Date Started: Aug 1984
**
-** $Id: procsignal.cpp,v 1.20 2001/01/12 21:53:27 kevin Exp $
+** This is part of the CTSim program
+** Copyright (c) 1983-2009 Kevin Rosenberg
**
** 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 "../src/dlgezplot.h"
+#include "nographics.h"
#endif
// FilterMethod ID/Names
const int ProcessSignal::FILTER_METHOD_FFTW = 4;
const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
#endif
-const char* ProcessSignal::s_aszFilterMethodName[] = {
- {"convolution"},
- {"fourier"},
- {"fouier_table"},
- {"fft"},
+const char* const ProcessSignal::s_aszFilterMethodName[] = {
+ "convolution",
+ "fourier",
+ "fouier-table",
+ "fft",
#if HAVE_FFTW
- {"fftw"},
- {"rfftw"},
+ "fftw",
+ "rfftw",
#endif
};
-const char* ProcessSignal::s_aszFilterMethodTitle[] = {
- {"Convolution"},
- {"Fourier"},
- {"Fouier Trigometric Table"},
- {"FFT"},
+const char* const ProcessSignal::s_aszFilterMethodTitle[] = {
+ "Convolution",
+ "Fourier",
+ "Fouier Trigometric Table",
+ "FFT",
#if HAVE_FFTW
- {"FFTW"},
- {"Real/Half-Complex FFTW"},
+ "FFTW",
+ "Real/Half-Complex FFTW",
#endif
};
const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
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[] = {
- {"direct"},
- {"inverse_fourier"},
+const char* const ProcessSignal::s_aszFilterGenerationName[] = {
+ "direct",
+ "inverse-fourier",
};
-const char* ProcessSignal::s_aszFilterGenerationTitle[] = {
- {"Direct"},
- {"Inverse Fourier"},
+const char* const ProcessSignal::s_aszFilterGenerationTitle[] = {
+ "Direct",
+ "Inverse Fourier",
};
const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
// 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,
- int iGeometry, double dFocalLength, SGP* pSGP)
+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)
{
m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
m_failMessage += szDomainName;
return;
}
-
- init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
- m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
+
+ 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,
- double dFocalLength, SGP* pSGP)
+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;
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;
return;
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 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) {
#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);
double* adFrequencyFilter = new double [m_nFilterPoints];
filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
#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
+#endif
Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
#ifdef HAVE_SGP
if (g_bRunningWXWindows && m_traceLevel > 0) {
#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 = 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;
-#if defined(DEBUG) || defined(_DEBUG)
- if (m_traceLevel >= Trace::TRACE_CONSOLE)
- sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
-#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);
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;
}
#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
+ // 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;
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))
dlgEZPlot.ShowModal();
}
#endif
-
-#define PRE_JAN_2001 1
-#ifdef PRE_JAN_2001
+
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
for (i = 0; i < nSpatialPoints; i++)
adSpatialFilter[i] *= 0.5;
#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];
finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
dlgEZPlot.ShowModal();
}
#endif
-
-#else
- for (i = nSpatialPoints; i < m_nFilterPoints; i++)
- adSpatialFilter[i] = 0;
-
- std::complex<double>* acInverseFilter = new std::complex<double> [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]);
- 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);
- 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) || 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
}
}
-
+
// 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;
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);
- m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
- m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
- for (i = 0; i < m_nFilterPoints; i++)
+ m_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+ m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+ m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE);
+ m_adRealFftSignal = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
+ m_adRealFftBackwardOutput = static_cast<double*>(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++)
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_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
- m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
- for (i = 0; i < m_nFilterPoints; i++)
- m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
- for (i = 0; i < m_nOutputPoints; i++)
- m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
+ m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+ m_adComplexFftOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+ m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD, FFTW_ESTIMATE);
+ m_adComplexFftSignal = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
+ m_adComplexFftBackwardOutput = static_cast<fftw_complex*>(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++)
+ m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0;
+ for (i = 0; i < m_nOutputPoints; i++)
+ m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0;
}
#endif
-
+
}
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);
fftw_destroy_plan(m_complexPlanBackward);
- delete [] m_adComplexFftInput;
- delete [] m_adComplexFftSignal;
+ fftw_free (m_adComplexFftInput);
+ fftw_free (m_adComplexFftOutput);
+ fftw_free (m_adComplexFftSignal);
+ fftw_free (m_adComplexFftBackwardOutput);
}
if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
- rfftw_destroy_plan(m_realPlanForward);
- rfftw_destroy_plan(m_realPlanBackward);
- delete [] m_adRealFftInput;
- delete [] m_adRealFftSignal;
+ fftw_destroy_plan(m_realPlanForward);
+ fftw_destroy_plan(m_realPlanBackward);
+ fftw_free (m_adRealFftInput);
+ fftw_free (m_adRealFftOutput);
+ fftw_free (m_adRealFftSignal);
+ fftw_free (m_adRealFftBackwardOutput);
}
#endif
}
ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
{
int fmID = FILTER_METHOD_INVALID;
-
- for (int i = 0; i < s_iFilterMethodCount; i++)
+
+ for (int i = 0; i < s_iFilterMethodCount; i++) {
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 (name);
}
ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
{
static const char *title = "";
-
+
if (fmID >= 0 && fmID < s_iFilterMethodCount)
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++)
+
+ for (int i = 0; i < s_iFilterGenerationCount; i++) {
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 (name);
}
ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
{
static const char *name = "";
-
+
if (fgID >= 0 && fgID < s_iFilterGenerationCount)
return (s_aszFilterGenerationTitle [fgID]);
-
+
return (name);
}
{
double* input = new double [m_nSignalPoints];
int i;
+
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
for (i = 0; i < m_nSignalPoints; i++)
input[i] = constInput[i];
-
+
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
for (int i = 0; i < m_nSignalPoints; i++) {
int iDetFromCenter = i - (m_nSignalPoints / 2);
input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
}
} else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
for (int i = 0; i < m_nSignalPoints; i++) {
int iDetFromCenter = i - (m_nSignalPoints / 2);
input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
}
}
if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
for (i = 0; i < m_nSignalPoints; i++)
output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
} else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
delete inputSignal;
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
for (i = 0; i < m_nFilterPoints; i++)
fftSignal[i] *= m_adFilter[i];
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) {
std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
finiteFourierTransform (inputSignal, fftSignal, FORWARD);
delete inputSignal;
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
for (i = 0; i < m_nFilterPoints; i++)
fftSignal[i] *= m_adFilter[i];
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;
}
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);
+
+ fftw_execute (m_realPlanForward);
for (i = 0; i < m_nFilterPoints; i++)
- m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
- delete [] fftOutput;
+ m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
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);
+ m_adRealFftSignal[i] = 0;
+
+ fftw_execute (m_realPlanBackward);
for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
- output[i] = ifftOutput[i];
- delete [] ifftOutput;
+ output[i] = m_adRealFftBackwardOutput[i];
} 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);
+ m_adComplexFftInput[i][0] = input[i];
+
+ fftw_execute (m_complexPlanForward);
+#if HAVE_OPENMP
+ #pragma omp parallel for
+#endif
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;
+ m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
+ m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
}
- 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;
- delete [] ifftOutput;
+ fftw_execute (m_complexPlanBackward);
+ for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
+ output[i] = m_adComplexFftBackwardOutput[i][0];
}
#endif
delete input;
/* 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)];
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)];
for (int i = 0; i < np; i++)
sum += *func++ * *f2--;
#endif
-
+
return (sum * dx);
}
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();
{
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;
{
if (direction < 0)
direction = -1;
- else
+ else
direction = 1;
-
+
double angleIncrement = direction * 2 * PI / n;
for (int i = 0; i < n; i++) {
std::complex<double> sum (0,0);
{
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;
{
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++) {
{
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];
{
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];
}
}
}
}
+
+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