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.20 2001/01/12 21:53:27 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 doesn't work: Need to add filtering for divergent geometries & Frequency/Direct filtering
 301       // Jan 2001: Direct seems to work for equilinear and equiangular
 302       // however, inverse_fourier doesn't work for equiangular on all versions of CTSim tested
 303       // Scaling is done with data in frequency space, natural order
 304       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 305         for (i = 0; i < m_nFilterPoints; i++)
 306           m_adFilter[i] *= 0.5;
 307       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 308         for (i = 0; i < m_nFilterPoints; i++) {
 309           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 310           double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
 311           double dScale = 0.5 * sinScale * sinScale;
 312           m_adFilter[i] *= dScale;
 313         }
 314       }
 315 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 316       if (g_bRunningWXWindows && m_traceLevel > 0) {
 317         EZPlotDialog dlgEZPlot;
 318         dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order");
 319         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 320         dlgEZPlot.ShowModal();
 321       }
 322 #endif
 323       Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
 324 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 325       if (g_bRunningWXWindows && m_traceLevel > 0) {
 326         EZPlotDialog dlgEZPlot;
 327         dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order");
 328         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 329         dlgEZPlot.ShowModal();
 330       }
 331 #endif
 332
 333       // FILTERING:  FREQUENCY - INVERSE FOURIER
 334
 335     } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
 336       // calculate number of filter points with zeropadding
 337       int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
 338       m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
 339       m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
 340       m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
 341       m_nFilterPoints = nSpatialPoints;
 342       if (m_iZeropad > 0) {
 343         double logBase2 = log(nSpatialPoints) / log(2);
 344         int nextPowerOf2 = static_cast<int>(floor(logBase2));
 345         if (logBase2 != floor(logBase2))
 346           nextPowerOf2++;
 347         nextPowerOf2 += (m_iZeropad - 1);
 348         m_nFilterPoints = 1 << nextPowerOf2;
 349       }
 350       m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
 351 #if defined(DEBUG) || defined(_DEBUG)
 352       if (m_traceLevel >= Trace::TRACE_CONSOLE)
 353         sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
 354 #endif
 355       double* adSpatialFilter = new double [m_nFilterPoints];
 356       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth,
 357         m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
 358       filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
 359 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 360       if (g_bRunningWXWindows && m_traceLevel > 0) {
 361         EZPlotDialog dlgEZPlot;;
 362         dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order");
 363         dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
 364         dlgEZPlot.ShowModal();
 365       }
 366 #endif
 367
 368 #define PRE_JAN_2001 1
 369 #ifdef PRE_JAN_2001
 370       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 371         for (i = 0; i < nSpatialPoints; i++)
 372           adSpatialFilter[i] *= 0.5;
 373       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 374         for (i = 0; i < nSpatialPoints; i++) {
 375           int iDetFromZero = i - ((nSpatialPoints - 1) / 2);
 376           double sinScale = sin (iDetFromZero * m_dSignalInc);
 377           if (fabs(sinScale) < 1E-7)
 378             sinScale = 1;
 379           else
 380             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
 381           double dScale = 0.5 * sinScale * sinScale;
 382           adSpatialFilter[i] *= dScale;
 383         }
 384       }
 385 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 386       if (g_bRunningWXWindows && m_traceLevel > 0) {
 387         EZPlotDialog dlgEZPlot;;
 388         dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order");
 389         dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints);
 390         dlgEZPlot.ShowModal();
 391       }
 392 #endif
 393       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 394         adSpatialFilter[i] = 0;
 395
 396       m_adFilter = new double [m_nFilterPoints];
 397       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
 398       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
 399       delete adSpatialFilter;
 400       for (i = 0; i < m_nFilterPoints; i++)
 401         m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc;
 402       delete acInverseFilter;
 403 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 404       if (g_bRunningWXWindows && m_traceLevel > 0) {
 405         EZPlotDialog dlgEZPlot;
 406         dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order");
 407         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 408         dlgEZPlot.ShowModal();
 409       }
 410 #endif
 411
 412 #else
 413       for (i = nSpatialPoints; i < m_nFilterPoints; i++)
 414         adSpatialFilter[i] = 0;
 415
 416       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];
 417       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, FORWARD);
 418       delete adSpatialFilter;
 419       m_adFilter = new double [m_nFilterPoints];
 420       for (i = 0; i < m_nFilterPoints; i++)
 421         m_adFilter[i] = std::abs(acInverseFilter[i]);
 422       delete acInverseFilter;
 423
 424 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 425       if (g_bRunningWXWindows && m_traceLevel > 0) {
 426         EZPlotDialog dlgEZPlot;
 427         dlgEZPlot.getEZPlot()->ezset ("title Inverse Spatial Filter: Fourier order");
 428         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 429         dlgEZPlot.ShowModal();
 430       }
 431 #endif
 432       Fourier::shuffleFourierToNaturalOrder(m_adFilter, m_nFilterPoints);
 433 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 434       if (g_bRunningWXWindows && m_traceLevel > 0) {
 435         EZPlotDialog dlgEZPlot;
 436         dlgEZPlot.getEZPlot()->ezset ("title Inverse Spatial Filter: Natural order");
 437         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 438         dlgEZPlot.ShowModal();
 439       }
 440 #endif
 441       if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 442         for (i = 0; i < m_nFilterPoints; i++)
 443           m_adFilter[i] *= 0.5;
 444       } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 445         for (i = 0; i < m_nFilterPoints; i++) {
 446           int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
 447           if (abs(iDetFromZero) < m_nSignalPoints) {
 448             double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc);
 449             double dScale = 0.5 * sinScale * sinScale;
 450             m_adFilter[i] *= dScale;
 451           }
 452         }
 453       }
 454 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 455       if (g_bRunningWXWindows && m_traceLevel > 0) {
 456         EZPlotDialog dlgEZPlot;
 457         dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Spatial Filter: Natural order");
 458         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 459         dlgEZPlot.ShowModal();
 460       }
 461 #endif
 462       Fourier::shuffleNaturalToFourierOrder(m_adFilter, m_nFilterPoints);
 463 #endif
 464
 465 #if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG))
 466       if (g_bRunningWXWindows && m_traceLevel > 0) {
 467         EZPlotDialog dlgEZPlot;
 468         dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter Inverse Post Geometry Filtering");
 469         dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints);
 470         dlgEZPlot.ShowModal();
 471       }
 472 #endif
 473     }
 474   }
 475
 476   // precalculate sin and cosine tables for fourier transform
 477   if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 478     int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1;
 479     double angleIncrement = (2. * PI) / m_nFilterPoints;
 480     m_adFourierCosTable = new double[ nFourier ];
 481     m_adFourierSinTable = new double[ nFourier ];
 482     double angle = 0;
 483     for (i = 0; i < nFourier; i++) {
 484       m_adFourierCosTable[i] = cos (angle);
 485       m_adFourierSinTable[i] = sin (angle);
 486       angle += angleIncrement;
 487     }
 488   }
 489
 490 #if HAVE_FFTW
 491   if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
 492     for (i = 0; i < m_nFilterPoints; i++)  //fftw uses unnormalized fft
 493       m_adFilter[i] /= m_nFilterPoints;
 494   }
 495
 496   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 497     m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
 498     m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
 499     m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
 500     m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
 501     for (i = 0; i < m_nFilterPoints; i++)
 502       m_adRealFftInput[i] = 0;
 503   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 504     m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
 505     m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
 506     m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
 507     m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
 508     for (i = 0; i < m_nFilterPoints; i++)
 509       m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
 510     for (i = 0; i < m_nOutputPoints; i++)
 511       m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
 512   }
 513 #endif
 514
 515 }
 516
 517 ProcessSignal::~ProcessSignal (void)
 518 {
 519   delete [] m_adFourierSinTable;
 520   delete [] m_adFourierCosTable;
 521   delete [] m_adFilter;
 522
 523 #if HAVE_FFTW
 524   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 525     fftw_destroy_plan(m_complexPlanForward);
 526     fftw_destroy_plan(m_complexPlanBackward);
 527     delete [] m_adComplexFftInput;
 528     delete [] m_adComplexFftSignal;
 529   }
 530   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 531     rfftw_destroy_plan(m_realPlanForward);
 532     rfftw_destroy_plan(m_realPlanBackward);
 533     delete [] m_adRealFftInput;
 534     delete [] m_adRealFftSignal;
 535   }
 536 #endif
 537 }
 538
 539 int
 540 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 541 {
 542   int fmID = FILTER_METHOD_INVALID;
 543
 544   for (int i = 0; i < s_iFilterMethodCount; i++)
 545     if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
 546       fmID = i;
 547       break;
 548     }
 549
 550     return (fmID);
 551 }
 552
 553 const char *
 554 ProcessSignal::convertFilterMethodIDToName (const int fmID)
 555 {
 556   static const char *name = "";
 557
 558   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 559     return (s_aszFilterMethodName [fmID]);
 560
 561   return (name);
 562 }
 563
 564 const char *
 565 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
 566 {
 567   static const char *title = "";
 568
 569   if (fmID >= 0 && fmID < s_iFilterMethodCount)
 570     return (s_aszFilterMethodTitle [fmID]);
 571
 572   return (title);
 573 }
 574
 575
 576 int
 577 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 578 {
 579   int fgID = FILTER_GENERATION_INVALID;
 580
 581   for (int i = 0; i < s_iFilterGenerationCount; i++)
 582     if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
 583       fgID = i;
 584       break;
 585     }
 586
 587     return (fgID);
 588 }
 589
 590 const char *
 591 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
 592 {
 593   static const char *name = "";
 594
 595   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 596     return (s_aszFilterGenerationName [fgID]);
 597
 598   return (name);
 599 }
 600
 601 const char *
 602 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 603 {
 604   static const char *name = "";
 605
 606   if (fgID >= 0 && fgID < s_iFilterGenerationCount)
 607     return (s_aszFilterGenerationTitle [fgID]);
 608
 609   return (name);
 610 }
 611
 612 void
 613 ProcessSignal::filterSignal (const float constInput[], double output[]) const
 614 {
 615   double* input = new double [m_nSignalPoints];
 616   int i;
 617   for (i = 0; i < m_nSignalPoints; i++)
 618     input[i] = constInput[i];
 619
 620   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
 621     for (int i = 0; i < m_nSignalPoints; i++) {
 622       int iDetFromCenter = i - (m_nSignalPoints / 2);
 623       input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
 624     }
 625   } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
 626     for (int i = 0; i < m_nSignalPoints; i++) {
 627       int iDetFromCenter = i - (m_nSignalPoints / 2);
 628       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
 629     }
 630   }
 631   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
 632     for (i = 0; i < m_nSignalPoints; i++)
 633       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
 634   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
 635     double* inputSignal = new double [m_nFilterPoints];
 636     for (i = 0; i < m_nSignalPoints; i++)
 637       inputSignal[i] = input[i];
 638     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 639       inputSignal[i] = 0;  // zeropad
 640     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 641     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
 642     delete inputSignal;
 643     for (i = 0; i < m_nFilterPoints; i++)
 644       fftSignal[i] *= m_adFilter[i];
 645     double* inverseFourier = new double [m_nFilterPoints];
 646     finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD);
 647     delete fftSignal;
 648     for (i = 0; i < m_nSignalPoints; i++)
 649       output[i] = inverseFourier[i];
 650     delete inverseFourier;
 651   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
 652     double* inputSignal = new double [m_nFilterPoints];
 653     for (i = 0; i < m_nSignalPoints; i++)
 654       inputSignal[i] = input[i];
 655     for (i = m_nSignalPoints; i < m_nFilterPoints; i++)
 656       inputSignal[i] = 0;  // zeropad
 657     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 658     finiteFourierTransform (inputSignal, fftSignal, FORWARD);
 659     delete inputSignal;
 660     for (i = 0; i < m_nFilterPoints; i++)
 661       fftSignal[i] *= m_adFilter[i];
 662     double* inverseFourier = new double [m_nFilterPoints];
 663     finiteFourierTransform (fftSignal, inverseFourier, BACKWARD);
 664     delete fftSignal;
 665     for (i = 0; i < m_nSignalPoints; i++)
 666       output[i] = inverseFourier[i];
 667     delete inverseFourier;
 668   }
 669 #if HAVE_FFTW
 670   else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
 671     for (i = 0; i < m_nSignalPoints; i++)
 672       m_adRealFftInput[i] = input[i];
 673
 674     fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
 675     rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
 676     for (i = 0; i < m_nFilterPoints; i++)
 677       m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
 678     delete [] fftOutput;
 679     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
 680              m_adRealFftSignal[i] = 0;
 681
 682     fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
 683     rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
 684     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 685       output[i] = ifftOutput[i];
 686     delete [] ifftOutput;
 687   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
 688     for (i = 0; i < m_nSignalPoints; i++)
 689       m_adComplexFftInput[i].re = input[i];
 690
 691     fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
 692     fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
 693     for (i = 0; i < m_nFilterPoints; i++) {
 694       m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
 695       m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
 696     }
 697     delete [] fftOutput;
 698     fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
 699     fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
 700     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
 701       output[i] = ifftOutput[i].re;
 702     delete [] ifftOutput;
 703   }
 704 #endif
 705   delete input;
 706 }
 707
 708
 709 /* NAME
 710 *    convolve                   Discrete convolution of two functions
 711 *
 712 * SYNOPSIS
 713 *    r = convolve (f1, f2, dx, n, np, func_type)
 714 *    double r                   Convolved result
 715 *    double f1[], f2[]          Functions to be convolved
 716 *    double dx                  Difference between successive x values
 717 *    int n                      Array index to center convolution about
 718 *    int np                     Number of points in f1 array
 719 *    int func_type              EVEN or ODD or EVEN_AND_ODD function f2
 720 *
 721 * NOTES
 722 *    f1 is the projection data, its indices range from 0 to np - 1.
 723 *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 724 *    There are 3 ways to handle the negative vertices of f2:
 725 *       1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 726 *          All filters used in reconstruction are even.
 727 *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 728 *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 729 *          for negative indices.  Since f2 must range from -(np-1) to (np-1),
 730 *          if we add (np - 1) to f2's array index, then f2's index will
 731 *          range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 732 *          stored at f2[np-1].
 733 */
 734
 735 double
 736 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
 737 {
 738   double sum = 0.0;
 739
 740 #if UNOPTIMIZED_CONVOLUTION
 741   for (int i = 0; i < np; i++)
 742     sum += func[i] * m_adFilter[n - i + (np - 1)];
 743 #else
 744   double* f2 = m_adFilter + n + (np - 1);
 745   for (int i = 0; i < np; i++)
 746     sum += *func++ * *f2--;
 747 #endif
 748
 749   return (sum * dx);
 750 }
 751
 752
 753 double
 754 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
 755 {
 756   double sum = 0.0;
 757
 758 #if UNOPTIMIZED_CONVOLUTION
 759   for (int i = 0; i < np; i++)
 760     sum += func[i] * m_adFilter[n - i + (np - 1)];
 761 #else
 762   double* f2 = m_adFilter + n + (np - 1);
 763   for (int i = 0; i < np; i++)
 764     sum += *func++ * *f2--;
 765 #endif
 766
 767   return (sum * dx);
 768 }
 769
 770
 771 void
 772 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
 773 {
 774   std::complex<double>* complexOutput = new std::complex<double> [n];
 775
 776   finiteFourierTransform (input, complexOutput, n, direction);
 777   for (int i = 0; i < n; i++)
 778     output[i] = complexOutput[i].real();
 779   delete [] complexOutput;
 780 }
 781
 782 void
 783 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], const int n, int direction)
 784 {
 785   if (direction < 0)
 786     direction = -1;
 787   else
 788     direction = 1;
 789
 790   double angleIncrement = direction * 2 * PI / n;
 791   for (int i = 0; i < n; i++) {
 792     double sumReal = 0;
 793     double sumImag = 0;
 794     for (int j = 0; j < n; j++) {
 795       double angle = i * j * angleIncrement;
 796       sumReal += input[j] * cos(angle);
 797       sumImag += input[j] * sin(angle);
 798     }
 799     if (direction < 0) {
 800       sumReal /= n;
 801       sumImag /= n;
 802     }
 803     output[i] = std::complex<double> (sumReal, sumImag);
 804   }
 805 }
 806
 807
 808 void
 809 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, int direction)
 810 {
 811   if (direction < 0)
 812     direction = -1;
 813   else
 814     direction = 1;
 815
 816   double angleIncrement = direction * 2 * PI / n;
 817   for (int i = 0; i < n; i++) {
 818     std::complex<double> sum (0,0);
 819     for (int j = 0; j < n; j++) {
 820       double angle = i * j * angleIncrement;
 821       std::complex<double> exponentTerm (cos(angle), sin(angle));
 822       sum += input[j] * exponentTerm;
 823     }
 824     if (direction < 0) {
 825       sum /= n;
 826     }
 827     output[i] = sum;
 828   }
 829 }
 830
 831 void
 832 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], const int n, int direction)
 833 {
 834   if (direction < 0)
 835     direction = -1;
 836   else
 837     direction = 1;
 838
 839   double angleIncrement = direction * 2 * PI / n;
 840   for (int i = 0; i < n; i++) {
 841     double sumReal = 0;
 842     for (int j = 0; j < n; j++) {
 843       double angle = i * j * angleIncrement;
 844       sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
 845     }
 846     if (direction < 0) {
 847       sumReal /= n;
 848     }
 849     output[i] = sumReal;
 850   }
 851 }
 852
 853 // Table-based routines
 854
 855 void
 856 ProcessSignal::finiteFourierTransform (const double input[], std::complex<double> output[], int direction) const
 857 {
 858   if (direction < 0)
 859     direction = -1;
 860   else
 861     direction = 1;
 862
 863   for (int i = 0; i < m_nFilterPoints; i++) {
 864     double sumReal = 0, sumImag = 0;
 865     for (int j = 0; j < m_nFilterPoints; j++) {
 866       int tableIndex = i * j;
 867       if (direction > 0) {
 868         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 869         sumImag += input[j] * m_adFourierSinTable[tableIndex];
 870       } else {
 871         sumReal += input[j] * m_adFourierCosTable[tableIndex];
 872         sumImag -= input[j] * m_adFourierSinTable[tableIndex];
 873       }
 874     }
 875     if (direction < 0) {
 876       sumReal /= m_nFilterPoints;
 877       sumImag /= m_nFilterPoints;
 878     }
 879     output[i] = std::complex<double> (sumReal, sumImag);
 880   }
 881 }
 882
 883 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 884 void
 885 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], int direction) const
 886 {
 887   if (direction < 0)
 888     direction = -1;
 889   else
 890     direction = 1;
 891
 892   for (int i = 0; i < m_nFilterPoints; i++) {
 893     double sumReal = 0, sumImag = 0;
 894     for (int j = 0; j < m_nFilterPoints; j++) {
 895       int tableIndex = i * j;
 896       if (direction > 0) {
 897         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 898           - input[j].imag() * m_adFourierSinTable[tableIndex];
 899         sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
 900           + input[j].imag() * m_adFourierCosTable[tableIndex];
 901       } else {
 902         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 903           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 904         sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
 905           + input[j].imag() * m_adFourierCosTable[tableIndex];
 906       }
 907     }
 908     if (direction < 0) {
 909       sumReal /= m_nFilterPoints;
 910       sumImag /= m_nFilterPoints;
 911     }
 912     output[i] = std::complex<double> (sumReal, sumImag);
 913   }
 914 }
 915
 916 void
 917 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], double output[], int direction) const
 918 {
 919   if (direction < 0)
 920     direction = -1;
 921   else
 922     direction = 1;
 923
 924   for (int i = 0; i < m_nFilterPoints; i++) {
 925     double sumReal = 0;
 926     for (int j = 0; j < m_nFilterPoints; j++) {
 927       int tableIndex = i * j;
 928       if (direction > 0) {
 929         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 930           - input[j].imag() * m_adFourierSinTable[tableIndex];
 931       } else {
 932         sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
 933           - input[j].imag() * -m_adFourierSinTable[tableIndex];
 934       }
 935     }
 936     if (direction < 0) {
 937       sumReal /= m_nFilterPoints;
 938     }
 939     output[i] = sumReal;
 940   }
 941 }
 942