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.13 2001/01/02 05:34:57 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 ();
 186         pEZPlot->ezset ("title Filter Response: Natural Order");
 187         pEZPlot->ezset ("ylength 0.25");
 188         pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints);
 189         pEZPlot->plot (pSGP);
 190       }
 191 #endif
 192       Fourier::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 (pSGP);
 200       }
 201 #endif
 202       ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD);
 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 (pSGP);
 211       }
 212 #endif
 213       Fourier::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 (pSGP);
 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           std::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;
 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 (pSGP);
 305       }
 306 #endif
 307       Fourier::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 (pSGP);
 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         std::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;
 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 (pSGP);
 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       std::complex<double>* acInverseFilter = new std::complex<double> [m_nFilterPoints];\r
 373       finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD);
 374       delete adSpatialFilter;\r
 375       for (i = 0; i < m_nFilterPoints; i++)
 376         m_adFilter[i] = std::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 (pSGP);
 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     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 556     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);\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, BACKWARD);\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     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
 573     finiteFourierTransform (inputSignal, fftSignal, FORWARD);\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, BACKWARD);\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   std::complex<double>* complexOutput = new std::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[], std::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] = std::complex<double> (sumReal, sumImag);
 719   }
 720 }
 721
 722
 723 void
 724 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::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     std::complex<double> sum (0,0);
 734     for (int j = 0; j < n; j++) {
 735       double angle = i * j * angleIncrement;
 736       std::complex<double> exponentTerm (cos(angle), sin(angle));\r
 737       sum += input[j] * exponentTerm;\r
 738     }
 739     if (direction < 0) {
 740       sum /= n;
 741     }
 742     output[i] = sum;
 743   }
 744 }
 745
 746 void
 747 ProcessSignal::finiteFourierTransform (const std::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[], std::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] = std::complex<double> (sumReal, sumImag);
 795   }
 796 }
 797
 798 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
 799 void
 800 ProcessSignal::finiteFourierTransform (const std::complex<double> input[], std::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] = std::complex<double> (sumReal, sumImag);
 828   }
 829 }
 830
 831 void
 832 ProcessSignal::finiteFourierTransform (const std::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