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-2001 Kevin Rosenberg
12 ** $Id: procsignal.cpp,v 1.27 2001/03/01 07:30:49 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 "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* const ProcessSignal::s_aszFilterMethodName[] = {
54 const char* const 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* const ProcessSignal::s_aszFilterGenerationName[] = {
74 const char* const 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, double dSourceDetectorLength, 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,
121 dSourceDetectorLength, pSGP);
126 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement,
127 int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration,
128 const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry,
129 double dFocalLength, double dSourceDetectorLength, SGP* pSGP)
132 m_idFilter = idFilter;
133 m_idDomain = idDomain;
134 m_idFilterMethod = idFilterMethod;
135 m_idFilterGeneration = idFilterGeneration;
136 m_idGeometry = iGeometry;
137 m_dFocalLength = dFocalLength;
138 m_dSourceDetectorLength = dSourceDetectorLength;
140 if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
144 m_traceLevel = iTraceLevel;
145 m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
146 m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
147 m_dBandwidth = dBandwidth;
148 m_nSignalPoints = nSignalPoints;
149 m_dSignalInc = dSignalIncrement;
150 m_dFilterParam = dFilterParam;
151 m_iZeropad = iZeropad;
152 m_iPreinterpolationFactor = iPreinterpolationFactor;
154 // scale signalInc/BW to adjust for imaginary detector through origin of phantom
155 // see Kak-Slaney Fig 3.22, for Collinear diagram
156 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
157 double dEquilinearScale = m_dSourceDetectorLength / m_dFocalLength;
158 m_dSignalInc /= dEquilinearScale;
159 m_dBandwidth *= dEquilinearScale;
162 if (m_idFilterMethod == FILTER_METHOD_FFT) {
164 m_idFilterMethod = FILTER_METHOD_RFFTW;
167 m_failMessage = "FFT not yet implemented";
172 bool m_bFrequencyFiltering = true;
173 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
174 m_bFrequencyFiltering = false;
176 // Spatial-based filtering
177 if (! m_bFrequencyFiltering) {
179 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
180 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
181 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
182 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
183 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
184 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
185 m_adFilter = new double[ m_nFilterPoints ];
186 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
187 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
188 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
189 m_dFilterMin = -1. / (2 * m_dSignalInc);
190 m_dFilterMax = 1. / (2 * m_dSignalInc);
191 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
192 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
193 m_adFilter = new double[ m_nFilterPoints ];
194 double* adFrequencyFilter = new double [m_nFilterPoints];
195 filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
196 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
197 if (g_bRunningWXWindows && m_traceLevel > 0) {
198 EZPlotDialog dlgEZPlot;
199 dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Natural Order");
200 dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints);
201 dlgEZPlot.ShowModal();
204 Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
206 if (g_bRunningWXWindows && m_traceLevel > 0) {
207 EZPlotDialog dlgEZPlot;
208 dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Fourier Order");
209 dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints);
210 dlgEZPlot.ShowModal();
213 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
214 delete adFrequencyFilter;
215 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
216 if (g_bRunningWXWindows && m_traceLevel > 0) {
217 EZPlotDialog dlgEZPlot;
218 dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Fourier Order");
219 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
220 dlgEZPlot.ShowModal();
223 Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
224 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
225 if (g_bRunningWXWindows && m_traceLevel > 0) {
226 EZPlotDialog dlgEZPlot;
227 dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Natural Order");
228 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
229 dlgEZPlot.ShowModal();
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 = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
243 double dScale = 0.5 * sinScale * sinScale;
244 m_adFilter[i] *= dScale;
246 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
247 if (g_bRunningWXWindows && m_traceLevel > 0) {
248 EZPlotDialog dlgEZPlot;
249 dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Fourier Frequency: Natural Order");
250 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
251 dlgEZPlot.ShowModal();
255 } // if (spatial filtering)
257 else if (m_bFrequencyFiltering) { // Frequency-based filtering
259 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
260 // calculate number of filter points with zeropadding
261 m_nFilterPoints = m_nSignalPoints;
262 if (m_iZeropad > 0) {
263 double logBase2 = log(m_nFilterPoints) / log(2);
264 int nextPowerOf2 = static_cast<int>(floor(logBase2));
265 if (logBase2 != floor(logBase2))
267 nextPowerOf2 += (m_iZeropad - 1);
268 m_nFilterPoints = 1 << nextPowerOf2;
269 #if defined(DEBUG) || defined(_DEBUG)
270 if (m_traceLevel >= Trace::TRACE_CONSOLE)
271 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
274 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
276 if (m_nFilterPoints % 2) { // Odd
277 m_dFilterMin = -1. / (2 * m_dSignalInc);
278 m_dFilterMax = 1. / (2 * m_dSignalInc);
279 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
281 m_dFilterMin = -1. / (2 * m_dSignalInc);
282 m_dFilterMax = 1. / (2 * m_dSignalInc);
283 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
284 m_dFilterMax -= m_dFilterInc;
287 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
288 m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
289 m_adFilter = new double [m_nFilterPoints];
290 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
292 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
293 if (g_bRunningWXWindows && m_traceLevel > 0) {
294 EZPlotDialog dlgEZPlot;
295 dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
296 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
297 dlgEZPlot.ShowModal();
301 // This works fairly well. I'm not sure why since scaling for geometries is done on
302 // frequency filter rather than spatial filter as it should be.
303 // It gives values slightly off than freq/inverse filtering
304 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
305 for (i = 0; i < m_nFilterPoints; i++)
306 m_adFilter[i] *= 0.5;
307 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
308 for (i = 0; i < m_nFilterPoints; i++) {
309 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
310 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
311 double dScale = 0.5 * sinScale * sinScale;
312 m_adFilter[i] *= dScale;
315 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
316 if (g_bRunningWXWindows && m_traceLevel > 0) {
317 EZPlotDialog dlgEZPlot;
318 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
319 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
320 dlgEZPlot.ShowModal();
323 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
324 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
325 if (g_bRunningWXWindows && m_traceLevel > 0) {
326 EZPlotDialog dlgEZPlot;
327 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
328 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
329 dlgEZPlot.ShowModal();
333 // FILTERING: FREQUENCY - INVERSE FOURIER
335 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
336 // calculate number of filter points with zeropadding
337 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
338 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
339 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
340 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
341 m_nFilterPoints = nSpatialPoints;
342 if (m_iZeropad > 0) {
343 double logBase2 = log(nSpatialPoints) / log(2);
344 int nextPowerOf2 = static_cast<int>(floor(logBase2));
345 if (logBase2 != floor(logBase2))
347 nextPowerOf2 += (m_iZeropad - 1);
348 m_nFilterPoints = 1 << nextPowerOf2;
350 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
351 #if defined(DEBUG) || defined(_DEBUG)
352 if (m_traceLevel >= Trace::TRACE_CONSOLE)
353 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
355 double* adSpatialFilter = new double [m_nFilterPoints];
356 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
357 m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
358 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
359 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
360 if (g_bRunningWXWindows && m_traceLevel > 0) {
361 EZPlotDialog dlgEZPlot;;
362 dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
363 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
364 dlgEZPlot.ShowModal();
368 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
369 for (i = 0; i < nSpatialPoints; i++)
370 adSpatialFilter[i] *= 0.5;
371 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
372 for (i = 0; i < nSpatialPoints; i++) {
373 int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
374 double sinScale = sin (iDetFromZero * m_dSignalInc);
375 if (fabs(sinScale) < 1E-7)
378 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
379 double dScale = 0.5 * sinScale * sinScale;
380 adSpatialFilter[i] *= dScale;
383 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
384 if (g_bRunningWXWindows && m_traceLevel > 0) {
385 EZPlotDialog dlgEZPlot;;
386 dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
387 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
388 dlgEZPlot.ShowModal();
391 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
392 adSpatialFilter[i] = 0;
394 m_adFilter = new double [m_nFilterPoints];
395 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
396 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
397 delete adSpatialFilter;
398 for (i = 0; i < m_nFilterPoints; i++)
399 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
400 delete acInverseFilter;
401 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
402 if (g_bRunningWXWindows && m_traceLevel > 0) {
403 EZPlotDialog dlgEZPlot;
404 dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
405 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
406 dlgEZPlot.ShowModal();
412 // precalculate sin and cosine tables for fourier transform
413 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
414 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
415 double angleIncrement = (2. * PI) / m_nFilterPoints;
416 m_adFourierCosTable = new double[ nFourier ];
417 m_adFourierSinTable = new double[ nFourier ];
419 for (i = 0; i < nFourier; i++) {
420 m_adFourierCosTable[i] = cos (angle);
421 m_adFourierSinTable[i] = sin (angle);
422 angle += angleIncrement;
427 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
428 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
429 m_adFilter[i] /= m_nFilterPoints;
432 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
433 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
434 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
435 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
436 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
437 for (i = 0; i < m_nFilterPoints; i++)
438 m_adRealFftInput[i] = 0;
439 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
440 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
441 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
442 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
443 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
444 for (i = 0; i < m_nFilterPoints; i++)
445 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
446 for (i = 0; i < m_nOutputPoints; i++)
447 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
453 ProcessSignal::~ProcessSignal (void)
455 delete [] m_adFourierSinTable;
456 delete [] m_adFourierCosTable;
457 delete [] m_adFilter;
460 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
461 fftw_destroy_plan(m_complexPlanForward);
462 fftw_destroy_plan(m_complexPlanBackward);
463 delete [] m_adComplexFftInput;
464 delete [] m_adComplexFftSignal;
466 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
467 rfftw_destroy_plan(m_realPlanForward);
468 rfftw_destroy_plan(m_realPlanBackward);
469 delete [] m_adRealFftInput;
470 delete [] m_adRealFftSignal;
476 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
478 int fmID = FILTER_METHOD_INVALID;
480 for (int i = 0; i < s_iFilterMethodCount; i++)
481 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
490 ProcessSignal::convertFilterMethodIDToName (const int fmID)
492 static const char *name = "";
494 if (fmID >= 0 && fmID < s_iFilterMethodCount)
495 return (s_aszFilterMethodName [fmID]);
501 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
503 static const char *title = "";
505 if (fmID >= 0 && fmID < s_iFilterMethodCount)
506 return (s_aszFilterMethodTitle [fmID]);
513 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
515 int fgID = FILTER_GENERATION_INVALID;
517 for (int i = 0; i < s_iFilterGenerationCount; i++)
518 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
527 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
529 static const char *name = "";
531 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
532 return (s_aszFilterGenerationName [fgID]);
538 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
540 static const char *name = "";
542 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
543 return (s_aszFilterGenerationTitle [fgID]);
549 ProcessSignal::filterSignal (const float constInput[], double output[]) const
551 double* input = new double [m_nSignalPoints];
553 for (i = 0; i < m_nSignalPoints; i++)
554 input[i] = constInput[i];
556 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
557 for (int i = 0; i < m_nSignalPoints; i++) {
558 int iDetFromCenter = i - (m_nSignalPoints / 2);
559 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
561 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
562 for (int i = 0; i < m_nSignalPoints; i++) {
563 int iDetFromCenter = i - (m_nSignalPoints / 2);
564 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
567 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
568 for (i = 0; i < m_nSignalPoints; i++)
569 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
570 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
571 double* inputSignal = new double [m_nFilterPoints];
572 for (i = 0; i < m_nSignalPoints; i++)
573 inputSignal[i] = input[i];
574 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
575 inputSignal[i] = 0; // zeropad
576 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
577 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
579 for (i = 0; i < m_nFilterPoints; i++)
580 fftSignal[i] *= m_adFilter[i];
581 double* inverseFourier = new double [m_nFilterPoints];
582 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
584 for (i = 0; i < m_nSignalPoints; i++)
585 output[i] = inverseFourier[i];
586 delete inverseFourier;
587 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
588 double* inputSignal = new double [m_nFilterPoints];
589 for (i = 0; i < m_nSignalPoints; i++)
590 inputSignal[i] = input[i];
591 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
592 inputSignal[i] = 0; // zeropad
593 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
594 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
596 for (i = 0; i < m_nFilterPoints; i++)
597 fftSignal[i] *= m_adFilter[i];
598 double* inverseFourier = new double [m_nFilterPoints];
599 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
601 for (i = 0; i < m_nSignalPoints; i++)
602 output[i] = inverseFourier[i];
603 delete inverseFourier;
606 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
607 for (i = 0; i < m_nSignalPoints; i++)
608 m_adRealFftInput[i] = input[i];
610 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
611 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
612 for (i = 0; i < m_nFilterPoints; i++)
613 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
615 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
616 m_adRealFftSignal[i] = 0;
618 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
619 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
620 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
621 output[i] = ifftOutput[i];
622 delete [] ifftOutput;
623 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
624 for (i = 0; i < m_nSignalPoints; i++)
625 m_adComplexFftInput[i].re = input[i];
627 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
628 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
629 for (i = 0; i < m_nFilterPoints; i++) {
630 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
631 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
634 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
635 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
636 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
637 output[i] = ifftOutput[i].re;
638 delete [] ifftOutput;
646 * convolve Discrete convolution of two functions
649 * r = convolve (f1, f2, dx, n, np, func_type)
650 * double r Convolved result
651 * double f1[], f2[] Functions to be convolved
652 * double dx Difference between successive x values
653 * int n Array index to center convolution about
654 * int np Number of points in f1 array
655 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
658 * f1 is the projection data, its indices range from 0 to np - 1.
659 * The index for f2, the filter, ranges from -(np-1) to (np-1).
660 * There are 3 ways to handle the negative vertices of f2:
661 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
662 * All filters used in reconstruction are even.
663 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
664 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
665 * for negative indices. Since f2 must range from -(np-1) to (np-1),
666 * if we add (np - 1) to f2's array index, then f2's index will
667 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
668 * stored at f2[np-1].
672 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
676 #if UNOPTIMIZED_CONVOLUTION
677 for (int i = 0; i < np; i++)
678 sum += func[i] * m_adFilter[n - i + (np - 1)];
680 double* f2 = m_adFilter + n + (np - 1);
681 for (int i = 0; i < np; i++)
682 sum += *func++ * *f2--;
690 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
694 #if UNOPTIMIZED_CONVOLUTION
695 for (int i = 0; i < np; i++)
696 sum += func[i] * m_adFilter[n - i + (np - 1)];
698 double* f2 = m_adFilter + n + (np - 1);
699 for (int i = 0; i < np; i++)
700 sum += *func++ * *f2--;
708 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
710 std::complex<double>* complexOutput = new std::complex<double> [n];
712 finiteFourierTransform (input, complexOutput, n, direction);
713 for (int i = 0; i < n; i++)
714 output[i] = complexOutput[i].real();
715 delete [] complexOutput;
719 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
726 double angleIncrement = direction * 2 * PI / n;
727 for (int i = 0; i < n; i++) {
730 for (int j = 0; j < n; j++) {
731 double angle = i * j * angleIncrement;
732 sumReal += input[j] * cos(angle);
733 sumImag += input[j] * sin(angle);
739 output[i] = std::complex<double> (sumReal, sumImag);
745 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
752 double angleIncrement = direction * 2 * PI / n;
753 for (int i = 0; i < n; i++) {
754 std::complex<double> sum (0,0);
755 for (int j = 0; j < n; j++) {
756 double angle = i * j * angleIncrement;
757 std::complex<double> exponentTerm (cos(angle), sin(angle));
758 sum += input[j] * exponentTerm;
768 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
775 double angleIncrement = direction * 2 * PI / n;
776 for (int i = 0; i < n; i++) {
778 for (int j = 0; j < n; j++) {
779 double angle = i * j * angleIncrement;
780 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
789 // Table-based routines
792 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
799 for (int i = 0; i < m_nFilterPoints; i++) {
800 double sumReal = 0, sumImag = 0;
801 for (int j = 0; j < m_nFilterPoints; j++) {
802 int tableIndex = i * j;
804 sumReal += input[j] * m_adFourierCosTable[tableIndex];
805 sumImag += input[j] * m_adFourierSinTable[tableIndex];
807 sumReal += input[j] * m_adFourierCosTable[tableIndex];
808 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
812 sumReal /= m_nFilterPoints;
813 sumImag /= m_nFilterPoints;
815 output[i] = std::complex<double> (sumReal, sumImag);
819 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
821 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
828 for (int i = 0; i < m_nFilterPoints; i++) {
829 double sumReal = 0, sumImag = 0;
830 for (int j = 0; j < m_nFilterPoints; j++) {
831 int tableIndex = i * j;
833 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
834 - input[j].imag() * m_adFourierSinTable[tableIndex];
835 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
836 + input[j].imag() * m_adFourierCosTable[tableIndex];
838 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
839 - input[j].imag() * -m_adFourierSinTable[tableIndex];
840 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
841 + input[j].imag() * m_adFourierCosTable[tableIndex];
845 sumReal /= m_nFilterPoints;
846 sumImag /= m_nFilterPoints;
848 output[i] = std::complex<double> (sumReal, sumImag);
853 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
860 for (int i = 0; i < m_nFilterPoints; i++) {
862 for (int j = 0; j < m_nFilterPoints; j++) {
863 int tableIndex = i * j;
865 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
866 - input[j].imag() * m_adFourierSinTable[tableIndex];
868 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
869 - input[j].imag() * -m_adFourierSinTable[tableIndex];
873 sumReal /= m_nFilterPoints;