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.13 2001/01/02 05:34:57 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)
120 m_idFilter = idFilter;
121 m_idDomain = idDomain;
122 m_idFilterMethod = idFilterMethod;
123 m_idFilterGeneration = idFilterGeneration;
124 m_idGeometry = iGeometry;
125 m_dFocalLength = dFocalLength;
127 if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
131 m_traceLevel = iTraceLevel;
132 m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
133 m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
134 m_dBandwidth = dBandwidth;
135 m_nSignalPoints = nSignalPoints;
136 m_dSignalInc = dSignalIncrement;
137 m_dFilterParam = dFilterParam;
138 m_iZeropad = iZeropad;
139 m_iPreinterpolationFactor = iPreinterpolationFactor;
141 // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
142 // through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear
143 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
148 if (m_idFilterMethod == FILTER_METHOD_FFT) {
150 m_idFilterMethod = FILTER_METHOD_RFFTW;
153 m_failMessage = "FFT not yet implemented";
158 bool m_bFrequencyFiltering = true;
159 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
160 m_bFrequencyFiltering = false;
162 // Spatial-based filtering
163 if (! m_bFrequencyFiltering) {
165 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
166 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
167 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
168 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
169 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
170 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
171 m_adFilter = new double[ m_nFilterPoints ];
172 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
173 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
174 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
175 m_dFilterMin = -1. / (2 * m_dSignalInc);
176 m_dFilterMax = 1. / (2 * m_dSignalInc);
177 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
178 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
179 m_adFilter = new double[ m_nFilterPoints ];
180 double* adFrequencyFilter = new double [m_nFilterPoints];
181 filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
183 EZPlot* pEZPlot = NULL;
184 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
185 pEZPlot = new EZPlot ();
186 pEZPlot->ezset ("title Filter Response: Natural Order");
187 pEZPlot->ezset ("ylength 0.25");
188 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
189 pEZPlot->plot (pSGP);
192 Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
194 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
195 pEZPlot->ezset ("title Filter Response: Fourier Order");
196 pEZPlot->ezset ("ylength 0.25");
197 pEZPlot->ezset ("yporigin 0.25");
198 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
199 pEZPlot->plot (pSGP);
202 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
203 delete adFrequencyFilter;
\r
205 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
206 pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
207 pEZPlot->ezset ("ylength 0.25");
208 pEZPlot->ezset ("yporigin 0.50");
209 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
210 pEZPlot->plot (pSGP);
213 Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
215 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
216 pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
217 pEZPlot->ezset ("ylength 0.25");
218 pEZPlot->ezset ("yporigin 0.75");
219 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
220 pEZPlot->plot (pSGP);
224 for (i = 0; i < m_nFilterPoints; i++) {
225 m_adFilter[i] /= m_dSignalInc;
228 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
229 for (i = 0; i < m_nFilterPoints; i++)
230 m_adFilter[i] *= 0.5;
231 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
232 for (i = 0; i < m_nFilterPoints; i++) {
233 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
234 double sinScale = sin (iDetFromZero * m_dSignalInc);
235 if (fabs(sinScale) < 1E-7)
238 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
239 double dScale = 0.5 * sinScale * sinScale;
240 m_adFilter[i] *= dScale;
243 } // if (spatial filtering)
245 else if (m_bFrequencyFiltering) { // Frequency-based filtering
247 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
248 // calculate number of filter points with zeropadding
249 m_nFilterPoints = m_nSignalPoints;
250 if (m_iZeropad > 0) {
251 double logBase2 = log(m_nFilterPoints) / log(2);
252 int nextPowerOf2 = static_cast<int>(floor(logBase2));
253 if (logBase2 != floor(logBase2))
255 nextPowerOf2 += (m_iZeropad - 1);
256 m_nFilterPoints = 1 << nextPowerOf2;
258 if (m_traceLevel >= Trace::TRACE_CONSOLE)
259 std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
262 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
264 if (m_nFilterPoints % 2) { // Odd
265 m_dFilterMin = -1. / (2 * m_dSignalInc);
266 m_dFilterMax = 1. / (2 * m_dSignalInc);
267 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
269 m_dFilterMin = -1. / (2 * m_dSignalInc);
270 m_dFilterMax = 1. / (2 * m_dSignalInc);
271 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
272 m_dFilterMax -= m_dFilterInc;
275 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
276 m_adFilter = new double [m_nFilterPoints];
277 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
279 // This doesn't work!
280 // Need to add filtering for divergent geometries & Frequency/Direct filtering
281 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
282 for (i = 0; i < m_nFilterPoints; i++)
283 m_adFilter[i] *= 0.5;
284 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
285 for (i = 0; i < m_nFilterPoints; i++) {
286 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
287 double sinScale = sin (iDetFromZero * m_dSignalInc);
288 if (fabs(sinScale) < 1E-7)
291 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
292 double dScale = 0.5 * sinScale * sinScale;
293 m_adFilter[i] *= dScale;
297 EZPlot* pEZPlot = NULL;
298 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
299 pEZPlot = new EZPlot;
300 pEZPlot->ezset ("title Filter Filter: Natural Order");
301 pEZPlot->ezset ("ylength 0.50");
302 pEZPlot->ezset ("yporigin 0.00");
303 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
304 pEZPlot->plot (pSGP);
307 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
309 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
310 pEZPlot->ezset ("title Filter Filter: Fourier Order");
311 pEZPlot->ezset ("ylength 0.50");
312 pEZPlot->ezset ("yporigin 0.50");
313 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
314 pEZPlot->plot (pSGP);
318 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
319 // calculate number of filter points with zeropadding
320 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
321 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
322 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
323 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
324 m_nFilterPoints = nSpatialPoints;
325 if (m_iZeropad > 0) {
326 double logBase2 = log(nSpatialPoints) / log(2);
327 int nextPowerOf2 = static_cast<int>(floor(logBase2));
328 if (logBase2 != floor(logBase2))
330 nextPowerOf2 += (m_iZeropad - 1);
331 m_nFilterPoints = 1 << nextPowerOf2;
333 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
335 if (m_traceLevel >= Trace::TRACE_CONSOLE)
336 std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
338 double* adSpatialFilter = new double [m_nFilterPoints];
\r
339 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
340 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
342 EZPlot* pEZPlot = NULL;
343 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
344 pEZPlot = new EZPlot;
345 pEZPlot->ezset ("title Spatial Filter: Natural Order");
346 pEZPlot->ezset ("ylength 0.50");
347 pEZPlot->ezset ("yporigin 0.00");
348 pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
349 pEZPlot->plot (pSGP);
353 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
354 for (i = 0; i < m_nFilterPoints; i++)
355 adSpatialFilter[i] *= 0.5;
356 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
357 for (i = 0; i < m_nFilterPoints; i++) {
358 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
359 double sinScale = sin (iDetFromZero * m_dSignalInc);
360 if (fabs(sinScale) < 1E-7)
363 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
364 double dScale = 0.5 * sinScale * sinScale;
365 adSpatialFilter[i] *= dScale;
368 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
369 adSpatialFilter[i] = 0;
371 m_adFilter = new double [m_nFilterPoints];
372 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
\r
373 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
374 delete adSpatialFilter;
\r
375 for (i = 0; i < m_nFilterPoints; i++)
376 m_adFilter[i] = std::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);
384 pEZPlot->plot (pSGP);
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 (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 (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 (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 (i = 0; i < m_nFilterPoints; i++)
424 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
425 for (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 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
556 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
\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, BACKWARD);
\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 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
573 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
\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, BACKWARD);
\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 (i = 0; i < m_nSignalPoints; i++)
587 m_adRealFftInput[i] = input[i];
589 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
590 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
591 for (i = 0; i < m_nFilterPoints; i++)
592 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
\r
594 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
595 m_adRealFftSignal[i] = 0;
597 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
598 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
599 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
600 output[i] = ifftOutput[i];
\r
601 delete [] ifftOutput;
602 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
603 for (i = 0; i < m_nSignalPoints; i++)
604 m_adComplexFftInput[i].re = input[i];
606 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
607 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
608 for (i = 0; i < m_nFilterPoints; i++) {
609 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
610 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
613 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
614 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
615 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
616 output[i] = ifftOutput[i].re;
\r
617 delete [] ifftOutput;
625 * convolve Discrete convolution of two functions
628 * r = convolve (f1, f2, dx, n, np, func_type)
629 * double r Convolved result
630 * double f1[], f2[] Functions to be convolved
631 * double dx Difference between successive x values
632 * int n Array index to center convolution about
633 * int np Number of points in f1 array
634 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
637 * f1 is the projection data, its indices range from 0 to np - 1.
638 * The index for f2, the filter, ranges from -(np-1) to (np-1).
639 * There are 3 ways to handle the negative vertices of f2:
640 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
641 * All filters used in reconstruction are even.
642 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
643 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
644 * for negative indices. Since f2 must range from -(np-1) to (np-1),
645 * if we add (np - 1) to f2's array index, then f2's index will
646 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
647 * stored at f2[np-1].
651 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
655 #if UNOPTIMIZED_CONVOLUTION
656 for (int i = 0; i < np; i++)
657 sum += func[i] * m_adFilter[n - i + (np - 1)];
659 double* f2 = m_adFilter + n + (np - 1);
660 for (int i = 0; i < np; i++)
661 sum += *func++ * *f2--;
669 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
673 #if UNOPTIMIZED_CONVOLUTION
674 for (int i = 0; i < np; i++)
675 sum += func[i] * m_adFilter[n - i + (np - 1)];
677 double* f2 = m_adFilter + n + (np - 1);
678 for (int i = 0; i < np; i++)
679 sum += *func++ * *f2--;
687 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
689 std::complex<double>* complexOutput = new std::complex<double> [n];
691 finiteFourierTransform (input, complexOutput, n, direction);
692 for (int i = 0; i < n; i++)
693 output[i] = complexOutput[i].real();
\r
694 delete [] complexOutput;
698 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
705 double angleIncrement = direction * 2 * PI / n;
706 for (int i = 0; i < n; i++) {
709 for (int j = 0; j < n; j++) {
710 double angle = i * j * angleIncrement;
711 sumReal += input[j] * cos(angle);
712 sumImag += input[j] * sin(angle);
718 output[i] = std::complex<double> (sumReal, sumImag);
724 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
731 double angleIncrement = direction * 2 * PI / n;
732 for (int i = 0; i < n; i++) {
733 std::complex<double> sum (0,0);
734 for (int j = 0; j < n; j++) {
735 double angle = i * j * angleIncrement;
736 std::complex<double> exponentTerm (cos(angle), sin(angle));
\r
737 sum += input[j] * exponentTerm;
\r
747 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
754 double angleIncrement = direction * 2 * PI / n;
755 for (int i = 0; i < n; i++) {
757 for (int j = 0; j < n; j++) {
758 double angle = i * j * angleIncrement;
759 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
768 // Table-based routines
771 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
778 for (int i = 0; i < m_nFilterPoints; i++) {
779 double sumReal = 0, sumImag = 0;
780 for (int j = 0; j < m_nFilterPoints; j++) {
781 int tableIndex = i * j;
783 sumReal += input[j] * m_adFourierCosTable[tableIndex];
784 sumImag += input[j] * m_adFourierSinTable[tableIndex];
786 sumReal += input[j] * m_adFourierCosTable[tableIndex];
787 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
791 sumReal /= m_nFilterPoints;
792 sumImag /= m_nFilterPoints;
794 output[i] = std::complex<double> (sumReal, sumImag);
798 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
800 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
807 for (int i = 0; i < m_nFilterPoints; i++) {
808 double sumReal = 0, sumImag = 0;
809 for (int j = 0; j < m_nFilterPoints; j++) {
810 int tableIndex = i * j;
812 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
813 - input[j].imag() * m_adFourierSinTable[tableIndex];
814 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
815 + input[j].imag() * m_adFourierCosTable[tableIndex];
817 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
818 - input[j].imag() * -m_adFourierSinTable[tableIndex];
819 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
820 + input[j].imag() * m_adFourierCosTable[tableIndex];
824 sumReal /= m_nFilterPoints;
825 sumImag /= m_nFilterPoints;
827 output[i] = std::complex<double> (sumReal, sumImag);
832 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
839 for (int i = 0; i < m_nFilterPoints; i++) {
841 for (int j = 0; j < m_nFilterPoints; j++) {
842 int tableIndex = i * j;
844 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
845 - input[j].imag() * m_adFourierSinTable[tableIndex];
847 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
848 - input[j].imag() * -m_adFourierSinTable[tableIndex];
852 sumReal /= m_nFilterPoints;