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];
550 #pragma omp parallel for
552 for (i = 0; i < m_nSignalPoints; i++)
553 input[i] = constInput[i];
555 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
557 #pragma omp parallel for
559 for (int i = 0; i < m_nSignalPoints; i++) {
560 int iDetFromCenter = i - (m_nSignalPoints / 2);
561 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
563 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
565 #pragma omp parallel for
567 for (int i = 0; i < m_nSignalPoints; i++) {
568 int iDetFromCenter = i - (m_nSignalPoints / 2);
569 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
572 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
574 #pragma omp parallel for
576 for (i = 0; i < m_nSignalPoints; i++)
577 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
578 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
579 double* inputSignal = new double [m_nFilterPoints];
580 for (i = 0; i < m_nSignalPoints; i++)
581 inputSignal[i] = input[i];
582 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
583 inputSignal[i] = 0; // zeropad
584 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
585 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
588 #pragma omp parallel for
590 for (i = 0; i < m_nFilterPoints; i++)
591 fftSignal[i] *= m_adFilter[i];
592 double* inverseFourier = new double [m_nFilterPoints];
593 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
595 for (i = 0; i < m_nSignalPoints; i++)
596 output[i] = inverseFourier[i];
597 delete inverseFourier;
598 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
599 double* inputSignal = new double [m_nFilterPoints];
600 for (i = 0; i < m_nSignalPoints; i++)
601 inputSignal[i] = input[i];
602 for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
603 inputSignal[i] = 0; // zeropad
604 std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
605 finiteFourierTransform (inputSignal, fftSignal, FORWARD);
608 #pragma omp parallel for
610 for (i = 0; i < m_nFilterPoints; i++)
611 fftSignal[i] *= m_adFilter[i];
612 double* inverseFourier = new double [m_nFilterPoints];
613 finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
615 for (i = 0; i < m_nSignalPoints; i++)
616 output[i] = inverseFourier[i];
617 delete inverseFourier;
620 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
621 for (i = 0; i < m_nSignalPoints; i++)
622 m_adRealFftInput[i] = input[i];
624 fftw_execute (m_realPlanForward);
625 for (i = 0; i < m_nFilterPoints; i++)
626 m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
627 for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
628 m_adRealFftSignal[i] = 0;
630 fftw_execute (m_realPlanBackward);
631 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
632 output[i] = m_adRealFftBackwardOutput[i];
633 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
634 for (i = 0; i < m_nSignalPoints; i++)
635 m_adComplexFftInput[i][0] = input[i];
637 fftw_execute (m_complexPlanForward);
639 #pragma omp parallel for
641 for (i = 0; i < m_nFilterPoints; i++) {
642 m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
643 m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
645 fftw_execute (m_complexPlanBackward);
646 for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
647 output[i] = m_adComplexFftBackwardOutput[i][0];
655 * convolve Discrete convolution of two functions
658 * r = convolve (f1, f2, dx, n, np, func_type)
659 * double r Convolved result
660 * double f1[], f2[] Functions to be convolved
661 * double dx Difference between successive x values
662 * int n Array index to center convolution about
663 * int np Number of points in f1 array
664 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
667 * f1 is the projection data, its indices range from 0 to np - 1.
668 * The index for f2, the filter, ranges from -(np-1) to (np-1).
669 * There are 3 ways to handle the negative vertices of f2:
670 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
671 * All filters used in reconstruction are even.
672 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
673 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
674 * for negative indices. Since f2 must range from -(np-1) to (np-1),
675 * if we add (np - 1) to f2's array index, then f2's index will
676 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
677 * stored at f2[np-1].
681 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
685 #if UNOPTIMIZED_CONVOLUTION
686 for (int i = 0; i < np; i++)
687 sum += func[i] * m_adFilter[n - i + (np - 1)];
689 double* f2 = m_adFilter + n + (np - 1);
690 for (int i = 0; i < np; i++)
691 sum += *func++ * *f2--;
699 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
703 #if UNOPTIMIZED_CONVOLUTION
704 for (int i = 0; i < np; i++)
705 sum += func[i] * m_adFilter[n - i + (np - 1)];
707 double* f2 = m_adFilter + n + (np - 1);
708 for (int i = 0; i < np; i++)
709 sum += *func++ * *f2--;
717 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
719 std::complex<double>* complexOutput = new std::complex<double> [n];
721 finiteFourierTransform (input, complexOutput, n, direction);
722 for (int i = 0; i < n; i++)
723 output[i] = complexOutput[i].real();
724 delete [] complexOutput;
728 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
735 double angleIncrement = direction * 2 * PI / n;
736 for (int i = 0; i < n; i++) {
739 for (int j = 0; j < n; j++) {
740 double angle = i * j * angleIncrement;
741 sumReal += input[j] * cos(angle);
742 sumImag += input[j] * sin(angle);
748 output[i] = std::complex<double> (sumReal, sumImag);
754 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
761 double angleIncrement = direction * 2 * PI / n;
762 for (int i = 0; i < n; i++) {
763 std::complex<double> sum (0,0);
764 for (int j = 0; j < n; j++) {
765 double angle = i * j * angleIncrement;
766 std::complex<double> exponentTerm (cos(angle), sin(angle));
767 sum += input[j] * exponentTerm;
777 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
784 double angleIncrement = direction * 2 * PI / n;
785 for (int i = 0; i < n; i++) {
787 for (int j = 0; j < n; j++) {
788 double angle = i * j * angleIncrement;
789 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
798 // Table-based routines
801 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
808 for (int i = 0; i < m_nFilterPoints; i++) {
809 double sumReal = 0, sumImag = 0;
810 for (int j = 0; j < m_nFilterPoints; j++) {
811 int tableIndex = i * j;
813 sumReal += input[j] * m_adFourierCosTable[tableIndex];
814 sumImag += input[j] * m_adFourierSinTable[tableIndex];
816 sumReal += input[j] * m_adFourierCosTable[tableIndex];
817 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
821 sumReal /= m_nFilterPoints;
822 sumImag /= m_nFilterPoints;
824 output[i] = std::complex<double> (sumReal, sumImag);
828 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
830 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
837 for (int i = 0; i < m_nFilterPoints; i++) {
838 double sumReal = 0, sumImag = 0;
839 for (int j = 0; j < m_nFilterPoints; j++) {
840 int tableIndex = i * j;
842 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
843 - input[j].imag() * m_adFourierSinTable[tableIndex];
844 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
845 + input[j].imag() * m_adFourierCosTable[tableIndex];
847 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
848 - input[j].imag() * -m_adFourierSinTable[tableIndex];
849 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
850 + input[j].imag() * m_adFourierCosTable[tableIndex];
854 sumReal /= m_nFilterPoints;
855 sumImag /= m_nFilterPoints;
857 output[i] = std::complex<double> (sumReal, sumImag);
862 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
869 for (int i = 0; i < m_nFilterPoints; i++) {
871 for (int j = 0; j < m_nFilterPoints; j++) {
872 int tableIndex = i * j;
874 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
875 - input[j].imag() * m_adFourierSinTable[tableIndex];
877 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
878 - input[j].imag() * -m_adFourierSinTable[tableIndex];
882 sumReal /= m_nFilterPoints;
890 ProcessSignal::addZeropadFactor (int n, int iZeropad)
893 double dLogBase2 = log(n) / log(2);
894 int iLogBase2 = static_cast<int>(floor (dLogBase2));
895 int iPaddedN = 1 << (iLogBase2 + iZeropad);
897 sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);