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.23 2001/01/13 04:52:01 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* 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 if (g_bRunningWXWindows && m_traceLevel > 0) {
196 EZPlotDialog dlgEZPlot;
197 dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Natural Order");
198 dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints);
199 dlgEZPlot.ShowModal();
202 Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
204 if (g_bRunningWXWindows && m_traceLevel > 0) {
205 EZPlotDialog dlgEZPlot;
206 dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Fourier Order");
207 dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints);
208 dlgEZPlot.ShowModal();
211 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
212 delete adFrequencyFilter;
213 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
214 if (g_bRunningWXWindows && m_traceLevel > 0) {
215 EZPlotDialog dlgEZPlot;
216 dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Fourier Order");
217 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
218 dlgEZPlot.ShowModal();
221 Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
222 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
223 if (g_bRunningWXWindows && m_traceLevel > 0) {
224 EZPlotDialog dlgEZPlot;
225 dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Natural Order");
226 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
227 dlgEZPlot.ShowModal();
230 for (i = 0; i < m_nFilterPoints; i++) {
231 m_adFilter[i] /= m_dSignalInc;
234 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
235 for (i = 0; i < m_nFilterPoints; i++)
236 m_adFilter[i] *= 0.5;
237 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
238 for (i = 0; i < m_nFilterPoints; i++) {
239 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
240 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
241 double dScale = 0.5 * sinScale * sinScale;
242 m_adFilter[i] *= dScale;
244 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
245 if (g_bRunningWXWindows && m_traceLevel > 0) {
246 EZPlotDialog dlgEZPlot;
247 dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Fourier Frequency: Natural Order");
248 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
249 dlgEZPlot.ShowModal();
253 } // if (spatial filtering)
255 else if (m_bFrequencyFiltering) { // Frequency-based filtering
257 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
258 // calculate number of filter points with zeropadding
259 m_nFilterPoints = m_nSignalPoints;
260 if (m_iZeropad > 0) {
261 double logBase2 = log(m_nFilterPoints) / log(2);
262 int nextPowerOf2 = static_cast<int>(floor(logBase2));
263 if (logBase2 != floor(logBase2))
265 nextPowerOf2 += (m_iZeropad - 1);
266 m_nFilterPoints = 1 << nextPowerOf2;
267 #if defined(DEBUG) || defined(_DEBUG)
268 if (m_traceLevel >= Trace::TRACE_CONSOLE)
269 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
272 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
274 if (m_nFilterPoints % 2) { // Odd
275 m_dFilterMin = -1. / (2 * m_dSignalInc);
276 m_dFilterMax = 1. / (2 * m_dSignalInc);
277 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
279 m_dFilterMin = -1. / (2 * m_dSignalInc);
280 m_dFilterMax = 1. / (2 * m_dSignalInc);
281 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
282 m_dFilterMax -= m_dFilterInc;
285 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
286 m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
287 m_adFilter = new double [m_nFilterPoints];
288 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
290 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
291 if (g_bRunningWXWindows && m_traceLevel > 0) {
292 EZPlotDialog dlgEZPlot;
293 dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
294 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
295 dlgEZPlot.ShowModal();
299 // This works fairly well. I'm not sure why since scaling for geometries is done on
300 // frequency filter rather than spatial filter as it should be.
301 // It gives values slightly off than freq/inverse filtering
302 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
303 for (i = 0; i < m_nFilterPoints; i++)
304 m_adFilter[i] *= 0.5;
305 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
306 for (i = 0; i < m_nFilterPoints; i++) {
307 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
308 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
309 double dScale = 0.5 * sinScale * sinScale;
310 m_adFilter[i] *= dScale;
313 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
314 if (g_bRunningWXWindows && m_traceLevel > 0) {
315 EZPlotDialog dlgEZPlot;
316 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
317 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
318 dlgEZPlot.ShowModal();
321 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
322 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
323 if (g_bRunningWXWindows && m_traceLevel > 0) {
324 EZPlotDialog dlgEZPlot;
325 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
326 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
327 dlgEZPlot.ShowModal();
331 // FILTERING: FREQUENCY - INVERSE FOURIER
333 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
334 // calculate number of filter points with zeropadding
335 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
336 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
337 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
338 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
339 m_nFilterPoints = nSpatialPoints;
340 if (m_iZeropad > 0) {
341 double logBase2 = log(nSpatialPoints) / log(2);
342 int nextPowerOf2 = static_cast<int>(floor(logBase2));
343 if (logBase2 != floor(logBase2))
345 nextPowerOf2 += (m_iZeropad - 1);
346 m_nFilterPoints = 1 << nextPowerOf2;
348 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
349 #if defined(DEBUG) || defined(_DEBUG)
350 if (m_traceLevel >= Trace::TRACE_CONSOLE)
351 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
353 double* adSpatialFilter = new double [m_nFilterPoints];
354 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
355 m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
356 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
357 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
358 if (g_bRunningWXWindows && m_traceLevel > 0) {
359 EZPlotDialog dlgEZPlot;;
360 dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
361 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
362 dlgEZPlot.ShowModal();
366 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
367 for (i = 0; i < nSpatialPoints; i++)
368 adSpatialFilter[i] *= 0.5;
369 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
370 for (i = 0; i < nSpatialPoints; i++) {
371 int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
372 double sinScale = sin (iDetFromZero * m_dSignalInc);
373 if (fabs(sinScale) < 1E-7)
376 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
377 double dScale = 0.5 * sinScale * sinScale;
378 adSpatialFilter[i] *= dScale;
381 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
382 if (g_bRunningWXWindows && m_traceLevel > 0) {
383 EZPlotDialog dlgEZPlot;;
384 dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
385 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
386 dlgEZPlot.ShowModal();
389 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
390 adSpatialFilter[i] = 0;
392 m_adFilter = new double [m_nFilterPoints];
393 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
394 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
395 delete adSpatialFilter;
396 for (i = 0; i < m_nFilterPoints; i++)
397 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
398 delete acInverseFilter;
399 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
400 if (g_bRunningWXWindows && m_traceLevel > 0) {
401 EZPlotDialog dlgEZPlot;
402 dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
403 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
404 dlgEZPlot.ShowModal();
410 // precalculate sin and cosine tables for fourier transform
411 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
412 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
413 double angleIncrement = (2. * PI) / m_nFilterPoints;
414 m_adFourierCosTable = new double[ nFourier ];
415 m_adFourierSinTable = new double[ nFourier ];
417 for (i = 0; i < nFourier; i++) {
418 m_adFourierCosTable[i] = cos (angle);
419 m_adFourierSinTable[i] = sin (angle);
420 angle += angleIncrement;
425 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
426 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
427 m_adFilter[i] /= m_nFilterPoints;
430 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
431 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
432 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
433 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
434 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
435 for (i = 0; i < m_nFilterPoints; i++)
436 m_adRealFftInput[i] = 0;
437 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
438 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
439 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
440 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
441 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
442 for (i = 0; i < m_nFilterPoints; i++)
443 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
444 for (i = 0; i < m_nOutputPoints; i++)
445 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
451 ProcessSignal::~ProcessSignal (void)
453 delete [] m_adFourierSinTable;
454 delete [] m_adFourierCosTable;
455 delete [] m_adFilter;
458 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
459 fftw_destroy_plan(m_complexPlanForward);
460 fftw_destroy_plan(m_complexPlanBackward);
461 delete [] m_adComplexFftInput;
462 delete [] m_adComplexFftSignal;
464 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
465 rfftw_destroy_plan(m_realPlanForward);
466 rfftw_destroy_plan(m_realPlanBackward);
467 delete [] m_adRealFftInput;
468 delete [] m_adRealFftSignal;
474 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
476 int fmID = FILTER_METHOD_INVALID;
478 for (int i = 0; i < s_iFilterMethodCount; i++)
479 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
488 ProcessSignal::convertFilterMethodIDToName (const int fmID)
490 static const char *name = "";
492 if (fmID >= 0 && fmID < s_iFilterMethodCount)
493 return (s_aszFilterMethodName [fmID]);
499 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
501 static const char *title = "";
503 if (fmID >= 0 && fmID < s_iFilterMethodCount)
504 return (s_aszFilterMethodTitle [fmID]);
511 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
513 int fgID = FILTER_GENERATION_INVALID;
515 for (int i = 0; i < s_iFilterGenerationCount; i++)
516 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
525 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
527 static const char *name = "";
529 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
530 return (s_aszFilterGenerationName [fgID]);
536 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
538 static const char *name = "";
540 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
541 return (s_aszFilterGenerationTitle [fgID]);
547 ProcessSignal::filterSignal (const float constInput[], double output[]) const
549 double* input = new double [m_nSignalPoints];
551 for (i = 0; i < m_nSignalPoints; i++)
552 input[i] = constInput[i];
554 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
555 for (int i = 0; i < m_nSignalPoints; i++) {
556 int iDetFromCenter = i - (m_nSignalPoints / 2);
557 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
559 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
560 for (int i = 0; i < m_nSignalPoints; i++) {
561 int iDetFromCenter = i - (m_nSignalPoints / 2);
562 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
565 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
566 for (i = 0; i < m_nSignalPoints; i++)
567 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
568 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
569 double* inputSignal = new double [m_nFilterPoints];
570 for (i = 0; i < m_nSignalPoints; i++)
571 inputSignal[i] = input[i];
572 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
573 inputSignal[i] = 0; // zeropad
574 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
575 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
577 for (i = 0; i < m_nFilterPoints; i++)
578 fftSignal[i] *= m_adFilter[i];
579 double* inverseFourier = new double [m_nFilterPoints];
580 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
582 for (i = 0; i < m_nSignalPoints; i++)
583 output[i] = inverseFourier[i];
584 delete inverseFourier;
585 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
586 double* inputSignal = new double [m_nFilterPoints];
587 for (i = 0; i < m_nSignalPoints; i++)
588 inputSignal[i] = input[i];
589 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
590 inputSignal[i] = 0; // zeropad
591 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
592 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
594 for (i = 0; i < m_nFilterPoints; i++)
595 fftSignal[i] *= m_adFilter[i];
596 double* inverseFourier = new double [m_nFilterPoints];
597 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
599 for (i = 0; i < m_nSignalPoints; i++)
600 output[i] = inverseFourier[i];
601 delete inverseFourier;
604 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
605 for (i = 0; i < m_nSignalPoints; i++)
606 m_adRealFftInput[i] = input[i];
608 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
609 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
610 for (i = 0; i < m_nFilterPoints; i++)
611 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
613 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
614 m_adRealFftSignal[i] = 0;
616 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
617 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
618 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
619 output[i] = ifftOutput[i];
620 delete [] ifftOutput;
621 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
622 for (i = 0; i < m_nSignalPoints; i++)
623 m_adComplexFftInput[i].re = input[i];
625 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
626 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
627 for (i = 0; i < m_nFilterPoints; i++) {
628 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
629 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
632 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
633 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
634 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
635 output[i] = ifftOutput[i].re;
636 delete [] ifftOutput;
644 * convolve Discrete convolution of two functions
647 * r = convolve (f1, f2, dx, n, np, func_type)
648 * double r Convolved result
649 * double f1[], f2[] Functions to be convolved
650 * double dx Difference between successive x values
651 * int n Array index to center convolution about
652 * int np Number of points in f1 array
653 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
656 * f1 is the projection data, its indices range from 0 to np - 1.
657 * The index for f2, the filter, ranges from -(np-1) to (np-1).
658 * There are 3 ways to handle the negative vertices of f2:
659 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
660 * All filters used in reconstruction are even.
661 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
662 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
663 * for negative indices. Since f2 must range from -(np-1) to (np-1),
664 * if we add (np - 1) to f2's array index, then f2's index will
665 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
666 * stored at f2[np-1].
670 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
674 #if UNOPTIMIZED_CONVOLUTION
675 for (int i = 0; i < np; i++)
676 sum += func[i] * m_adFilter[n - i + (np - 1)];
678 double* f2 = m_adFilter + n + (np - 1);
679 for (int i = 0; i < np; i++)
680 sum += *func++ * *f2--;
688 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
692 #if UNOPTIMIZED_CONVOLUTION
693 for (int i = 0; i < np; i++)
694 sum += func[i] * m_adFilter[n - i + (np - 1)];
696 double* f2 = m_adFilter + n + (np - 1);
697 for (int i = 0; i < np; i++)
698 sum += *func++ * *f2--;
706 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
708 std::complex<double>* complexOutput = new std::complex<double> [n];
710 finiteFourierTransform (input, complexOutput, n, direction);
711 for (int i = 0; i < n; i++)
712 output[i] = complexOutput[i].real();
713 delete [] complexOutput;
717 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
724 double angleIncrement = direction * 2 * PI / n;
725 for (int i = 0; i < n; i++) {
728 for (int j = 0; j < n; j++) {
729 double angle = i * j * angleIncrement;
730 sumReal += input[j] * cos(angle);
731 sumImag += input[j] * sin(angle);
737 output[i] = std::complex<double> (sumReal, sumImag);
743 ProcessSignal::finiteFourierTransform (const std::complex<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++) {
752 std::complex<double> sum (0,0);
753 for (int j = 0; j < n; j++) {
754 double angle = i * j * angleIncrement;
755 std::complex<double> exponentTerm (cos(angle), sin(angle));
756 sum += input[j] * exponentTerm;
766 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
773 double angleIncrement = direction * 2 * PI / n;
774 for (int i = 0; i < n; i++) {
776 for (int j = 0; j < n; j++) {
777 double angle = i * j * angleIncrement;
778 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
787 // Table-based routines
790 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
797 for (int i = 0; i < m_nFilterPoints; i++) {
798 double sumReal = 0, sumImag = 0;
799 for (int j = 0; j < m_nFilterPoints; j++) {
800 int tableIndex = i * j;
802 sumReal += input[j] * m_adFourierCosTable[tableIndex];
803 sumImag += input[j] * m_adFourierSinTable[tableIndex];
805 sumReal += input[j] * m_adFourierCosTable[tableIndex];
806 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
810 sumReal /= m_nFilterPoints;
811 sumImag /= m_nFilterPoints;
813 output[i] = std::complex<double> (sumReal, sumImag);
817 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
819 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
826 for (int i = 0; i < m_nFilterPoints; i++) {
827 double sumReal = 0, sumImag = 0;
828 for (int j = 0; j < m_nFilterPoints; j++) {
829 int tableIndex = i * j;
831 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
832 - input[j].imag() * m_adFourierSinTable[tableIndex];
833 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
834 + input[j].imag() * m_adFourierCosTable[tableIndex];
836 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
837 - input[j].imag() * -m_adFourierSinTable[tableIndex];
838 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
839 + input[j].imag() * m_adFourierCosTable[tableIndex];
843 sumReal /= m_nFilterPoints;
844 sumImag /= m_nFilterPoints;
846 output[i] = std::complex<double> (sumReal, sumImag);
851 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
858 for (int i = 0; i < m_nFilterPoints; i++) {
860 for (int j = 0; j < m_nFilterPoints; j++) {
861 int tableIndex = i * j;
863 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
864 - input[j].imag() * m_adFourierSinTable[tableIndex];
866 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
867 - input[j].imag() * -m_adFourierSinTable[tableIndex];
871 sumReal /= m_nFilterPoints;