libctsim/procsignal.cpp

   1 /*****************************************************************************
   2 ** File IDENTIFICATION
   3 **
   4 **     Name:            procsignal.cpp
   5 **     Purpose:         Routines for processing signals and projections
   6 **     Progammer:           Kevin Rosenberg
   7 **     Date Started:    Aug 1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (c) 1983-2001 Kevin Rosenberg
  11 **
  12 **  $Id: procsignal.cpp,v 1.32 2001/03/30 19:17:32 kevin Exp $
  13 **
  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.
  17 **
  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.
  22 **
  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 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30 #ifdef HAVE_WXWINDOWS
  31 #include "nographics.h"
  32 #endif
  33
  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;
  40 #if HAVE_FFTW
  41 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
  42 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
  43 #endif
  44 const char* const ProcessSignal::s_aszFilterMethodName[] = {
  45   {"convolution"},
  46   {"fourier"},
  47   {"fouier-table"},
  48   {"fft"},
  49 #if HAVE_FFTW
  50   {"fftw"},
  51   {"rfftw"},
  52 #endif
  53 };
  54 const char* const ProcessSignal::s_aszFilterMethodTitle[] = {
  55   {"Convolution"},
  56   {"Fourier"},
  57   {"Fouier Trigometric Table"},
  58   {"FFT"},
  59 #if HAVE_FFTW
  60   {"FFTW"},
  61   {"Real/Half-Complex FFTW"},
  62 #endif
  63 };
  64 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
  65
  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[] = {
  71   {"direct"},
  72   {"inverse-fourier"},
  73 };
  74 const char* const ProcessSignal::s_aszFilterGenerationTitle[] = {
  75   {"Direct"},
  76   {"Inverse Fourier"},
  77 };
  78 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
  79
  80
  81 // CLASS IDENTIFICATION
  82 //   ProcessSignal
  83 //
  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)
  89 {
  90   m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
  91   if (m_idFilterMethod == FILTER_METHOD_INVALID) {
  92     m_fail = true;
  93     m_failMessage = "Invalid filter method name ";
  94     m_failMessage += szFilterMethodName;
  95     return;
  96   }
  97   m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
  98   if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
  99     m_fail = true;
 100     m_failMessage = "Invalid frequency filter name ";
 101     m_failMessage += szFilterGenerationName;
 102     return;
 103   }
 104   m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
 105   if (m_idFilter == SignalFilter::FILTER_INVALID) {
 106     m_fail = true;
 107     m_failMessage = "Invalid Filter name ";
 108     m_failMessage += szFilterName;
 109     return;
 110   }
 111   m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
 112   if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
 113     m_fail = true;
 114     m_failMessage = "Invalid domain name ";
 115     m_failMessage += szDomainName;
 116     return;
 117   }
 118
 119   init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
 120     m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength,
 121     dSourceDetectorLength, pSGP);
 122 }
 123
 124
 125 void
 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)
 130 {
 131   int i;
 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;
 139
 140   if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
 141     m_fail = true;
 142     return;
 143   }
 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;
 153
 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;
 160   }
 161
 162   if (m_idFilterMethod == FILTER_METHOD_FFT) {
 163 #if HAVE_FFTW
 164     m_idFilterMethod = FILTER_METHOD_RFFTW;
 165 #else
 166     m_fail = true;
 167     m_failMessage = "FFT not yet implemented";
 168     return;
 169 #endif
 170   }
 171
 172   bool m_bFrequencyFiltering = true;
 173   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
 174     m_bFrequencyFiltering = false;
 175
 176   // Spatial-based filtering
 177   if (! m_bFrequencyFiltering) {
 178
 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();
 202       }
 203 #endif
 204       Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
 205 #ifdef HAVE_SGP
 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();
 211       }
 212 #endif
 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();
 221       }
 222 #endif
 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();
 230       }
 231 #endif
 232       for (i = 0; i < m_nFilterPoints; i++) {
 233         m_adFilter[i] /= m_dSignalInc;
 234       }
 235     }
 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;
 245       }
 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();
 252       }
 253 #endif
 254     } // if (geometry)
 255   } // if (spatial filtering)
 256
 257   else if (m_bFrequencyFiltering) {  // Frequency-based filtering
 258
 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;
 263
 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);
 268       } else { // Even
 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;
 273       }
 274
 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);
 279
 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();
 286       }
 287 #endif
 288
 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;
 301         }
 302       }
 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();
 309       }
 310 #endif
 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();
 318       }
 319 #endif
 320
 321       // FILTERING:  FREQUENCY - INVERSE FOURIER
 322
 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))
 334           nextPowerOf2++;
 335         nextPowerOf2 += (m_iZeropad - 1);
 336         m_nFilterPoints = 1 << nextPowerOf2;
 337       }
 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);
 342 #endif
 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();
 353       }
 354 #endif
 355
 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)
 364             sinScale = 1;
 365           else
 366             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 367           double dScale = 0.5 * sinScale * sinScale;
 368           adSpatialFilter[i] *= dScale;
 369         }
 370       }
 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();
 377       }
 378 #endif
 379       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 380         adSpatialFilter[i] = 0;
 381
 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();
 395       }
 396 #endif
 397     }
 398   }
 399
 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 ];
 406     double angle = 0;
 407     for (i = 0; i < nFourier; i++) {
 408       m_adFourierCosTable[i] = cos (angle);
 409       m_adFourierSinTable[i] = sin (angle);
 410       angle += angleIncrement;
 411     }
 412   }
 413
 414 #if HAVE_FFTW
 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;
 418   }
 419
 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;
 436   }
 437 #endif
 438
 439 }
 440
 441 ProcessSignal::~ProcessSignal (void)
 442 {
 443   delete [] m_adFourierSinTable;
 444   delete [] m_adFourierCosTable;
 445   delete [] m_adFilter;
 446
 447 #if HAVE_FFTW
 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;
 453   }
 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;
 459   }
 460 #endif
 461 }
 462
 463 int
 464 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 465 {
 466   int fmID = FILTER_METHOD_INVALID;
 467
 468   for (int i = 0; i < s_iFilterMethodCount; i++)
 469     if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 470       fmID = i;
 471       break;
 472     }
 473
 474     return (fmID);
 475 }
 476
 477 const char *
 478 ProcessSignal::convertFilterMethodIDToName (const int fmID)
 479 {
 480   static const char *name = "";
 481
 482   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 483     return (s_aszFilterMethodName [fmID]);
 484
 485   return (name);
 486 }
 487
 488 const char *
 489 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 490 {
 491   static const char *title = "";
 492
 493   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 494     return (s_aszFilterMethodTitle [fmID]);
 495
 496   return (title);
 497 }
 498
 499
 500 int
 501 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 502 {
 503   int fgID = FILTER_GENERATION_INVALID;
 504
 505   for (int i = 0; i < s_iFilterGenerationCount; i++)
 506     if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 507       fgID = i;
 508       break;
 509     }
 510
 511     return (fgID);
 512 }
 513
 514 const char *
 515 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 516 {
 517   static const char *name = "";
 518
 519   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 520     return (s_aszFilterGenerationName [fgID]);
 521
 522   return (name);
 523 }
 524
 525 const char *
 526 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 527 {
 528   static const char *name = "";
 529
 530   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 531     return (s_aszFilterGenerationTitle [fgID]);
 532
 533   return (name);
 534 }
 535
 536 void
 537 ProcessSignal::filterSignal (const float constInput[], double output[]) const
 538 {
 539   double* input = new double [m_nSignalPoints];
 540   int i;
 541   for (i = 0; i < m_nSignalPoints; i++)
 542     input[i] = constInput[i];
 543
 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);
 548     }
 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);
 553     }
 554   }
 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);
 566     delete inputSignal;
 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);
 571     delete fftSignal;
 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);
 583     delete inputSignal;
 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);
 588     delete fftSignal;
 589     for (i = 0; i < m_nSignalPoints; i++)
 590       output[i] = inverseFourier[i];
 591     delete inverseFourier;
 592   }
 593 #if HAVE_FFTW
 594   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 595     for (i = 0; i < m_nSignalPoints; i++)
 596       m_adRealFftInput[i] = input[i];
 597
 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];
 602     delete [] fftOutput;
 603     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
 604              m_adRealFftSignal[i] = 0;
 605
 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];
 614
 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;
 620     }
 621     delete [] fftOutput;
 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;
 627   }
 628 #endif
 629   delete input;
 630 }
 631
 632
 633 /* NAME
 634 *    convolve                   Discrete convolution of two functions
 635 *
 636 * SYNOPSIS
 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
 644 *
 645 * NOTES
 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].
 657 */
 658
 659 double
 660 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 661 {
 662   double sum = 0.0;
 663
 664 #if UNOPTIMIZED_CONVOLUTION
 665   for (int i = 0; i < np; i++)
 666     sum += func[i] * m_adFilter[n - i + (np - 1)];
 667 #else
 668   double* f2 = m_adFilter + n + (np - 1);
 669   for (int i = 0; i < np; i++)
 670     sum += *func++ * *f2--;
 671 #endif
 672
 673   return (sum * dx);
 674 }
 675
 676
 677 double
 678 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 679 {
 680   double sum = 0.0;
 681
 682 #if UNOPTIMIZED_CONVOLUTION
 683   for (int i = 0; i < np; i++)
 684     sum += func[i] * m_adFilter[n - i + (np - 1)];
 685 #else
 686   double* f2 = m_adFilter + n + (np - 1);
 687   for (int i = 0; i < np; i++)
 688     sum += *func++ * *f2--;
 689 #endif
 690
 691   return (sum * dx);
 692 }
 693
 694
 695 void
 696 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 697 {
 698   std::complex<double>* complexOutput = new std::complex<double> [n];
 699
 700   finiteFourierTransform (input, complexOutput, n, direction);
 701   for (int i = 0; i < n; i++)
 702     output[i] = complexOutput[i].real();
 703   delete [] complexOutput;
 704 }
 705
 706 void
 707 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
 708 {
 709   if (direction < 0)
 710     direction = -1;
 711   else
 712     direction = 1;
 713
 714   double angleIncrement = direction * 2 * PI / n;
 715   for (int i = 0; i < n; i++) {
 716     double sumReal = 0;
 717     double sumImag = 0;
 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);
 722     }
 723     if (direction < 0) {
 724       sumReal /= n;
 725       sumImag /= n;
 726     }
 727     output[i] = std::complex<double> (sumReal, sumImag);
 728   }
 729 }
 730
 731
 732 void
 733 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
 734 {
 735   if (direction < 0)
 736     direction = -1;
 737   else
 738     direction = 1;
 739
 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;
 747     }
 748     if (direction < 0) {
 749       sum /= n;
 750     }
 751     output[i] = sum;
 752   }
 753 }
 754
 755 void
 756 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
 757 {
 758   if (direction < 0)
 759     direction = -1;
 760   else
 761     direction = 1;
 762
 763   double angleIncrement = direction * 2 * PI / n;
 764   for (int i = 0; i < n; i++) {
 765     double sumReal = 0;
 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);
 769     }
 770     if (direction < 0) {
 771       sumReal /= n;
 772     }
 773     output[i] = sumReal;
 774   }
 775 }
 776
 777 // Table-based routines
 778
 779 void
 780 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
 781 {
 782   if (direction < 0)
 783     direction = -1;
 784   else
 785     direction = 1;
 786
 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;
 791       if (direction > 0) {
 792         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 793         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 794       } else {
 795         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 796         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 797       }
 798     }
 799     if (direction < 0) {
 800       sumReal /= m_nFilterPoints;
 801       sumImag /= m_nFilterPoints;
 802     }
 803     output[i] = std::complex<double> (sumReal, sumImag);
 804   }
 805 }
 806
 807 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 808 void
 809 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
 810 {
 811   if (direction < 0)
 812     direction = -1;
 813   else
 814     direction = 1;
 815
 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;
 820       if (direction > 0) {
 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];
 825       } else {
 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];
 830       }
 831     }
 832     if (direction < 0) {
 833       sumReal /= m_nFilterPoints;
 834       sumImag /= m_nFilterPoints;
 835     }
 836     output[i] = std::complex<double> (sumReal, sumImag);
 837   }
 838 }
 839
 840 void
 841 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
 842 {
 843   if (direction < 0)
 844     direction = -1;
 845   else
 846     direction = 1;
 847
 848   for (int i = 0; i < m_nFilterPoints; i++) {
 849     double sumReal = 0;
 850     for (int j = 0; j < m_nFilterPoints; j++) {
 851       int tableIndex = i * j;
 852       if (direction > 0) {
 853         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 854           - input[j].imag() * m_adFourierSinTable[tableIndex];
 855       } else {
 856         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 857           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 858       }
 859     }
 860     if (direction < 0) {
 861       sumReal /= m_nFilterPoints;
 862     }
 863     output[i] = sumReal;
 864   }
 865 }
 866
 867
 868 int
 869 ProcessSignal::addZeropadFactor (int n, int iZeropad)
 870 {
 871   if (iZeropad > 0) {
 872     double dLogBase2 = log(n) / log(2);
 873     int iLogBase2 = static_cast<int>(floor (dLogBase2));
 874     int iPaddedN = 1 << (iLogBase2 + iZeropad);
 875 #ifdef DEBUG
 876     sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);
 877 #endif
 878     return iPaddedN;
 879   }
 880
 881   return n;
 882 }