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