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