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.6 2000/09/02 05:10:39 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);
181 EZPlot* pEZPlot = NULL;
182 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
183 pEZPlot = new EZPlot (*pSGP);
184 pEZPlot->ezset ("title Filter Response: Natural Order");
185 pEZPlot->ezset ("ylength 0.25");
186 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
190 shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
191 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
192 pEZPlot->ezset ("title Filter Response: Fourier Order");
193 pEZPlot->ezset ("ylength 0.25");
194 pEZPlot->ezset ("yporigin 0.25");
195 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
198 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
199 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
200 pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
201 pEZPlot->ezset ("ylength 0.25");
202 pEZPlot->ezset ("yporigin 0.50");
203 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
206 shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
207 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
208 pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
209 pEZPlot->ezset ("ylength 0.25");
210 pEZPlot->ezset ("yporigin 0.75");
211 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
215 for (int i = 0; i < m_nFilterPoints; i++) {
216 m_adFilter[i] /= m_dSignalInc;
219 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
220 for (int i = 0; i < m_nFilterPoints; i++)
221 m_adFilter[i] *= 0.5;
222 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
223 for (int i = 0; i < m_nFilterPoints; i++) {
224 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
225 double sinScale = sin (iDetFromZero * m_dSignalInc);
226 if (fabs(sinScale) < 1E-7)
229 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
230 double dScale = 0.5 * sinScale * sinScale;
231 m_adFilter[i] *= dScale;
234 } // if (spatial filtering)
236 else if (m_bFrequencyFiltering) { // Frequency-based filtering
238 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
239 // calculate number of filter points with zeropadding
240 m_nFilterPoints = m_nSignalPoints;
241 if (m_iZeropad > 0) {
242 double logBase2 = log(m_nFilterPoints) / log(2);
243 int nextPowerOf2 = static_cast<int>(floor(logBase2));
244 if (logBase2 != floor(logBase2))
246 nextPowerOf2 += (m_iZeropad - 1);
247 m_nFilterPoints = 1 << nextPowerOf2;
249 if (m_traceLevel >= Trace::TRACE_CONSOLE)
250 cout << "nFilterPoints = " << m_nFilterPoints << endl;
253 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
255 if (m_nFilterPoints % 2) { // Odd
256 m_dFilterMin = -1. / (2 * m_dSignalInc);
257 m_dFilterMax = 1. / (2 * m_dSignalInc);
258 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
260 m_dFilterMin = -1. / (2 * m_dSignalInc);
261 m_dFilterMax = 1. / (2 * m_dSignalInc);
262 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
263 m_dFilterMax -= m_dFilterInc;
266 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
267 m_adFilter = new double [m_nFilterPoints];
268 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
270 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
271 for (int i = 0; i < m_nFilterPoints; i++)
272 m_adFilter[i] *= 0.5;
273 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
274 for (int i = 0; i < m_nFilterPoints; i++) {
275 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
276 double sinScale = sin (iDetFromZero * m_dSignalInc);
277 if (fabs(sinScale) < 1E-7)
280 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
281 double dScale = 0.5 * sinScale * sinScale;
282 // m_adFilter[i] *= dScale;
285 EZPlot* pEZPlot = NULL;
286 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
287 pEZPlot = new EZPlot (*pSGP);
288 pEZPlot->ezset ("title Filter Filter: Natural Order");
289 pEZPlot->ezset ("ylength 0.50");
290 pEZPlot->ezset ("yporigin 0.00");
291 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
294 shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
295 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
296 pEZPlot->ezset ("title Filter Filter: Fourier Order");
297 pEZPlot->ezset ("ylength 0.50");
298 pEZPlot->ezset ("yporigin 0.50");
299 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
303 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
304 // calculate number of filter points with zeropadding
305 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
306 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
307 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
308 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
309 m_nFilterPoints = nSpatialPoints;
310 if (m_iZeropad > 0) {
311 double logBase2 = log(nSpatialPoints) / log(2);
312 int nextPowerOf2 = static_cast<int>(floor(logBase2));
313 if (logBase2 != floor(logBase2))
315 nextPowerOf2 += (m_iZeropad - 1);
316 m_nFilterPoints = 1 << nextPowerOf2;
318 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
320 if (m_traceLevel >= Trace::TRACE_CONSOLE)
321 cout << "nFilterPoints = " << m_nFilterPoints << endl;
323 double adSpatialFilter [m_nFilterPoints];
324 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
325 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
326 EZPlot* pEZPlot = NULL;
327 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
328 pEZPlot = new EZPlot (*pSGP);
329 pEZPlot->ezset ("title Spatial Filter: Natural Order");
330 pEZPlot->ezset ("ylength 0.50");
331 pEZPlot->ezset ("yporigin 0.00");
332 pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
336 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
337 for (int i = 0; i < m_nFilterPoints; i++)
338 adSpatialFilter[i] *= 0.5;
339 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
340 for (int i = 0; i < m_nFilterPoints; i++) {
341 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
342 double sinScale = sin (iDetFromZero * m_dSignalInc);
343 if (fabs(sinScale) < 1E-7)
346 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
347 double dScale = 0.5 * sinScale * sinScale;
348 adSpatialFilter[i] *= dScale;
351 for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
352 adSpatialFilter[i] = 0;
354 m_adFilter = new double [m_nFilterPoints];
355 complex<double> acInverseFilter [m_nFilterPoints];
356 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
357 for (int i = 0; i < m_nFilterPoints; i++)
358 m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
359 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
360 pEZPlot->ezset ("title Spatial Filter: Inverse");
361 pEZPlot->ezset ("ylength 0.50");
362 pEZPlot->ezset ("yporigin 0.50");
363 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
370 // precalculate sin and cosine tables for fourier transform
371 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
372 int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
373 double angleIncrement = (2. * PI) / m_nFilterPoints;
374 m_adFourierCosTable = new double[ nFourier ];
375 m_adFourierSinTable = new double[ nFourier ];
377 for (int i = 0; i < nFourier; i++) {
378 m_adFourierCosTable[i] = cos (angle);
379 m_adFourierSinTable[i] = sin (angle);
380 angle += angleIncrement;
385 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
386 for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
387 m_adFilter[i] /= m_nFilterPoints;
390 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
391 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
392 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
393 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
394 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
395 for (int i = 0; i < m_nFilterPoints; i++)
396 m_adRealFftInput[i] = 0;
397 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
398 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
399 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
400 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
401 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
402 for (int i = 0; i < m_nFilterPoints; i++)
403 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
404 for (int i = 0; i < m_nOutputPoints; i++)
405 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
411 ProcessSignal::~ProcessSignal (void)
413 delete [] m_adFourierSinTable;
414 delete [] m_adFourierCosTable;
415 delete [] m_adFilter;
418 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
419 fftw_destroy_plan(m_complexPlanForward);
420 fftw_destroy_plan(m_complexPlanBackward);
421 delete [] m_adComplexFftInput;
422 delete [] m_adComplexFftSignal;
424 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
425 rfftw_destroy_plan(m_realPlanForward);
426 rfftw_destroy_plan(m_realPlanBackward);
427 delete [] m_adRealFftInput;
428 delete [] m_adRealFftSignal;
434 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
436 int fmID = FILTER_METHOD_INVALID;
438 for (int i = 0; i < s_iFilterMethodCount; i++)
439 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
448 ProcessSignal::convertFilterMethodIDToName (const int fmID)
450 static const char *name = "";
452 if (fmID >= 0 && fmID < s_iFilterMethodCount)
453 return (s_aszFilterMethodName [fmID]);
459 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
461 static const char *title = "";
463 if (fmID >= 0 && fmID < s_iFilterMethodCount)
464 return (s_aszFilterMethodTitle [fmID]);
471 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
473 int fgID = FILTER_GENERATION_INVALID;
475 for (int i = 0; i < s_iFilterGenerationCount; i++)
476 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
485 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
487 static const char *name = "";
489 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
490 return (s_aszFilterGenerationName [fgID]);
496 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
498 static const char *name = "";
500 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
501 return (s_aszFilterGenerationTitle [fgID]);
507 ProcessSignal::filterSignal (const float constInput[], double output[]) const
509 double input [m_nSignalPoints];
510 for (int i = 0; i < m_nSignalPoints; i++)
511 input[i] = constInput[i];
513 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
514 for (int i = 0; i < m_nSignalPoints; i++) {
515 int iDetFromCenter = i - (m_nSignalPoints / 2);
516 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
518 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
519 for (int i = 0; i < m_nSignalPoints; i++) {
520 int iDetFromCenter = i - (m_nSignalPoints / 2);
521 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
524 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
525 for (int i = 0; i < m_nSignalPoints; i++)
526 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
527 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
528 double inputSignal[m_nFilterPoints];
529 for (int i = 0; i < m_nSignalPoints; i++)
530 inputSignal[i] = input[i];
531 for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
532 inputSignal[i] = 0; // zeropad
533 complex<double> fftSignal[m_nFilterPoints];
534 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
535 for (int i = 0; i < m_nFilterPoints; i++)
536 fftSignal[i] *= m_adFilter[i];
537 double inverseFourier[m_nFilterPoints];
538 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
539 for (int i = 0; i < m_nSignalPoints; i++)
540 output[i] = inverseFourier[i];
541 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
542 double inputSignal[m_nFilterPoints];
543 for (int i = 0; i < m_nSignalPoints; i++)
544 inputSignal[i] = input[i];
545 for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
546 inputSignal[i] = 0; // zeropad
547 complex<double> fftSignal[m_nFilterPoints];
548 finiteFourierTransform (inputSignal, fftSignal, -1);
549 for (int i = 0; i < m_nFilterPoints; i++)
550 fftSignal[i] *= m_adFilter[i];
551 double inverseFourier[m_nFilterPoints];
552 finiteFourierTransform (fftSignal, inverseFourier, 1);
553 for (int i = 0; i < m_nSignalPoints; i++)
554 output[i] = inverseFourier[i];
557 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
558 for (int i = 0; i < m_nSignalPoints; i++)
559 m_adRealFftInput[i] = input[i];
561 fftw_real fftOutput [ m_nFilterPoints ];
562 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
563 for (int i = 0; i < m_nFilterPoints; i++)
564 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
565 for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
566 m_adRealFftSignal[i] = 0;
568 fftw_real ifftOutput [ m_nOutputPoints ];
569 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
570 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
571 output[i] = ifftOutput[i];
572 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
573 for (int i = 0; i < m_nSignalPoints; i++)
574 m_adComplexFftInput[i].re = input[i];
576 fftw_complex fftOutput [ m_nFilterPoints ];
577 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
578 for (int i = 0; i < m_nFilterPoints; i++) {
579 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
580 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
582 fftw_complex ifftOutput [ m_nOutputPoints ];
583 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
584 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
585 output[i] = ifftOutput[i].re;
592 * convolve Discrete convolution of two functions
595 * r = convolve (f1, f2, dx, n, np, func_type)
596 * double r Convolved result
597 * double f1[], f2[] Functions to be convolved
598 * double dx Difference between successive x values
599 * int n Array index to center convolution about
600 * int np Number of points in f1 array
601 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
604 * f1 is the projection data, its indices range from 0 to np - 1.
605 * The index for f2, the filter, ranges from -(np-1) to (np-1).
606 * There are 3 ways to handle the negative vertices of f2:
607 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
608 * All filters used in reconstruction are even.
609 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
610 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
611 * for negative indices. Since f2 must range from -(np-1) to (np-1),
612 * if we add (np - 1) to f2's array index, then f2's index will
613 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
614 * stored at f2[np-1].
618 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
622 #if UNOPTIMIZED_CONVOLUTION
623 for (int i = 0; i < np; i++)
624 sum += func[i] * m_adFilter[n - i + (np - 1)];
626 double* f2 = m_adFilter + n + (np - 1);
627 for (int i = 0; i < np; i++)
628 sum += *func++ * *f2--;
636 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
640 #if UNOPTIMIZED_CONVOLUTION
641 for (int i = 0; i < np; i++)
642 sum += func[i] * m_adFilter[n - i + (np - 1)];
644 double* f2 = m_adFilter + n + (np - 1);
645 for (int i = 0; i < np; i++)
646 sum += *func++ * *f2--;
654 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
656 complex<double> complexOutput[n];
658 finiteFourierTransform (input, complexOutput, n, direction);
659 for (int i = 0; i < n; i++)
660 output[i] = complexOutput[i].real();
664 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
671 double angleIncrement = direction * 2 * PI / n;
672 for (int i = 0; i < n; i++) {
675 for (int j = 0; j < n; j++) {
676 double angle = i * j * angleIncrement;
677 sumReal += input[j] * cos(angle);
678 sumImag += input[j] * sin(angle);
684 output[i] = complex<double> (sumReal, sumImag);
690 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
697 double angleIncrement = direction * 2 * PI / n;
698 for (int i = 0; i < n; i++) {
699 complex<double> sum (0,0);
700 for (int j = 0; j < n; j++) {
701 double angle = i * j * angleIncrement;
702 complex<double> exponentTerm (cos(angle), sin(angle));
703 sum += input[j] * exponentTerm;
713 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
720 double angleIncrement = direction * 2 * PI / n;
721 for (int i = 0; i < n; i++) {
723 for (int j = 0; j < n; j++) {
724 double angle = i * j * angleIncrement;
725 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
734 // Table-based routines
737 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
744 for (int i = 0; i < m_nFilterPoints; i++) {
745 double sumReal = 0, sumImag = 0;
746 for (int j = 0; j < m_nFilterPoints; j++) {
747 int tableIndex = i * j;
749 sumReal += input[j] * m_adFourierCosTable[tableIndex];
750 sumImag += input[j] * m_adFourierSinTable[tableIndex];
752 sumReal += input[j] * m_adFourierCosTable[tableIndex];
753 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
757 sumReal /= m_nFilterPoints;
758 sumImag /= m_nFilterPoints;
760 output[i] = complex<double> (sumReal, sumImag);
764 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
766 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
773 for (int i = 0; i < m_nFilterPoints; i++) {
774 double sumReal = 0, sumImag = 0;
775 for (int j = 0; j < m_nFilterPoints; j++) {
776 int tableIndex = i * j;
778 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
779 - input[j].imag() * m_adFourierSinTable[tableIndex];
780 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
781 + input[j].imag() * m_adFourierCosTable[tableIndex];
783 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
784 - input[j].imag() * -m_adFourierSinTable[tableIndex];
785 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
786 + input[j].imag() * m_adFourierCosTable[tableIndex];
790 sumReal /= m_nFilterPoints;
791 sumImag /= m_nFilterPoints;
793 output[i] = complex<double> (sumReal, sumImag);
798 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
805 for (int i = 0; i < m_nFilterPoints; i++) {
807 for (int j = 0; j < m_nFilterPoints; j++) {
808 int tableIndex = i * j;
810 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
811 - input[j].imag() * m_adFourierSinTable[tableIndex];
813 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
814 - input[j].imag() * -m_adFourierSinTable[tableIndex];
818 sumReal /= m_nFilterPoints;
824 // Odd Number of Points
825 // Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
826 // Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
827 // Even Number of Points
828 // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
829 // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
832 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
834 double* pdTemp = new double [n];
836 int iHalfN = (n - 1) / 2;
838 pdTemp[0] = pdVector[iHalfN];
839 for (int i = 0; i < iHalfN; i++)
840 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
841 for (int i = 0; i < iHalfN; i++)
842 pdTemp[i + iHalfN + 1] = pdVector[i];
845 pdTemp[0] = pdVector[iHalfN];
846 for (int i = 0; i < iHalfN; i++)
847 pdTemp[i + 1] = pdVector[i + iHalfN];
848 for (int i = 0; i < iHalfN - 1; i++)
849 pdTemp[i + iHalfN + 1] = pdVector[i];
852 for (int i = 0; i < n; i++)
853 pdVector[i] = pdTemp[i];
859 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
861 double* pdTemp = new double [n];
863 int iHalfN = (n - 1) / 2;
865 pdTemp[iHalfN] = pdVector[0];
866 for (int i = 0; i < iHalfN; i++)
867 pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
868 for (int i = 0; i < iHalfN; i++)
869 pdTemp[i] = pdVector[i + iHalfN + 1];
872 pdTemp[iHalfN] = pdVector[0];
873 for (int i = 0; i < iHalfN; i++)
874 pdTemp[i] = pdVector[i + iHalfN];
875 for (int i = 0; i < iHalfN - 1; i++)
876 pdTemp[i + iHalfN + 1] = pdVector[i+1];
879 for (int i = 0; i < n; i++)
880 pdVector[i] = pdTemp[i];