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