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