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.20 2001/01/12 21:53:27 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 doesn't work: Need to add filtering for divergent geometries & Frequency/Direct filtering
301 // Jan 2001: Direct seems to work for equilinear and equiangular
302 // however, inverse_fourier doesn't work for equiangular on all versions of CTSim tested
303 // Scaling is done with data in frequency space, natural order
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 #define PRE_JAN_2001 1
370 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
371 for (i = 0; i < nSpatialPoints; i++)
372 adSpatialFilter[i] *= 0.5;
373 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
374 for (i = 0; i < nSpatialPoints; i++) {
375 int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
376 double sinScale = sin (iDetFromZero * m_dSignalInc);
377 if (fabs(sinScale) < 1E-7)
380 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
381 double dScale = 0.5 * sinScale * sinScale;
382 adSpatialFilter[i] *= dScale;
385 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
386 if (g_bRunningWXWindows && m_traceLevel > 0) {
387 EZPlotDialog dlgEZPlot;;
388 dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
389 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
390 dlgEZPlot.ShowModal();
393 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
394 adSpatialFilter[i] = 0;
396 m_adFilter = new double [m_nFilterPoints];
397 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
398 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
399 delete adSpatialFilter;
400 for (i = 0; i < m_nFilterPoints; i++)
401 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
402 delete acInverseFilter;
403 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
404 if (g_bRunningWXWindows && m_traceLevel > 0) {
405 EZPlotDialog dlgEZPlot;
406 dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
407 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
408 dlgEZPlot.ShowModal();
413 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
414 adSpatialFilter[i] = 0;
416 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
417 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, FORWARD);
418 delete adSpatialFilter;
419 m_adFilter = new double [m_nFilterPoints];
420 for (i = 0; i < m_nFilterPoints; i++)
421 m_adFilter[i] = std::abs(acInverseFilter[i]);
422 delete acInverseFilter;
424 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
425 if (g_bRunningWXWindows && m_traceLevel > 0) {
426 EZPlotDialog dlgEZPlot;
427 dlgEZPlot.getEZPlot()->ezset ("title Inverse Spatial Filter: Fourier order");
428 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
429 dlgEZPlot.ShowModal();
432 Fourier::shuffleFourierToNaturalOrder(m_adFilter, m_nFilterPoints);
433 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
434 if (g_bRunningWXWindows && m_traceLevel > 0) {
435 EZPlotDialog dlgEZPlot;
436 dlgEZPlot.getEZPlot()->ezset ("title Inverse Spatial Filter: Natural order");
437 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
438 dlgEZPlot.ShowModal();
441 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
442 for (i = 0; i < m_nFilterPoints; i++)
443 m_adFilter[i] *= 0.5;
444 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
445 for (i = 0; i < m_nFilterPoints; i++) {
446 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
447 if (abs(iDetFromZero) < m_nSignalPoints) {
448 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
449 double dScale = 0.5 * sinScale * sinScale;
450 m_adFilter[i] *= dScale;
454 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
455 if (g_bRunningWXWindows && m_traceLevel > 0) {
456 EZPlotDialog dlgEZPlot;
457 dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Spatial Filter: Natural order");
458 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
459 dlgEZPlot.ShowModal();
462 Fourier::shuffleNaturalToFourierOrder(m_adFilter, m_nFilterPoints);
465 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
466 if (g_bRunningWXWindows && m_traceLevel > 0) {
467 EZPlotDialog dlgEZPlot;
468 dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter Inverse Post Geometry Filtering");
469 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
470 dlgEZPlot.ShowModal();
476 // precalculate sin and cosine tables for fourier transform
477 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
478 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
479 double angleIncrement = (2. * PI) / m_nFilterPoints;
480 m_adFourierCosTable = new double[ nFourier ];
481 m_adFourierSinTable = new double[ nFourier ];
483 for (i = 0; i < nFourier; i++) {
484 m_adFourierCosTable[i] = cos (angle);
485 m_adFourierSinTable[i] = sin (angle);
486 angle += angleIncrement;
491 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
492 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
493 m_adFilter[i] /= m_nFilterPoints;
496 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
497 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
498 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
499 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
500 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
501 for (i = 0; i < m_nFilterPoints; i++)
502 m_adRealFftInput[i] = 0;
503 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
504 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
505 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
506 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
507 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
508 for (i = 0; i < m_nFilterPoints; i++)
509 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
510 for (i = 0; i < m_nOutputPoints; i++)
511 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
517 ProcessSignal::~ProcessSignal (void)
519 delete [] m_adFourierSinTable;
520 delete [] m_adFourierCosTable;
521 delete [] m_adFilter;
524 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
525 fftw_destroy_plan(m_complexPlanForward);
526 fftw_destroy_plan(m_complexPlanBackward);
527 delete [] m_adComplexFftInput;
528 delete [] m_adComplexFftSignal;
530 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
531 rfftw_destroy_plan(m_realPlanForward);
532 rfftw_destroy_plan(m_realPlanBackward);
533 delete [] m_adRealFftInput;
534 delete [] m_adRealFftSignal;
540 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
542 int fmID = FILTER_METHOD_INVALID;
544 for (int i = 0; i < s_iFilterMethodCount; i++)
545 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
554 ProcessSignal::convertFilterMethodIDToName (const int fmID)
556 static const char *name = "";
558 if (fmID >= 0 && fmID < s_iFilterMethodCount)
559 return (s_aszFilterMethodName [fmID]);
565 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
567 static const char *title = "";
569 if (fmID >= 0 && fmID < s_iFilterMethodCount)
570 return (s_aszFilterMethodTitle [fmID]);
577 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
579 int fgID = FILTER_GENERATION_INVALID;
581 for (int i = 0; i < s_iFilterGenerationCount; i++)
582 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
591 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
593 static const char *name = "";
595 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
596 return (s_aszFilterGenerationName [fgID]);
602 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
604 static const char *name = "";
606 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
607 return (s_aszFilterGenerationTitle [fgID]);
613 ProcessSignal::filterSignal (const float constInput[], double output[]) const
615 double* input = new double [m_nSignalPoints];
617 for (i = 0; i < m_nSignalPoints; i++)
618 input[i] = constInput[i];
620 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
621 for (int i = 0; i < m_nSignalPoints; i++) {
622 int iDetFromCenter = i - (m_nSignalPoints / 2);
623 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
625 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
626 for (int i = 0; i < m_nSignalPoints; i++) {
627 int iDetFromCenter = i - (m_nSignalPoints / 2);
628 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
631 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
632 for (i = 0; i < m_nSignalPoints; i++)
633 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
634 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
635 double* inputSignal = new double [m_nFilterPoints];
636 for (i = 0; i < m_nSignalPoints; i++)
637 inputSignal[i] = input[i];
638 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
639 inputSignal[i] = 0; // zeropad
640 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
641 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
643 for (i = 0; i < m_nFilterPoints; i++)
644 fftSignal[i] *= m_adFilter[i];
645 double* inverseFourier = new double [m_nFilterPoints];
646 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
648 for (i = 0; i < m_nSignalPoints; i++)
649 output[i] = inverseFourier[i];
650 delete inverseFourier;
651 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
652 double* inputSignal = new double [m_nFilterPoints];
653 for (i = 0; i < m_nSignalPoints; i++)
654 inputSignal[i] = input[i];
655 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
656 inputSignal[i] = 0; // zeropad
657 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
658 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
660 for (i = 0; i < m_nFilterPoints; i++)
661 fftSignal[i] *= m_adFilter[i];
662 double* inverseFourier = new double [m_nFilterPoints];
663 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
665 for (i = 0; i < m_nSignalPoints; i++)
666 output[i] = inverseFourier[i];
667 delete inverseFourier;
670 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
671 for (i = 0; i < m_nSignalPoints; i++)
672 m_adRealFftInput[i] = input[i];
674 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
675 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
676 for (i = 0; i < m_nFilterPoints; i++)
677 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
679 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
680 m_adRealFftSignal[i] = 0;
682 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
683 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
684 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
685 output[i] = ifftOutput[i];
686 delete [] ifftOutput;
687 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
688 for (i = 0; i < m_nSignalPoints; i++)
689 m_adComplexFftInput[i].re = input[i];
691 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
692 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
693 for (i = 0; i < m_nFilterPoints; i++) {
694 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
695 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
698 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
699 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
700 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
701 output[i] = ifftOutput[i].re;
702 delete [] ifftOutput;
710 * convolve Discrete convolution of two functions
713 * r = convolve (f1, f2, dx, n, np, func_type)
714 * double r Convolved result
715 * double f1[], f2[] Functions to be convolved
716 * double dx Difference between successive x values
717 * int n Array index to center convolution about
718 * int np Number of points in f1 array
719 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
722 * f1 is the projection data, its indices range from 0 to np - 1.
723 * The index for f2, the filter, ranges from -(np-1) to (np-1).
724 * There are 3 ways to handle the negative vertices of f2:
725 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
726 * All filters used in reconstruction are even.
727 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
728 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
729 * for negative indices. Since f2 must range from -(np-1) to (np-1),
730 * if we add (np - 1) to f2's array index, then f2's index will
731 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
732 * stored at f2[np-1].
736 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
740 #if UNOPTIMIZED_CONVOLUTION
741 for (int i = 0; i < np; i++)
742 sum += func[i] * m_adFilter[n - i + (np - 1)];
744 double* f2 = m_adFilter + n + (np - 1);
745 for (int i = 0; i < np; i++)
746 sum += *func++ * *f2--;
754 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
758 #if UNOPTIMIZED_CONVOLUTION
759 for (int i = 0; i < np; i++)
760 sum += func[i] * m_adFilter[n - i + (np - 1)];
762 double* f2 = m_adFilter + n + (np - 1);
763 for (int i = 0; i < np; i++)
764 sum += *func++ * *f2--;
772 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
774 std::complex<double>* complexOutput = new std::complex<double> [n];
776 finiteFourierTransform (input, complexOutput, n, direction);
777 for (int i = 0; i < n; i++)
778 output[i] = complexOutput[i].real();
779 delete [] complexOutput;
783 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
790 double angleIncrement = direction * 2 * PI / n;
791 for (int i = 0; i < n; i++) {
794 for (int j = 0; j < n; j++) {
795 double angle = i * j * angleIncrement;
796 sumReal += input[j] * cos(angle);
797 sumImag += input[j] * sin(angle);
803 output[i] = std::complex<double> (sumReal, sumImag);
809 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
816 double angleIncrement = direction * 2 * PI / n;
817 for (int i = 0; i < n; i++) {
818 std::complex<double> sum (0,0);
819 for (int j = 0; j < n; j++) {
820 double angle = i * j * angleIncrement;
821 std::complex<double> exponentTerm (cos(angle), sin(angle));
822 sum += input[j] * exponentTerm;
832 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
839 double angleIncrement = direction * 2 * PI / n;
840 for (int i = 0; i < n; i++) {
842 for (int j = 0; j < n; j++) {
843 double angle = i * j * angleIncrement;
844 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
853 // Table-based routines
856 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
863 for (int i = 0; i < m_nFilterPoints; i++) {
864 double sumReal = 0, sumImag = 0;
865 for (int j = 0; j < m_nFilterPoints; j++) {
866 int tableIndex = i * j;
868 sumReal += input[j] * m_adFourierCosTable[tableIndex];
869 sumImag += input[j] * m_adFourierSinTable[tableIndex];
871 sumReal += input[j] * m_adFourierCosTable[tableIndex];
872 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
876 sumReal /= m_nFilterPoints;
877 sumImag /= m_nFilterPoints;
879 output[i] = std::complex<double> (sumReal, sumImag);
883 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
885 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
892 for (int i = 0; i < m_nFilterPoints; i++) {
893 double sumReal = 0, sumImag = 0;
894 for (int j = 0; j < m_nFilterPoints; j++) {
895 int tableIndex = i * j;
897 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
898 - input[j].imag() * m_adFourierSinTable[tableIndex];
899 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
900 + input[j].imag() * m_adFourierCosTable[tableIndex];
902 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
903 - input[j].imag() * -m_adFourierSinTable[tableIndex];
904 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
905 + input[j].imag() * m_adFourierCosTable[tableIndex];
909 sumReal /= m_nFilterPoints;
910 sumImag /= m_nFilterPoints;
912 output[i] = std::complex<double> (sumReal, sumImag);
917 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
924 for (int i = 0; i < m_nFilterPoints; i++) {
926 for (int j = 0; j < m_nFilterPoints; j++) {
927 int tableIndex = i * j;
929 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
930 - input[j].imag() * m_adFourierSinTable[tableIndex];
932 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
933 - input[j].imag() * -m_adFourierSinTable[tableIndex];
937 sumReal /= m_nFilterPoints;