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