libctsim/procsignal.cpp

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