home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / 3dscale_geom / filt.c next >
Encoding:
C/C++ Source or Header  |  1991-08-30  |  11.1 KB  |  476 lines
  1. /*
  2.  * filt: package of 1-d signal filters, both FIR and IIR
  3.  *
  4.  * Paul Heckbert, ph@miro.berkeley.edu    23 Oct 1986, 10 Sept 1988
  5.  *
  6.  * Copyright (c) 1989  Paul S. Heckbert
  7.  * This source may be used for peaceful, nonprofit purposes only, unless
  8.  * under licence from the author. This notice should remain in the source.
  9.  *
  10.  *
  11.  * Documentation:
  12.  *  To use these routines,
  13.  *    #include <filt.h>
  14.  *  Then call filt_find to select the desired filter, e.g.
  15.  *    Filt *f;
  16.  *    f = filt_find("catrom");
  17.  *  This filter function (impulse response) is then evaluated by calling
  18.  *  filt_func(f, x).  Each filter is nonzero between -f->supp and f->supp.
  19.  *  Typically one will use the filter something like this:
  20.  *    double phase, x, weight;
  21.  *    for (x=phase-f->supp; x<f->supp; x+=deltax) {
  22.  *        weight = filt_func(f, x);
  23.  *        # do something with weight
  24.  *    }
  25.  *
  26.  *  Example of windowing an IIR filter:
  27.  *    f = filt_find("sinc");
  28.  *    # if you wanted to vary sinc width, set f->supp = desired_width here
  29.  *    # e.g. f->supp = 6;
  30.  *    f = filt_window(f, "blackman");
  31.  *  You can then use f as above.
  32.  */
  33.  
  34. /*
  35.  * references:
  36.  *
  37.  * A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975
  38.  *
  39.  * R.W. Hamming, Digital Filters, Prentice-Hall, Englewood Cliffs, NJ, 1983
  40.  *
  41.  * W.K. Pratt, Digital Image Processing, John Wiley and Sons, 1978
  42.  *
  43.  * H.S. Hou, H.C. Andrews, "Cubic Splines for Image Interpolation and
  44.  *    Digital Filtering", IEEE Trans. Acoustics, Speech, and Signal Proc.,
  45.  *    vol. ASSP-26, no. 6, Dec. 1978, pp. 508-517
  46.  */
  47.  
  48. static char rcsid[] = "$Header: filt.c,v 2.2 88/12/30 21:32:17 ph Locked $";
  49. #include <math.h>
  50.  
  51. #include "simple.h"
  52. #include "filt.h"
  53.  
  54. #define EPSILON 1e-7
  55.  
  56. typedef struct {    /* data for parameterized Mitchell filter */
  57.     double p0, p2, p3;
  58.     double q0, q1, q2, q3;
  59. } mitchell_data;
  60.  
  61. typedef struct {    /* data for parameterized Kaiser window */
  62.     double a;        /* = w*(N-1)/2 in Oppenheim&Schafer notation */
  63.     double i0a;
  64.     /*
  65.      * typically 4<a<9
  66.      * param a trades off main lobe width (sharpness)
  67.      * for side lobe amplitude (ringing)
  68.      */
  69. } kaiser_data;
  70.  
  71. typedef struct {    /* data for windowed function compound filter */
  72.     Filt *filter;
  73.     Filt *window;
  74. } window_data;
  75.  
  76. void mitchell_init(), mitchell_print();
  77. void kaiser_init(), kaiser_print();
  78. void window_print();
  79. double window_func();
  80. static mitchell_data md;
  81. static kaiser_data kd;
  82.  
  83. #define NFILTMAX 30
  84. static int nfilt = 0;
  85.  
  86. /*
  87.  * note: for the IIR (infinite impulse response) filters,
  88.  * gaussian, sinc, etc, the support values given are arbitrary,
  89.  * convenient cutoff points, while for the FIR (finite impulse response)
  90.  * filters the support is finite and absolute.
  91.  */
  92.  
  93. static Filt filt[NFILTMAX] = {
  94. /*  NAME    FUNC        SUPP    BLUR WIN CARD UNIT  OPT... */
  95.    {"point",    filt_box,    0.0,    1.0,  0,  1,  1},
  96.    {"box",    filt_box,       0.5,    1.0,  0,  1,  1},
  97.    {"triangle",    filt_triangle,  1.0,    1.0,  0,  1,  1},
  98.    {"quadratic",filt_quadratic, 1.5,    1.0,  0,  0,  1},
  99.    {"cubic",    filt_cubic,     2.0,    1.0,  0,  0,  1},
  100.  
  101.    {"catrom",    filt_catrom,    2.0,    1.0,  0,  1,  0},
  102.    {"mitchell",    filt_mitchell,    2.0,    1.0,  0,  0,  0,
  103.     mitchell_init, mitchell_print, (char *)&md},
  104.  
  105.    {"gaussian",    filt_gaussian,  1.25,    1.0,  0,  0,  1},
  106.    {"sinc",    filt_sinc,      4.0,    1.0,  1,  1,  0},
  107.    {"bessel",    filt_bessel,    3.2383,    1.0,  1,  0,  0},
  108.  
  109.    {"hanning",    filt_hanning,    1.0,    1.0,  0,  1,  1},
  110.    {"hamming",    filt_hamming,    1.0,    1.0,  0,  1,  1},
  111.    {"blackman",    filt_blackman,    1.0,    1.0,  0,  1,  1},
  112.    {"kaiser",    filt_kaiser,    1.0,    1.0,  0,  1,  1,
  113.     kaiser_init, kaiser_print, (char *)&kd},
  114.    {0}
  115. };
  116.  
  117. /*-------------------- general filter routines --------------------*/
  118.  
  119. static filt_init()
  120. {
  121.     /* count the filters to set nfilt */
  122.     for (nfilt=0; nfilt<NFILTMAX && filt[nfilt].name; nfilt++);
  123.     /* better way?? */
  124.     mitchell_init(1./3., 1./3., (char *)&md);
  125.     kaiser_init(6.5, 0., (char *)&kd);
  126. }
  127.  
  128. /*
  129.  * filt_find: return ptr to filter descriptor given filter name
  130.  */
  131.  
  132. Filt *filt_find(name)
  133. char *name;
  134. {
  135.     int i;
  136.  
  137.     if (nfilt==0) filt_init();
  138.     for (i=0; i<nfilt; i++)
  139.     if (str_eq(name, filt[i].name))
  140.         return &filt[i];
  141.     return 0;
  142. }
  143.  
  144. /*
  145.  * filt_insert: insert filter f in filter collection
  146.  */
  147.  
  148. void filt_insert(f)
  149. Filt *f;
  150. {
  151.     if (nfilt==0) filt_init();
  152.     if (filt_find(f->name)!=0) {
  153.     fprintf(stderr, "filt_insert: there's already a filter called %s\n",
  154.         f->name);
  155.     return;
  156.     }
  157.     if (nfilt>=NFILTMAX) {
  158.     fprintf(stderr, "filt_insert: too many filters: %d\n", nfilt+1);
  159.     return;
  160.     }
  161.     filt[nfilt++] = *f;
  162. }
  163.  
  164. /*
  165.  * filt_catalog: print a filter catalog to stdout
  166.  */
  167.  
  168. void filt_catalog()
  169. {
  170.     int i;
  171.  
  172.     if (nfilt==0) filt_init();
  173.     for (i=0; i<nfilt; i++)
  174.     filt_print(&filt[i]);
  175. }
  176.  
  177. /*
  178.  * filt_print: print info about a filter to stdout
  179.  */
  180.  
  181. void filt_print(f)
  182. Filt *f;
  183. {
  184.     fprintf(stderr,"%-9s\t%4.2f%s",
  185.     f->name, f->supp, f->windowme ? " (windowed by default)" : "");
  186.     if (f->printproc) {
  187.     fprintf(stderr,"\n    ");
  188.     filt_print_client(f);
  189.     }
  190.     fprintf(stderr,"\n");
  191. }
  192.  
  193. /*-------------------- windowing a filter --------------------*/
  194.  
  195. /*
  196.  * filt_window: given an IIR filter f and the name of a window function,
  197.  * create a compound filter that is the product of the two:
  198.  * wf->func(x) = f->func(x) * w->func(x/s)
  199.  *
  200.  * note: allocates memory that is (probably) never freed
  201.  */
  202.  
  203. Filt *filt_window(f, windowname)
  204. Filt *f;
  205. char *windowname;
  206. {
  207.     Filt *w, *wf;
  208.     window_data *d;
  209.  
  210.     if (str_eq(windowname, "box")) return f;    /* window with box is NOP */
  211.     w = filt_find(windowname);
  212.     ALLOC(wf, Filt, 1);
  213.     *wf = *f;
  214.     ALLOC(wf->name, char, 50);
  215.     sprintf(wf->name, "%s*%s", f->name, w->name);
  216.     wf->func = window_func;
  217.     wf->initproc = 0;
  218.     if (f->printproc || w->printproc) wf->printproc = window_print;
  219.     else wf->printproc = 0;
  220.     ALLOC(d, window_data, 1);
  221.     d->filter = f;
  222.     d->window = w;
  223.     wf->clientdata = (char *)d;
  224.     return wf;
  225. }
  226.  
  227. static double window_func(x, d)
  228. double x;
  229. char *d;
  230. {
  231.     register window_data *w;
  232.  
  233.     w = (window_data *)d;
  234. #   ifdef DEBUG
  235.     fprintf(stderr, "%s*%s(%g) = %g*%g = %g\n",
  236.     w->filter->name, w->window->name, x);
  237.     filt_func(w->filter, x), filt_func(w->window, x/w->filter->supp),
  238.     filt_func(w->filter, x) * filt_func(w->window, x/w->filter->supp));
  239. #   endif
  240.     return filt_func(w->filter, x) * filt_func(w->window, x/w->filter->supp);
  241. }
  242.  
  243. static void window_print(d)
  244. char *d;
  245. {
  246.     window_data *w;
  247.  
  248.     w = (window_data *)d;
  249.     if (w->filter->printproc) filt_print_client(w->filter);
  250.     if (w->window->printproc) filt_print_client(w->window);
  251. }
  252.  
  253. /*--------------- unit-area filters for unit-spaced samples ---------------*/
  254.  
  255. /* all filters centered on 0 */
  256.  
  257. double filt_box(x, d)        /* box, pulse, Fourier window, */
  258. double x;            /* 1st order (constant) b-spline */
  259. char *d;
  260. {
  261.     if (x<-.5) return 0.;
  262.     if (x<.5) return 1.;
  263.     return 0.;
  264. }
  265.  
  266. double filt_triangle(x, d)    /* triangle, Bartlett window, */
  267. double x;            /* 2nd order (linear) b-spline */
  268. char *d;
  269. {
  270.     if (x<-1.) return 0.;
  271.     if (x<0.) return 1.+x;
  272.     if (x<1.) return 1.-x;
  273.     return 0.;
  274. }
  275.  
  276. double filt_quadratic(x, d)    /* 3rd order (quadratic) b-spline */
  277. double x;
  278. char *d;
  279. {
  280.     double t;
  281.  
  282.     if (x<-1.5) return 0.;
  283.     if (x<-.5) {t = x+1.5; return .5*t*t;}
  284.     if (x<.5) return .75-x*x;
  285.     if (x<1.5) {t = x-1.5; return .5*t*t;}
  286.     return 0.;
  287. }
  288.  
  289. double filt_cubic(x, d)        /* 4th order (cubic) b-spline */
  290. double x;
  291. char *d;
  292. {
  293.     double t;
  294.  
  295.     if (x<-2.) return 0.;
  296.     if (x<-1.) {t = 2.+x; return t*t*t/6.;}
  297.     if (x<0.) return (4.+x*x*(-6.+x*-3.))/6.;
  298.     if (x<1.) return (4.+x*x*(-6.+x*3.))/6.;
  299.     if (x<2.) {t = 2.-x; return t*t*t/6.;}
  300.     return 0.;
  301. }
  302.  
  303. double filt_catrom(x, d)    /* Catmull-Rom spline, Overhauser spline */
  304. double x;
  305. char *d;
  306. {
  307.     if (x<-2.) return 0.;
  308.     if (x<-1.) return .5*(4.+x*(8.+x*(5.+x)));
  309.     if (x<0.) return .5*(2.+x*x*(-5.+x*-3.));
  310.     if (x<1.) return .5*(2.+x*x*(-5.+x*3.));
  311.     if (x<2.) return .5*(4.+x*(-8.+x*(5.-x)));
  312.     return 0.;
  313. }
  314.  
  315. double filt_gaussian(x, d)    /* Gaussian (infinite) */
  316. double x;
  317. char *d;
  318. {
  319.     return exp(-2.*x*x)*sqrt(2./PI);
  320. }
  321.  
  322. double filt_sinc(x, d)        /* Sinc, perfect lowpass filter (infinite) */
  323. double x;
  324. char *d;
  325. {
  326.     return x==0. ? 1. : sin(PI*x)/(PI*x);
  327. }
  328.  
  329. double filt_bessel(x, d)    /* Bessel (for circularly symm. 2-d filt, inf)*/
  330. double x;
  331. char *d;
  332. {
  333.     /*
  334.      * See Pratt "Digital Image Processing" p. 97 for Bessel functions
  335.      * zeros are at approx x=1.2197, 2.2331, 3.2383, 4.2411, 5.2428, 6.2439,
  336.      * 7.2448, 8.2454
  337.      */
  338.     return x==0. ? PI/4. : j1(PI*x)/(2.*x);
  339. }
  340.  
  341. /*-------------------- parameterized filters --------------------*/
  342.  
  343. double filt_mitchell(x, d)    /* Mitchell & Netravali's two-param cubic */
  344. double x;
  345. char *d;
  346. {
  347.     register mitchell_data *m;
  348.  
  349.     /*
  350.      * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
  351.      * SIGGRAPH 88
  352.      */
  353.     m = (mitchell_data *)d;
  354.     if (x<-2.) return 0.;
  355.     if (x<-1.) return m->q0-x*(m->q1-x*(m->q2-x*m->q3));
  356.     if (x<0.) return m->p0+x*x*(m->p2-x*m->p3);
  357.     if (x<1.) return m->p0+x*x*(m->p2+x*m->p3);
  358.     if (x<2.) return m->q0+x*(m->q1+x*(m->q2+x*m->q3));
  359.     return 0.;
  360. }
  361.  
  362. static void mitchell_init(b, c, d)
  363. double b, c;
  364. char *d;
  365. {
  366.     mitchell_data *m;
  367.  
  368.     m = (mitchell_data *)d;
  369.     m->p0 = (  6. -  2.*b        ) / 6.;
  370.     m->p2 = (-18. + 12.*b +  6.*c) / 6.;
  371.     m->p3 = ( 12. -  9.*b -  6.*c) / 6.;
  372.     m->q0 = (         8.*b + 24.*c) / 6.;
  373.     m->q1 = (      - 12.*b - 48.*c) / 6.;
  374.     m->q2 = (         6.*b + 30.*c) / 6.;
  375.     m->q3 = (     -     b -  6.*c) / 6.;
  376. }
  377.  
  378. static void mitchell_print(d)
  379. char *d;
  380. {
  381.     mitchell_data *m;
  382.  
  383.     m = (mitchell_data *)d;
  384.     fprintf(stderr,"mitchell: p0=%g p2=%g p3=%g q0=%g q1=%g q2=%g q3=%g\n",
  385.     m->p0, m->p2, m->p3, m->q0, m->q1, m->q2, m->q3);
  386. }
  387.  
  388. /*-------------------- window functions --------------------*/
  389.  
  390. double filt_hanning(x, d)    /* Hanning window */
  391. double x;
  392. char *d;
  393. {
  394.     return .5+.5*cos(PI*x);
  395. }
  396.  
  397. double filt_hamming(x, d)    /* Hamming window */
  398. double x;
  399. char *d;
  400. {
  401.     return .54+.46*cos(PI*x);
  402. }
  403.  
  404. double filt_blackman(x, d)    /* Blackman window */
  405. double x;
  406. char *d;
  407. {
  408.     return .42+.50*cos(PI*x)+.08*cos(2.*PI*x);
  409. }
  410.  
  411. /*-------------------- parameterized windows --------------------*/
  412.  
  413. double filt_kaiser(x, d)    /* parameterized Kaiser window */
  414. double x;
  415. char *d;
  416. {
  417.     /* from Oppenheim & Schafer, Hamming */
  418.     kaiser_data *k;
  419.  
  420.     k = (kaiser_data *)d;
  421.     return bessel_i0(k->a*sqrt(1.-x*x))*k->i0a;
  422. }
  423.  
  424. static void kaiser_init(a, b, d)
  425. double a, b;
  426. char *d;
  427. {
  428.     kaiser_data *k;
  429.  
  430.     k = (kaiser_data *)d;
  431.     k->a = a;
  432.     k->i0a = 1./bessel_i0(a);
  433. }
  434.  
  435. static void kaiser_print(d)
  436. char *d;
  437. {
  438.     kaiser_data *k;
  439.  
  440.     k = (kaiser_data *)d;
  441.     fprintf(stderr,"kaiser: a=%g i0a=%g\n", k->a, k->i0a);
  442. }
  443.  
  444. double bessel_i0(x)
  445. double x;
  446. {
  447.     /*
  448.      * modified zeroth order Bessel function of the first kind.
  449.      * Are there better ways to compute this than the power series?
  450.      */
  451.     register int i;
  452.     double sum, y, t;
  453.  
  454.     sum = 1.;
  455.     y = x*x/4.;
  456.     t = y;
  457.     for (i=2; t>EPSILON; i++) {
  458.     sum += t;
  459.     t *= (double)y/(i*i);
  460.     }
  461.     return sum;
  462. }
  463.  
  464. /*--------------- filters for non-unit spaced samples ---------------*/
  465.  
  466. double filt_normal(x, d)    /* normal distribution (infinite) */
  467. double x;
  468. char *d;
  469. {
  470.     /*
  471.      * normal distribution: has unit area, but it's not for unit spaced samples
  472.      * Normal(x) = Gaussian(x/2)/2
  473.      */
  474.     return exp(-x*x/2.)/sqrt(2.*PI);
  475. }
  476.