libctsim/filter.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: filter.cpp,v 1.2 2000/06/20 17:54:51 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
  31 /* NAME
  32  *   filter_generate                            Generate a filter
  33  *
  34  * SYNOPSIS
  35  *   f = filter_generate (filt_type, bw, xmin, xmax, n, param, domain, analytic)
  36  *   double f                           Generated filter vector
  37  *   int filt_type                      Type of filter wanted
  38  *   double bw                          Bandwidth of filter
  39  *   double xmin, xmax                  Filter limits
  40  *   int n                              Number of points in filter
  41  *   double param                       General input parameter to filters
  42  *   int domain                         FREQ or SPATIAL domain wanted
  43  *   int numint                         Number if intervals for calculating
  44  *                                      discrete inverse fourier xform
  45  *                                      for spatial domain filters.  For
  46  *                                      ANALYTIC solutions, use numint = 0
  47  */
  48
  49 SignalFilter::SignalFilter (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint)
  50 {
  51   m_vecFilter = new double [n];
  52   m_filterType = filt_type;
  53   m_bw = bw;
  54
  55   double xinc = (xmax - xmin) / (n - 1);
  56
  57   if (m_filterType == FILTER_SHEPP) {
  58     double a = 2 * m_bw;
  59     double c = - 4. / (a * a);
  60     int center = (n - 1) / 2;
  61     int sidelen = center;
  62     m_vecFilter[center] = 4. / (a * a);
  63
  64     for (int i = 1; i <= sidelen; i++ )
  65       m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
  66   } else if (domain == D_FREQ) {
  67     double x;
  68     int i;
  69     for (x = xmin, i = 0; i < n; x += xinc, i++)
  70       m_vecFilter[i] = frequencyResponse (x, param);
  71   } else if (domain == D_SPATIAL) {
  72     double x;
  73     int i;
  74     for (x = xmin, i = 0; i < n; x += xinc, i++)
  75       if (numint == 0)
  76         m_vecFilter[i] = spatialResponseAnalytic (x, param);
  77       else
  78         m_vecFilter[i] = spatialResponseCalc (x, param, numint);
  79   } else {
  80     sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain);
  81   }
  82 }
  83
  84 SignalFilter::~SignalFilter (void)
  85 {
  86   delete m_vecFilter;
  87 }
  88
  89
  90 /* NAME
  91  *   filter_spatial_response_calc       Calculate filter by discrete inverse fourier
  92  *                                      transform of filters's frequency
  93  *                                      response
  94  *
  95  * SYNOPSIS
  96  *   y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
  97  *   double y                   Filter's response in spatial domain
  98  *   int filt_type              Type of filter (definitions in ct.h)
  99  *   double x                   Spatial position to evaluate filter
 100  *   double m_bw                        Bandwidth of window
 101  *   double param               General parameter for various filters
 102  *   int n                      Number of points to calculate integrations
 103  */
 104
 105 double
 106 SignalFilter::spatialResponseCalc (double x, double param, int n) const
 107 {
 108   return (spatialResponseCalc (m_filterType, m_bw, x, param, n));
 109 }
 110
 111 double
 112 SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double param, int n)
 113 {
 114   double zmin, zmax;
 115
 116   if (fType == FILTER_TRIANGLE) {
 117     zmin = 0;
 118     zmax = bw;
 119   } else {
 120     zmin = 0;
 121     zmax = bw / 2;
 122   }
 123   double zinc = (zmax - zmin) / (n - 1);
 124
 125   double z = zmin;
 126   double q [n];
 127   for (int i = 0; i < n; i++, z += zinc)
 128     q[i] = frequencyResponse (fType, bw, z, param) * cos (TWOPI * z * x);
 129
 130   double y = 2 * integrateSimpson (zmin, zmax, q, n);
 131
 132   return (y);
 133 }
 134
 135
 136 /* NAME
 137  *    filter_frequency_response                 Return filter frequency response
 138  *
 139  * SYNOPSIS
 140  *    h = filter_frequency_response (filt_type, u, m_bw, param)
 141  *    double h                  Filters frequency response at u
 142  *    int filt_type             Type of filter
 143  *    double u                  Frequency to evaluate filter at
 144  *    double m_bw                       Bandwidth of filter
 145  *    double param              General input parameter for various filters
 146  */
 147
 148 double
 149 SignalFilter::frequencyResponse (double u, double param) const
 150 {
 151   return frequencyResponse (m_filterType, m_bw, u, param);
 152 }
 153
 154
 155 double
 156 SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double param)
 157 {
 158   double q;
 159   double au = fabs (u);
 160
 161   switch (fType) {
 162   case FILTER_BANDLIMIT:
 163     if (au >= bw / 2)
 164       q = 0.;
 165     else
 166       q = 1;
 167     break;
 168   case FILTER_ABS_BANDLIMIT:
 169     if (au >= bw / 2)
 170       q = 0.;
 171     else
 172       q = au;
 173     break;
 174   case FILTER_TRIANGLE:
 175     if (au >= bw)
 176       q = 0;
 177     else
 178       q = 1 - au / bw;
 179     break;
 180   case FILTER_COSINE:
 181     if (au >= bw / 2)
 182       q = 0;
 183     else
 184       q = cos(PI * u / bw);
 185     break;
 186   case FILTER_ABS_COSINE:
 187     if (au >= bw / 2)
 188       q = 0;
 189     else
 190       q = au * cos(PI * u / bw);
 191     break;
 192   case FILTER_SINC:
 193     q = bw * sinc (PI * bw * u, 1.);
 194     break;
 195   case FILTER_ABS_SINC:
 196     q = au * bw * sinc (PI * bw * u, 1.);
 197     break;
 198   case FILTER_G_HAMMING:
 199     if (au >= bw / 2)
 200       q = 0;
 201     else
 202       q = param + (1 - param) * cos (TWOPI * u / bw);
 203     break;
 204   case FILTER_ABS_G_HAMMING:
 205     if (au >= bw / 2)
 206       q = 0;
 207     else
 208       q = au * (param + (1 - param) * cos(TWOPI * u / bw));
 209     break;
 210   default:
 211     q = 0;
 212     sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", fType);
 213     break;
 214   }
 215   return (q);
 216 }
 217
 218
 219
 220 /* NAME
 221  *   filter_spatial_response_analytic                   Calculate filter by analytic inverse fourier
 222  *                              transform of filters's frequency
 223  *                              response
 224  *
 225  * SYNOPSIS
 226  *   y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
 227  *   double y                   Filter's response in spatial domain
 228  *   int filt_type              Type of filter (definitions in ct.h)
 229  *   double x                   Spatial position to evaluate filter
 230  *   double m_bw                        Bandwidth of window
 231  *   double param               General parameter for various filters
 232  */
 233
 234 double
 235 SignalFilter::spatialResponseAnalytic (double x, double param) const
 236 {
 237   return spatialResponseAnalytic (m_filterType, m_bw, x, param);
 238 }
 239
 240 double
 241 SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, double param)
 242 {
 243   double q, temp;
 244   double u = TWOPI * x;
 245   double w = bw / 2;
 246   double b = PI / bw;
 247   double b2 = TWOPI / bw;
 248
 249   switch (fType) {
 250   case FILTER_BANDLIMIT:
 251     q = bw * sinc(u * w, 1.0);
 252     break;
 253   case FILTER_TRIANGLE:
 254     temp = sinc (u * w, 1.0);
 255     q = bw * temp * temp;
 256     break;
 257   case FILTER_COSINE:
 258     q = sinc(b-u,w) + sinc(b+u,w);
 259     break;
 260   case FILTER_G_HAMMING:
 261     q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
 262     break;
 263   case FILTER_ABS_BANDLIMIT:
 264     q = 2 * integral_abscos (u, w);
 265     break;
 266   case FILTER_ABS_COSINE:
 267     q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
 268     break;
 269   case FILTER_ABS_G_HAMMING:
 270     q = 2 * param * integral_abscos(u,w) +
 271       (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
 272     break;
 273   case FILTER_SHEPP:
 274     if (fabs (u) < 1E-7)
 275       q = 4. / (PI * bw * bw);
 276     else
 277       q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
 278     break;
 279   case FILTER_SINC:
 280     if (fabs (x) < bw / 2)
 281       q = 1.;
 282     else
 283       q = 0.;
 284     break;
 285   case FILTER_ABS_SINC:
 286   default:
 287     sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", fType);
 288     q = 0;
 289     break;
 290   }
 291
 292   return (q);
 293 }
 294
 295
 296 /* NAME
 297  *   sinc                       Return sin(x)/x function
 298  *
 299  * SYNOPSIS
 300  *   v = sinc (x, mult)
 301  *   double v                   sinc value
 302  *   double x, mult
 303  *
 304  * DESCRIPTION
 305  *   v = sin(x * mult) / x;
 306  */
 307
 308
 309 /* NAME
 310  *   integral_abscos                    Returns integral of u*cos(u)
 311  *
 312  * SYNOPSIS
 313  *   q = integral_abscos (u, w)
 314  *   double q                   Integral value
 315  *   double u                   Integration variable
 316  *   double w                   Upper integration boundary
 317  *
 318  * DESCRIPTION
 319  *   Returns the value of integral of u*cos(u)*dV for V = 0 to w
 320  */
 321
 322 double
 323 SignalFilter::integral_abscos (double u, double w)
 324 {
 325   if (fabs (u) > F_EPSILON)
 326     return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
 327   else
 328     return (w * w / 2);
 329 }
 330
 331
 332 /* NAME
 333  *    convolve                  Discrete convolution of two functions
 334  *
 335  * SYNOPSIS
 336  *    r = convolve (f1, f2, dx, n, np, func_type)
 337  *    double r                  Convolved result
 338  *    double f1[], f2[]         Functions to be convolved
 339  *    double dx                 Difference between successive x values
 340  *    int n                     Array index to center convolution about
 341  *    int np                    Number of points in f1 array
 342  *    int func_type             EVEN or ODD or EVEN_AND_ODD function f2
 343  *
 344  * NOTES
 345  *    f1 is the projection data, its indices range from 0 to np - 1.
 346  *    The index for f2, the filter, ranges from -(np-1) to (np-1).
 347  *    There are 3 ways to handle the negative vertices of f2:
 348  *      1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
 349  *         All filters used in reconstruction are even.
 350  *      2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
 351  *      3. If f2 is both ODD AND EVEN, then we must store the value of f2
 352  *         for negative indices.  Since f2 must range from -(np-1) to (np-1),
 353  *         if we add (np - 1) to f2's array index, then f2's index will
 354  *         range from 0 to 2 * (np - 1), and the origin, x = 0, will be
 355  *         stored at f2[np-1].
 356  */
 357
 358 double
 359 SignalFilter::convolve (const double func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
 360 {
 361   double sum = 0.0;
 362
 363   if (func_type == FUNC_BOTH) {
 364 #if UNOPTIMIZED_CONVOLUTION
 365     for (int i = 0; i < np; i++)
 366       sum += func[i] * m_vecFilter[n - i + (np - 1)];
 367 #else
 368     double* f2 = m_vecFilter + n + (np - 1);
 369     for (int i = 0; i < np; i++)
 370       sum += *func++ * *f2--;
 371 #endif
 372   }  else if (func_type == FUNC_EVEN) {
 373     for (int i = 0; i < np; i++) {
 374       int k = abs (n - i);
 375       sum += func[i] * m_vecFilter[k];
 376     }
 377   } else if (func_type == FUNC_ODD) {
 378     for (int i = 0; i < np; i++) {
 379       int k = n - i;
 380       if (k < 0)
 381         sum -= func[i] * m_vecFilter[k];
 382       else
 383         sum += func[i] * m_vecFilter[k];
 384     }
 385   } else
 386     sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
 387
 388   return (sum * dx);
 389 }
 390
 391
 392 double
 393 SignalFilter::convolve (const float func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
 394 {
 395   double sum = 0.0;
 396
 397   if (func_type == FUNC_BOTH) {
 398 #if UNOPTIMIZED_CONVOLUTION
 399     for (int i = 0; i < np; i++)
 400       sum += func[i] * m_vecFilter[n - i + (np - 1)];
 401 #else
 402     double* f2 = m_vecFilter + n + (np - 1);
 403     for (int i = 0; i < np; i++)
 404       sum += *func++ * *f2--;
 405 #endif
 406   }  else if (func_type == FUNC_EVEN) {
 407     for (int i = 0; i < np; i++) {
 408       int k = abs (n - i);
 409       sum += func[i] * m_vecFilter[k];
 410     }
 411   } else if (func_type == FUNC_ODD) {
 412     for (int i = 0; i < np; i++) {
 413       int k = n - i;
 414       if (k < 0)
 415         sum -= func[i] * m_vecFilter[k];
 416       else
 417         sum += func[i] * m_vecFilter[k];
 418     }
 419   } else
 420     sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
 421
 422   return (sum * dx);
 423 }
 424