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
 549 #if HAVE_OPENMP
 550   #pragma omp parallel for
 551 #endif
 552   for (i = 0; i < m_nSignalPoints; i++)
 553     input[i] = constInput[i];
 554
 555   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 556 #if HAVE_OPENMP
 557     #pragma omp parallel for
 558 #endif
 559     for (int i = 0; i < m_nSignalPoints; i++) {
 560       int iDetFromCenter = i - (m_nSignalPoints / 2);
 561       input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
 562     }
 563   } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 564 #if HAVE_OPENMP
 565     #pragma omp parallel for
 566 #endif
 567     for (int i = 0; i < m_nSignalPoints; i++) {
 568       int iDetFromCenter = i - (m_nSignalPoints / 2);
 569       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
 570     }
 571   }
 572   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 573 #if HAVE_OPENMP
 574     #pragma omp parallel for
 575 #endif
 576     for (i = 0; i < m_nSignalPoints; i++)
 577       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 578   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 579     double* inputSignal = new double [m_nFilterPoints];
 580     for (i = 0; i < m_nSignalPoints; i++)
 581       inputSignal[i] = input[i];
 582     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 583       inputSignal[i] = 0;  // zeropad
 584     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 585     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
 586     delete inputSignal;
 587 #if HAVE_OPENMP
 588     #pragma omp parallel for
 589 #endif
 590     for (i = 0; i < m_nFilterPoints; i++)
 591       fftSignal[i] *= m_adFilter[i];
 592     double* inverseFourier = new double [m_nFilterPoints];
 593     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
 594     delete fftSignal;
 595     for (i = 0; i < m_nSignalPoints; i++)
 596       output[i] = inverseFourier[i];
 597     delete inverseFourier;
 598   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 599     double* inputSignal = new double [m_nFilterPoints];
 600     for (i = 0; i < m_nSignalPoints; i++)
 601       inputSignal[i] = input[i];
 602     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 603       inputSignal[i] = 0;  // zeropad
 604     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 605     finiteFourierTransform (inputSignal, fftSignal, FORWARD);
 606     delete inputSignal;
 607 #if HAVE_OPENMP
 608     #pragma omp parallel for
 609 #endif
 610     for (i = 0; i < m_nFilterPoints; i++)
 611       fftSignal[i] *= m_adFilter[i];
 612     double* inverseFourier = new double [m_nFilterPoints];
 613     finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
 614     delete fftSignal;
 615     for (i = 0; i < m_nSignalPoints; i++)
 616       output[i] = inverseFourier[i];
 617     delete inverseFourier;
 618   }
 619 #if HAVE_FFTW
 620   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 621     for (i = 0; i < m_nSignalPoints; i++)
 622       m_adRealFftInput[i] = input[i];
 623
 624     fftw_execute (m_realPlanForward);
 625     for (i = 0; i < m_nFilterPoints; i++)
 626       m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
 627     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
 628              m_adRealFftSignal[i] = 0;
 629
 630     fftw_execute (m_realPlanBackward);
 631     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 632       output[i] = m_adRealFftBackwardOutput[i];
 633   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 634     for (i = 0; i < m_nSignalPoints; i++)
 635       m_adComplexFftInput[i][0] = input[i];
 636
 637     fftw_execute (m_complexPlanForward);
 638 #if HAVE_OPENMP
 639     #pragma omp parallel for
 640 #endif
 641     for (i = 0; i < m_nFilterPoints; i++) {
 642       m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
 643       m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
 644     }
 645     fftw_execute (m_complexPlanBackward);
 646     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 647       output[i] = m_adComplexFftBackwardOutput[i][0];
 648   }
 649 #endif
 650   delete input;
 651 }
 652
 653
 654 /* NAME
 655 *    convolve                   Discrete convolution of two functions
 656 *
 657 * SYNOPSIS
 658 *    r = convolve (f1, f2, dx, n, np, func_type)
 659 *    double r                   Convolved result
 660 *    double f1[], f2[]          Functions to be convolved
 661 *    double dx                  Difference between successive x values
 662 *    int n                      Array index to center convolution about
 663 *    int np                     Number of points in f1 array
 664 *    int func_type              EVEN or ODD or EVEN_AND_ODD function f2
 665 *
 666 * NOTES
 667 *    f1 is the projection data, its indices range from 0 to np - 1.
 668 *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 669 *    There are 3 ways to handle the negative vertices of f2:
 670 *       1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 671 *          All filters used in reconstruction are even.
 672 *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 673 *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 674 *          for negative indices.  Since f2 must range from -(np-1) to (np-1),
 675 *          if we add (np - 1) to f2's array index, then f2's index will
 676 *          range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 677 *          stored at f2[np-1].
 678 */
 679
 680 double
 681 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 682 {
 683   double sum = 0.0;
 684
 685 #if UNOPTIMIZED_CONVOLUTION
 686   for (int i = 0; i < np; i++)
 687     sum += func[i] * m_adFilter[n - i + (np - 1)];
 688 #else
 689   double* f2 = m_adFilter + n + (np - 1);
 690   for (int i = 0; i < np; i++)
 691     sum += *func++ * *f2--;
 692 #endif
 693
 694   return (sum * dx);
 695 }
 696
 697
 698 double
 699 ProcessSignal::convolve (const float 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 void
 717 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 718 {
 719   std::complex<double>* complexOutput = new std::complex<double> [n];
 720
 721   finiteFourierTransform (input, complexOutput, n, direction);
 722   for (int i = 0; i < n; i++)
 723     output[i] = complexOutput[i].real();
 724   delete [] complexOutput;
 725 }
 726
 727 void
 728 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
 729 {
 730   if (direction < 0)
 731     direction = -1;
 732   else
 733     direction = 1;
 734
 735   double angleIncrement = direction * 2 * PI / n;
 736   for (int i = 0; i < n; i++) {
 737     double sumReal = 0;
 738     double sumImag = 0;
 739     for (int j = 0; j < n; j++) {
 740       double angle = i * j * angleIncrement;
 741       sumReal += input[j] * cos(angle);
 742       sumImag += input[j] * sin(angle);
 743     }
 744     if (direction < 0) {
 745       sumReal /= n;
 746       sumImag /= n;
 747     }
 748     output[i] = std::complex<double> (sumReal, sumImag);
 749   }
 750 }
 751
 752
 753 void
 754 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
 755 {
 756   if (direction < 0)
 757     direction = -1;
 758   else
 759     direction = 1;
 760
 761   double angleIncrement = direction * 2 * PI / n;
 762   for (int i = 0; i < n; i++) {
 763     std::complex<double> sum (0,0);
 764     for (int j = 0; j < n; j++) {
 765       double angle = i * j * angleIncrement;
 766       std::complex<double> exponentTerm (cos(angle), sin(angle));
 767       sum += input[j] * exponentTerm;
 768     }
 769     if (direction < 0) {
 770       sum /= n;
 771     }
 772     output[i] = sum;
 773   }
 774 }
 775
 776 void
 777 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
 778 {
 779   if (direction < 0)
 780     direction = -1;
 781   else
 782     direction = 1;
 783
 784   double angleIncrement = direction * 2 * PI / n;
 785   for (int i = 0; i < n; i++) {
 786     double sumReal = 0;
 787     for (int j = 0; j < n; j++) {
 788       double angle = i * j * angleIncrement;
 789       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 790     }
 791     if (direction < 0) {
 792       sumReal /= n;
 793     }
 794     output[i] = sumReal;
 795   }
 796 }
 797
 798 // Table-based routines
 799
 800 void
 801 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
 802 {
 803   if (direction < 0)
 804     direction = -1;
 805   else
 806     direction = 1;
 807
 808   for (int i = 0; i < m_nFilterPoints; i++) {
 809     double sumReal = 0, sumImag = 0;
 810     for (int j = 0; j < m_nFilterPoints; j++) {
 811       int tableIndex = i * j;
 812       if (direction > 0) {
 813         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 814         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 815       } else {
 816         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 817         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 818       }
 819     }
 820     if (direction < 0) {
 821       sumReal /= m_nFilterPoints;
 822       sumImag /= m_nFilterPoints;
 823     }
 824     output[i] = std::complex<double> (sumReal, sumImag);
 825   }
 826 }
 827
 828 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 829 void
 830 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
 831 {
 832   if (direction < 0)
 833     direction = -1;
 834   else
 835     direction = 1;
 836
 837   for (int i = 0; i < m_nFilterPoints; i++) {
 838     double sumReal = 0, sumImag = 0;
 839     for (int j = 0; j < m_nFilterPoints; j++) {
 840       int tableIndex = i * j;
 841       if (direction > 0) {
 842         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 843           - input[j].imag() * m_adFourierSinTable[tableIndex];
 844         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 845           + input[j].imag() * m_adFourierCosTable[tableIndex];
 846       } else {
 847         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 848           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 849         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 850           + input[j].imag() * m_adFourierCosTable[tableIndex];
 851       }
 852     }
 853     if (direction < 0) {
 854       sumReal /= m_nFilterPoints;
 855       sumImag /= m_nFilterPoints;
 856     }
 857     output[i] = std::complex<double> (sumReal, sumImag);
 858   }
 859 }
 860
 861 void
 862 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
 863 {
 864   if (direction < 0)
 865     direction = -1;
 866   else
 867     direction = 1;
 868
 869   for (int i = 0; i < m_nFilterPoints; i++) {
 870     double sumReal = 0;
 871     for (int j = 0; j < m_nFilterPoints; j++) {
 872       int tableIndex = i * j;
 873       if (direction > 0) {
 874         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 875           - input[j].imag() * m_adFourierSinTable[tableIndex];
 876       } else {
 877         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 878           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 879       }
 880     }
 881     if (direction < 0) {
 882       sumReal /= m_nFilterPoints;
 883     }
 884     output[i] = sumReal;
 885   }
 886 }
 887
 888
 889 int
 890 ProcessSignal::addZeropadFactor (int n, int iZeropad)
 891 {
 892   if (iZeropad > 0) {
 893     double dLogBase2 = log(n) / log(2);
 894     int iLogBase2 = static_cast<int>(floor (dLogBase2));
 895     int iPaddedN = 1 << (iLogBase2 + iZeropad);
 896 #ifdef DEBUG
 897     sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);
 898 #endif
 899     return iPaddedN;
 900   }
 901
 902   return n;
 903 }