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.21 2001/01/13 02:55:14 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 ******************************************************************************/
31 #include "../src/dlgezplot.h"
34 // FilterMethod ID/Names
35 const int ProcessSignal::FILTER_METHOD_INVALID = -1;
36 const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
37 const int ProcessSignal::FILTER_METHOD_FOURIER = 1;
38 const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2;
39 const int ProcessSignal::FILTER_METHOD_FFT = 3;
41 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
42 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
44 const char* ProcessSignal::s_aszFilterMethodName[] = {
54 const char* ProcessSignal::s_aszFilterMethodTitle[] = {
57 {"Fouier Trigometric Table"},
61 {"Real/Half-Complex FFTW"},
64 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
66 // FilterGeneration ID/Names
67 const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
68 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
69 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
70 const char* ProcessSignal::s_aszFilterGenerationName[] = {
74 const char* ProcessSignal::s_aszFilterGenerationTitle[] = {
78 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
81 // CLASS IDENTIFICATION
84 ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth,
85 double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName,
86 const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel,
87 int iGeometry, double dFocalLength, SGP* pSGP)
88 : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
90 m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
91 if (m_idFilterMethod == FILTER_METHOD_INVALID) {
93 m_failMessage = "Invalid filter method name ";
94 m_failMessage += szFilterMethodName;
97 m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
98 if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
100 m_failMessage = "Invalid frequency filter name ";
101 m_failMessage += szFilterGenerationName;
104 m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
105 if (m_idFilter == SignalFilter::FILTER_INVALID) {
107 m_failMessage = "Invalid Filter name ";
108 m_failMessage += szFilterName;
111 m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
112 if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
114 m_failMessage = "Invalid domain name ";
115 m_failMessage += szDomainName;
119 init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
120 m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
125 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement,
126 int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration,
127 const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry,
128 double dFocalLength, SGP* pSGP)
131 m_idFilter = idFilter;
132 m_idDomain = idDomain;
133 m_idFilterMethod = idFilterMethod;
134 m_idFilterGeneration = idFilterGeneration;
135 m_idGeometry = iGeometry;
136 m_dFocalLength = dFocalLength;
138 if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
142 m_traceLevel = iTraceLevel;
143 m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
144 m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
145 m_dBandwidth = dBandwidth;
146 m_nSignalPoints = nSignalPoints;
147 m_dSignalInc = dSignalIncrement;
148 m_dFilterParam = dFilterParam;
149 m_iZeropad = iZeropad;
150 m_iPreinterpolationFactor = iPreinterpolationFactor;
152 // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
153 // through origin of phantom rather than 2 times distance to detector,
154 // see Kak-Slaney Fig 3.22, for Collinear diagram
155 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
160 if (m_idFilterMethod == FILTER_METHOD_FFT) {
162 m_idFilterMethod = FILTER_METHOD_RFFTW;
165 m_failMessage = "FFT not yet implemented";
170 bool m_bFrequencyFiltering = true;
171 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
172 m_bFrequencyFiltering = false;
174 // Spatial-based filtering
175 if (! m_bFrequencyFiltering) {
177 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
178 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
179 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
180 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
181 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
182 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
183 m_adFilter = new double[ m_nFilterPoints ];
184 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
185 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
186 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
187 m_dFilterMin = -1. / (2 * m_dSignalInc);
188 m_dFilterMax = 1. / (2 * m_dSignalInc);
189 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
190 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
191 m_adFilter = new double[ m_nFilterPoints ];
192 double* adFrequencyFilter = new double [m_nFilterPoints];
193 filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
194 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
195 EZPlotDialog* pEZPlotDlg = NULL;
196 if (g_bRunningWXWindows && m_traceLevel > 0) {
197 EZPlotDialog dlgEZPlot;
198 dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Natural Order");
199 dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints);
200 dlgEZPlot.ShowModal();
203 Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
205 if (g_bRunningWXWindows && m_traceLevel > 0) {
206 EZPlotDialog dlgEZPlot;
207 dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Fourier Order");
208 dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints);
209 dlgEZPlot.ShowModal();
212 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
213 delete adFrequencyFilter;
214 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
215 if (g_bRunningWXWindows && m_traceLevel > 0) {
216 EZPlotDialog dlgEZPlot;
217 dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Fourier Order");
218 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
219 dlgEZPlot.ShowModal();
222 Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
223 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
224 if (g_bRunningWXWindows && m_traceLevel > 0) {
225 EZPlotDialog dlgEZPlot;
226 dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Natural Order");
227 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
228 dlgEZPlot.ShowModal();
231 for (i = 0; i < m_nFilterPoints; i++) {
232 m_adFilter[i] /= m_dSignalInc;
235 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
236 for (i = 0; i < m_nFilterPoints; i++)
237 m_adFilter[i] *= 0.5;
238 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
239 for (i = 0; i < m_nFilterPoints; i++) {
240 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
241 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
242 double dScale = 0.5 * sinScale * sinScale;
243 m_adFilter[i] *= dScale;
245 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
246 if (g_bRunningWXWindows && m_traceLevel > 0) {
247 EZPlotDialog dlgEZPlot;
248 dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Fourier Frequency: Natural Order");
249 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
250 dlgEZPlot.ShowModal();
254 } // if (spatial filtering)
256 else if (m_bFrequencyFiltering) { // Frequency-based filtering
258 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
259 // calculate number of filter points with zeropadding
260 m_nFilterPoints = m_nSignalPoints;
261 if (m_iZeropad > 0) {
262 double logBase2 = log(m_nFilterPoints) / log(2);
263 int nextPowerOf2 = static_cast<int>(floor(logBase2));
264 if (logBase2 != floor(logBase2))
266 nextPowerOf2 += (m_iZeropad - 1);
267 m_nFilterPoints = 1 << nextPowerOf2;
268 #if defined(DEBUG) || defined(_DEBUG)
269 if (m_traceLevel >= Trace::TRACE_CONSOLE)
270 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
273 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
275 if (m_nFilterPoints % 2) { // Odd
276 m_dFilterMin = -1. / (2 * m_dSignalInc);
277 m_dFilterMax = 1. / (2 * m_dSignalInc);
278 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
280 m_dFilterMin = -1. / (2 * m_dSignalInc);
281 m_dFilterMax = 1. / (2 * m_dSignalInc);
282 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
283 m_dFilterMax -= m_dFilterInc;
286 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
287 m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
288 m_adFilter = new double [m_nFilterPoints];
289 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
291 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
292 if (g_bRunningWXWindows && m_traceLevel > 0) {
293 EZPlotDialog dlgEZPlot;
294 dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
295 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
296 dlgEZPlot.ShowModal();
300 // This works fairly well. I'm not sure why since scaling for geometries is done on
301 // frequency filter rather than spatial filter as it should be.
302 // It gives values slightly off than freq/inverse filtering
303 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
304 for (i = 0; i < m_nFilterPoints; i++)
305 m_adFilter[i] *= 0.5;
306 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
307 for (i = 0; i < m_nFilterPoints; i++) {
308 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
309 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
310 double dScale = 0.5 * sinScale * sinScale;
311 m_adFilter[i] *= dScale;
314 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
315 if (g_bRunningWXWindows && m_traceLevel > 0) {
316 EZPlotDialog dlgEZPlot;
317 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
318 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
319 dlgEZPlot.ShowModal();
322 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
323 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
324 if (g_bRunningWXWindows && m_traceLevel > 0) {
325 EZPlotDialog dlgEZPlot;
326 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
327 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
328 dlgEZPlot.ShowModal();
332 // FILTERING: FREQUENCY - INVERSE FOURIER
334 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
335 // calculate number of filter points with zeropadding
336 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
337 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
338 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
339 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
340 m_nFilterPoints = nSpatialPoints;
341 if (m_iZeropad > 0) {
342 double logBase2 = log(nSpatialPoints) / log(2);
343 int nextPowerOf2 = static_cast<int>(floor(logBase2));
344 if (logBase2 != floor(logBase2))
346 nextPowerOf2 += (m_iZeropad - 1);
347 m_nFilterPoints = 1 << nextPowerOf2;
349 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
350 #if defined(DEBUG) || defined(_DEBUG)
351 if (m_traceLevel >= Trace::TRACE_CONSOLE)
352 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
354 double* adSpatialFilter = new double [m_nFilterPoints];
355 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
356 m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
357 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
358 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
359 if (g_bRunningWXWindows && m_traceLevel > 0) {
360 EZPlotDialog dlgEZPlot;;
361 dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
362 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
363 dlgEZPlot.ShowModal();
367 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
368 for (i = 0; i < nSpatialPoints; i++)
369 adSpatialFilter[i] *= 0.5;
370 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
371 for (i = 0; i < nSpatialPoints; i++) {
372 int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
373 double sinScale = sin (iDetFromZero * m_dSignalInc);
374 if (fabs(sinScale) < 1E-7)
377 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
378 double dScale = 0.5 * sinScale * sinScale;
379 adSpatialFilter[i] *= dScale;
382 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
383 if (g_bRunningWXWindows && m_traceLevel > 0) {
384 EZPlotDialog dlgEZPlot;;
385 dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
386 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
387 dlgEZPlot.ShowModal();
390 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
391 adSpatialFilter[i] = 0;
393 m_adFilter = new double [m_nFilterPoints];
394 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
395 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
396 delete adSpatialFilter;
397 for (i = 0; i < m_nFilterPoints; i++)
398 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
399 delete acInverseFilter;
400 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
401 if (g_bRunningWXWindows && m_traceLevel > 0) {
402 EZPlotDialog dlgEZPlot;
403 dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
404 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
405 dlgEZPlot.ShowModal();
411 // precalculate sin and cosine tables for fourier transform
412 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
413 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
414 double angleIncrement = (2. * PI) / m_nFilterPoints;
415 m_adFourierCosTable = new double[ nFourier ];
416 m_adFourierSinTable = new double[ nFourier ];
418 for (i = 0; i < nFourier; i++) {
419 m_adFourierCosTable[i] = cos (angle);
420 m_adFourierSinTable[i] = sin (angle);
421 angle += angleIncrement;
426 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
427 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
428 m_adFilter[i] /= m_nFilterPoints;
431 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
432 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
433 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
434 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
435 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
436 for (i = 0; i < m_nFilterPoints; i++)
437 m_adRealFftInput[i] = 0;
438 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
439 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
440 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
441 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
442 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
443 for (i = 0; i < m_nFilterPoints; i++)
444 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
445 for (i = 0; i < m_nOutputPoints; i++)
446 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
452 ProcessSignal::~ProcessSignal (void)
454 delete [] m_adFourierSinTable;
455 delete [] m_adFourierCosTable;
456 delete [] m_adFilter;
459 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
460 fftw_destroy_plan(m_complexPlanForward);
461 fftw_destroy_plan(m_complexPlanBackward);
462 delete [] m_adComplexFftInput;
463 delete [] m_adComplexFftSignal;
465 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
466 rfftw_destroy_plan(m_realPlanForward);
467 rfftw_destroy_plan(m_realPlanBackward);
468 delete [] m_adRealFftInput;
469 delete [] m_adRealFftSignal;
475 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
477 int fmID = FILTER_METHOD_INVALID;
479 for (int i = 0; i < s_iFilterMethodCount; i++)
480 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
489 ProcessSignal::convertFilterMethodIDToName (const int fmID)
491 static const char *name = "";
493 if (fmID >= 0 && fmID < s_iFilterMethodCount)
494 return (s_aszFilterMethodName [fmID]);
500 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
502 static const char *title = "";
504 if (fmID >= 0 && fmID < s_iFilterMethodCount)
505 return (s_aszFilterMethodTitle [fmID]);
512 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
514 int fgID = FILTER_GENERATION_INVALID;
516 for (int i = 0; i < s_iFilterGenerationCount; i++)
517 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
526 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
528 static const char *name = "";
530 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
531 return (s_aszFilterGenerationName [fgID]);
537 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
539 static const char *name = "";
541 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
542 return (s_aszFilterGenerationTitle [fgID]);
548 ProcessSignal::filterSignal (const float constInput[], double output[]) const
550 double* input = new double [m_nSignalPoints];
552 for (i = 0; i < m_nSignalPoints; i++)
553 input[i] = constInput[i];
555 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
556 for (int i = 0; i < m_nSignalPoints; i++) {
557 int iDetFromCenter = i - (m_nSignalPoints / 2);
558 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
560 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
561 for (int i = 0; i < m_nSignalPoints; i++) {
562 int iDetFromCenter = i - (m_nSignalPoints / 2);
563 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
566 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
567 for (i = 0; i < m_nSignalPoints; i++)
568 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
569 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
570 double* inputSignal = new double [m_nFilterPoints];
571 for (i = 0; i < m_nSignalPoints; i++)
572 inputSignal[i] = input[i];
573 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
574 inputSignal[i] = 0; // zeropad
575 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
576 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
578 for (i = 0; i < m_nFilterPoints; i++)
579 fftSignal[i] *= m_adFilter[i];
580 double* inverseFourier = new double [m_nFilterPoints];
581 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
583 for (i = 0; i < m_nSignalPoints; i++)
584 output[i] = inverseFourier[i];
585 delete inverseFourier;
586 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
587 double* inputSignal = new double [m_nFilterPoints];
588 for (i = 0; i < m_nSignalPoints; i++)
589 inputSignal[i] = input[i];
590 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
591 inputSignal[i] = 0; // zeropad
592 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
593 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
595 for (i = 0; i < m_nFilterPoints; i++)
596 fftSignal[i] *= m_adFilter[i];
597 double* inverseFourier = new double [m_nFilterPoints];
598 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
600 for (i = 0; i < m_nSignalPoints; i++)
601 output[i] = inverseFourier[i];
602 delete inverseFourier;
605 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
606 for (i = 0; i < m_nSignalPoints; i++)
607 m_adRealFftInput[i] = input[i];
609 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
610 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
611 for (i = 0; i < m_nFilterPoints; i++)
612 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
614 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
615 m_adRealFftSignal[i] = 0;
617 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
618 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
619 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
620 output[i] = ifftOutput[i];
621 delete [] ifftOutput;
622 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
623 for (i = 0; i < m_nSignalPoints; i++)
624 m_adComplexFftInput[i].re = input[i];
626 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
627 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
628 for (i = 0; i < m_nFilterPoints; i++) {
629 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
630 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
633 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
634 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
635 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
636 output[i] = ifftOutput[i].re;
637 delete [] ifftOutput;
645 * convolve Discrete convolution of two functions
648 * r = convolve (f1, f2, dx, n, np, func_type)
649 * double r Convolved result
650 * double f1[], f2[] Functions to be convolved
651 * double dx Difference between successive x values
652 * int n Array index to center convolution about
653 * int np Number of points in f1 array
654 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
657 * f1 is the projection data, its indices range from 0 to np - 1.
658 * The index for f2, the filter, ranges from -(np-1) to (np-1).
659 * There are 3 ways to handle the negative vertices of f2:
660 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
661 * All filters used in reconstruction are even.
662 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
663 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
664 * for negative indices. Since f2 must range from -(np-1) to (np-1),
665 * if we add (np - 1) to f2's array index, then f2's index will
666 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
667 * stored at f2[np-1].
671 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
675 #if UNOPTIMIZED_CONVOLUTION
676 for (int i = 0; i < np; i++)
677 sum += func[i] * m_adFilter[n - i + (np - 1)];
679 double* f2 = m_adFilter + n + (np - 1);
680 for (int i = 0; i < np; i++)
681 sum += *func++ * *f2--;
689 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
693 #if UNOPTIMIZED_CONVOLUTION
694 for (int i = 0; i < np; i++)
695 sum += func[i] * m_adFilter[n - i + (np - 1)];
697 double* f2 = m_adFilter + n + (np - 1);
698 for (int i = 0; i < np; i++)
699 sum += *func++ * *f2--;
707 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
709 std::complex<double>* complexOutput = new std::complex<double> [n];
711 finiteFourierTransform (input, complexOutput, n, direction);
712 for (int i = 0; i < n; i++)
713 output[i] = complexOutput[i].real();
714 delete [] complexOutput;
718 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
725 double angleIncrement = direction * 2 * PI / n;
726 for (int i = 0; i < n; i++) {
729 for (int j = 0; j < n; j++) {
730 double angle = i * j * angleIncrement;
731 sumReal += input[j] * cos(angle);
732 sumImag += input[j] * sin(angle);
738 output[i] = std::complex<double> (sumReal, sumImag);
744 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
751 double angleIncrement = direction * 2 * PI / n;
752 for (int i = 0; i < n; i++) {
753 std::complex<double> sum (0,0);
754 for (int j = 0; j < n; j++) {
755 double angle = i * j * angleIncrement;
756 std::complex<double> exponentTerm (cos(angle), sin(angle));
757 sum += input[j] * exponentTerm;
767 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
774 double angleIncrement = direction * 2 * PI / n;
775 for (int i = 0; i < n; i++) {
777 for (int j = 0; j < n; j++) {
778 double angle = i * j * angleIncrement;
779 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
788 // Table-based routines
791 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
798 for (int i = 0; i < m_nFilterPoints; i++) {
799 double sumReal = 0, sumImag = 0;
800 for (int j = 0; j < m_nFilterPoints; j++) {
801 int tableIndex = i * j;
803 sumReal += input[j] * m_adFourierCosTable[tableIndex];
804 sumImag += input[j] * m_adFourierSinTable[tableIndex];
806 sumReal += input[j] * m_adFourierCosTable[tableIndex];
807 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
811 sumReal /= m_nFilterPoints;
812 sumImag /= m_nFilterPoints;
814 output[i] = std::complex<double> (sumReal, sumImag);
818 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
820 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
827 for (int i = 0; i < m_nFilterPoints; i++) {
828 double sumReal = 0, sumImag = 0;
829 for (int j = 0; j < m_nFilterPoints; j++) {
830 int tableIndex = i * j;
832 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
833 - input[j].imag() * m_adFourierSinTable[tableIndex];
834 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
835 + input[j].imag() * m_adFourierCosTable[tableIndex];
837 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
838 - input[j].imag() * -m_adFourierSinTable[tableIndex];
839 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
840 + input[j].imag() * m_adFourierCosTable[tableIndex];
844 sumReal /= m_nFilterPoints;
845 sumImag /= m_nFilterPoints;
847 output[i] = std::complex<double> (sumReal, sumImag);
852 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
859 for (int i = 0; i < m_nFilterPoints; i++) {
861 for (int j = 0; j < m_nFilterPoints; j++) {
862 int tableIndex = i * j;
864 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
865 - input[j].imag() * m_adFourierSinTable[tableIndex];
867 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
868 - input[j].imag() * -m_adFourierSinTable[tableIndex];
872 sumReal /= m_nFilterPoints;