1 /*****************************************************************************
5 ** Purpose: Routines for signal-procesing filters
6 ** Progammer: Kevin Rosenberg
7 ** Date Started: Aug 1984
9 ** This is part of the CTSim program
10 ** Copyright (C) 1983-2000 Kevin Rosenberg
12 ** $Id: procsignal.cpp,v 1.8 2000/12/06 01:46:43 kevin Exp $
14 ** This program is free software; you can redistribute it and/or modify
15 ** it under the terms of the GNU General Public License (version 2) as
16 ** published by the Free Software Foundation.
18 ** This program is distributed in the hope that it will be useful,
19 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 ** GNU General Public License for more details.
23 ** You should have received a copy of the GNU General Public License
24 ** along with this program; if not, write to the Free Software
25 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 ******************************************************************************/
30 // FilterMethod ID/Names
31 const int ProcessSignal::FILTER_METHOD_INVALID = -1;
32 const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
33 const int ProcessSignal::FILTER_METHOD_FOURIER = 1;
34 const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2;
35 const int ProcessSignal::FILTER_METHOD_FFT = 3;
37 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
38 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
40 const char* ProcessSignal::s_aszFilterMethodName[] = {
50 const char* ProcessSignal::s_aszFilterMethodTitle[] = {
53 {"Fouier Trigometric Table"},
57 {"Real/Half-Complex FFTW"},
60 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
62 // FilterGeneration ID/Names
63 const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
64 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
65 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
66 const char* ProcessSignal::s_aszFilterGenerationName[] = {
70 const char* ProcessSignal::s_aszFilterGenerationTitle[] = {
74 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
77 // CLASS IDENTIFICATION
80 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)
81 : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
83 m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
84 if (m_idFilterMethod == FILTER_METHOD_INVALID) {
86 m_failMessage = "Invalid filter method name ";
87 m_failMessage += szFilterMethodName;
90 m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
91 if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
93 m_failMessage = "Invalid frequency filter name ";
94 m_failMessage += szFilterGenerationName;
97 m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
98 if (m_idFilter == SignalFilter::FILTER_INVALID) {
100 m_failMessage = "Invalid Filter name ";
101 m_failMessage += szFilterName;
104 m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
105 if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
107 m_failMessage = "Invalid domain name ";
108 m_failMessage += szDomainName;
112 init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
117 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)
119 m_idFilter = idFilter;
120 m_idDomain = idDomain;
121 m_idFilterMethod = idFilterMethod;
122 m_idFilterGeneration = idFilterGeneration;
123 m_idGeometry = iGeometry;
124 m_dFocalLength = dFocalLength;
126 if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
130 m_traceLevel = iTraceLevel;
131 m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
132 m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
133 m_dBandwidth = dBandwidth;
134 m_nSignalPoints = nSignalPoints;
135 m_dSignalInc = dSignalIncrement;
136 m_dFilterParam = dFilterParam;
137 m_iZeropad = iZeropad;
138 m_iPreinterpolationFactor = iPreinterpolationFactor;
140 // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
141 // through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear
142 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
147 if (m_idFilterMethod == FILTER_METHOD_FFT) {
149 m_idFilterMethod = FILTER_METHOD_RFFTW;
152 m_failMessage = "FFT not yet implemented";
157 bool m_bFrequencyFiltering = true;
158 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
159 m_bFrequencyFiltering = false;
161 // Spatial-based filtering
162 if (! m_bFrequencyFiltering) {
164 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
165 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
166 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
167 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
168 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
169 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
170 m_adFilter = new double[ m_nFilterPoints ];
171 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
172 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
173 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
174 m_dFilterMin = -1. / (2 * m_dSignalInc);
175 m_dFilterMax = 1. / (2 * m_dSignalInc);
176 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
177 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
178 m_adFilter = new double[ m_nFilterPoints ];
179 double* adFrequencyFilter = new double [m_nFilterPoints];
180 filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
182 EZPlot* pEZPlot = NULL;
183 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
184 pEZPlot = new EZPlot (*pSGP);
185 pEZPlot->ezset ("title Filter Response: Natural Order");
186 pEZPlot->ezset ("ylength 0.25");
187 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
191 shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
193 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
194 pEZPlot->ezset ("title Filter Response: Fourier Order");
195 pEZPlot->ezset ("ylength 0.25");
196 pEZPlot->ezset ("yporigin 0.25");
197 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
201 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
202 delete adFrequencyFilter;
\r
204 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
205 pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
206 pEZPlot->ezset ("ylength 0.25");
207 pEZPlot->ezset ("yporigin 0.50");
208 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
212 shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
214 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
215 pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
216 pEZPlot->ezset ("ylength 0.25");
217 pEZPlot->ezset ("yporigin 0.75");
218 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
223 for (int i = 0; i < m_nFilterPoints; i++) {
224 m_adFilter[i] /= m_dSignalInc;
227 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
228 for (int i = 0; i < m_nFilterPoints; i++)
229 m_adFilter[i] *= 0.5;
230 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
231 for (int i = 0; i < m_nFilterPoints; i++) {
232 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
233 double sinScale = sin (iDetFromZero * m_dSignalInc);
234 if (fabs(sinScale) < 1E-7)
237 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
238 double dScale = 0.5 * sinScale * sinScale;
239 m_adFilter[i] *= dScale;
242 } // if (spatial filtering)
244 else if (m_bFrequencyFiltering) { // Frequency-based filtering
246 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
247 // calculate number of filter points with zeropadding
248 m_nFilterPoints = m_nSignalPoints;
249 if (m_iZeropad > 0) {
250 double logBase2 = log(m_nFilterPoints) / log(2);
251 int nextPowerOf2 = static_cast<int>(floor(logBase2));
252 if (logBase2 != floor(logBase2))
254 nextPowerOf2 += (m_iZeropad - 1);
255 m_nFilterPoints = 1 << nextPowerOf2;
257 if (m_traceLevel >= Trace::TRACE_CONSOLE)
258 cout << "nFilterPoints = " << m_nFilterPoints << endl;
261 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
263 if (m_nFilterPoints % 2) { // Odd
264 m_dFilterMin = -1. / (2 * m_dSignalInc);
265 m_dFilterMax = 1. / (2 * m_dSignalInc);
266 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
268 m_dFilterMin = -1. / (2 * m_dSignalInc);
269 m_dFilterMax = 1. / (2 * m_dSignalInc);
270 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
271 m_dFilterMax -= m_dFilterInc;
274 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
275 m_adFilter = new double [m_nFilterPoints];
276 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
278 // This doesn't work!
279 // Need to add filtering for divergent geometries & Frequency/Direct filtering
280 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
281 for (int i = 0; i < m_nFilterPoints; i++)
282 m_adFilter[i] *= 0.5;
283 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
284 for (int i = 0; i < m_nFilterPoints; i++) {
285 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
286 double sinScale = sin (iDetFromZero * m_dSignalInc);
287 if (fabs(sinScale) < 1E-7)
290 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
291 double dScale = 0.5 * sinScale * sinScale;
292 m_adFilter[i] *= dScale;
296 EZPlot* pEZPlot = NULL;
297 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
298 pEZPlot = new EZPlot (*pSGP);
299 pEZPlot->ezset ("title Filter Filter: Natural Order");
300 pEZPlot->ezset ("ylength 0.50");
301 pEZPlot->ezset ("yporigin 0.00");
302 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
306 shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
308 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
309 pEZPlot->ezset ("title Filter Filter: Fourier Order");
310 pEZPlot->ezset ("ylength 0.50");
311 pEZPlot->ezset ("yporigin 0.50");
312 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
317 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
318 // calculate number of filter points with zeropadding
319 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
320 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
321 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
322 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
323 m_nFilterPoints = nSpatialPoints;
324 if (m_iZeropad > 0) {
325 double logBase2 = log(nSpatialPoints) / log(2);
326 int nextPowerOf2 = static_cast<int>(floor(logBase2));
327 if (logBase2 != floor(logBase2))
329 nextPowerOf2 += (m_iZeropad - 1);
330 m_nFilterPoints = 1 << nextPowerOf2;
332 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
334 if (m_traceLevel >= Trace::TRACE_CONSOLE)
335 cout << "nFilterPoints = " << m_nFilterPoints << endl;
337 double* adSpatialFilter = new double [m_nFilterPoints];
\r
338 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
339 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
341 EZPlot* pEZPlot = NULL;
342 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
343 pEZPlot = new EZPlot (*pSGP);
344 pEZPlot->ezset ("title Spatial Filter: Natural Order");
345 pEZPlot->ezset ("ylength 0.50");
346 pEZPlot->ezset ("yporigin 0.00");
347 pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
352 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
353 for (int i = 0; i < m_nFilterPoints; i++)
354 adSpatialFilter[i] *= 0.5;
355 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
356 for (int i = 0; i < m_nFilterPoints; i++) {
357 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
358 double sinScale = sin (iDetFromZero * m_dSignalInc);
359 if (fabs(sinScale) < 1E-7)
362 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
363 double dScale = 0.5 * sinScale * sinScale;
364 adSpatialFilter[i] *= dScale;
368 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
369 adSpatialFilter[i] = 0;
371 m_adFilter = new double [m_nFilterPoints];
372 complex<double>* acInverseFilter = new complex<double> [m_nFilterPoints];
\r
373 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
374 delete adSpatialFilter;
\r
375 for (i = 0; i < m_nFilterPoints; i++)
376 m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
377 delete acInverseFilter;
\r
379 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
380 pEZPlot->ezset ("title Spatial Filter: Inverse");
381 pEZPlot->ezset ("ylength 0.50");
382 pEZPlot->ezset ("yporigin 0.50");
383 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
391 // precalculate sin and cosine tables for fourier transform
392 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
393 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
394 double angleIncrement = (2. * PI) / m_nFilterPoints;
395 m_adFourierCosTable = new double[ nFourier ];
396 m_adFourierSinTable = new double[ nFourier ];
398 for (int i = 0; i < nFourier; i++) {
399 m_adFourierCosTable[i] = cos (angle);
400 m_adFourierSinTable[i] = sin (angle);
401 angle += angleIncrement;
406 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
407 for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
408 m_adFilter[i] /= m_nFilterPoints;
411 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
412 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
413 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
414 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
415 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
416 for (int i = 0; i < m_nFilterPoints; i++)
417 m_adRealFftInput[i] = 0;
418 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
419 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
420 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
421 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
422 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
423 for (int i = 0; i < m_nFilterPoints; i++)
424 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
425 for (int i = 0; i < m_nOutputPoints; i++)
426 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
432 ProcessSignal::~ProcessSignal (void)
434 delete [] m_adFourierSinTable;
435 delete [] m_adFourierCosTable;
436 delete [] m_adFilter;
439 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
440 fftw_destroy_plan(m_complexPlanForward);
441 fftw_destroy_plan(m_complexPlanBackward);
442 delete [] m_adComplexFftInput;
443 delete [] m_adComplexFftSignal;
445 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
446 rfftw_destroy_plan(m_realPlanForward);
447 rfftw_destroy_plan(m_realPlanBackward);
448 delete [] m_adRealFftInput;
449 delete [] m_adRealFftSignal;
455 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
457 int fmID = FILTER_METHOD_INVALID;
459 for (int i = 0; i < s_iFilterMethodCount; i++)
460 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
469 ProcessSignal::convertFilterMethodIDToName (const int fmID)
471 static const char *name = "";
473 if (fmID >= 0 && fmID < s_iFilterMethodCount)
474 return (s_aszFilterMethodName [fmID]);
480 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
482 static const char *title = "";
484 if (fmID >= 0 && fmID < s_iFilterMethodCount)
485 return (s_aszFilterMethodTitle [fmID]);
492 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
494 int fgID = FILTER_GENERATION_INVALID;
496 for (int i = 0; i < s_iFilterGenerationCount; i++)
497 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
506 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
508 static const char *name = "";
510 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
511 return (s_aszFilterGenerationName [fgID]);
517 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
519 static const char *name = "";
521 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
522 return (s_aszFilterGenerationTitle [fgID]);
528 ProcessSignal::filterSignal (const float constInput[], double output[]) const
530 double* input = new double [m_nSignalPoints];
532 for (i = 0; i < m_nSignalPoints; i++)
533 input[i] = constInput[i];
535 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
536 for (int i = 0; i < m_nSignalPoints; i++) {
537 int iDetFromCenter = i - (m_nSignalPoints / 2);
538 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
540 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
541 for (int i = 0; i < m_nSignalPoints; i++) {
542 int iDetFromCenter = i - (m_nSignalPoints / 2);
543 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
546 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
547 for (i = 0; i < m_nSignalPoints; i++)
548 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
549 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
550 double* inputSignal = new double [m_nFilterPoints];
551 for (i = 0; i < m_nSignalPoints; i++)
552 inputSignal[i] = input[i];
553 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
554 inputSignal[i] = 0; // zeropad
555 complex<double>* fftSignal = new complex<double> [m_nFilterPoints];
556 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
\r
558 for (i = 0; i < m_nFilterPoints; i++)
559 fftSignal[i] *= m_adFilter[i];
560 double* inverseFourier = new double [m_nFilterPoints];
561 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
\r
563 for (i = 0; i < m_nSignalPoints; i++)
564 output[i] = inverseFourier[i];
\r
565 delete inverseFourier;
566 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
567 double* inputSignal = new double [m_nFilterPoints];
568 for (i = 0; i < m_nSignalPoints; i++)
569 inputSignal[i] = input[i];
570 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
571 inputSignal[i] = 0; // zeropad
572 complex<double>* fftSignal = new complex<double> [m_nFilterPoints];
573 finiteFourierTransform (inputSignal, fftSignal, -1);
\r
575 for (i = 0; i < m_nFilterPoints; i++)
576 fftSignal[i] *= m_adFilter[i];
577 double* inverseFourier = new double [m_nFilterPoints];
578 finiteFourierTransform (fftSignal, inverseFourier, 1);
\r
580 for (i = 0; i < m_nSignalPoints; i++)
581 output[i] = inverseFourier[i];
\r
582 delete inverseFourier;
585 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
586 for (int i = 0; i < m_nSignalPoints; i++)
587 m_adRealFftInput[i] = input[i];
589 fftw_real fftOutput [ m_nFilterPoints ];
590 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
591 for (int i = 0; i < m_nFilterPoints; i++)
592 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
593 for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
594 m_adRealFftSignal[i] = 0;
596 fftw_real ifftOutput [ m_nOutputPoints ];
597 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
598 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
599 output[i] = ifftOutput[i];
600 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
601 for (int i = 0; i < m_nSignalPoints; i++)
602 m_adComplexFftInput[i].re = input[i];
604 fftw_complex fftOutput [ m_nFilterPoints ];
605 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
606 for (int i = 0; i < m_nFilterPoints; i++) {
607 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
608 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
610 fftw_complex ifftOutput [ m_nOutputPoints ];
611 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
612 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
613 output[i] = ifftOutput[i].re;
621 * convolve Discrete convolution of two functions
624 * r = convolve (f1, f2, dx, n, np, func_type)
625 * double r Convolved result
626 * double f1[], f2[] Functions to be convolved
627 * double dx Difference between successive x values
628 * int n Array index to center convolution about
629 * int np Number of points in f1 array
630 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
633 * f1 is the projection data, its indices range from 0 to np - 1.
634 * The index for f2, the filter, ranges from -(np-1) to (np-1).
635 * There are 3 ways to handle the negative vertices of f2:
636 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
637 * All filters used in reconstruction are even.
638 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
639 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
640 * for negative indices. Since f2 must range from -(np-1) to (np-1),
641 * if we add (np - 1) to f2's array index, then f2's index will
642 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
643 * stored at f2[np-1].
647 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
651 #if UNOPTIMIZED_CONVOLUTION
652 for (int i = 0; i < np; i++)
653 sum += func[i] * m_adFilter[n - i + (np - 1)];
655 double* f2 = m_adFilter + n + (np - 1);
656 for (int i = 0; i < np; i++)
657 sum += *func++ * *f2--;
665 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
669 #if UNOPTIMIZED_CONVOLUTION
670 for (int i = 0; i < np; i++)
671 sum += func[i] * m_adFilter[n - i + (np - 1)];
673 double* f2 = m_adFilter + n + (np - 1);
674 for (int i = 0; i < np; i++)
675 sum += *func++ * *f2--;
683 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
685 complex<double>* complexOutput = new complex<double> [n];
687 finiteFourierTransform (input, complexOutput, n, direction);
688 for (int i = 0; i < n; i++)
689 output[i] = complexOutput[i].real();
\r
690 delete [] complexOutput;
694 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
701 double angleIncrement = direction * 2 * PI / n;
702 for (int i = 0; i < n; i++) {
705 for (int j = 0; j < n; j++) {
706 double angle = i * j * angleIncrement;
707 sumReal += input[j] * cos(angle);
708 sumImag += input[j] * sin(angle);
714 output[i] = complex<double> (sumReal, sumImag);
720 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
727 double angleIncrement = direction * 2 * PI / n;
728 for (int i = 0; i < n; i++) {
729 complex<double> sum (0,0);
730 for (int j = 0; j < n; j++) {
731 double angle = i * j * angleIncrement;
732 complex<double> exponentTerm (cos(angle), sin(angle));
733 sum += input[j] * exponentTerm;
743 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
750 double angleIncrement = direction * 2 * PI / n;
751 for (int i = 0; i < n; i++) {
753 for (int j = 0; j < n; j++) {
754 double angle = i * j * angleIncrement;
755 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
764 // Table-based routines
767 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
774 for (int i = 0; i < m_nFilterPoints; i++) {
775 double sumReal = 0, sumImag = 0;
776 for (int j = 0; j < m_nFilterPoints; j++) {
777 int tableIndex = i * j;
779 sumReal += input[j] * m_adFourierCosTable[tableIndex];
780 sumImag += input[j] * m_adFourierSinTable[tableIndex];
782 sumReal += input[j] * m_adFourierCosTable[tableIndex];
783 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
787 sumReal /= m_nFilterPoints;
788 sumImag /= m_nFilterPoints;
790 output[i] = complex<double> (sumReal, sumImag);
794 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
796 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
803 for (int i = 0; i < m_nFilterPoints; i++) {
804 double sumReal = 0, sumImag = 0;
805 for (int j = 0; j < m_nFilterPoints; j++) {
806 int tableIndex = i * j;
808 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
809 - input[j].imag() * m_adFourierSinTable[tableIndex];
810 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
811 + input[j].imag() * m_adFourierCosTable[tableIndex];
813 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
814 - input[j].imag() * -m_adFourierSinTable[tableIndex];
815 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
816 + input[j].imag() * m_adFourierCosTable[tableIndex];
820 sumReal /= m_nFilterPoints;
821 sumImag /= m_nFilterPoints;
823 output[i] = complex<double> (sumReal, sumImag);
828 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
835 for (int i = 0; i < m_nFilterPoints; i++) {
837 for (int j = 0; j < m_nFilterPoints; j++) {
838 int tableIndex = i * j;
840 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
841 - input[j].imag() * m_adFourierSinTable[tableIndex];
843 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
844 - input[j].imag() * -m_adFourierSinTable[tableIndex];
848 sumReal /= m_nFilterPoints;
854 // Odd Number of Points
855 // Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
856 // Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
857 // Even Number of Points
858 // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
859 // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
862 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
864 double* pdTemp = new double [n];
\r
867 int iHalfN = (n - 1) / 2;
869 pdTemp[0] = pdVector[iHalfN];
870 for (i = 0; i < iHalfN; i++)
871 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
872 for (i = 0; i < iHalfN; i++)
873 pdTemp[i + iHalfN + 1] = pdVector[i];
876 pdTemp[0] = pdVector[iHalfN];
877 for (i = 0; i < iHalfN; i++)
878 pdTemp[i + 1] = pdVector[i + iHalfN];
879 for (i = 0; i < iHalfN - 1; i++)
880 pdTemp[i + iHalfN + 1] = pdVector[i];
883 for (i = 0; i < n; i++)
884 pdVector[i] = pdTemp[i];
890 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
892 double* pdTemp = new double [n];
\r
895 int iHalfN = (n - 1) / 2;
897 pdTemp[iHalfN] = pdVector[0];
\r
898 for (i = 0; i < iHalfN; i++)
899 pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
900 for (i = 0; i < iHalfN; i++)
901 pdTemp[i] = pdVector[i + iHalfN + 1];
904 pdTemp[iHalfN] = pdVector[0];
905 for (i = 0; i < iHalfN; i++)
906 pdTemp[i] = pdVector[i + iHalfN];
907 for (i = 0; i < iHalfN - 1; i++)
908 pdTemp[i + iHalfN + 1] = pdVector[i+1];
911 for (i = 0; i < n; i++)
912 pdVector[i] = pdTemp[i];