libctsim/procsignal.cpp

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