**
** 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
**
** 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 (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,
+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)
{
int iGeometry, double dFocalLength, double dSourceDetectorLength, SGP* pSGP)
: m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
{
-
- init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
- m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength,
+
+ init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
+ m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength,
-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,
+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,
m_idGeometry = iGeometry;
m_dFocalLength = dFocalLength;
m_dSourceDetectorLength = dSourceDetectorLength;
m_idGeometry = iGeometry;
m_dFocalLength = dFocalLength;
m_dSourceDetectorLength = dSourceDetectorLength;
m_dBandwidth = dBandwidth;
m_nSignalPoints = nSignalPoints;
m_dSignalInc = dSignalIncrement;
m_dBandwidth = dBandwidth;
m_nSignalPoints = nSignalPoints;
m_dSignalInc = dSignalIncrement;
// see Kak-Slaney Fig 3.22, for Collinear diagram
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
double dEquilinearScale = m_dSourceDetectorLength / m_dFocalLength;
m_dSignalInc /= dEquilinearScale;
m_dBandwidth *= dEquilinearScale;
}
// see Kak-Slaney Fig 3.22, for Collinear diagram
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
double dEquilinearScale = m_dSourceDetectorLength / m_dFocalLength;
m_dSignalInc /= dEquilinearScale;
m_dBandwidth *= dEquilinearScale;
}
bool m_bFrequencyFiltering = true;
if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
m_bFrequencyFiltering = false;
bool m_bFrequencyFiltering = true;
if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
m_bFrequencyFiltering = false;
if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
#ifdef HAVE_SGP
if (g_bRunningWXWindows && m_traceLevel > 0) {
Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
#ifdef HAVE_SGP
if (g_bRunningWXWindows && m_traceLevel > 0) {
if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
// calculate number of filter points with zeropadding
m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad);
m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
// calculate number of filter points with zeropadding
m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad);
m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
if (isOdd (m_nFilterPoints)) { // Odd
m_dFilterMin = -1. / (2 * m_dSignalInc);
m_dFilterMax = 1. / (2 * m_dSignalInc);
if (isOdd (m_nFilterPoints)) { // Odd
m_dFilterMin = -1. / (2 * m_dSignalInc);
m_dFilterMax = 1. / (2 * m_dSignalInc);
-
- 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);
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;
#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
if (g_bRunningWXWindows && m_traceLevel > 0) {
EZPlotDialog dlgEZPlot;
// 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
// 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
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
// calculate number of filter points with zeropadding
int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
// calculate number of filter points with zeropadding
int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
#endif
double* adSpatialFilter = new double [m_nFilterPoints];
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))
m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
for (i = 0; i < nSpatialPoints; i++)
adSpatialFilter[i] *= 0.5;
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
for (i = 0; i < nSpatialPoints; i++)
adSpatialFilter[i] *= 0.5;
m_adFilter = new double [m_nFilterPoints];
std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
m_adFilter = new double [m_nFilterPoints];
std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
// 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;
// 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;
#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 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_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
m_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
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);
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);
m_adRealFftInput[i] = 0;
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
m_adRealFftInput[i] = 0;
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
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);
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);
ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
{
int fmID = FILTER_METHOD_INVALID;
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) {
fmID = i;
break;
}
for (int i = 0; i < s_iFilterMethodCount; i++)
if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
fmID = i;
break;
}
ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
{
int fgID = FILTER_GENERATION_INVALID;
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) {
fgID = i;
break;
}
for (int i = 0; i < s_iFilterGenerationCount; i++)
if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
fgID = i;
break;
}
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
for (int i = 0; i < m_nSignalPoints; i++) {
int iDetFromCenter = i - (m_nSignalPoints / 2);
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
for (int i = 0; i < m_nSignalPoints; i++) {
int iDetFromCenter = i - (m_nSignalPoints / 2);
double* inverseFourier = new double [m_nFilterPoints];
finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
delete fftSignal;
double* inverseFourier = new double [m_nFilterPoints];
finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
delete fftSignal;
double* inverseFourier = new double [m_nFilterPoints];
finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
delete fftSignal;
double* inverseFourier = new double [m_nFilterPoints];
finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
delete fftSignal;
else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
for (i = 0; i < m_nSignalPoints; i++)
m_adRealFftInput[i] = input[i];
else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
for (i = 0; i < m_nSignalPoints; i++)
m_adRealFftInput[i] = input[i];
fftw_execute (m_realPlanForward);
for (i = 0; i < m_nFilterPoints; i++)
m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
fftw_execute (m_realPlanForward);
for (i = 0; i < m_nFilterPoints; i++)
m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
fftw_execute (m_realPlanBackward);
for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
output[i] = m_adRealFftBackwardOutput[i];
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
for (i = 0; i < m_nSignalPoints; i++)
m_adComplexFftInput[i][0] = input[i];
fftw_execute (m_realPlanBackward);
for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
output[i] = m_adRealFftBackwardOutput[i];
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
for (i = 0; i < m_nSignalPoints; i++)
m_adComplexFftInput[i][0] = input[i];
fftw_execute (m_complexPlanForward);
for (i = 0; i < m_nFilterPoints; i++) {
m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
}
fftw_execute (m_complexPlanBackward);
fftw_execute (m_complexPlanForward);
for (i = 0; i < m_nFilterPoints; i++) {
m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
}
fftw_execute (m_complexPlanBackward);
*
* SYNOPSIS
* r = convolve (f1, f2, dx, n, np, func_type)
*
* 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:
*
* 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]
-* 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].
ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
{
double sum = 0.0;
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)];
#if UNOPTIMIZED_CONVOLUTION
for (int i = 0; i < np; i++)
sum += func[i] * m_adFilter[n - i + (np - 1)];
ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
{
double sum = 0.0;
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)];
#if UNOPTIMIZED_CONVOLUTION
for (int i = 0; i < np; i++)
sum += func[i] * m_adFilter[n - i + (np - 1)];
ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
{
std::complex<double>* complexOutput = new std::complex<double> [n];
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();
finiteFourierTransform (input, complexOutput, n, direction);
for (int i = 0; i < n; i++)
output[i] = complexOutput[i].real();
double angleIncrement = direction * 2 * PI / n;
for (int i = 0; i < n; i++) {
std::complex<double> sum (0,0);
double angleIncrement = direction * 2 * PI / n;
for (int i = 0; i < n; i++) {
std::complex<double> sum (0,0);
for (int i = 0; i < m_nFilterPoints; i++) {
double sumReal = 0, sumImag = 0;
for (int j = 0; j < m_nFilterPoints; j++) {
for (int i = 0; i < m_nFilterPoints; i++) {
double sumReal = 0, sumImag = 0;
for (int j = 0; j < m_nFilterPoints; j++) {
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) {
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) {
- input[j].imag() * m_adFourierSinTable[tableIndex];
sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
+ input[j].imag() * m_adFourierCosTable[tableIndex];
} else {
- input[j].imag() * m_adFourierSinTable[tableIndex];
sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
+ input[j].imag() * m_adFourierCosTable[tableIndex];
} else {
- input[j].imag() * -m_adFourierSinTable[tableIndex];
sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
+ input[j].imag() * m_adFourierCosTable[tableIndex];
- input[j].imag() * -m_adFourierSinTable[tableIndex];
sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
+ input[j].imag() * m_adFourierCosTable[tableIndex];