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