1 /*****************************************************************************
4 ** Name: procsignal.cpp
5 ** Purpose: Routines for processing signals and projections
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.32 2001/03/30 19:17:32 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 "nographics.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 = addZeropadFactor (m_nSignalPoints, m_iZeropad);
262 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
264 if (isOdd (m_nFilterPoints)) { // Odd
265 m_dFilterMin = -1. / (2 * m_dSignalInc);
266 m_dFilterMax = 1. / (2 * m_dSignalInc);
267 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
269 m_dFilterMin = -1. / (2 * m_dSignalInc);
270 m_dFilterMax = 1. / (2 * m_dSignalInc);
271 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
272 m_dFilterMax -= m_dFilterInc;
275 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
276 m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
277 m_adFilter = new double [m_nFilterPoints];
278 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
280 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
281 if (g_bRunningWXWindows && m_traceLevel > 0) {
282 EZPlotDialog dlgEZPlot;
283 dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
284 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
285 dlgEZPlot.ShowModal();
289 // This works fairly well. I'm not sure why since scaling for geometries is done on
290 // frequency filter rather than spatial filter as it should be.
291 // It gives values slightly off than freq/inverse filtering
292 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
293 for (i = 0; i < m_nFilterPoints; i++)
294 m_adFilter[i] *= 0.5;
295 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
296 for (i = 0; i < m_nFilterPoints; i++) {
297 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
298 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
299 double dScale = 0.5 * sinScale * sinScale;
300 m_adFilter[i] *= dScale;
303 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
304 if (g_bRunningWXWindows && m_traceLevel > 0) {
305 EZPlotDialog dlgEZPlot;
306 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
307 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
308 dlgEZPlot.ShowModal();
311 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
312 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
313 if (g_bRunningWXWindows && m_traceLevel > 0) {
314 EZPlotDialog dlgEZPlot;
315 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
316 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
317 dlgEZPlot.ShowModal();
321 // FILTERING: FREQUENCY - INVERSE FOURIER
323 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
324 // calculate number of filter points with zeropadding
325 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
326 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
327 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
328 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
329 m_nFilterPoints = nSpatialPoints;
330 if (m_iZeropad > 0) {
331 double logBase2 = log(nSpatialPoints) / log(2);
332 int nextPowerOf2 = static_cast<int>(floor(logBase2));
333 if (logBase2 != floor(logBase2))
335 nextPowerOf2 += (m_iZeropad - 1);
336 m_nFilterPoints = 1 << nextPowerOf2;
338 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
339 #if defined(DEBUG) || defined(_DEBUG)
340 if (m_traceLevel >= Trace::TRACE_CONSOLE)
341 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
343 double* adSpatialFilter = new double [m_nFilterPoints];
344 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
345 m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
346 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
347 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
348 if (g_bRunningWXWindows && m_traceLevel > 0) {
349 EZPlotDialog dlgEZPlot;;
350 dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
351 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
352 dlgEZPlot.ShowModal();
356 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
357 for (i = 0; i < nSpatialPoints; i++)
358 adSpatialFilter[i] *= 0.5;
359 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
360 for (i = 0; i < nSpatialPoints; i++) {
361 int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
362 double sinScale = sin (iDetFromZero * m_dSignalInc);
363 if (fabs(sinScale) < 1E-7)
366 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
367 double dScale = 0.5 * sinScale * sinScale;
368 adSpatialFilter[i] *= dScale;
371 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
372 if (g_bRunningWXWindows && m_traceLevel > 0) {
373 EZPlotDialog dlgEZPlot;;
374 dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
375 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
376 dlgEZPlot.ShowModal();
379 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
380 adSpatialFilter[i] = 0;
382 m_adFilter = new double [m_nFilterPoints];
383 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
384 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
385 delete adSpatialFilter;
386 for (i = 0; i < m_nFilterPoints; i++)
387 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
388 delete acInverseFilter;
389 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
390 if (g_bRunningWXWindows && m_traceLevel > 0) {
391 EZPlotDialog dlgEZPlot;
392 dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
393 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
394 dlgEZPlot.ShowModal();
400 // precalculate sin and cosine tables for fourier transform
401 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
402 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
403 double angleIncrement = (2. * PI) / m_nFilterPoints;
404 m_adFourierCosTable = new double[ nFourier ];
405 m_adFourierSinTable = new double[ nFourier ];
407 for (i = 0; i < nFourier; i++) {
408 m_adFourierCosTable[i] = cos (angle);
409 m_adFourierSinTable[i] = sin (angle);
410 angle += angleIncrement;
415 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
416 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
417 m_adFilter[i] /= m_nFilterPoints;
420 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
421 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_USE_WISDOM);
422 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_USE_WISDOM);
423 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
424 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
425 for (i = 0; i < m_nFilterPoints; i++)
426 m_adRealFftInput[i] = 0;
427 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
428 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM);
429 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM);
430 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
431 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
432 for (i = 0; i < m_nFilterPoints; i++)
433 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
434 for (i = 0; i < m_nOutputPoints; i++)
435 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
441 ProcessSignal::~ProcessSignal (void)
443 delete [] m_adFourierSinTable;
444 delete [] m_adFourierCosTable;
445 delete [] m_adFilter;
448 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
449 fftw_destroy_plan(m_complexPlanForward);
450 fftw_destroy_plan(m_complexPlanBackward);
451 delete [] m_adComplexFftInput;
452 delete [] m_adComplexFftSignal;
454 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
455 rfftw_destroy_plan(m_realPlanForward);
456 rfftw_destroy_plan(m_realPlanBackward);
457 delete [] m_adRealFftInput;
458 delete [] m_adRealFftSignal;
464 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
466 int fmID = FILTER_METHOD_INVALID;
468 for (int i = 0; i < s_iFilterMethodCount; i++)
469 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
478 ProcessSignal::convertFilterMethodIDToName (const int fmID)
480 static const char *name = "";
482 if (fmID >= 0 && fmID < s_iFilterMethodCount)
483 return (s_aszFilterMethodName [fmID]);
489 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
491 static const char *title = "";
493 if (fmID >= 0 && fmID < s_iFilterMethodCount)
494 return (s_aszFilterMethodTitle [fmID]);
501 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
503 int fgID = FILTER_GENERATION_INVALID;
505 for (int i = 0; i < s_iFilterGenerationCount; i++)
506 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
515 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
517 static const char *name = "";
519 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
520 return (s_aszFilterGenerationName [fgID]);
526 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
528 static const char *name = "";
530 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
531 return (s_aszFilterGenerationTitle [fgID]);
537 ProcessSignal::filterSignal (const float constInput[], double output[]) const
539 double* input = new double [m_nSignalPoints];
541 for (i = 0; i < m_nSignalPoints; i++)
542 input[i] = constInput[i];
544 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
545 for (int i = 0; i < m_nSignalPoints; i++) {
546 int iDetFromCenter = i - (m_nSignalPoints / 2);
547 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
549 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
550 for (int i = 0; i < m_nSignalPoints; i++) {
551 int iDetFromCenter = i - (m_nSignalPoints / 2);
552 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
555 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
556 for (i = 0; i < m_nSignalPoints; i++)
557 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
558 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
559 double* inputSignal = new double [m_nFilterPoints];
560 for (i = 0; i < m_nSignalPoints; i++)
561 inputSignal[i] = input[i];
562 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
563 inputSignal[i] = 0; // zeropad
564 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
565 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
567 for (i = 0; i < m_nFilterPoints; i++)
568 fftSignal[i] *= m_adFilter[i];
569 double* inverseFourier = new double [m_nFilterPoints];
570 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
572 for (i = 0; i < m_nSignalPoints; i++)
573 output[i] = inverseFourier[i];
574 delete inverseFourier;
575 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
576 double* inputSignal = new double [m_nFilterPoints];
577 for (i = 0; i < m_nSignalPoints; i++)
578 inputSignal[i] = input[i];
579 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
580 inputSignal[i] = 0; // zeropad
581 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
582 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
584 for (i = 0; i < m_nFilterPoints; i++)
585 fftSignal[i] *= m_adFilter[i];
586 double* inverseFourier = new double [m_nFilterPoints];
587 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
589 for (i = 0; i < m_nSignalPoints; i++)
590 output[i] = inverseFourier[i];
591 delete inverseFourier;
594 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
595 for (i = 0; i < m_nSignalPoints; i++)
596 m_adRealFftInput[i] = input[i];
598 fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
599 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
600 for (i = 0; i < m_nFilterPoints; i++)
601 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
603 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
604 m_adRealFftSignal[i] = 0;
606 fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
607 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
608 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
609 output[i] = ifftOutput[i];
610 delete [] ifftOutput;
611 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
612 for (i = 0; i < m_nSignalPoints; i++)
613 m_adComplexFftInput[i].re = input[i];
615 fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
616 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
617 for (i = 0; i < m_nFilterPoints; i++) {
618 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
619 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
622 fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
623 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
624 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
625 output[i] = ifftOutput[i].re;
626 delete [] ifftOutput;
634 * convolve Discrete convolution of two functions
637 * r = convolve (f1, f2, dx, n, np, func_type)
638 * double r Convolved result
639 * double f1[], f2[] Functions to be convolved
640 * double dx Difference between successive x values
641 * int n Array index to center convolution about
642 * int np Number of points in f1 array
643 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
646 * f1 is the projection data, its indices range from 0 to np - 1.
647 * The index for f2, the filter, ranges from -(np-1) to (np-1).
648 * There are 3 ways to handle the negative vertices of f2:
649 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
650 * All filters used in reconstruction are even.
651 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
652 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
653 * for negative indices. Since f2 must range from -(np-1) to (np-1),
654 * if we add (np - 1) to f2's array index, then f2's index will
655 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
656 * stored at f2[np-1].
660 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
664 #if UNOPTIMIZED_CONVOLUTION
665 for (int i = 0; i < np; i++)
666 sum += func[i] * m_adFilter[n - i + (np - 1)];
668 double* f2 = m_adFilter + n + (np - 1);
669 for (int i = 0; i < np; i++)
670 sum += *func++ * *f2--;
678 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
682 #if UNOPTIMIZED_CONVOLUTION
683 for (int i = 0; i < np; i++)
684 sum += func[i] * m_adFilter[n - i + (np - 1)];
686 double* f2 = m_adFilter + n + (np - 1);
687 for (int i = 0; i < np; i++)
688 sum += *func++ * *f2--;
696 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
698 std::complex<double>* complexOutput = new std::complex<double> [n];
700 finiteFourierTransform (input, complexOutput, n, direction);
701 for (int i = 0; i < n; i++)
702 output[i] = complexOutput[i].real();
703 delete [] complexOutput;
707 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
714 double angleIncrement = direction * 2 * PI / n;
715 for (int i = 0; i < n; i++) {
718 for (int j = 0; j < n; j++) {
719 double angle = i * j * angleIncrement;
720 sumReal += input[j] * cos(angle);
721 sumImag += input[j] * sin(angle);
727 output[i] = std::complex<double> (sumReal, sumImag);
733 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
740 double angleIncrement = direction * 2 * PI / n;
741 for (int i = 0; i < n; i++) {
742 std::complex<double> sum (0,0);
743 for (int j = 0; j < n; j++) {
744 double angle = i * j * angleIncrement;
745 std::complex<double> exponentTerm (cos(angle), sin(angle));
746 sum += input[j] * exponentTerm;
756 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
763 double angleIncrement = direction * 2 * PI / n;
764 for (int i = 0; i < n; i++) {
766 for (int j = 0; j < n; j++) {
767 double angle = i * j * angleIncrement;
768 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
777 // Table-based routines
780 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
787 for (int i = 0; i < m_nFilterPoints; i++) {
788 double sumReal = 0, sumImag = 0;
789 for (int j = 0; j < m_nFilterPoints; j++) {
790 int tableIndex = i * j;
792 sumReal += input[j] * m_adFourierCosTable[tableIndex];
793 sumImag += input[j] * m_adFourierSinTable[tableIndex];
795 sumReal += input[j] * m_adFourierCosTable[tableIndex];
796 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
800 sumReal /= m_nFilterPoints;
801 sumImag /= m_nFilterPoints;
803 output[i] = std::complex<double> (sumReal, sumImag);
807 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
809 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
816 for (int i = 0; i < m_nFilterPoints; i++) {
817 double sumReal = 0, sumImag = 0;
818 for (int j = 0; j < m_nFilterPoints; j++) {
819 int tableIndex = i * j;
821 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
822 - input[j].imag() * m_adFourierSinTable[tableIndex];
823 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
824 + input[j].imag() * m_adFourierCosTable[tableIndex];
826 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
827 - input[j].imag() * -m_adFourierSinTable[tableIndex];
828 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
829 + input[j].imag() * m_adFourierCosTable[tableIndex];
833 sumReal /= m_nFilterPoints;
834 sumImag /= m_nFilterPoints;
836 output[i] = std::complex<double> (sumReal, sumImag);
841 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
848 for (int i = 0; i < m_nFilterPoints; i++) {
850 for (int j = 0; j < m_nFilterPoints; j++) {
851 int tableIndex = i * j;
853 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
854 - input[j].imag() * m_adFourierSinTable[tableIndex];
856 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
857 - input[j].imag() * -m_adFourierSinTable[tableIndex];
861 sumReal /= m_nFilterPoints;
869 ProcessSignal::addZeropadFactor (int n, int iZeropad)
872 double dLogBase2 = log(n) / log(2);
873 int iLogBase2 = static_cast<int>(floor (dLogBase2));
874 int iPaddedN = 1 << (iLogBase2 + iZeropad);
876 sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);