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
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_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
422 m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
423 m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE);
424 m_adRealFftSignal = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
425 m_adRealFftBackwardOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
426 m_realPlanBackward = fftw_plan_r2r_1d (m_nOutputPoints, m_adRealFftSignal, m_adRealFftBackwardOutput, FFTW_HC2R, FFTW_ESTIMATE);
427 for (i = 0; i < m_nFilterPoints; i++)
428 m_adRealFftInput[i] = 0;
429 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
430 m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
431 m_adComplexFftOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
432 m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD, FFTW_ESTIMATE);
433 m_adComplexFftSignal = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
434 m_adComplexFftBackwardOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
435 m_complexPlanBackward = fftw_plan_dft_1d (m_nOutputPoints, m_adComplexFftSignal, m_adComplexFftBackwardOutput, FFTW_BACKWARD, FFTW_ESTIMATE);
437 for (i = 0; i < m_nFilterPoints; i++)
438 m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0;
439 for (i = 0; i < m_nOutputPoints; i++)
440 m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0;
446 ProcessSignal::~ProcessSignal (void)
448 delete [] m_adFourierSinTable;
449 delete [] m_adFourierCosTable;
450 delete [] m_adFilter;
453 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
454 fftw_destroy_plan(m_complexPlanForward);
455 fftw_destroy_plan(m_complexPlanBackward);
456 fftw_free (m_adComplexFftInput);
457 fftw_free (m_adComplexFftOutput);
458 fftw_free (m_adComplexFftSignal);
459 fftw_free (m_adComplexFftBackwardOutput);
461 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
462 fftw_destroy_plan(m_realPlanForward);
463 fftw_destroy_plan(m_realPlanBackward);
464 fftw_free (m_adRealFftInput);
465 fftw_free (m_adRealFftOutput);
466 fftw_free (m_adRealFftSignal);
467 fftw_free (m_adRealFftBackwardOutput);
473 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
475 int fmID = FILTER_METHOD_INVALID;
477 for (int i = 0; i < s_iFilterMethodCount; i++)
478 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
487 ProcessSignal::convertFilterMethodIDToName (const int fmID)
489 static const char *name = "";
491 if (fmID >= 0 && fmID < s_iFilterMethodCount)
492 return (s_aszFilterMethodName [fmID]);
498 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
500 static const char *title = "";
502 if (fmID >= 0 && fmID < s_iFilterMethodCount)
503 return (s_aszFilterMethodTitle [fmID]);
510 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
512 int fgID = FILTER_GENERATION_INVALID;
514 for (int i = 0; i < s_iFilterGenerationCount; i++)
515 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
524 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
526 static const char *name = "";
528 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
529 return (s_aszFilterGenerationName [fgID]);
535 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
537 static const char *name = "";
539 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
540 return (s_aszFilterGenerationTitle [fgID]);
546 ProcessSignal::filterSignal (const float constInput[], double output[]) const
548 double* input = new double [m_nSignalPoints];
550 for (i = 0; i < m_nSignalPoints; i++)
551 input[i] = constInput[i];
553 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
554 for (int i = 0; i < m_nSignalPoints; i++) {
555 int iDetFromCenter = i - (m_nSignalPoints / 2);
556 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
558 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
559 for (int i = 0; i < m_nSignalPoints; i++) {
560 int iDetFromCenter = i - (m_nSignalPoints / 2);
561 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
564 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
565 for (i = 0; i < m_nSignalPoints; i++)
566 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
567 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
568 double* inputSignal = new double [m_nFilterPoints];
569 for (i = 0; i < m_nSignalPoints; i++)
570 inputSignal[i] = input[i];
571 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
572 inputSignal[i] = 0; // zeropad
573 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
574 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
576 for (i = 0; i < m_nFilterPoints; i++)
577 fftSignal[i] *= m_adFilter[i];
578 double* inverseFourier = new double [m_nFilterPoints];
579 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
581 for (i = 0; i < m_nSignalPoints; i++)
582 output[i] = inverseFourier[i];
583 delete inverseFourier;
584 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
585 double* inputSignal = new double [m_nFilterPoints];
586 for (i = 0; i < m_nSignalPoints; i++)
587 inputSignal[i] = input[i];
588 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
589 inputSignal[i] = 0; // zeropad
590 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
591 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
593 for (i = 0; i < m_nFilterPoints; i++)
594 fftSignal[i] *= m_adFilter[i];
595 double* inverseFourier = new double [m_nFilterPoints];
596 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
598 for (i = 0; i < m_nSignalPoints; i++)
599 output[i] = inverseFourier[i];
600 delete inverseFourier;
603 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
604 for (i = 0; i < m_nSignalPoints; i++)
605 m_adRealFftInput[i] = input[i];
607 fftw_execute (m_realPlanForward);
608 for (i = 0; i < m_nFilterPoints; i++)
609 m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
610 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
611 m_adRealFftSignal[i] = 0;
613 fftw_execute (m_realPlanBackward);
614 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
615 output[i] = m_adRealFftBackwardOutput[i];
616 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
617 for (i = 0; i < m_nSignalPoints; i++)
618 m_adComplexFftInput[i][0] = input[i];
620 fftw_execute (m_complexPlanForward);
621 for (i = 0; i < m_nFilterPoints; i++) {
622 m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
623 m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
625 fftw_execute (m_complexPlanBackward);
626 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
627 output[i] = m_adComplexFftBackwardOutput[i][0];
635 * convolve Discrete convolution of two functions
638 * r = convolve (f1, f2, dx, n, np, func_type)
639 * double r Convolved result
640 * double f1[], f2[] Functions to be convolved
641 * double dx Difference between successive x values
642 * int n Array index to center convolution about
643 * int np Number of points in f1 array
644 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
647 * f1 is the projection data, its indices range from 0 to np - 1.
648 * The index for f2, the filter, ranges from -(np-1) to (np-1).
649 * There are 3 ways to handle the negative vertices of f2:
650 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
651 * All filters used in reconstruction are even.
652 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
653 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
654 * for negative indices. Since f2 must range from -(np-1) to (np-1),
655 * if we add (np - 1) to f2's array index, then f2's index will
656 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
657 * stored at f2[np-1].
661 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
665 #if UNOPTIMIZED_CONVOLUTION
666 for (int i = 0; i < np; i++)
667 sum += func[i] * m_adFilter[n - i + (np - 1)];
669 double* f2 = m_adFilter + n + (np - 1);
670 for (int i = 0; i < np; i++)
671 sum += *func++ * *f2--;
679 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
683 #if UNOPTIMIZED_CONVOLUTION
684 for (int i = 0; i < np; i++)
685 sum += func[i] * m_adFilter[n - i + (np - 1)];
687 double* f2 = m_adFilter + n + (np - 1);
688 for (int i = 0; i < np; i++)
689 sum += *func++ * *f2--;
697 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
699 std::complex<double>* complexOutput = new std::complex<double> [n];
701 finiteFourierTransform (input, complexOutput, n, direction);
702 for (int i = 0; i < n; i++)
703 output[i] = complexOutput[i].real();
704 delete [] complexOutput;
708 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
715 double angleIncrement = direction * 2 * PI / n;
716 for (int i = 0; i < n; i++) {
719 for (int j = 0; j < n; j++) {
720 double angle = i * j * angleIncrement;
721 sumReal += input[j] * cos(angle);
722 sumImag += input[j] * sin(angle);
728 output[i] = std::complex<double> (sumReal, sumImag);
734 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
741 double angleIncrement = direction * 2 * PI / n;
742 for (int i = 0; i < n; i++) {
743 std::complex<double> sum (0,0);
744 for (int j = 0; j < n; j++) {
745 double angle = i * j * angleIncrement;
746 std::complex<double> exponentTerm (cos(angle), sin(angle));
747 sum += input[j] * exponentTerm;
757 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
764 double angleIncrement = direction * 2 * PI / n;
765 for (int i = 0; i < n; i++) {
767 for (int j = 0; j < n; j++) {
768 double angle = i * j * angleIncrement;
769 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
778 // Table-based routines
781 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
788 for (int i = 0; i < m_nFilterPoints; i++) {
789 double sumReal = 0, sumImag = 0;
790 for (int j = 0; j < m_nFilterPoints; j++) {
791 int tableIndex = i * j;
793 sumReal += input[j] * m_adFourierCosTable[tableIndex];
794 sumImag += input[j] * m_adFourierSinTable[tableIndex];
796 sumReal += input[j] * m_adFourierCosTable[tableIndex];
797 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
801 sumReal /= m_nFilterPoints;
802 sumImag /= m_nFilterPoints;
804 output[i] = std::complex<double> (sumReal, sumImag);
808 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
810 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
817 for (int i = 0; i < m_nFilterPoints; i++) {
818 double sumReal = 0, sumImag = 0;
819 for (int j = 0; j < m_nFilterPoints; j++) {
820 int tableIndex = i * j;
822 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
823 - input[j].imag() * m_adFourierSinTable[tableIndex];
824 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
825 + input[j].imag() * m_adFourierCosTable[tableIndex];
827 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
828 - input[j].imag() * -m_adFourierSinTable[tableIndex];
829 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
830 + input[j].imag() * m_adFourierCosTable[tableIndex];
834 sumReal /= m_nFilterPoints;
835 sumImag /= m_nFilterPoints;
837 output[i] = std::complex<double> (sumReal, sumImag);
842 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
849 for (int i = 0; i < m_nFilterPoints; i++) {
851 for (int j = 0; j < m_nFilterPoints; j++) {
852 int tableIndex = i * j;
854 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
855 - input[j].imag() * m_adFourierSinTable[tableIndex];
857 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
858 - input[j].imag() * -m_adFourierSinTable[tableIndex];
862 sumReal /= m_nFilterPoints;
870 ProcessSignal::addZeropadFactor (int n, int iZeropad)
873 double dLogBase2 = log(n) / log(2);
874 int iLogBase2 = static_cast<int>(floor (dLogBase2));
875 int iPaddedN = 1 << (iLogBase2 + iZeropad);
877 sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);