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.7 2000/09/07 14:29:05 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 [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);
203 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
204 pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
205 pEZPlot->ezset ("ylength 0.25");
206 pEZPlot->ezset ("yporigin 0.50");
207 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
211 shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
213 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
214 pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
215 pEZPlot->ezset ("ylength 0.25");
216 pEZPlot->ezset ("yporigin 0.75");
217 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
222 for (int i = 0; i < m_nFilterPoints; i++) {
223 m_adFilter[i] /= m_dSignalInc;
226 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
227 for (int i = 0; i < m_nFilterPoints; i++)
228 m_adFilter[i] *= 0.5;
229 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
230 for (int i = 0; i < m_nFilterPoints; i++) {
231 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
232 double sinScale = sin (iDetFromZero * m_dSignalInc);
233 if (fabs(sinScale) < 1E-7)
236 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
237 double dScale = 0.5 * sinScale * sinScale;
238 m_adFilter[i] *= dScale;
241 } // if (spatial filtering)
243 else if (m_bFrequencyFiltering) { // Frequency-based filtering
245 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
246 // calculate number of filter points with zeropadding
247 m_nFilterPoints = m_nSignalPoints;
248 if (m_iZeropad > 0) {
249 double logBase2 = log(m_nFilterPoints) / log(2);
250 int nextPowerOf2 = static_cast<int>(floor(logBase2));
251 if (logBase2 != floor(logBase2))
253 nextPowerOf2 += (m_iZeropad - 1);
254 m_nFilterPoints = 1 << nextPowerOf2;
256 if (m_traceLevel >= Trace::TRACE_CONSOLE)
257 cout << "nFilterPoints = " << m_nFilterPoints << endl;
260 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
262 if (m_nFilterPoints % 2) { // Odd
263 m_dFilterMin = -1. / (2 * m_dSignalInc);
264 m_dFilterMax = 1. / (2 * m_dSignalInc);
265 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
267 m_dFilterMin = -1. / (2 * m_dSignalInc);
268 m_dFilterMax = 1. / (2 * m_dSignalInc);
269 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
270 m_dFilterMax -= m_dFilterInc;
273 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
274 m_adFilter = new double [m_nFilterPoints];
275 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
277 // This doesn't work!
278 // Need to add filtering for divergent geometries & Frequency/Direct filtering
279 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
280 for (int i = 0; i < m_nFilterPoints; i++)
281 m_adFilter[i] *= 0.5;
282 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
283 for (int i = 0; i < m_nFilterPoints; i++) {
284 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
285 double sinScale = sin (iDetFromZero * m_dSignalInc);
286 if (fabs(sinScale) < 1E-7)
289 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
290 double dScale = 0.5 * sinScale * sinScale;
291 m_adFilter[i] *= dScale;
295 EZPlot* pEZPlot = NULL;
296 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
297 pEZPlot = new EZPlot (*pSGP);
298 pEZPlot->ezset ("title Filter Filter: Natural Order");
299 pEZPlot->ezset ("ylength 0.50");
300 pEZPlot->ezset ("yporigin 0.00");
301 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
305 shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
307 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
308 pEZPlot->ezset ("title Filter Filter: Fourier Order");
309 pEZPlot->ezset ("ylength 0.50");
310 pEZPlot->ezset ("yporigin 0.50");
311 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
316 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
317 // calculate number of filter points with zeropadding
318 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
319 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
320 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
321 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
322 m_nFilterPoints = nSpatialPoints;
323 if (m_iZeropad > 0) {
324 double logBase2 = log(nSpatialPoints) / log(2);
325 int nextPowerOf2 = static_cast<int>(floor(logBase2));
326 if (logBase2 != floor(logBase2))
328 nextPowerOf2 += (m_iZeropad - 1);
329 m_nFilterPoints = 1 << nextPowerOf2;
331 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
333 if (m_traceLevel >= Trace::TRACE_CONSOLE)
334 cout << "nFilterPoints = " << m_nFilterPoints << endl;
336 double adSpatialFilter [m_nFilterPoints];
337 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
338 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
340 EZPlot* pEZPlot = NULL;
341 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
342 pEZPlot = new EZPlot (*pSGP);
343 pEZPlot->ezset ("title Spatial Filter: Natural Order");
344 pEZPlot->ezset ("ylength 0.50");
345 pEZPlot->ezset ("yporigin 0.00");
346 pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
351 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
352 for (int i = 0; i < m_nFilterPoints; i++)
353 adSpatialFilter[i] *= 0.5;
354 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
355 for (int i = 0; i < m_nFilterPoints; i++) {
356 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
357 double sinScale = sin (iDetFromZero * m_dSignalInc);
358 if (fabs(sinScale) < 1E-7)
361 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
362 double dScale = 0.5 * sinScale * sinScale;
363 adSpatialFilter[i] *= dScale;
366 for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
367 adSpatialFilter[i] = 0;
369 m_adFilter = new double [m_nFilterPoints];
370 complex<double> acInverseFilter [m_nFilterPoints];
371 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
372 for (int i = 0; i < m_nFilterPoints; i++)
373 m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
375 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
376 pEZPlot->ezset ("title Spatial Filter: Inverse");
377 pEZPlot->ezset ("ylength 0.50");
378 pEZPlot->ezset ("yporigin 0.50");
379 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
387 // precalculate sin and cosine tables for fourier transform
388 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
389 int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
390 double angleIncrement = (2. * PI) / m_nFilterPoints;
391 m_adFourierCosTable = new double[ nFourier ];
392 m_adFourierSinTable = new double[ nFourier ];
394 for (int i = 0; i < nFourier; i++) {
395 m_adFourierCosTable[i] = cos (angle);
396 m_adFourierSinTable[i] = sin (angle);
397 angle += angleIncrement;
402 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
403 for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
404 m_adFilter[i] /= m_nFilterPoints;
407 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
408 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
409 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
410 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
411 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
412 for (int i = 0; i < m_nFilterPoints; i++)
413 m_adRealFftInput[i] = 0;
414 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
415 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
416 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
417 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
418 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
419 for (int i = 0; i < m_nFilterPoints; i++)
420 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
421 for (int i = 0; i < m_nOutputPoints; i++)
422 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
428 ProcessSignal::~ProcessSignal (void)
430 delete [] m_adFourierSinTable;
431 delete [] m_adFourierCosTable;
432 delete [] m_adFilter;
435 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
436 fftw_destroy_plan(m_complexPlanForward);
437 fftw_destroy_plan(m_complexPlanBackward);
438 delete [] m_adComplexFftInput;
439 delete [] m_adComplexFftSignal;
441 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
442 rfftw_destroy_plan(m_realPlanForward);
443 rfftw_destroy_plan(m_realPlanBackward);
444 delete [] m_adRealFftInput;
445 delete [] m_adRealFftSignal;
451 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
453 int fmID = FILTER_METHOD_INVALID;
455 for (int i = 0; i < s_iFilterMethodCount; i++)
456 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
465 ProcessSignal::convertFilterMethodIDToName (const int fmID)
467 static const char *name = "";
469 if (fmID >= 0 && fmID < s_iFilterMethodCount)
470 return (s_aszFilterMethodName [fmID]);
476 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
478 static const char *title = "";
480 if (fmID >= 0 && fmID < s_iFilterMethodCount)
481 return (s_aszFilterMethodTitle [fmID]);
488 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
490 int fgID = FILTER_GENERATION_INVALID;
492 for (int i = 0; i < s_iFilterGenerationCount; i++)
493 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
502 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
504 static const char *name = "";
506 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
507 return (s_aszFilterGenerationName [fgID]);
513 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
515 static const char *name = "";
517 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
518 return (s_aszFilterGenerationTitle [fgID]);
524 ProcessSignal::filterSignal (const float constInput[], double output[]) const
526 double input [m_nSignalPoints];
527 for (int i = 0; i < m_nSignalPoints; i++)
528 input[i] = constInput[i];
530 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
531 for (int i = 0; i < m_nSignalPoints; i++) {
532 int iDetFromCenter = i - (m_nSignalPoints / 2);
533 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
535 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
536 for (int i = 0; i < m_nSignalPoints; i++) {
537 int iDetFromCenter = i - (m_nSignalPoints / 2);
538 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
541 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
542 for (int i = 0; i < m_nSignalPoints; i++)
543 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
544 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
545 double inputSignal[m_nFilterPoints];
546 for (int i = 0; i < m_nSignalPoints; i++)
547 inputSignal[i] = input[i];
548 for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
549 inputSignal[i] = 0; // zeropad
550 complex<double> fftSignal[m_nFilterPoints];
551 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
552 for (int i = 0; i < m_nFilterPoints; i++)
553 fftSignal[i] *= m_adFilter[i];
554 double inverseFourier[m_nFilterPoints];
555 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
556 for (int i = 0; i < m_nSignalPoints; i++)
557 output[i] = inverseFourier[i];
558 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
559 double inputSignal[m_nFilterPoints];
560 for (int i = 0; i < m_nSignalPoints; i++)
561 inputSignal[i] = input[i];
562 for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
563 inputSignal[i] = 0; // zeropad
564 complex<double> fftSignal[m_nFilterPoints];
565 finiteFourierTransform (inputSignal, fftSignal, -1);
566 for (int i = 0; i < m_nFilterPoints; i++)
567 fftSignal[i] *= m_adFilter[i];
568 double inverseFourier[m_nFilterPoints];
569 finiteFourierTransform (fftSignal, inverseFourier, 1);
570 for (int i = 0; i < m_nSignalPoints; i++)
571 output[i] = inverseFourier[i];
574 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
575 for (int i = 0; i < m_nSignalPoints; i++)
576 m_adRealFftInput[i] = input[i];
578 fftw_real fftOutput [ m_nFilterPoints ];
579 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
580 for (int i = 0; i < m_nFilterPoints; i++)
581 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
582 for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
583 m_adRealFftSignal[i] = 0;
585 fftw_real ifftOutput [ m_nOutputPoints ];
586 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
587 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
588 output[i] = ifftOutput[i];
589 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
590 for (int i = 0; i < m_nSignalPoints; i++)
591 m_adComplexFftInput[i].re = input[i];
593 fftw_complex fftOutput [ m_nFilterPoints ];
594 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
595 for (int i = 0; i < m_nFilterPoints; i++) {
596 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
597 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
599 fftw_complex ifftOutput [ m_nOutputPoints ];
600 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
601 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
602 output[i] = ifftOutput[i].re;
609 * convolve Discrete convolution of two functions
612 * r = convolve (f1, f2, dx, n, np, func_type)
613 * double r Convolved result
614 * double f1[], f2[] Functions to be convolved
615 * double dx Difference between successive x values
616 * int n Array index to center convolution about
617 * int np Number of points in f1 array
618 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
621 * f1 is the projection data, its indices range from 0 to np - 1.
622 * The index for f2, the filter, ranges from -(np-1) to (np-1).
623 * There are 3 ways to handle the negative vertices of f2:
624 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
625 * All filters used in reconstruction are even.
626 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
627 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
628 * for negative indices. Since f2 must range from -(np-1) to (np-1),
629 * if we add (np - 1) to f2's array index, then f2's index will
630 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
631 * stored at f2[np-1].
635 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
639 #if UNOPTIMIZED_CONVOLUTION
640 for (int i = 0; i < np; i++)
641 sum += func[i] * m_adFilter[n - i + (np - 1)];
643 double* f2 = m_adFilter + n + (np - 1);
644 for (int i = 0; i < np; i++)
645 sum += *func++ * *f2--;
653 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
657 #if UNOPTIMIZED_CONVOLUTION
658 for (int i = 0; i < np; i++)
659 sum += func[i] * m_adFilter[n - i + (np - 1)];
661 double* f2 = m_adFilter + n + (np - 1);
662 for (int i = 0; i < np; i++)
663 sum += *func++ * *f2--;
671 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
673 complex<double> complexOutput[n];
675 finiteFourierTransform (input, complexOutput, n, direction);
676 for (int i = 0; i < n; i++)
677 output[i] = complexOutput[i].real();
681 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
688 double angleIncrement = direction * 2 * PI / n;
689 for (int i = 0; i < n; i++) {
692 for (int j = 0; j < n; j++) {
693 double angle = i * j * angleIncrement;
694 sumReal += input[j] * cos(angle);
695 sumImag += input[j] * sin(angle);
701 output[i] = complex<double> (sumReal, sumImag);
707 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
714 double angleIncrement = direction * 2 * PI / n;
715 for (int i = 0; i < n; i++) {
716 complex<double> sum (0,0);
717 for (int j = 0; j < n; j++) {
718 double angle = i * j * angleIncrement;
719 complex<double> exponentTerm (cos(angle), sin(angle));
720 sum += input[j] * exponentTerm;
730 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
737 double angleIncrement = direction * 2 * PI / n;
738 for (int i = 0; i < n; i++) {
740 for (int j = 0; j < n; j++) {
741 double angle = i * j * angleIncrement;
742 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
751 // Table-based routines
754 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
761 for (int i = 0; i < m_nFilterPoints; i++) {
762 double sumReal = 0, sumImag = 0;
763 for (int j = 0; j < m_nFilterPoints; j++) {
764 int tableIndex = i * j;
766 sumReal += input[j] * m_adFourierCosTable[tableIndex];
767 sumImag += input[j] * m_adFourierSinTable[tableIndex];
769 sumReal += input[j] * m_adFourierCosTable[tableIndex];
770 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
774 sumReal /= m_nFilterPoints;
775 sumImag /= m_nFilterPoints;
777 output[i] = complex<double> (sumReal, sumImag);
781 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
783 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
790 for (int i = 0; i < m_nFilterPoints; i++) {
791 double sumReal = 0, sumImag = 0;
792 for (int j = 0; j < m_nFilterPoints; j++) {
793 int tableIndex = i * j;
795 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
796 - input[j].imag() * m_adFourierSinTable[tableIndex];
797 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
798 + input[j].imag() * m_adFourierCosTable[tableIndex];
800 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
801 - input[j].imag() * -m_adFourierSinTable[tableIndex];
802 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
803 + input[j].imag() * m_adFourierCosTable[tableIndex];
807 sumReal /= m_nFilterPoints;
808 sumImag /= m_nFilterPoints;
810 output[i] = complex<double> (sumReal, sumImag);
815 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
822 for (int i = 0; i < m_nFilterPoints; i++) {
824 for (int j = 0; j < m_nFilterPoints; j++) {
825 int tableIndex = i * j;
827 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
828 - input[j].imag() * m_adFourierSinTable[tableIndex];
830 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
831 - input[j].imag() * -m_adFourierSinTable[tableIndex];
835 sumReal /= m_nFilterPoints;
841 // Odd Number of Points
842 // Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
843 // Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
844 // Even Number of Points
845 // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
846 // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
849 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
851 double* pdTemp = new double [n];
853 int iHalfN = (n - 1) / 2;
855 pdTemp[0] = pdVector[iHalfN];
856 for (int i = 0; i < iHalfN; i++)
857 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
858 for (int i = 0; i < iHalfN; i++)
859 pdTemp[i + iHalfN + 1] = pdVector[i];
862 pdTemp[0] = pdVector[iHalfN];
863 for (int i = 0; i < iHalfN; i++)
864 pdTemp[i + 1] = pdVector[i + iHalfN];
865 for (int i = 0; i < iHalfN - 1; i++)
866 pdTemp[i + iHalfN + 1] = pdVector[i];
869 for (int i = 0; i < n; i++)
870 pdVector[i] = pdTemp[i];
876 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
878 double* pdTemp = new double [n];
880 int iHalfN = (n - 1) / 2;
882 pdTemp[iHalfN] = pdVector[0];
883 for (int i = 0; i < iHalfN; i++)
884 pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
885 for (int i = 0; i < iHalfN; i++)
886 pdTemp[i] = pdVector[i + iHalfN + 1];
889 pdTemp[iHalfN] = pdVector[0];
890 for (int i = 0; i < iHalfN; i++)
891 pdTemp[i] = pdVector[i + iHalfN];
892 for (int i = 0; i < iHalfN - 1; i++)
893 pdTemp[i + iHalfN + 1] = pdVector[i+1];
896 for (int i = 0; i < n; i++)
897 pdVector[i] = pdTemp[i];