libctsim/procsignal.cpp

   1 /*****************************************************************************
   2 ** File IDENTIFICATION
   3 **
   4 **     Name:                   filter.cpp
   5 **     Purpose:                Routines for signal-procesing filters
   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.27 2001/03/01 07:30:49 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 "dlgezplot.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 = m_nSignalPoints;
 262       if (m_iZeropad > 0) {
 263         double logBase2 = log(m_nFilterPoints) / log(2);
 264         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 265         if (logBase2 != floor(logBase2))
 266           nextPowerOf2++;
 267         nextPowerOf2 += (m_iZeropad - 1);
 268         m_nFilterPoints = 1 << nextPowerOf2;
 269 #if defined(DEBUG) || defined(_DEBUG)
 270         if (m_traceLevel >= Trace::TRACE_CONSOLE)
 271           sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
 272 #endif
 273       }
 274       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 275
 276       if (m_nFilterPoints % 2) { // Odd
 277         m_dFilterMin = -1. / (2 * m_dSignalInc);
 278         m_dFilterMax = 1. / (2 * m_dSignalInc);
 279         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 280       } else { // Even
 281         m_dFilterMin = -1. / (2 * m_dSignalInc);
 282         m_dFilterMax = 1. / (2 * m_dSignalInc);
 283         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
 284         m_dFilterMax -= m_dFilterInc;
 285       }
 286
 287       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
 288         m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 289       m_adFilter = new double [m_nFilterPoints];
 290       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 291
 292 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 293       if (g_bRunningWXWindows && m_traceLevel > 0) {
 294         EZPlotDialog dlgEZPlot;
 295         dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
 296         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 297         dlgEZPlot.ShowModal();
 298       }
 299 #endif
 300
 301       // This works fairly well. I'm not sure why since scaling for geometries is done on
 302       // frequency filter rather than spatial filter as it should be.
 303       // It gives values slightly off than freq/inverse filtering
 304       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 305         for (i = 0; i < m_nFilterPoints; i++)
 306           m_adFilter[i] *= 0.5;
 307       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 308         for (i = 0; i < m_nFilterPoints; i++) {
 309           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 310           double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
 311           double dScale = 0.5 * sinScale * sinScale;
 312           m_adFilter[i] *= dScale;
 313         }
 314       }
 315 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 316       if (g_bRunningWXWindows && m_traceLevel > 0) {
 317         EZPlotDialog dlgEZPlot;
 318         dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
 319         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 320         dlgEZPlot.ShowModal();
 321       }
 322 #endif
 323       Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 324 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 325       if (g_bRunningWXWindows && m_traceLevel > 0) {
 326         EZPlotDialog dlgEZPlot;
 327         dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
 328         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 329         dlgEZPlot.ShowModal();
 330       }
 331 #endif
 332
 333       // FILTERING:  FREQUENCY - INVERSE FOURIER
 334
 335     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 336       // calculate number of filter points with zeropadding
 337       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
 338       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 339       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 340       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
 341       m_nFilterPoints = nSpatialPoints;
 342       if (m_iZeropad > 0) {
 343         double logBase2 = log(nSpatialPoints) / log(2);
 344         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 345         if (logBase2 != floor(logBase2))
 346           nextPowerOf2++;
 347         nextPowerOf2 += (m_iZeropad - 1);
 348         m_nFilterPoints = 1 << nextPowerOf2;
 349       }
 350       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 351 #if defined(DEBUG) || defined(_DEBUG)
 352       if (m_traceLevel >= Trace::TRACE_CONSOLE)
 353         sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
 354 #endif
 355       double* adSpatialFilter = new double [m_nFilterPoints];
 356       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
 357         m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 358       filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
 359 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 360       if (g_bRunningWXWindows && m_traceLevel > 0) {
 361         EZPlotDialog dlgEZPlot;;
 362         dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
 363         dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
 364         dlgEZPlot.ShowModal();
 365       }
 366 #endif
 367
 368       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 369         for (i = 0; i < nSpatialPoints; i++)
 370           adSpatialFilter[i] *= 0.5;
 371       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 372         for (i = 0; i < nSpatialPoints; i++) {
 373           int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
 374           double sinScale = sin (iDetFromZero * m_dSignalInc);
 375           if (fabs(sinScale) < 1E-7)
 376             sinScale = 1;
 377           else
 378             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 379           double dScale = 0.5 * sinScale * sinScale;
 380           adSpatialFilter[i] *= dScale;
 381         }
 382       }
 383 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 384       if (g_bRunningWXWindows && m_traceLevel > 0) {
 385         EZPlotDialog dlgEZPlot;;
 386         dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
 387         dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
 388         dlgEZPlot.ShowModal();
 389       }
 390 #endif
 391       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 392         adSpatialFilter[i] = 0;
 393
 394       m_adFilter = new double [m_nFilterPoints];
 395       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
 396       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
 397       delete adSpatialFilter;
 398       for (i = 0; i < m_nFilterPoints; i++)
 399         m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
 400       delete acInverseFilter;
 401 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 402       if (g_bRunningWXWindows && m_traceLevel > 0) {
 403         EZPlotDialog dlgEZPlot;
 404         dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
 405         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 406         dlgEZPlot.ShowModal();
 407       }
 408 #endif
 409     }
 410   }
 411
 412   // precalculate sin and cosine tables for fourier transform
 413   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 414     int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
 415     double angleIncrement = (2. * PI) / m_nFilterPoints;
 416     m_adFourierCosTable = new double[ nFourier ];
 417     m_adFourierSinTable = new double[ nFourier ];
 418     double angle = 0;
 419     for (i = 0; i < nFourier; i++) {
 420       m_adFourierCosTable[i] = cos (angle);
 421       m_adFourierSinTable[i] = sin (angle);
 422       angle += angleIncrement;
 423     }
 424   }
 425
 426 #if HAVE_FFTW
 427   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 428     for (i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 429       m_adFilter[i] /= m_nFilterPoints;
 430   }
 431
 432   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 433     m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 434     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 435     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
 436     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
 437     for (i = 0; i < m_nFilterPoints; i++)
 438       m_adRealFftInput[i] = 0;
 439   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 440     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 441     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 442     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 443     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 444     for (i = 0; i < m_nFilterPoints; i++)
 445       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
 446     for (i = 0; i < m_nOutputPoints; i++)
 447       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
 448   }
 449 #endif
 450
 451 }
 452
 453 ProcessSignal::~ProcessSignal (void)
 454 {
 455   delete [] m_adFourierSinTable;
 456   delete [] m_adFourierCosTable;
 457   delete [] m_adFilter;
 458
 459 #if HAVE_FFTW
 460   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 461     fftw_destroy_plan(m_complexPlanForward);
 462     fftw_destroy_plan(m_complexPlanBackward);
 463     delete [] m_adComplexFftInput;
 464     delete [] m_adComplexFftSignal;
 465   }
 466   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 467     rfftw_destroy_plan(m_realPlanForward);
 468     rfftw_destroy_plan(m_realPlanBackward);
 469     delete [] m_adRealFftInput;
 470     delete [] m_adRealFftSignal;
 471   }
 472 #endif
 473 }
 474
 475 int
 476 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 477 {
 478   int fmID = FILTER_METHOD_INVALID;
 479
 480   for (int i = 0; i < s_iFilterMethodCount; i++)
 481     if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 482       fmID = i;
 483       break;
 484     }
 485
 486     return (fmID);
 487 }
 488
 489 const char *
 490 ProcessSignal::convertFilterMethodIDToName (const int fmID)
 491 {
 492   static const char *name = "";
 493
 494   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 495     return (s_aszFilterMethodName [fmID]);
 496
 497   return (name);
 498 }
 499
 500 const char *
 501 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 502 {
 503   static const char *title = "";
 504
 505   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 506     return (s_aszFilterMethodTitle [fmID]);
 507
 508   return (title);
 509 }
 510
 511
 512 int
 513 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 514 {
 515   int fgID = FILTER_GENERATION_INVALID;
 516
 517   for (int i = 0; i < s_iFilterGenerationCount; i++)
 518     if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 519       fgID = i;
 520       break;
 521     }
 522
 523     return (fgID);
 524 }
 525
 526 const char *
 527 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 528 {
 529   static const char *name = "";
 530
 531   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 532     return (s_aszFilterGenerationName [fgID]);
 533
 534   return (name);
 535 }
 536
 537 const char *
 538 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 539 {
 540   static const char *name = "";
 541
 542   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 543     return (s_aszFilterGenerationTitle [fgID]);
 544
 545   return (name);
 546 }
 547
 548 void
 549 ProcessSignal::filterSignal (const float constInput[], double output[]) const
 550 {
 551   double* input = new double [m_nSignalPoints];
 552   int i;
 553   for (i = 0; i < m_nSignalPoints; i++)
 554     input[i] = constInput[i];
 555
 556   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 557     for (int i = 0; i < m_nSignalPoints; i++) {
 558       int iDetFromCenter = i - (m_nSignalPoints / 2);
 559       input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
 560     }
 561   } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 562     for (int i = 0; i < m_nSignalPoints; i++) {
 563       int iDetFromCenter = i - (m_nSignalPoints / 2);
 564       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
 565     }
 566   }
 567   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 568     for (i = 0; i < m_nSignalPoints; i++)
 569       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 570   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 571     double* inputSignal = new double [m_nFilterPoints];
 572     for (i = 0; i < m_nSignalPoints; i++)
 573       inputSignal[i] = input[i];
 574     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 575       inputSignal[i] = 0;  // zeropad
 576     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 577     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
 578     delete inputSignal;
 579     for (i = 0; i < m_nFilterPoints; i++)
 580       fftSignal[i] *= m_adFilter[i];
 581     double* inverseFourier = new double [m_nFilterPoints];
 582     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
 583     delete fftSignal;
 584     for (i = 0; i < m_nSignalPoints; i++)
 585       output[i] = inverseFourier[i];
 586     delete inverseFourier;
 587   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 588     double* inputSignal = new double [m_nFilterPoints];
 589     for (i = 0; i < m_nSignalPoints; i++)
 590       inputSignal[i] = input[i];
 591     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 592       inputSignal[i] = 0;  // zeropad
 593     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 594     finiteFourierTransform (inputSignal, fftSignal, FORWARD);
 595     delete inputSignal;
 596     for (i = 0; i < m_nFilterPoints; i++)
 597       fftSignal[i] *= m_adFilter[i];
 598     double* inverseFourier = new double [m_nFilterPoints];
 599     finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
 600     delete fftSignal;
 601     for (i = 0; i < m_nSignalPoints; i++)
 602       output[i] = inverseFourier[i];
 603     delete inverseFourier;
 604   }
 605 #if HAVE_FFTW
 606   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 607     for (i = 0; i < m_nSignalPoints; i++)
 608       m_adRealFftInput[i] = input[i];
 609
 610     fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
 611     rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
 612     for (i = 0; i < m_nFilterPoints; i++)
 613       m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
 614     delete [] fftOutput;
 615     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
 616              m_adRealFftSignal[i] = 0;
 617
 618     fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
 619     rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
 620     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 621       output[i] = ifftOutput[i];
 622     delete [] ifftOutput;
 623   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 624     for (i = 0; i < m_nSignalPoints; i++)
 625       m_adComplexFftInput[i].re = input[i];
 626
 627     fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
 628     fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
 629     for (i = 0; i < m_nFilterPoints; i++) {
 630       m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
 631       m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
 632     }
 633     delete [] fftOutput;
 634     fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
 635     fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
 636     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 637       output[i] = ifftOutput[i].re;
 638     delete [] ifftOutput;
 639   }
 640 #endif
 641   delete input;
 642 }
 643
 644
 645 /* NAME
 646 *    convolve                   Discrete convolution of two functions
 647 *
 648 * SYNOPSIS
 649 *    r = convolve (f1, f2, dx, n, np, func_type)
 650 *    double r                   Convolved result
 651 *    double f1[], f2[]          Functions to be convolved
 652 *    double dx                  Difference between successive x values
 653 *    int n                      Array index to center convolution about
 654 *    int np                     Number of points in f1 array
 655 *    int func_type              EVEN or ODD or EVEN_AND_ODD function f2
 656 *
 657 * NOTES
 658 *    f1 is the projection data, its indices range from 0 to np - 1.
 659 *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 660 *    There are 3 ways to handle the negative vertices of f2:
 661 *       1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 662 *          All filters used in reconstruction are even.
 663 *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 664 *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 665 *          for negative indices.  Since f2 must range from -(np-1) to (np-1),
 666 *          if we add (np - 1) to f2's array index, then f2's index will
 667 *          range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 668 *          stored at f2[np-1].
 669 */
 670
 671 double
 672 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 673 {
 674   double sum = 0.0;
 675
 676 #if UNOPTIMIZED_CONVOLUTION
 677   for (int i = 0; i < np; i++)
 678     sum += func[i] * m_adFilter[n - i + (np - 1)];
 679 #else
 680   double* f2 = m_adFilter + n + (np - 1);
 681   for (int i = 0; i < np; i++)
 682     sum += *func++ * *f2--;
 683 #endif
 684
 685   return (sum * dx);
 686 }
 687
 688
 689 double
 690 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 691 {
 692   double sum = 0.0;
 693
 694 #if UNOPTIMIZED_CONVOLUTION
 695   for (int i = 0; i < np; i++)
 696     sum += func[i] * m_adFilter[n - i + (np - 1)];
 697 #else
 698   double* f2 = m_adFilter + n + (np - 1);
 699   for (int i = 0; i < np; i++)
 700     sum += *func++ * *f2--;
 701 #endif
 702
 703   return (sum * dx);
 704 }
 705
 706
 707 void
 708 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 709 {
 710   std::complex<double>* complexOutput = new std::complex<double> [n];
 711
 712   finiteFourierTransform (input, complexOutput, n, direction);
 713   for (int i = 0; i < n; i++)
 714     output[i] = complexOutput[i].real();
 715   delete [] complexOutput;
 716 }
 717
 718 void
 719 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
 720 {
 721   if (direction < 0)
 722     direction = -1;
 723   else
 724     direction = 1;
 725
 726   double angleIncrement = direction * 2 * PI / n;
 727   for (int i = 0; i < n; i++) {
 728     double sumReal = 0;
 729     double sumImag = 0;
 730     for (int j = 0; j < n; j++) {
 731       double angle = i * j * angleIncrement;
 732       sumReal += input[j] * cos(angle);
 733       sumImag += input[j] * sin(angle);
 734     }
 735     if (direction < 0) {
 736       sumReal /= n;
 737       sumImag /= n;
 738     }
 739     output[i] = std::complex<double> (sumReal, sumImag);
 740   }
 741 }
 742
 743
 744 void
 745 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
 746 {
 747   if (direction < 0)
 748     direction = -1;
 749   else
 750     direction = 1;
 751
 752   double angleIncrement = direction * 2 * PI / n;
 753   for (int i = 0; i < n; i++) {
 754     std::complex<double> sum (0,0);
 755     for (int j = 0; j < n; j++) {
 756       double angle = i * j * angleIncrement;
 757       std::complex<double> exponentTerm (cos(angle), sin(angle));
 758       sum += input[j] * exponentTerm;
 759     }
 760     if (direction < 0) {
 761       sum /= n;
 762     }
 763     output[i] = sum;
 764   }
 765 }
 766
 767 void
 768 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
 769 {
 770   if (direction < 0)
 771     direction = -1;
 772   else
 773     direction = 1;
 774
 775   double angleIncrement = direction * 2 * PI / n;
 776   for (int i = 0; i < n; i++) {
 777     double sumReal = 0;
 778     for (int j = 0; j < n; j++) {
 779       double angle = i * j * angleIncrement;
 780       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 781     }
 782     if (direction < 0) {
 783       sumReal /= n;
 784     }
 785     output[i] = sumReal;
 786   }
 787 }
 788
 789 // Table-based routines
 790
 791 void
 792 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
 793 {
 794   if (direction < 0)
 795     direction = -1;
 796   else
 797     direction = 1;
 798
 799   for (int i = 0; i < m_nFilterPoints; i++) {
 800     double sumReal = 0, sumImag = 0;
 801     for (int j = 0; j < m_nFilterPoints; j++) {
 802       int tableIndex = i * j;
 803       if (direction > 0) {
 804         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 805         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 806       } else {
 807         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 808         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 809       }
 810     }
 811     if (direction < 0) {
 812       sumReal /= m_nFilterPoints;
 813       sumImag /= m_nFilterPoints;
 814     }
 815     output[i] = std::complex<double> (sumReal, sumImag);
 816   }
 817 }
 818
 819 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 820 void
 821 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
 822 {
 823   if (direction < 0)
 824     direction = -1;
 825   else
 826     direction = 1;
 827
 828   for (int i = 0; i < m_nFilterPoints; i++) {
 829     double sumReal = 0, sumImag = 0;
 830     for (int j = 0; j < m_nFilterPoints; j++) {
 831       int tableIndex = i * j;
 832       if (direction > 0) {
 833         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 834           - input[j].imag() * m_adFourierSinTable[tableIndex];
 835         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 836           + input[j].imag() * m_adFourierCosTable[tableIndex];
 837       } else {
 838         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 839           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 840         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 841           + input[j].imag() * m_adFourierCosTable[tableIndex];
 842       }
 843     }
 844     if (direction < 0) {
 845       sumReal /= m_nFilterPoints;
 846       sumImag /= m_nFilterPoints;
 847     }
 848     output[i] = std::complex<double> (sumReal, sumImag);
 849   }
 850 }
 851
 852 void
 853 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
 854 {
 855   if (direction < 0)
 856     direction = -1;
 857   else
 858     direction = 1;
 859
 860   for (int i = 0; i < m_nFilterPoints; i++) {
 861     double sumReal = 0;
 862     for (int j = 0; j < m_nFilterPoints; j++) {
 863       int tableIndex = i * j;
 864       if (direction > 0) {
 865         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 866           - input[j].imag() * m_adFourierSinTable[tableIndex];
 867       } else {
 868         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 869           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 870       }
 871     }
 872     if (direction < 0) {
 873       sumReal /= m_nFilterPoints;
 874     }
 875     output[i] = sumReal;
 876   }
 877 }
 878