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-2009 Kevin Rosenberg
12 ** This program is free software; you can redistribute it and/or modify
13 ** it under the terms of the GNU General Public License (version 2) as
14 ** published by the Free Software Foundation.
16 ** This program is distributed in the hope that it will be useful,
17 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ** GNU General Public License for more details.
21 ** You should have received a copy of the GNU General Public License
22 ** along with this program; if not, write to the Free Software
23 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 ******************************************************************************/
29 #include "nographics.h"
32 // FilterMethod ID/Names
33 const int ProcessSignal::FILTER_METHOD_INVALID = -1;
34 const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
35 const int ProcessSignal::FILTER_METHOD_FOURIER = 1;
36 const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2;
37 const int ProcessSignal::FILTER_METHOD_FFT = 3;
39 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
40 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
42 const char* const ProcessSignal::s_aszFilterMethodName[] = {
52 const char* const ProcessSignal::s_aszFilterMethodTitle[] = {
55 "Fouier Trigometric Table",
59 "Real/Half-Complex FFTW",
62 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
64 // FilterGeneration ID/Names
65 const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
66 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
67 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
68 const char* const ProcessSignal::s_aszFilterGenerationName[] = {
72 const char* const ProcessSignal::s_aszFilterGenerationTitle[] = {
76 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
79 // CLASS IDENTIFICATION
82 ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth,
83 double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName,
84 const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel,
85 int iGeometry, double dFocalLength, double dSourceDetectorLength, SGP* pSGP)
86 : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
88 m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
89 if (m_idFilterMethod == FILTER_METHOD_INVALID) {
91 m_failMessage = "Invalid filter method name ";
92 m_failMessage += szFilterMethodName;
95 m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
96 if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
98 m_failMessage = "Invalid frequency filter name ";
99 m_failMessage += szFilterGenerationName;
102 m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
103 if (m_idFilter == SignalFilter::FILTER_INVALID) {
105 m_failMessage = "Invalid Filter name ";
106 m_failMessage += szFilterName;
109 m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
110 if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
112 m_failMessage = "Invalid domain name ";
113 m_failMessage += szDomainName;
117 init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
118 m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength,
119 dSourceDetectorLength, pSGP);
124 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement,
125 int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration,
126 const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry,
127 double dFocalLength, double dSourceDetectorLength, SGP* pSGP)
130 m_idFilter = idFilter;
131 m_idDomain = idDomain;
132 m_idFilterMethod = idFilterMethod;
133 m_idFilterGeneration = idFilterGeneration;
134 m_idGeometry = iGeometry;
135 m_dFocalLength = dFocalLength;
136 m_dSourceDetectorLength = dSourceDetectorLength;
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 adjust for imaginary detector through origin of phantom
153 // see Kak-Slaney Fig 3.22, for Collinear diagram
154 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
155 double dEquilinearScale = m_dSourceDetectorLength / m_dFocalLength;
156 m_dSignalInc /= dEquilinearScale;
157 m_dBandwidth *= dEquilinearScale;
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 = addZeropadFactor (m_nSignalPoints, m_iZeropad);
260 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
262 if (isOdd (m_nFilterPoints)) { // Odd
263 m_dFilterMin = -1. / (2 * m_dSignalInc);
264 m_dFilterMax = 1. / (2 * m_dSignalInc);
265 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
267 m_dFilterMin = -1. / (2 * m_dSignalInc);
268 m_dFilterMax = 1. / (2 * m_dSignalInc);
269 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
270 m_dFilterMax -= m_dFilterInc;
273 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
274 m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
275 m_adFilter = new double [m_nFilterPoints];
276 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
278 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
279 if (g_bRunningWXWindows && m_traceLevel > 0) {
280 EZPlotDialog dlgEZPlot;
281 dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
282 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
283 dlgEZPlot.ShowModal();
287 // This works fairly well. I'm not sure why since scaling for geometries is done on
288 // frequency filter rather than spatial filter as it should be.
289 // It gives values slightly off than freq/inverse filtering
290 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
291 for (i = 0; i < m_nFilterPoints; i++)
292 m_adFilter[i] *= 0.5;
293 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
294 for (i = 0; i < m_nFilterPoints; i++) {
295 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
296 double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
297 double dScale = 0.5 * sinScale * sinScale;
298 m_adFilter[i] *= dScale;
301 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
302 if (g_bRunningWXWindows && m_traceLevel > 0) {
303 EZPlotDialog dlgEZPlot;
304 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
305 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
306 dlgEZPlot.ShowModal();
309 Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
310 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
311 if (g_bRunningWXWindows && m_traceLevel > 0) {
312 EZPlotDialog dlgEZPlot;
313 dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
314 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
315 dlgEZPlot.ShowModal();
319 // FILTERING: FREQUENCY - INVERSE FOURIER
321 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
322 // calculate number of filter points with zeropadding
323 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
324 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
325 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
326 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
327 m_nFilterPoints = nSpatialPoints;
328 if (m_iZeropad > 0) {
329 double logBase2 = log(nSpatialPoints) / log(2);
330 int nextPowerOf2 = static_cast<int>(floor(logBase2));
331 if (logBase2 != floor(logBase2))
333 nextPowerOf2 += (m_iZeropad - 1);
334 m_nFilterPoints = 1 << nextPowerOf2;
336 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
337 #if defined(DEBUG) || defined(_DEBUG)
338 if (m_traceLevel >= Trace::TRACE_CONSOLE)
339 sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
341 double* adSpatialFilter = new double [m_nFilterPoints];
342 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
343 m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
344 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
345 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
346 if (g_bRunningWXWindows && m_traceLevel > 0) {
347 EZPlotDialog dlgEZPlot;;
348 dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
349 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
350 dlgEZPlot.ShowModal();
354 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
355 for (i = 0; i < nSpatialPoints; i++)
356 adSpatialFilter[i] *= 0.5;
357 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
358 for (i = 0; i < nSpatialPoints; i++) {
359 int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
360 double sinScale = sin (iDetFromZero * m_dSignalInc);
361 if (fabs(sinScale) < 1E-7)
364 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
365 double dScale = 0.5 * sinScale * sinScale;
366 adSpatialFilter[i] *= dScale;
369 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
370 if (g_bRunningWXWindows && m_traceLevel > 0) {
371 EZPlotDialog dlgEZPlot;;
372 dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
373 dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
374 dlgEZPlot.ShowModal();
377 for (i = nSpatialPoints; i < m_nFilterPoints; i++)
378 adSpatialFilter[i] = 0;
380 m_adFilter = new double [m_nFilterPoints];
381 std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
382 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
383 delete adSpatialFilter;
384 for (i = 0; i < m_nFilterPoints; i++)
385 m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
386 delete acInverseFilter;
387 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
388 if (g_bRunningWXWindows && m_traceLevel > 0) {
389 EZPlotDialog dlgEZPlot;
390 dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
391 dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
392 dlgEZPlot.ShowModal();
398 // precalculate sin and cosine tables for fourier transform
399 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
400 int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
401 double angleIncrement = (2. * PI) / m_nFilterPoints;
402 m_adFourierCosTable = new double[ nFourier ];
403 m_adFourierSinTable = new double[ nFourier ];
405 for (i = 0; i < nFourier; i++) {
406 m_adFourierCosTable[i] = cos (angle);
407 m_adFourierSinTable[i] = sin (angle);
408 angle += angleIncrement;
413 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
414 for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
415 m_adFilter[i] /= m_nFilterPoints;
418 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
419 m_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
420 m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
421 m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE);
422 m_adRealFftSignal = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
423 m_adRealFftBackwardOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
424 m_realPlanBackward = fftw_plan_r2r_1d (m_nOutputPoints, m_adRealFftSignal, m_adRealFftBackwardOutput, FFTW_HC2R, FFTW_ESTIMATE);
425 for (i = 0; i < m_nFilterPoints; i++)
426 m_adRealFftInput[i] = 0;
427 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
428 m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
429 m_adComplexFftOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
430 m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD, FFTW_ESTIMATE);
431 m_adComplexFftSignal = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
432 m_adComplexFftBackwardOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
433 m_complexPlanBackward = fftw_plan_dft_1d (m_nOutputPoints, m_adComplexFftSignal, m_adComplexFftBackwardOutput, FFTW_BACKWARD, FFTW_ESTIMATE);
435 for (i = 0; i < m_nFilterPoints; i++)
436 m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0;
437 for (i = 0; i < m_nOutputPoints; i++)
438 m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0;
444 ProcessSignal::~ProcessSignal (void)
446 delete [] m_adFourierSinTable;
447 delete [] m_adFourierCosTable;
448 delete [] m_adFilter;
451 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
452 fftw_destroy_plan(m_complexPlanForward);
453 fftw_destroy_plan(m_complexPlanBackward);
454 fftw_free (m_adComplexFftInput);
455 fftw_free (m_adComplexFftOutput);
456 fftw_free (m_adComplexFftSignal);
457 fftw_free (m_adComplexFftBackwardOutput);
459 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
460 fftw_destroy_plan(m_realPlanForward);
461 fftw_destroy_plan(m_realPlanBackward);
462 fftw_free (m_adRealFftInput);
463 fftw_free (m_adRealFftOutput);
464 fftw_free (m_adRealFftSignal);
465 fftw_free (m_adRealFftBackwardOutput);
471 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
473 int fmID = FILTER_METHOD_INVALID;
475 for (int i = 0; i < s_iFilterMethodCount; i++)
476 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
485 ProcessSignal::convertFilterMethodIDToName (const int fmID)
487 static const char *name = "";
489 if (fmID >= 0 && fmID < s_iFilterMethodCount)
490 return (s_aszFilterMethodName [fmID]);
496 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
498 static const char *title = "";
500 if (fmID >= 0 && fmID < s_iFilterMethodCount)
501 return (s_aszFilterMethodTitle [fmID]);
508 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
510 int fgID = FILTER_GENERATION_INVALID;
512 for (int i = 0; i < s_iFilterGenerationCount; i++)
513 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
522 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
524 static const char *name = "";
526 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
527 return (s_aszFilterGenerationName [fgID]);
533 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
535 static const char *name = "";
537 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
538 return (s_aszFilterGenerationTitle [fgID]);
544 ProcessSignal::filterSignal (const float constInput[], double output[]) const
546 double* input = new double [m_nSignalPoints];
548 for (i = 0; i < m_nSignalPoints; i++)
549 input[i] = constInput[i];
551 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
552 for (int i = 0; i < m_nSignalPoints; i++) {
553 int iDetFromCenter = i - (m_nSignalPoints / 2);
554 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
556 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
557 for (int i = 0; i < m_nSignalPoints; i++) {
558 int iDetFromCenter = i - (m_nSignalPoints / 2);
559 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
562 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
563 for (i = 0; i < m_nSignalPoints; i++)
564 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
565 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
566 double* inputSignal = new double [m_nFilterPoints];
567 for (i = 0; i < m_nSignalPoints; i++)
568 inputSignal[i] = input[i];
569 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
570 inputSignal[i] = 0; // zeropad
571 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
572 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
574 for (i = 0; i < m_nFilterPoints; i++)
575 fftSignal[i] *= m_adFilter[i];
576 double* inverseFourier = new double [m_nFilterPoints];
577 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
579 for (i = 0; i < m_nSignalPoints; i++)
580 output[i] = inverseFourier[i];
581 delete inverseFourier;
582 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
583 double* inputSignal = new double [m_nFilterPoints];
584 for (i = 0; i < m_nSignalPoints; i++)
585 inputSignal[i] = input[i];
586 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
587 inputSignal[i] = 0; // zeropad
588 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
589 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
591 for (i = 0; i < m_nFilterPoints; i++)
592 fftSignal[i] *= m_adFilter[i];
593 double* inverseFourier = new double [m_nFilterPoints];
594 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
596 for (i = 0; i < m_nSignalPoints; i++)
597 output[i] = inverseFourier[i];
598 delete inverseFourier;
601 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
602 for (i = 0; i < m_nSignalPoints; i++)
603 m_adRealFftInput[i] = input[i];
605 fftw_execute (m_realPlanForward);
606 for (i = 0; i < m_nFilterPoints; i++)
607 m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
608 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
609 m_adRealFftSignal[i] = 0;
611 fftw_execute (m_realPlanBackward);
612 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
613 output[i] = m_adRealFftBackwardOutput[i];
614 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
615 for (i = 0; i < m_nSignalPoints; i++)
616 m_adComplexFftInput[i][0] = input[i];
618 fftw_execute (m_complexPlanForward);
619 for (i = 0; i < m_nFilterPoints; i++) {
620 m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
621 m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
623 fftw_execute (m_complexPlanBackward);
624 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
625 output[i] = m_adComplexFftBackwardOutput[i][0];
633 * convolve Discrete convolution of two functions
636 * r = convolve (f1, f2, dx, n, np, func_type)
637 * double r Convolved result
638 * double f1[], f2[] Functions to be convolved
639 * double dx Difference between successive x values
640 * int n Array index to center convolution about
641 * int np Number of points in f1 array
642 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
645 * f1 is the projection data, its indices range from 0 to np - 1.
646 * The index for f2, the filter, ranges from -(np-1) to (np-1).
647 * There are 3 ways to handle the negative vertices of f2:
648 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
649 * All filters used in reconstruction are even.
650 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
651 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
652 * for negative indices. Since f2 must range from -(np-1) to (np-1),
653 * if we add (np - 1) to f2's array index, then f2's index will
654 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
655 * stored at f2[np-1].
659 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
663 #if UNOPTIMIZED_CONVOLUTION
664 for (int i = 0; i < np; i++)
665 sum += func[i] * m_adFilter[n - i + (np - 1)];
667 double* f2 = m_adFilter + n + (np - 1);
668 for (int i = 0; i < np; i++)
669 sum += *func++ * *f2--;
677 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
681 #if UNOPTIMIZED_CONVOLUTION
682 for (int i = 0; i < np; i++)
683 sum += func[i] * m_adFilter[n - i + (np - 1)];
685 double* f2 = m_adFilter + n + (np - 1);
686 for (int i = 0; i < np; i++)
687 sum += *func++ * *f2--;
695 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
697 std::complex<double>* complexOutput = new std::complex<double> [n];
699 finiteFourierTransform (input, complexOutput, n, direction);
700 for (int i = 0; i < n; i++)
701 output[i] = complexOutput[i].real();
702 delete [] complexOutput;
706 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
713 double angleIncrement = direction * 2 * PI / n;
714 for (int i = 0; i < n; i++) {
717 for (int j = 0; j < n; j++) {
718 double angle = i * j * angleIncrement;
719 sumReal += input[j] * cos(angle);
720 sumImag += input[j] * sin(angle);
726 output[i] = std::complex<double> (sumReal, sumImag);
732 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
739 double angleIncrement = direction * 2 * PI / n;
740 for (int i = 0; i < n; i++) {
741 std::complex<double> sum (0,0);
742 for (int j = 0; j < n; j++) {
743 double angle = i * j * angleIncrement;
744 std::complex<double> exponentTerm (cos(angle), sin(angle));
745 sum += input[j] * exponentTerm;
755 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
762 double angleIncrement = direction * 2 * PI / n;
763 for (int i = 0; i < n; i++) {
765 for (int j = 0; j < n; j++) {
766 double angle = i * j * angleIncrement;
767 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
776 // Table-based routines
779 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
786 for (int i = 0; i < m_nFilterPoints; i++) {
787 double sumReal = 0, sumImag = 0;
788 for (int j = 0; j < m_nFilterPoints; j++) {
789 int tableIndex = i * j;
791 sumReal += input[j] * m_adFourierCosTable[tableIndex];
792 sumImag += input[j] * m_adFourierSinTable[tableIndex];
794 sumReal += input[j] * m_adFourierCosTable[tableIndex];
795 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
799 sumReal /= m_nFilterPoints;
800 sumImag /= m_nFilterPoints;
802 output[i] = std::complex<double> (sumReal, sumImag);
806 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
808 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
815 for (int i = 0; i < m_nFilterPoints; i++) {
816 double sumReal = 0, sumImag = 0;
817 for (int j = 0; j < m_nFilterPoints; j++) {
818 int tableIndex = i * j;
820 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
821 - input[j].imag() * m_adFourierSinTable[tableIndex];
822 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
823 + input[j].imag() * m_adFourierCosTable[tableIndex];
825 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
826 - input[j].imag() * -m_adFourierSinTable[tableIndex];
827 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
828 + input[j].imag() * m_adFourierCosTable[tableIndex];
832 sumReal /= m_nFilterPoints;
833 sumImag /= m_nFilterPoints;
835 output[i] = std::complex<double> (sumReal, sumImag);
840 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
847 for (int i = 0; i < m_nFilterPoints; i++) {
849 for (int j = 0; j < m_nFilterPoints; j++) {
850 int tableIndex = i * j;
852 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
853 - input[j].imag() * m_adFourierSinTable[tableIndex];
855 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
856 - input[j].imag() * -m_adFourierSinTable[tableIndex];
860 sumReal /= m_nFilterPoints;
868 ProcessSignal::addZeropadFactor (int n, int iZeropad)
871 double dLogBase2 = log(n) / log(2);
872 int iLogBase2 = static_cast<int>(floor (dLogBase2));
873 int iPaddedN = 1 << (iLogBase2 + iZeropad);
875 sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);