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