/*****************************************************************************
** File IDENTIFICATION
-**
+**
** Name: procsignal.cpp
** Purpose: Routines for processing signals and projections
-** Progammer: Kevin Rosenberg
+** Progammer: Kevin Rosenberg
** Date Started: Aug 1984
**
** This is part of the CTSim program
// 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,
+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_failMessage += szDomainName;
return;
}
-
- 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,
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,
+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_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 adjust for imaginary detector through origin of phantom
+
+ // 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) {
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);
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 = 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);
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;
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
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;
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
-
+
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);
#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_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);
- for (i = 0; i < m_nFilterPoints; i++)
+ for (i = 0; i < m_nFilterPoints; i++)
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);
- for (i = 0; i < m_nFilterPoints; i++)
+ for (i = 0; i < m_nFilterPoints; i++)
m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0;
- for (i = 0; i < m_nOutputPoints; i++)
+ 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);
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;
}
-
+
return (fmID);
}
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++)
if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
fgID = i;
break;
}
-
+
return (fgID);
}
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);
}
int i;
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);
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) {
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_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++)
- m_adRealFftSignal[i] = 0;
-
+ m_adRealFftSignal[i] = 0;
+
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);
- for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
+ for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
output[i] = m_adComplexFftBackwardOutput[i][0];
}
#endif
/* 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];
}
}
projects / ctsim.git / blobdiff