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