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-2000 Kevin Rosenberg
  11 **
  12 **  $Id: procsignal.cpp,v 1.17 2001/01/12 14:14:58 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 // FilterMethod ID/Names
  31 const int ProcessSignal::FILTER_METHOD_INVALID = -1;
  32 const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
  33 const int ProcessSignal::FILTER_METHOD_FOURIER = 1;
  34 const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2;
  35 const int ProcessSignal::FILTER_METHOD_FFT = 3;
  36 #if HAVE_FFTW
  37 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
  38 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
  39 #endif
  40 const char* ProcessSignal::s_aszFilterMethodName[] = {
  41   {"convolution"},
  42   {"fourier"},
  43   {"fouier_table"},
  44   {"fft"},
  45 #if HAVE_FFTW
  46   {"fftw"},
  47   {"rfftw"},
  48 #endif
  49 };
  50 const char* ProcessSignal::s_aszFilterMethodTitle[] = {
  51   {"Convolution"},
  52   {"Fourier"},
  53   {"Fouier Trigometric Table"},
  54   {"FFT"},
  55 #if HAVE_FFTW
  56   {"FFTW"},
  57   {"Real/Half-Complex FFTW"},
  58 #endif
  59 };
  60 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
  61
  62 // FilterGeneration ID/Names
  63 const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
  64 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
  65 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
  66 const char* ProcessSignal::s_aszFilterGenerationName[] = {
  67   {"direct"},
  68   {"inverse_fourier"},
  69 };
  70 const char* ProcessSignal::s_aszFilterGenerationTitle[] = {
  71   {"Direct"},
  72   {"Inverse Fourier"},
  73 };
  74 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
  75
  76
  77 // CLASS IDENTIFICATION
  78 //   ProcessSignal
  79 //
  80 ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth,
  81                                                           double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName,
  82                                                           const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel,
  83                                                           int iGeometry, double dFocalLength, SGP* pSGP)
  84 : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
  85 {
  86   m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
  87   if (m_idFilterMethod == FILTER_METHOD_INVALID) {
  88     m_fail = true;
  89     m_failMessage = "Invalid filter method name ";
  90     m_failMessage += szFilterMethodName;
  91     return;
  92   }
  93   m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
  94   if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
  95     m_fail = true;
  96     m_failMessage = "Invalid frequency filter name ";
  97     m_failMessage += szFilterGenerationName;
  98     return;
  99   }
 100   m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
 101   if (m_idFilter == SignalFilter::FILTER_INVALID) {
 102     m_fail = true;
 103     m_failMessage = "Invalid Filter name ";
 104     m_failMessage += szFilterName;
 105     return;
 106   }
 107   m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
 108   if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
 109     m_fail = true;
 110     m_failMessage = "Invalid domain name ";
 111     m_failMessage += szDomainName;
 112     return;
 113   }
 114
 115   init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain,
 116           m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP);
 117 }
 118
 119
 120 void
 121 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement,
 122                                          int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration,
 123                                          const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry,
 124                                          double dFocalLength, SGP* pSGP)
 125 {
 126   int i;
 127   m_idFilter = idFilter;
 128   m_idDomain = idDomain;
 129   m_idFilterMethod = idFilterMethod;
 130   m_idFilterGeneration = idFilterGeneration;
 131   m_idGeometry = iGeometry;
 132   m_dFocalLength = dFocalLength;
 133
 134   if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
 135     m_fail = true;
 136     return;
 137   }
 138   m_traceLevel = iTraceLevel;
 139   m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
 140   m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
 141   m_dBandwidth = dBandwidth;
 142   m_nSignalPoints = nSignalPoints;
 143   m_dSignalInc = dSignalIncrement;
 144   m_dFilterParam = dFilterParam;
 145   m_iZeropad = iZeropad;
 146   m_iPreinterpolationFactor = iPreinterpolationFactor;
 147
 148   // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
 149   // through origin of phantom rather than 2 times distance to detector,
 150   // see Kak-Slaney Fig 3.22, for Collinear diagram
 151   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 152     m_dSignalInc /= 2;
 153     m_dBandwidth *= 2;
 154   }
 155
 156   if (m_idFilterMethod == FILTER_METHOD_FFT) {
 157 #if HAVE_FFTW
 158     m_idFilterMethod = FILTER_METHOD_RFFTW;
 159 #else
 160     m_fail = true;
 161     m_failMessage = "FFT not yet implemented";
 162     return;
 163 #endif
 164   }
 165
 166   bool m_bFrequencyFiltering = true;
 167   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
 168     m_bFrequencyFiltering = false;
 169
 170   // Spatial-based filtering
 171   if (! m_bFrequencyFiltering) {
 172
 173     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
 174       m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
 175       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 176       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 177       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 178       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 179       m_adFilter = new double[ m_nFilterPoints ];
 180       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 181     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 182       m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
 183       m_dFilterMin = -1. / (2 * m_dSignalInc);
 184       m_dFilterMax = 1. / (2 * m_dSignalInc);
 185       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 186       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 187       m_adFilter = new double[ m_nFilterPoints ];
 188       double* adFrequencyFilter = new double [m_nFilterPoints];
 189       filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
 190 #ifdef HAVE_SGP
 191       EZPlot* pEZPlot = NULL;
 192       if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
 193         pEZPlot = new EZPlot ();
 194         pEZPlot->ezset ("title Filter Response: Natural Order");
 195         pEZPlot->ezset ("ylength 0.25");
 196         pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
 197         pEZPlot->plot (pSGP);
 198       }
 199 #endif
 200       Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
 201 #ifdef HAVE_SGP
 202       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
 203         pEZPlot->ezset ("title Filter Response: Fourier Order");
 204         pEZPlot->ezset ("ylength 0.25");
 205         pEZPlot->ezset ("yporigin 0.25");
 206         pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
 207         pEZPlot->plot (pSGP);
 208       }
 209 #endif
 210       ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
 211       delete adFrequencyFilter;
 212 #ifdef HAVE_SGP
 213       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
 214         pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
 215         pEZPlot->ezset ("ylength 0.25");
 216         pEZPlot->ezset ("yporigin 0.50");
 217         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
 218         pEZPlot->plot (pSGP);
 219       }
 220 #endif
 221       Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
 222 #ifdef HAVE_SGP
 223       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
 224         pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order");
 225         pEZPlot->ezset ("ylength 0.25");
 226         pEZPlot->ezset ("yporigin 0.75");
 227         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
 228         pEZPlot->plot (pSGP);
 229         delete pEZPlot;
 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 = sin (iDetFromZero * m_dSignalInc);
 243         if (fabs(sinScale) < 1E-7)
 244           sinScale = 1;
 245         else
 246           sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 247         double dScale = 0.5 * sinScale * sinScale;
 248         m_adFilter[i] *= dScale;
 249       }
 250     } // if (geometry)
 251   } // if (spatial filtering)
 252
 253   else if (m_bFrequencyFiltering) {  // Frequency-based filtering
 254
 255     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
 256       // calculate number of filter points with zeropadding
 257       m_nFilterPoints = m_nSignalPoints;
 258       if (m_iZeropad > 0) {
 259         double logBase2 = log(m_nFilterPoints) / log(2);
 260         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 261         if (logBase2 != floor(logBase2))
 262           nextPowerOf2++;
 263         nextPowerOf2 += (m_iZeropad - 1);
 264         m_nFilterPoints = 1 << nextPowerOf2;
 265 #ifdef DEBUG
 266         if (m_traceLevel >= Trace::TRACE_CONSOLE)
 267           std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
 268 #endif
 269       }
 270       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 271
 272       if (m_nFilterPoints % 2) { // Odd
 273         m_dFilterMin = -1. / (2 * m_dSignalInc);
 274         m_dFilterMax = 1. / (2 * m_dSignalInc);
 275         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 276       } else { // Even
 277         m_dFilterMin = -1. / (2 * m_dSignalInc);
 278         m_dFilterMax = 1. / (2 * m_dSignalInc);
 279         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
 280         m_dFilterMax -= m_dFilterInc;
 281       }
 282
 283       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
 284                   m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 285       m_adFilter = new double [m_nFilterPoints];
 286       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 287
 288       // This doesn't work!
 289       // Need to add filtering for divergent geometries & Frequency/Direct filtering
 290           // Jan 2001: Direct seems to work for equilinear and equiangular
 291           // however, inverse_fourier doesn't work for equiangular on all versions of CTSim tested
 292
 293       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 294         for (i = 0; i < m_nFilterPoints; i++)
 295           m_adFilter[i] *= 0.5;
 296       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 297         for (i = 0; i < m_nFilterPoints; i++) {
 298           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 299           double sinScale = sin (iDetFromZero * m_dSignalInc);
 300           if (fabs(sinScale) < 1E-7)
 301             sinScale = 1;
 302           else
 303             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 304           double dScale = 0.5 * sinScale * sinScale;
 305           m_adFilter[i] *= dScale;
 306         }
 307       }
 308 #ifdef HAVE_SGP
 309       EZPlot* pEZPlot = NULL;
 310       if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
 311         pEZPlot = new EZPlot;
 312         pEZPlot->ezset ("title Filter Filter: Natural Order");
 313         pEZPlot->ezset ("ylength 0.50");
 314         pEZPlot->ezset ("yporigin 0.00");
 315         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
 316         pEZPlot->plot (pSGP);
 317       }
 318 #endif
 319       Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 320 #ifdef HAVE_SGP
 321       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
 322         pEZPlot->ezset ("title Filter Filter: Fourier Order");
 323         pEZPlot->ezset ("ylength 0.50");
 324         pEZPlot->ezset ("yporigin 0.50");
 325         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
 326         pEZPlot->plot (pSGP);
 327         delete pEZPlot;
 328       }
 329 #endif
 330     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 331       // calculate number of filter points with zeropadding
 332       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
 333       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 334       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 335       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
 336       m_nFilterPoints = nSpatialPoints;
 337       if (m_iZeropad > 0) {
 338         double logBase2 = log(nSpatialPoints) / log(2);
 339         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 340         if (logBase2 != floor(logBase2))
 341           nextPowerOf2++;
 342         nextPowerOf2 += (m_iZeropad - 1);
 343         m_nFilterPoints = 1 << nextPowerOf2;
 344       }
 345       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 346 #ifdef DEBUG
 347       if (m_traceLevel >= Trace::TRACE_CONSOLE)
 348         std::cout << "nFilterPoints = " << m_nFilterPoints << endl;
 349 #endif
 350       double* adSpatialFilter = new double [m_nFilterPoints];
 351       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
 352                   m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 353       filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
 354 #ifdef HAVE_SGP
 355       EZPlot* pEZPlot = NULL;
 356       if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) {
 357         pEZPlot = new EZPlot;
 358         pEZPlot->ezset ("title Spatial Filter: Natural Order");
 359         pEZPlot->ezset ("ylength 0.50");
 360         pEZPlot->ezset ("yporigin 0.00");
 361         pEZPlot->addCurve (adSpatialFilter, nSpatialPoints);
 362         pEZPlot->plot (pSGP);
 363         delete pEZPlot;
 364       }
 365 #endif
 366
 367 // #define PRE_JAN_2001 1
 368 #ifdef PRE_JAN_2001
 369       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 370         for (i = 0; i < m_nFilterPoints; i++)
 371           adSpatialFilter[i] *= 0.5;
 372       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 373         for (i = 0; i < m_nFilterPoints; i++) {
 374           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 375           double sinScale = sin (iDetFromZero * m_dSignalInc);
 376           if (fabs(sinScale) < 1E-7)
 377             sinScale = 1;
 378           else
 379             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 380           double dScale = 0.5 * sinScale * sinScale;
 381           adSpatialFilter[i] *= dScale;
 382         }
 383       }
 384       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 385         adSpatialFilter[i] = 0;
 386
 387       m_adFilter = new double [m_nFilterPoints];
 388       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
 389       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
 390       delete adSpatialFilter;
 391       for (i = 0; i < m_nFilterPoints; i++)
 392         m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
 393       delete acInverseFilter;
 394 #else
 395       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 396         adSpatialFilter[i] = 0;
 397
 398       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
 399       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, FORWARD);
 400       delete adSpatialFilter;
 401       m_adFilter = new double [m_nFilterPoints];
 402       for (i = 0; i < m_nFilterPoints; i++)
 403                   m_adFilter[i] = std::abs(acInverseFilter[i]);
 404       delete acInverseFilter;
 405
 406       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 407         for (i = 0; i < m_nFilterPoints; i++)
 408           m_adFilter[i] *= 0.5;
 409       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 410         for (i = 0; i < m_nFilterPoints; i++) {
 411           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 412           double sinScale = sin (iDetFromZero * m_dSignalInc);
 413           if (fabs(sinScale) < 1E-7)
 414             sinScale = 1;
 415           else
 416             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 417           double dScale = 0.5 * sinScale * sinScale;
 418           m_adFilter[i] *= dScale;
 419         }
 420       }
 421 #endif
 422
 423 #ifdef HAVE_SGP
 424       if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
 425         pEZPlot->ezset ("title Spatial Filter: Inverse");
 426         pEZPlot->ezset ("ylength 0.50");
 427         pEZPlot->ezset ("yporigin 0.50");
 428         pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
 429         pEZPlot->plot (pSGP);
 430         delete pEZPlot;
 431       }
 432 #endif
 433     }
 434   }
 435
 436   // precalculate sin and cosine tables for fourier transform
 437   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 438     int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
 439     double angleIncrement = (2. * PI) / m_nFilterPoints;
 440     m_adFourierCosTable = new double[ nFourier ];
 441     m_adFourierSinTable = new double[ nFourier ];
 442     double angle = 0;
 443     for (i = 0; i < nFourier; i++) {
 444       m_adFourierCosTable[i] = cos (angle);
 445       m_adFourierSinTable[i] = sin (angle);
 446       angle += angleIncrement;
 447     }
 448   }
 449
 450 #if HAVE_FFTW
 451   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 452     for (i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 453       m_adFilter[i] /= m_nFilterPoints;
 454   }
 455
 456   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 457     m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 458     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 459     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
 460     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
 461     for (i = 0; i < m_nFilterPoints; i++)
 462       m_adRealFftInput[i] = 0;
 463   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 464     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 465     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 466     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 467     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 468     for (i = 0; i < m_nFilterPoints; i++)
 469       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
 470     for (i = 0; i < m_nOutputPoints; i++)
 471       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
 472   }
 473 #endif
 474
 475 }
 476
 477 ProcessSignal::~ProcessSignal (void)
 478 {
 479   delete [] m_adFourierSinTable;
 480   delete [] m_adFourierCosTable;
 481   delete [] m_adFilter;
 482
 483 #if HAVE_FFTW
 484   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 485     fftw_destroy_plan(m_complexPlanForward);
 486     fftw_destroy_plan(m_complexPlanBackward);
 487     delete [] m_adComplexFftInput;
 488     delete [] m_adComplexFftSignal;
 489   }
 490   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 491     rfftw_destroy_plan(m_realPlanForward);
 492     rfftw_destroy_plan(m_realPlanBackward);
 493     delete [] m_adRealFftInput;
 494     delete [] m_adRealFftSignal;
 495   }
 496 #endif
 497 }
 498
 499 int
 500 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 501 {
 502   int fmID = FILTER_METHOD_INVALID;
 503
 504   for (int i = 0; i < s_iFilterMethodCount; i++)
 505     if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 506       fmID = i;
 507       break;
 508     }
 509
 510     return (fmID);
 511 }
 512
 513 const char *
 514 ProcessSignal::convertFilterMethodIDToName (const int fmID)
 515 {
 516   static const char *name = "";
 517
 518   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 519     return (s_aszFilterMethodName [fmID]);
 520
 521   return (name);
 522 }
 523
 524 const char *
 525 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 526 {
 527   static const char *title = "";
 528
 529   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 530     return (s_aszFilterMethodTitle [fmID]);
 531
 532   return (title);
 533 }
 534
 535
 536 int
 537 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 538 {
 539   int fgID = FILTER_GENERATION_INVALID;
 540
 541   for (int i = 0; i < s_iFilterGenerationCount; i++)
 542     if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 543       fgID = i;
 544       break;
 545     }
 546
 547     return (fgID);
 548 }
 549
 550 const char *
 551 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 552 {
 553   static const char *name = "";
 554
 555   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 556     return (s_aszFilterGenerationName [fgID]);
 557
 558   return (name);
 559 }
 560
 561 const char *
 562 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 563 {
 564   static const char *name = "";
 565
 566   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 567     return (s_aszFilterGenerationTitle [fgID]);
 568
 569   return (name);
 570 }
 571
 572 void
 573 ProcessSignal::filterSignal (const float constInput[], double output[]) const
 574 {
 575   double* input = new double [m_nSignalPoints];
 576   int i;
 577   for (i = 0; i < m_nSignalPoints; i++)
 578     input[i] = constInput[i];
 579
 580   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 581     for (int i = 0; i < m_nSignalPoints; i++) {
 582       int iDetFromCenter = i - (m_nSignalPoints / 2);
 583       input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
 584     }
 585   } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 586     for (int i = 0; i < m_nSignalPoints; i++) {
 587       int iDetFromCenter = i - (m_nSignalPoints / 2);
 588       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
 589     }
 590   }
 591   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 592     for (i = 0; i < m_nSignalPoints; i++)
 593       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 594   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 595     double* inputSignal = new double [m_nFilterPoints];
 596     for (i = 0; i < m_nSignalPoints; i++)
 597       inputSignal[i] = input[i];
 598     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 599       inputSignal[i] = 0;  // zeropad
 600     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 601     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
 602     delete inputSignal;
 603     for (i = 0; i < m_nFilterPoints; i++)
 604       fftSignal[i] *= m_adFilter[i];
 605     double* inverseFourier = new double [m_nFilterPoints];
 606     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
 607     delete fftSignal;
 608     for (i = 0; i < m_nSignalPoints; i++)
 609       output[i] = inverseFourier[i];
 610     delete inverseFourier;
 611   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 612     double* inputSignal = new double [m_nFilterPoints];
 613     for (i = 0; i < m_nSignalPoints; i++)
 614       inputSignal[i] = input[i];
 615     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 616       inputSignal[i] = 0;  // zeropad
 617     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 618     finiteFourierTransform (inputSignal, fftSignal, FORWARD);
 619     delete inputSignal;
 620     for (i = 0; i < m_nFilterPoints; i++)
 621       fftSignal[i] *= m_adFilter[i];
 622     double* inverseFourier = new double [m_nFilterPoints];
 623     finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
 624     delete fftSignal;
 625     for (i = 0; i < m_nSignalPoints; i++)
 626       output[i] = inverseFourier[i];
 627     delete inverseFourier;
 628   }
 629 #if HAVE_FFTW
 630   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 631     for (i = 0; i < m_nSignalPoints; i++)
 632       m_adRealFftInput[i] = input[i];
 633
 634     fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
 635     rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
 636     for (i = 0; i < m_nFilterPoints; i++)
 637       m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
 638     delete [] fftOutput;
 639     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
 640              m_adRealFftSignal[i] = 0;
 641
 642     fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
 643     rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
 644     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 645       output[i] = ifftOutput[i];
 646     delete [] ifftOutput;
 647   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 648     for (i = 0; i < m_nSignalPoints; i++)
 649       m_adComplexFftInput[i].re = input[i];
 650
 651     fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
 652     fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
 653     for (i = 0; i < m_nFilterPoints; i++) {
 654       m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
 655       m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
 656     }
 657     delete [] fftOutput;
 658     fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
 659     fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
 660     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 661       output[i] = ifftOutput[i].re;
 662     delete [] ifftOutput;
 663   }
 664 #endif
 665   delete input;
 666 }
 667
 668
 669 /* NAME
 670 *    convolve                   Discrete convolution of two functions
 671 *
 672 * SYNOPSIS
 673 *    r = convolve (f1, f2, dx, n, np, func_type)
 674 *    double r                   Convolved result
 675 *    double f1[], f2[]          Functions to be convolved
 676 *    double dx                  Difference between successive x values
 677 *    int n                      Array index to center convolution about
 678 *    int np                     Number of points in f1 array
 679 *    int func_type              EVEN or ODD or EVEN_AND_ODD function f2
 680 *
 681 * NOTES
 682 *    f1 is the projection data, its indices range from 0 to np - 1.
 683 *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 684 *    There are 3 ways to handle the negative vertices of f2:
 685 *       1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 686 *          All filters used in reconstruction are even.
 687 *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 688 *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 689 *          for negative indices.  Since f2 must range from -(np-1) to (np-1),
 690 *          if we add (np - 1) to f2's array index, then f2's index will
 691 *          range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 692 *          stored at f2[np-1].
 693 */
 694
 695 double
 696 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 697 {
 698   double sum = 0.0;
 699
 700 #if UNOPTIMIZED_CONVOLUTION
 701   for (int i = 0; i < np; i++)
 702     sum += func[i] * m_adFilter[n - i + (np - 1)];
 703 #else
 704   double* f2 = m_adFilter + n + (np - 1);
 705   for (int i = 0; i < np; i++)
 706     sum += *func++ * *f2--;
 707 #endif
 708
 709   return (sum * dx);
 710 }
 711
 712
 713 double
 714 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 715 {
 716   double sum = 0.0;
 717
 718 #if UNOPTIMIZED_CONVOLUTION
 719   for (int i = 0; i < np; i++)
 720     sum += func[i] * m_adFilter[n - i + (np - 1)];
 721 #else
 722   double* f2 = m_adFilter + n + (np - 1);
 723   for (int i = 0; i < np; i++)
 724     sum += *func++ * *f2--;
 725 #endif
 726
 727   return (sum * dx);
 728 }
 729
 730
 731 void
 732 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 733 {
 734   std::complex<double>* complexOutput = new std::complex<double> [n];
 735
 736   finiteFourierTransform (input, complexOutput, n, direction);
 737   for (int i = 0; i < n; i++)
 738     output[i] = complexOutput[i].real();
 739   delete [] complexOutput;
 740 }
 741
 742 void
 743 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
 744 {
 745   if (direction < 0)
 746     direction = -1;
 747   else
 748     direction = 1;
 749
 750   double angleIncrement = direction * 2 * PI / n;
 751   for (int i = 0; i < n; i++) {
 752     double sumReal = 0;
 753     double sumImag = 0;
 754     for (int j = 0; j < n; j++) {
 755       double angle = i * j * angleIncrement;
 756       sumReal += input[j] * cos(angle);
 757       sumImag += input[j] * sin(angle);
 758     }
 759     if (direction < 0) {
 760       sumReal /= n;
 761       sumImag /= n;
 762     }
 763     output[i] = std::complex<double> (sumReal, sumImag);
 764   }
 765 }
 766
 767
 768 void
 769 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
 770 {
 771   if (direction < 0)
 772     direction = -1;
 773   else
 774     direction = 1;
 775
 776   double angleIncrement = direction * 2 * PI / n;
 777   for (int i = 0; i < n; i++) {
 778     std::complex<double> sum (0,0);
 779     for (int j = 0; j < n; j++) {
 780       double angle = i * j * angleIncrement;
 781       std::complex<double> exponentTerm (cos(angle), sin(angle));
 782       sum += input[j] * exponentTerm;
 783     }
 784     if (direction < 0) {
 785       sum /= n;
 786     }
 787     output[i] = sum;
 788   }
 789 }
 790
 791 void
 792 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
 793 {
 794   if (direction < 0)
 795     direction = -1;
 796   else
 797     direction = 1;
 798
 799   double angleIncrement = direction * 2 * PI / n;
 800   for (int i = 0; i < n; i++) {
 801     double sumReal = 0;
 802     for (int j = 0; j < n; j++) {
 803       double angle = i * j * angleIncrement;
 804       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 805     }
 806     if (direction < 0) {
 807       sumReal /= n;
 808     }
 809     output[i] = sumReal;
 810   }
 811 }
 812
 813 // Table-based routines
 814
 815 void
 816 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
 817 {
 818   if (direction < 0)
 819     direction = -1;
 820   else
 821     direction = 1;
 822
 823   for (int i = 0; i < m_nFilterPoints; i++) {
 824     double sumReal = 0, sumImag = 0;
 825     for (int j = 0; j < m_nFilterPoints; j++) {
 826       int tableIndex = i * j;
 827       if (direction > 0) {
 828         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 829         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 830       } else {
 831         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 832         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 833       }
 834     }
 835     if (direction < 0) {
 836       sumReal /= m_nFilterPoints;
 837       sumImag /= m_nFilterPoints;
 838     }
 839     output[i] = std::complex<double> (sumReal, sumImag);
 840   }
 841 }
 842
 843 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 844 void
 845 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
 846 {
 847   if (direction < 0)
 848     direction = -1;
 849   else
 850     direction = 1;
 851
 852   for (int i = 0; i < m_nFilterPoints; i++) {
 853     double sumReal = 0, sumImag = 0;
 854     for (int j = 0; j < m_nFilterPoints; j++) {
 855       int tableIndex = i * j;
 856       if (direction > 0) {
 857         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 858           - input[j].imag() * m_adFourierSinTable[tableIndex];
 859         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 860           + input[j].imag() * m_adFourierCosTable[tableIndex];
 861       } else {
 862         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 863           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 864         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 865           + input[j].imag() * m_adFourierCosTable[tableIndex];
 866       }
 867     }
 868     if (direction < 0) {
 869       sumReal /= m_nFilterPoints;
 870       sumImag /= m_nFilterPoints;
 871     }
 872     output[i] = std::complex<double> (sumReal, sumImag);
 873   }
 874 }
 875
 876 void
 877 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
 878 {
 879   if (direction < 0)
 880     direction = -1;
 881   else
 882     direction = 1;
 883
 884   for (int i = 0; i < m_nFilterPoints; i++) {
 885     double sumReal = 0;
 886     for (int j = 0; j < m_nFilterPoints; j++) {
 887       int tableIndex = i * j;
 888       if (direction > 0) {
 889         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 890           - input[j].imag() * m_adFourierSinTable[tableIndex];
 891       } else {
 892         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 893           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 894       }
 895     }
 896     if (direction < 0) {
 897       sumReal /= m_nFilterPoints;
 898     }
 899     output[i] = sumReal;
 900   }
 901 }
 902