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