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