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.17 2001/01/12 14:14:58 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,
81 double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName,
82 const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel,
83 int iGeometry, double dFocalLength, SGP* pSGP)
84 : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
86 m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
87 if (m_idFilterMethod == FILTER_METHOD_INVALID) {
89 m_failMessage = "Invalid filter method name ";
90 m_failMessage += szFilterMethodName;
93 m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
94 if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
96 m_failMessage = "Invalid frequency filter name ";
97 m_failMessage += szFilterGenerationName;
100 m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
101 if (m_idFilter == SignalFilter::FILTER_INVALID) {
103 m_failMessage = "Invalid Filter name ";
104 m_failMessage += szFilterName;
107 m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
108 if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
110 m_failMessage = "Invalid domain name ";
111 m_failMessage += szDomainName;
115 init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
116 m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
121 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement,
122 int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration,
123 const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry,
124 double dFocalLength, SGP* pSGP)
127 m_idFilter = idFilter;
128 m_idDomain = idDomain;
129 m_idFilterMethod = idFilterMethod;
130 m_idFilterGeneration = idFilterGeneration;
131 m_idGeometry = iGeometry;
132 m_dFocalLength = dFocalLength;
134 if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
138 m_traceLevel = iTraceLevel;
139 m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
140 m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
141 m_dBandwidth = dBandwidth;
142 m_nSignalPoints = nSignalPoints;
143 m_dSignalInc = dSignalIncrement;
144 m_dFilterParam = dFilterParam;
145 m_iZeropad = iZeropad;
146 m_iPreinterpolationFactor = iPreinterpolationFactor;
148 // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
149 // through origin of phantom rather than 2 times distance to detector,
150 // see Kak-Slaney Fig 3.22, for Collinear diagram
151 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
156 if (m_idFilterMethod == FILTER_METHOD_FFT) {
158 m_idFilterMethod = FILTER_METHOD_RFFTW;
161 m_failMessage = "FFT not yet implemented";
166 bool m_bFrequencyFiltering = true;
167 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
168 m_bFrequencyFiltering = false;
170 // Spatial-based filtering
171 if (! m_bFrequencyFiltering) {
173 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
174 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
175 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
176 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
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_SPATIAL);
179 m_adFilter = new double[ m_nFilterPoints ];
180 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
181 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
182 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
183 m_dFilterMin = -1. / (2 * m_dSignalInc);
184 m_dFilterMax = 1. / (2 * m_dSignalInc);
185 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
186 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
187 m_adFilter = new double[ m_nFilterPoints ];
188 double* adFrequencyFilter = new double [m_nFilterPoints];
189 filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
191 EZPlot* pEZPlot = NULL;
192 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
193 pEZPlot = new EZPlot ();
194 pEZPlot->ezset ("title Filter Response: Natural Order");
195 pEZPlot->ezset ("ylength 0.25");
196 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
197 pEZPlot->plot (pSGP);
200 Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
202 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
203 pEZPlot->ezset ("title Filter Response: Fourier Order");
204 pEZPlot->ezset ("ylength 0.25");
205 pEZPlot->ezset ("yporigin 0.25");
206 pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
207 pEZPlot->plot (pSGP);
210 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
211 delete adFrequencyFilter;
213 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
214 pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
215 pEZPlot->ezset ("ylength 0.25");
216 pEZPlot->ezset ("yporigin 0.50");
217 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
218 pEZPlot->plot (pSGP);
221 Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
223 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
224 pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
225 pEZPlot->ezset ("ylength 0.25");
226 pEZPlot->ezset ("yporigin 0.75");
227 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
228 pEZPlot->plot (pSGP);
232 for (i = 0; i < m_nFilterPoints; i++) {
233 m_adFilter[i] /= m_dSignalInc;
236 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
237 for (i = 0; i < m_nFilterPoints; i++)
238 m_adFilter[i] *= 0.5;
239 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
240 for (i = 0; i < m_nFilterPoints; i++) {
241 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
242 double sinScale = sin (iDetFromZero * m_dSignalInc);
243 if (fabs(sinScale) < 1E-7)
246 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
247 double dScale = 0.5 * sinScale * sinScale;
248 m_adFilter[i] *= dScale;
251 } // if (spatial filtering)
253 else if (m_bFrequencyFiltering) { // Frequency-based filtering
255 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
256 // calculate number of filter points with zeropadding
257 m_nFilterPoints = m_nSignalPoints;
258 if (m_iZeropad > 0) {
259 double logBase2 = log(m_nFilterPoints) / log(2);
260 int nextPowerOf2 = static_cast<int>(floor(logBase2));
261 if (logBase2 != floor(logBase2))
263 nextPowerOf2 += (m_iZeropad - 1);
264 m_nFilterPoints = 1 << nextPowerOf2;
266 if (m_traceLevel >= Trace::TRACE_CONSOLE)
267 std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
270 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
272 if (m_nFilterPoints % 2) { // Odd
273 m_dFilterMin = -1. / (2 * m_dSignalInc);
274 m_dFilterMax = 1. / (2 * m_dSignalInc);
275 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
277 m_dFilterMin = -1. / (2 * m_dSignalInc);
278 m_dFilterMax = 1. / (2 * m_dSignalInc);
279 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
280 m_dFilterMax -= m_dFilterInc;
283 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
284 m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
285 m_adFilter = new double [m_nFilterPoints];
286 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
288 // This doesn't work!
289 // Need to add filtering for divergent geometries & Frequency/Direct filtering
290 // Jan 2001: Direct seems to work for equilinear and equiangular
291 // however, inverse_fourier doesn't work for equiangular on all versions of CTSim tested
293 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
294 for (i = 0; i < m_nFilterPoints; i++)
295 m_adFilter[i] *= 0.5;
296 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
297 for (i = 0; i < m_nFilterPoints; i++) {
298 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
299 double sinScale = sin (iDetFromZero * m_dSignalInc);
300 if (fabs(sinScale) < 1E-7)
303 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
304 double dScale = 0.5 * sinScale * sinScale;
305 m_adFilter[i] *= dScale;
309 EZPlot* pEZPlot = NULL;
310 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
311 pEZPlot = new EZPlot;
312 pEZPlot->ezset ("title Filter Filter: Natural Order");
313 pEZPlot->ezset ("ylength 0.50");
314 pEZPlot->ezset ("yporigin 0.00");
315 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
316 pEZPlot->plot (pSGP);
319 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
321 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
322 pEZPlot->ezset ("title Filter Filter: Fourier Order");
323 pEZPlot->ezset ("ylength 0.50");
324 pEZPlot->ezset ("yporigin 0.50");
325 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
326 pEZPlot->plot (pSGP);
330 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
331 // calculate number of filter points with zeropadding
332 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
333 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
334 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
335 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
336 m_nFilterPoints = nSpatialPoints;
337 if (m_iZeropad > 0) {
338 double logBase2 = log(nSpatialPoints) / log(2);
339 int nextPowerOf2 = static_cast<int>(floor(logBase2));
340 if (logBase2 != floor(logBase2))
342 nextPowerOf2 += (m_iZeropad - 1);
343 m_nFilterPoints = 1 << nextPowerOf2;
345 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
347 if (m_traceLevel >= Trace::TRACE_CONSOLE)
348 std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
350 double* adSpatialFilter = new double [m_nFilterPoints];
351 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
352 m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
353 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
355 EZPlot* pEZPlot = NULL;
356 if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
357 pEZPlot = new EZPlot;
358 pEZPlot->ezset ("title Spatial Filter: Natural Order");
359 pEZPlot->ezset ("ylength 0.50");
360 pEZPlot->ezset ("yporigin 0.00");
361 pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
362 pEZPlot->plot (pSGP);
367 // #define PRE_JAN_2001 1
369 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
370 for (i = 0; i < m_nFilterPoints; i++)
371 adSpatialFilter[i] *= 0.5;
372 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
373 for (i = 0; i < m_nFilterPoints; i++) {
374 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
375 double sinScale = sin (iDetFromZero * m_dSignalInc);
376 if (fabs(sinScale) < 1E-7)
379 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
380 double dScale = 0.5 * sinScale * sinScale;
381 adSpatialFilter[i] *= dScale;
384 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
385 adSpatialFilter[i] = 0;
387 m_adFilter = new double [m_nFilterPoints];
388 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
389 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
390 delete adSpatialFilter;
391 for (i = 0; i < m_nFilterPoints; i++)
392 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
393 delete acInverseFilter;
395 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
396 adSpatialFilter[i] = 0;
398 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
399 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, FORWARD);
400 delete adSpatialFilter;
401 m_adFilter = new double [m_nFilterPoints];
402 for (i = 0; i < m_nFilterPoints; i++)
403 m_adFilter[i] = std::abs(acInverseFilter[i]);
404 delete acInverseFilter;
406 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
407 for (i = 0; i < m_nFilterPoints; i++)
408 m_adFilter[i] *= 0.5;
409 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
410 for (i = 0; i < m_nFilterPoints; i++) {
411 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
412 double sinScale = sin (iDetFromZero * m_dSignalInc);
413 if (fabs(sinScale) < 1E-7)
416 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
417 double dScale = 0.5 * sinScale * sinScale;
418 m_adFilter[i] *= dScale;
424 if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
425 pEZPlot->ezset ("title Spatial Filter: Inverse");
426 pEZPlot->ezset ("ylength 0.50");
427 pEZPlot->ezset ("yporigin 0.50");
428 pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
429 pEZPlot->plot (pSGP);
436 // precalculate sin and cosine tables for fourier transform
437 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
438 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
439 double angleIncrement = (2. * PI) / m_nFilterPoints;
440 m_adFourierCosTable = new double[ nFourier ];
441 m_adFourierSinTable = new double[ nFourier ];
443 for (i = 0; i < nFourier; i++) {
444 m_adFourierCosTable[i] = cos (angle);
445 m_adFourierSinTable[i] = sin (angle);
446 angle += angleIncrement;
451 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
452 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
453 m_adFilter[i] /= m_nFilterPoints;
456 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
457 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
458 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
459 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
460 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
461 for (i = 0; i < m_nFilterPoints; i++)
462 m_adRealFftInput[i] = 0;
463 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
464 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
465 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
466 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
467 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
468 for (i = 0; i < m_nFilterPoints; i++)
469 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
470 for (i = 0; i < m_nOutputPoints; i++)
471 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
477 ProcessSignal::~ProcessSignal (void)
479 delete [] m_adFourierSinTable;
480 delete [] m_adFourierCosTable;
481 delete [] m_adFilter;
484 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
485 fftw_destroy_plan(m_complexPlanForward);
486 fftw_destroy_plan(m_complexPlanBackward);
487 delete [] m_adComplexFftInput;
488 delete [] m_adComplexFftSignal;
490 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
491 rfftw_destroy_plan(m_realPlanForward);
492 rfftw_destroy_plan(m_realPlanBackward);
493 delete [] m_adRealFftInput;
494 delete [] m_adRealFftSignal;
500 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
502 int fmID = FILTER_METHOD_INVALID;
504 for (int i = 0; i < s_iFilterMethodCount; i++)
505 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
514 ProcessSignal::convertFilterMethodIDToName (const int fmID)
516 static const char *name = "";
518 if (fmID >= 0 && fmID < s_iFilterMethodCount)
519 return (s_aszFilterMethodName [fmID]);
525 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
527 static const char *title = "";
529 if (fmID >= 0 && fmID < s_iFilterMethodCount)
530 return (s_aszFilterMethodTitle [fmID]);
537 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
539 int fgID = FILTER_GENERATION_INVALID;
541 for (int i = 0; i < s_iFilterGenerationCount; i++)
542 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
551 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
553 static const char *name = "";
555 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
556 return (s_aszFilterGenerationName [fgID]);
562 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
564 static const char *name = "";
566 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
567 return (s_aszFilterGenerationTitle [fgID]);
573 ProcessSignal::filterSignal (const float constInput[], double output[]) const
575 double* input = new double [m_nSignalPoints];
577 for (i = 0; i < m_nSignalPoints; i++)
578 input[i] = constInput[i];
580 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
581 for (int i = 0; i < m_nSignalPoints; i++) {
582 int iDetFromCenter = i - (m_nSignalPoints / 2);
583 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
585 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
586 for (int i = 0; i < m_nSignalPoints; i++) {
587 int iDetFromCenter = i - (m_nSignalPoints / 2);
588 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
591 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
592 for (i = 0; i < m_nSignalPoints; i++)
593 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
594 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
595 double* inputSignal = new double [m_nFilterPoints];
596 for (i = 0; i < m_nSignalPoints; i++)
597 inputSignal[i] = input[i];
598 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
599 inputSignal[i] = 0; // zeropad
600 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
601 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
603 for (i = 0; i < m_nFilterPoints; i++)
604 fftSignal[i] *= m_adFilter[i];
605 double* inverseFourier = new double [m_nFilterPoints];
606 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
608 for (i = 0; i < m_nSignalPoints; i++)
609 output[i] = inverseFourier[i];
610 delete inverseFourier;
611 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
612 double* inputSignal = new double [m_nFilterPoints];
613 for (i = 0; i < m_nSignalPoints; i++)
614 inputSignal[i] = input[i];
615 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
616 inputSignal[i] = 0; // zeropad
617 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
618 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
620 for (i = 0; i < m_nFilterPoints; i++)
621 fftSignal[i] *= m_adFilter[i];
622 double* inverseFourier = new double [m_nFilterPoints];
623 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
625 for (i = 0; i < m_nSignalPoints; i++)
626 output[i] = inverseFourier[i];
627 delete inverseFourier;
630 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
631 for (i = 0; i < m_nSignalPoints; i++)
632 m_adRealFftInput[i] = input[i];
634 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
635 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
636 for (i = 0; i < m_nFilterPoints; i++)
637 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
639 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
640 m_adRealFftSignal[i] = 0;
642 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
643 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
644 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
645 output[i] = ifftOutput[i];
646 delete [] ifftOutput;
647 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
648 for (i = 0; i < m_nSignalPoints; i++)
649 m_adComplexFftInput[i].re = input[i];
651 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
652 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
653 for (i = 0; i < m_nFilterPoints; i++) {
654 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
655 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
658 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
659 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
660 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
661 output[i] = ifftOutput[i].re;
662 delete [] ifftOutput;
670 * convolve Discrete convolution of two functions
673 * r = convolve (f1, f2, dx, n, np, func_type)
674 * double r Convolved result
675 * double f1[], f2[] Functions to be convolved
676 * double dx Difference between successive x values
677 * int n Array index to center convolution about
678 * int np Number of points in f1 array
679 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
682 * f1 is the projection data, its indices range from 0 to np - 1.
683 * The index for f2, the filter, ranges from -(np-1) to (np-1).
684 * There are 3 ways to handle the negative vertices of f2:
685 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
686 * All filters used in reconstruction are even.
687 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
688 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
689 * for negative indices. Since f2 must range from -(np-1) to (np-1),
690 * if we add (np - 1) to f2's array index, then f2's index will
691 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
692 * stored at f2[np-1].
696 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
700 #if UNOPTIMIZED_CONVOLUTION
701 for (int i = 0; i < np; i++)
702 sum += func[i] * m_adFilter[n - i + (np - 1)];
704 double* f2 = m_adFilter + n + (np - 1);
705 for (int i = 0; i < np; i++)
706 sum += *func++ * *f2--;
714 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
718 #if UNOPTIMIZED_CONVOLUTION
719 for (int i = 0; i < np; i++)
720 sum += func[i] * m_adFilter[n - i + (np - 1)];
722 double* f2 = m_adFilter + n + (np - 1);
723 for (int i = 0; i < np; i++)
724 sum += *func++ * *f2--;
732 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
734 std::complex<double>* complexOutput = new std::complex<double> [n];
736 finiteFourierTransform (input, complexOutput, n, direction);
737 for (int i = 0; i < n; i++)
738 output[i] = complexOutput[i].real();
739 delete [] complexOutput;
743 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
750 double angleIncrement = direction * 2 * PI / n;
751 for (int i = 0; i < n; i++) {
754 for (int j = 0; j < n; j++) {
755 double angle = i * j * angleIncrement;
756 sumReal += input[j] * cos(angle);
757 sumImag += input[j] * sin(angle);
763 output[i] = std::complex<double> (sumReal, sumImag);
769 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
776 double angleIncrement = direction * 2 * PI / n;
777 for (int i = 0; i < n; i++) {
778 std::complex<double> sum (0,0);
779 for (int j = 0; j < n; j++) {
780 double angle = i * j * angleIncrement;
781 std::complex<double> exponentTerm (cos(angle), sin(angle));
782 sum += input[j] * exponentTerm;
792 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
799 double angleIncrement = direction * 2 * PI / n;
800 for (int i = 0; i < n; i++) {
802 for (int j = 0; j < n; j++) {
803 double angle = i * j * angleIncrement;
804 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
813 // Table-based routines
816 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
823 for (int i = 0; i < m_nFilterPoints; i++) {
824 double sumReal = 0, sumImag = 0;
825 for (int j = 0; j < m_nFilterPoints; j++) {
826 int tableIndex = i * j;
828 sumReal += input[j] * m_adFourierCosTable[tableIndex];
829 sumImag += input[j] * m_adFourierSinTable[tableIndex];
831 sumReal += input[j] * m_adFourierCosTable[tableIndex];
832 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
836 sumReal /= m_nFilterPoints;
837 sumImag /= m_nFilterPoints;
839 output[i] = std::complex<double> (sumReal, sumImag);
843 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
845 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
852 for (int i = 0; i < m_nFilterPoints; i++) {
853 double sumReal = 0, sumImag = 0;
854 for (int j = 0; j < m_nFilterPoints; j++) {
855 int tableIndex = i * j;
857 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
858 - input[j].imag() * m_adFourierSinTable[tableIndex];
859 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
860 + input[j].imag() * m_adFourierCosTable[tableIndex];
862 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
863 - input[j].imag() * -m_adFourierSinTable[tableIndex];
864 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
865 + input[j].imag() * m_adFourierCosTable[tableIndex];
869 sumReal /= m_nFilterPoints;
870 sumImag /= m_nFilterPoints;
872 output[i] = std::complex<double> (sumReal, sumImag);
877 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
884 for (int i = 0; i < m_nFilterPoints; i++) {
886 for (int j = 0; j < m_nFilterPoints; j++) {
887 int tableIndex = i * j;
889 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
890 - input[j].imag() * m_adFourierSinTable[tableIndex];
892 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
893 - input[j].imag() * -m_adFourierSinTable[tableIndex];
897 sumReal /= m_nFilterPoints;