libctsim/procsignal.cpp

   1 /*****************************************************************************
   2 ** File IDENTIFICATION
   3 **
   4 **     Name:                   filter.cpp
   5 **     Purpose:                Routines for signal-procesing filters
   6 **     Progammer:              Kevin Rosenberg
   7 **     Date Started:           Aug 1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (c) 1983-2001 Kevin Rosenberg
  11 **
  12 **  $Id: procsignal.cpp,v 1.25 2001/02/11 04:56:37 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30 #ifdef HAVE_WXWINDOWS
  31 #include "dlgezplot.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* 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* 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* ProcessSignal::s_aszFilterGenerationName[] = {
  71   {"direct"},
  72   {"inverse-fourier"},
  73 };
  74 const char* 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, 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, pSGP);
 121 }
 122
 123
 124 void
 125 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement,
 126                      int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration,
 127                      const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry,
 128                      double dFocalLength, SGP* pSGP)
 129 {
 130   int i;
 131   m_idFilter = idFilter;
 132   m_idDomain = idDomain;
 133   m_idFilterMethod = idFilterMethod;
 134   m_idFilterGeneration = idFilterGeneration;
 135   m_idGeometry = iGeometry;
 136   m_dFocalLength = dFocalLength;
 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 signalInc/2 to adjust for imaginary detector
 153   // through origin of phantom rather than 2 times distance to detector,
 154   // see Kak-Slaney Fig 3.22, for Collinear diagram
 155   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 156     m_dSignalInc /= 2;
 157     m_dBandwidth *= 2;
 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 = m_nSignalPoints;
 260       if (m_iZeropad > 0) {
 261         double logBase2 = log(m_nFilterPoints) / log(2);
 262         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 263         if (logBase2 != floor(logBase2))
 264           nextPowerOf2++;
 265         nextPowerOf2 += (m_iZeropad - 1);
 266         m_nFilterPoints = 1 << nextPowerOf2;
 267 #if defined(DEBUG) || defined(_DEBUG)
 268         if (m_traceLevel >= Trace::TRACE_CONSOLE)
 269           sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
 270 #endif
 271       }
 272       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 273
 274       if (m_nFilterPoints % 2) { // Odd
 275         m_dFilterMin = -1. / (2 * m_dSignalInc);
 276         m_dFilterMax = 1. / (2 * m_dSignalInc);
 277         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
 278       } else { // Even
 279         m_dFilterMin = -1. / (2 * m_dSignalInc);
 280         m_dFilterMax = 1. / (2 * m_dSignalInc);
 281         m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
 282         m_dFilterMax -= m_dFilterInc;
 283       }
 284
 285       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth,
 286         m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
 287       m_adFilter = new double [m_nFilterPoints];
 288       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
 289
 290 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 291       if (g_bRunningWXWindows && m_traceLevel > 0) {
 292         EZPlotDialog dlgEZPlot;
 293         dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order");
 294         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 295         dlgEZPlot.ShowModal();
 296       }
 297 #endif
 298
 299       // This works fairly well. I'm not sure why since scaling for geometries is done on
 300       // frequency filter rather than spatial filter as it should be.
 301       // It gives values slightly off than freq/inverse filtering
 302       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 303         for (i = 0; i < m_nFilterPoints; i++)
 304           m_adFilter[i] *= 0.5;
 305       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 306         for (i = 0; i < m_nFilterPoints; i++) {
 307           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 308           double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
 309           double dScale = 0.5 * sinScale * sinScale;
 310           m_adFilter[i] *= dScale;
 311         }
 312       }
 313 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 314       if (g_bRunningWXWindows && m_traceLevel > 0) {
 315         EZPlotDialog dlgEZPlot;
 316         dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
 317         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 318         dlgEZPlot.ShowModal();
 319       }
 320 #endif
 321       Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 322 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 323       if (g_bRunningWXWindows && m_traceLevel > 0) {
 324         EZPlotDialog dlgEZPlot;
 325         dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
 326         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 327         dlgEZPlot.ShowModal();
 328       }
 329 #endif
 330
 331       // FILTERING:  FREQUENCY - INVERSE FOURIER
 332
 333     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 334       // calculate number of filter points with zeropadding
 335       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
 336       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 337       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 338       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
 339       m_nFilterPoints = nSpatialPoints;
 340       if (m_iZeropad > 0) {
 341         double logBase2 = log(nSpatialPoints) / log(2);
 342         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 343         if (logBase2 != floor(logBase2))
 344           nextPowerOf2++;
 345         nextPowerOf2 += (m_iZeropad - 1);
 346         m_nFilterPoints = 1 << nextPowerOf2;
 347       }
 348       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 349 #if defined(DEBUG) || defined(_DEBUG)
 350       if (m_traceLevel >= Trace::TRACE_CONSOLE)
 351         sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
 352 #endif
 353       double* adSpatialFilter = new double [m_nFilterPoints];
 354       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
 355         m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 356       filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
 357 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 358       if (g_bRunningWXWindows && m_traceLevel > 0) {
 359         EZPlotDialog dlgEZPlot;;
 360         dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
 361         dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
 362         dlgEZPlot.ShowModal();
 363       }
 364 #endif
 365
 366       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 367         for (i = 0; i < nSpatialPoints; i++)
 368           adSpatialFilter[i] *= 0.5;
 369       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 370         for (i = 0; i < nSpatialPoints; i++) {
 371           int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
 372           double sinScale = sin (iDetFromZero * m_dSignalInc);
 373           if (fabs(sinScale) < 1E-7)
 374             sinScale = 1;
 375           else
 376             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 377           double dScale = 0.5 * sinScale * sinScale;
 378           adSpatialFilter[i] *= dScale;
 379         }
 380       }
 381 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 382       if (g_bRunningWXWindows && m_traceLevel > 0) {
 383         EZPlotDialog dlgEZPlot;;
 384         dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
 385         dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
 386         dlgEZPlot.ShowModal();
 387       }
 388 #endif
 389       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 390         adSpatialFilter[i] = 0;
 391
 392       m_adFilter = new double [m_nFilterPoints];
 393       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
 394       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
 395       delete adSpatialFilter;
 396       for (i = 0; i < m_nFilterPoints; i++)
 397         m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
 398       delete acInverseFilter;
 399 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 400       if (g_bRunningWXWindows && m_traceLevel > 0) {
 401         EZPlotDialog dlgEZPlot;
 402         dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
 403         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 404         dlgEZPlot.ShowModal();
 405       }
 406 #endif
 407     }
 408   }
 409
 410   // precalculate sin and cosine tables for fourier transform
 411   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 412     int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
 413     double angleIncrement = (2. * PI) / m_nFilterPoints;
 414     m_adFourierCosTable = new double[ nFourier ];
 415     m_adFourierSinTable = new double[ nFourier ];
 416     double angle = 0;
 417     for (i = 0; i < nFourier; i++) {
 418       m_adFourierCosTable[i] = cos (angle);
 419       m_adFourierSinTable[i] = sin (angle);
 420       angle += angleIncrement;
 421     }
 422   }
 423
 424 #if HAVE_FFTW
 425   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 426     for (i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 427       m_adFilter[i] /= m_nFilterPoints;
 428   }
 429
 430   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 431     m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 432     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 433     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
 434     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
 435     for (i = 0; i < m_nFilterPoints; i++)
 436       m_adRealFftInput[i] = 0;
 437   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 438     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 439     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 440     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 441     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 442     for (i = 0; i < m_nFilterPoints; i++)
 443       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
 444     for (i = 0; i < m_nOutputPoints; i++)
 445       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
 446   }
 447 #endif
 448
 449 }
 450
 451 ProcessSignal::~ProcessSignal (void)
 452 {
 453   delete [] m_adFourierSinTable;
 454   delete [] m_adFourierCosTable;
 455   delete [] m_adFilter;
 456
 457 #if HAVE_FFTW
 458   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 459     fftw_destroy_plan(m_complexPlanForward);
 460     fftw_destroy_plan(m_complexPlanBackward);
 461     delete [] m_adComplexFftInput;
 462     delete [] m_adComplexFftSignal;
 463   }
 464   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 465     rfftw_destroy_plan(m_realPlanForward);
 466     rfftw_destroy_plan(m_realPlanBackward);
 467     delete [] m_adRealFftInput;
 468     delete [] m_adRealFftSignal;
 469   }
 470 #endif
 471 }
 472
 473 int
 474 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 475 {
 476   int fmID = FILTER_METHOD_INVALID;
 477
 478   for (int i = 0; i < s_iFilterMethodCount; i++)
 479     if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 480       fmID = i;
 481       break;
 482     }
 483
 484     return (fmID);
 485 }
 486
 487 const char *
 488 ProcessSignal::convertFilterMethodIDToName (const int fmID)
 489 {
 490   static const char *name = "";
 491
 492   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 493     return (s_aszFilterMethodName [fmID]);
 494
 495   return (name);
 496 }
 497
 498 const char *
 499 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 500 {
 501   static const char *title = "";
 502
 503   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 504     return (s_aszFilterMethodTitle [fmID]);
 505
 506   return (title);
 507 }
 508
 509
 510 int
 511 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 512 {
 513   int fgID = FILTER_GENERATION_INVALID;
 514
 515   for (int i = 0; i < s_iFilterGenerationCount; i++)
 516     if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 517       fgID = i;
 518       break;
 519     }
 520
 521     return (fgID);
 522 }
 523
 524 const char *
 525 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 526 {
 527   static const char *name = "";
 528
 529   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 530     return (s_aszFilterGenerationName [fgID]);
 531
 532   return (name);
 533 }
 534
 535 const char *
 536 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 537 {
 538   static const char *name = "";
 539
 540   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 541     return (s_aszFilterGenerationTitle [fgID]);
 542
 543   return (name);
 544 }
 545
 546 void
 547 ProcessSignal::filterSignal (const float constInput[], double output[]) const
 548 {
 549   double* input = new double [m_nSignalPoints];
 550   int i;
 551   for (i = 0; i < m_nSignalPoints; i++)
 552     input[i] = constInput[i];
 553
 554   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 555     for (int i = 0; i < m_nSignalPoints; i++) {
 556       int iDetFromCenter = i - (m_nSignalPoints / 2);
 557       input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
 558     }
 559   } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 560     for (int i = 0; i < m_nSignalPoints; i++) {
 561       int iDetFromCenter = i - (m_nSignalPoints / 2);
 562       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
 563     }
 564   }
 565   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 566     for (i = 0; i < m_nSignalPoints; i++)
 567       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 568   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 569     double* inputSignal = new double [m_nFilterPoints];
 570     for (i = 0; i < m_nSignalPoints; i++)
 571       inputSignal[i] = input[i];
 572     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 573       inputSignal[i] = 0;  // zeropad
 574     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 575     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
 576     delete inputSignal;
 577     for (i = 0; i < m_nFilterPoints; i++)
 578       fftSignal[i] *= m_adFilter[i];
 579     double* inverseFourier = new double [m_nFilterPoints];
 580     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
 581     delete fftSignal;
 582     for (i = 0; i < m_nSignalPoints; i++)
 583       output[i] = inverseFourier[i];
 584     delete inverseFourier;
 585   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 586     double* inputSignal = new double [m_nFilterPoints];
 587     for (i = 0; i < m_nSignalPoints; i++)
 588       inputSignal[i] = input[i];
 589     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 590       inputSignal[i] = 0;  // zeropad
 591     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 592     finiteFourierTransform (inputSignal, fftSignal, FORWARD);
 593     delete inputSignal;
 594     for (i = 0; i < m_nFilterPoints; i++)
 595       fftSignal[i] *= m_adFilter[i];
 596     double* inverseFourier = new double [m_nFilterPoints];
 597     finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
 598     delete fftSignal;
 599     for (i = 0; i < m_nSignalPoints; i++)
 600       output[i] = inverseFourier[i];
 601     delete inverseFourier;
 602   }
 603 #if HAVE_FFTW
 604   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 605     for (i = 0; i < m_nSignalPoints; i++)
 606       m_adRealFftInput[i] = input[i];
 607
 608     fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
 609     rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
 610     for (i = 0; i < m_nFilterPoints; i++)
 611       m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
 612     delete [] fftOutput;
 613     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
 614              m_adRealFftSignal[i] = 0;
 615
 616     fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
 617     rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
 618     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 619       output[i] = ifftOutput[i];
 620     delete [] ifftOutput;
 621   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 622     for (i = 0; i < m_nSignalPoints; i++)
 623       m_adComplexFftInput[i].re = input[i];
 624
 625     fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
 626     fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
 627     for (i = 0; i < m_nFilterPoints; i++) {
 628       m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
 629       m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
 630     }
 631     delete [] fftOutput;
 632     fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
 633     fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
 634     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 635       output[i] = ifftOutput[i].re;
 636     delete [] ifftOutput;
 637   }
 638 #endif
 639   delete input;
 640 }
 641
 642
 643 /* NAME
 644 *    convolve                   Discrete convolution of two functions
 645 *
 646 * SYNOPSIS
 647 *    r = convolve (f1, f2, dx, n, np, func_type)
 648 *    double r                   Convolved result
 649 *    double f1[], f2[]          Functions to be convolved
 650 *    double dx                  Difference between successive x values
 651 *    int n                      Array index to center convolution about
 652 *    int np                     Number of points in f1 array
 653 *    int func_type              EVEN or ODD or EVEN_AND_ODD function f2
 654 *
 655 * NOTES
 656 *    f1 is the projection data, its indices range from 0 to np - 1.
 657 *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 658 *    There are 3 ways to handle the negative vertices of f2:
 659 *       1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 660 *          All filters used in reconstruction are even.
 661 *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 662 *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 663 *          for negative indices.  Since f2 must range from -(np-1) to (np-1),
 664 *          if we add (np - 1) to f2's array index, then f2's index will
 665 *          range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 666 *          stored at f2[np-1].
 667 */
 668
 669 double
 670 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 671 {
 672   double sum = 0.0;
 673
 674 #if UNOPTIMIZED_CONVOLUTION
 675   for (int i = 0; i < np; i++)
 676     sum += func[i] * m_adFilter[n - i + (np - 1)];
 677 #else
 678   double* f2 = m_adFilter + n + (np - 1);
 679   for (int i = 0; i < np; i++)
 680     sum += *func++ * *f2--;
 681 #endif
 682
 683   return (sum * dx);
 684 }
 685
 686
 687 double
 688 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 689 {
 690   double sum = 0.0;
 691
 692 #if UNOPTIMIZED_CONVOLUTION
 693   for (int i = 0; i < np; i++)
 694     sum += func[i] * m_adFilter[n - i + (np - 1)];
 695 #else
 696   double* f2 = m_adFilter + n + (np - 1);
 697   for (int i = 0; i < np; i++)
 698     sum += *func++ * *f2--;
 699 #endif
 700
 701   return (sum * dx);
 702 }
 703
 704
 705 void
 706 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 707 {
 708   std::complex<double>* complexOutput = new std::complex<double> [n];
 709
 710   finiteFourierTransform (input, complexOutput, n, direction);
 711   for (int i = 0; i < n; i++)
 712     output[i] = complexOutput[i].real();
 713   delete [] complexOutput;
 714 }
 715
 716 void
 717 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
 718 {
 719   if (direction < 0)
 720     direction = -1;
 721   else
 722     direction = 1;
 723
 724   double angleIncrement = direction * 2 * PI / n;
 725   for (int i = 0; i < n; i++) {
 726     double sumReal = 0;
 727     double sumImag = 0;
 728     for (int j = 0; j < n; j++) {
 729       double angle = i * j * angleIncrement;
 730       sumReal += input[j] * cos(angle);
 731       sumImag += input[j] * sin(angle);
 732     }
 733     if (direction < 0) {
 734       sumReal /= n;
 735       sumImag /= n;
 736     }
 737     output[i] = std::complex<double> (sumReal, sumImag);
 738   }
 739 }
 740
 741
 742 void
 743 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
 744 {
 745   if (direction < 0)
 746     direction = -1;
 747   else
 748     direction = 1;
 749
 750   double angleIncrement = direction * 2 * PI / n;
 751   for (int i = 0; i < n; i++) {
 752     std::complex<double> sum (0,0);
 753     for (int j = 0; j < n; j++) {
 754       double angle = i * j * angleIncrement;
 755       std::complex<double> exponentTerm (cos(angle), sin(angle));
 756       sum += input[j] * exponentTerm;
 757     }
 758     if (direction < 0) {
 759       sum /= n;
 760     }
 761     output[i] = sum;
 762   }
 763 }
 764
 765 void
 766 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
 767 {
 768   if (direction < 0)
 769     direction = -1;
 770   else
 771     direction = 1;
 772
 773   double angleIncrement = direction * 2 * PI / n;
 774   for (int i = 0; i < n; i++) {
 775     double sumReal = 0;
 776     for (int j = 0; j < n; j++) {
 777       double angle = i * j * angleIncrement;
 778       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 779     }
 780     if (direction < 0) {
 781       sumReal /= n;
 782     }
 783     output[i] = sumReal;
 784   }
 785 }
 786
 787 // Table-based routines
 788
 789 void
 790 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
 791 {
 792   if (direction < 0)
 793     direction = -1;
 794   else
 795     direction = 1;
 796
 797   for (int i = 0; i < m_nFilterPoints; i++) {
 798     double sumReal = 0, sumImag = 0;
 799     for (int j = 0; j < m_nFilterPoints; j++) {
 800       int tableIndex = i * j;
 801       if (direction > 0) {
 802         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 803         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 804       } else {
 805         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 806         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 807       }
 808     }
 809     if (direction < 0) {
 810       sumReal /= m_nFilterPoints;
 811       sumImag /= m_nFilterPoints;
 812     }
 813     output[i] = std::complex<double> (sumReal, sumImag);
 814   }
 815 }
 816
 817 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 818 void
 819 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
 820 {
 821   if (direction < 0)
 822     direction = -1;
 823   else
 824     direction = 1;
 825
 826   for (int i = 0; i < m_nFilterPoints; i++) {
 827     double sumReal = 0, sumImag = 0;
 828     for (int j = 0; j < m_nFilterPoints; j++) {
 829       int tableIndex = i * j;
 830       if (direction > 0) {
 831         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 832           - input[j].imag() * m_adFourierSinTable[tableIndex];
 833         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 834           + input[j].imag() * m_adFourierCosTable[tableIndex];
 835       } else {
 836         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 837           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 838         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 839           + input[j].imag() * m_adFourierCosTable[tableIndex];
 840       }
 841     }
 842     if (direction < 0) {
 843       sumReal /= m_nFilterPoints;
 844       sumImag /= m_nFilterPoints;
 845     }
 846     output[i] = std::complex<double> (sumReal, sumImag);
 847   }
 848 }
 849
 850 void
 851 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
 852 {
 853   if (direction < 0)
 854     direction = -1;
 855   else
 856     direction = 1;
 857
 858   for (int i = 0; i < m_nFilterPoints; i++) {
 859     double sumReal = 0;
 860     for (int j = 0; j < m_nFilterPoints; j++) {
 861       int tableIndex = i * j;
 862       if (direction > 0) {
 863         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 864           - input[j].imag() * m_adFourierSinTable[tableIndex];
 865       } else {
 866         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 867           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 868       }
 869     }
 870     if (direction < 0) {
 871       sumReal /= m_nFilterPoints;
 872     }
 873     output[i] = sumReal;
 874   }
 875 }
 876