Path: utzoo!utgpu!water!watmath!clyde!cbosgd!ucbvax!degas.Berkeley.EDU!ph
From: ph@degas.Berkeley.EDU (Paul Heckbert)
Newsgroups: comp.graphics
Subject: Re: convolution kernel program wanted
Message-ID: <22210@ucbvax.BERKELEY.EDU>
Date: 17 Dec 87 11:41:59 GMT
Sender: usenet@ucbvax.BERKELEY.EDU
Reply-To: ph@degas.Berkeley.EDU (Paul Heckbert)
Organization: University of California, Berkeley
Lines: 248

# to unpack, cut here and run the following shell archive through sh
# contents: README filt1d.h filt1d.c demo.c Makefile
#
sed 's/^X//' <<'EOF10751' >README
X
X This is C source for a simple package of subroutines containing a few common
X 1-dimensional filter functions, both FIR (finite impulse response) and IIR
X (infinite impulse response).  These routines could form the foundation of
X an image processing library or a simple spline library.
X 
X The other standard IIR filters (Hamming, Hanning, Blackman, Kaiser, etc.)
X could easily be added to this package.
X
X All functions are double precision floating point.
X
X run "make" and "demo sinc" for a demo.
X
X
X Paul Heckbert, CS grad student
X 508-7 Evans Hall, UC Berkeley		UUCP: ucbvax!degas.berkeley.edu!ph
X Berkeley, CA 94720			ARPA: ph@degas.berkeley.edu
X
X Dec 87
EOF10751
sed 's/^X//' <<'EOF10752' >filt1d.h
X/* filt1d.h: header file for 1d filter package */
X
X#ifndef FILT1D_HDR
X#define FILT1D_HDR
X
Xtypedef struct {	/* FILTER DESCRIPTOR STRUCTURE */
X    char *name;		/* filter name */
X    double (*func)();	/* filter function */
X    double support;	/* radius of nonzero portion of function,
X			 * or for infinite filters, a convenient cutoff point */
X    int infinite;	/* is this filter infinite? */
X} Filt1d;
X
XFilt1d *Filt1dFind();
X
Xdouble Filt1dBox();
Xdouble Filt1dTriangle();
Xdouble Filt1dQuadratic();
Xdouble Filt1dCubic();
Xdouble Filt1dCatrom();
Xdouble Filt1dGaussian();
Xdouble Filt1dSinc();
Xdouble Filt1dBessel();
X
Xdouble Filt1dNormal();
X
X#endif
EOF10752
sed 's/^X//' <<'EOF10753' >filt1d.c
X/*
X * Filt1d: package of 1-d signal filters, both FIR and IIR
X *
X * Paul Heckbert	Dec 1987
X */
X/*
X * references:
X *
X * A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975
X *
X * W.K. Pratt, Digital Image Processing, John Wiley and Sons, 1978
X *
X * H.S. Hou, H.C. Andrews, "Cubic Splines for Image Interpolation and
X *	Digital Filtering", IEEE Trans. Acoustics, Speech, and Signal Proc.,
X *	vol. ASSP-26, no. 6, Dec. 1978, pp. 508-517
X */
X
X#include <math.h>
X#define PI 3.1415926535897932385
X#define STREQ(a, b) (strcmp(a, b)==0)
X#include "filt1d.h"
X
X/*
X * note: for the IIR (infinite impulse response) filters,
X * gaussian, sinc, etc, the support values given are arbitrary,
X * convenient cutoff points, while for the FIR (finite impulse response)
X * filters the support is finite and absolute.
X */
X
X/*  NAME	FILTERFUNC	SUPP   INF?  */
Xstatic Filt1d filt[] = {
X    "box",	Filt1dBox,       0.5,  0,
X    "triangle",	Filt1dTriangle,  1.0,  0,
X    "quadratic",Filt1dQuadratic, 1.5,  0,
X    "cubic",	Filt1dCubic,     2.0,  0,
X    "catrom",	Filt1dCatrom,    2.0,  0,
X    "gaussian",	Filt1dGaussian,  1.25, 1,
X    "sinc",	Filt1dSinc,      4.0,  1,
X    "bessel",	Filt1dBessel,    2.0,  1,
X
X    "normal",	Filt1dNormal,    2.5,  1,
X};
X#define NFILT (sizeof filt/sizeof filt[0])
X
X/*
X * Filt1dFind: return ptr to filter descriptor given filter name
X */
X
XFilt1d *Filt1dFind(name)
Xchar *name;
X{
X    int i;
X
X    for (i=0; i<NFILT; i++)
X	if (STREQ(name, filt[i].name))
X	    return &filt[i];
X    return 0;
X}
X
X/*----------------------------------------------------------------------*/
X
X/*
X * unit area filters for unit spaced samples
X */
X
Xdouble Filt1dBox(x)		/* box, pulse, Fourier window,
X					1st order (constant) b-spline */
Xdouble x;
X{
X    if (x<-.5) return 0.;
X    if (x<.5) return 1.;
X    return 0.;
X}
X
Xdouble Filt1dTriangle(x)	/* triangle, Bartlett window,
X					2nd order (linear) b-spline */
Xdouble x;
X{
X    if (x<-1.) return 0.;
X    if (x<0.) return 1.+x;
X    if (x<1.) return 1.-x;
X    return 0.;
X}
X
Xdouble Filt1dQuadratic(x)	/* 3rd order (quadratic) b-spline */
Xdouble x;
X{
X    double t;
X
X    if (x<-1.5) return 0.;
X    if (x<-.5) {t = x+1.5; return .5*t*t;}
X    if (x<.5) return .75-x*x;
X    if (x<1.5) {t = x-1.5; return .5*t*t;}
X    return 0.;
X}
X
Xdouble Filt1dCubic(x)		/* 4th order (cubic) b-spline */
Xdouble x;
X{
X    double t;
X
X    if (x<-2.) return 0.;
X    if (x<-1.) {t = 2.+x; return t*t*t/6.;}
X    if (x<0.) return (4.+x*x*(-6.+x*-3.))/6.;
X    if (x<1.) return (4.+x*x*(-6.+x*3.))/6.;
X    if (x<2.) {t = 2.-x; return t*t*t/6.;}
X    return 0.;
X}
X
Xdouble Filt1dCatrom(x)		/* Catmull-Rom spline, Overhauser spline */
Xdouble x;
X{
X    if (x<-2.) return 0.;
X    if (x<-1.) return .5*(4.+x*(8.+x*(5.+x)));
X    if (x<0.) return .5*(2.+x*x*(-5.+x*-3.));
X    if (x<1.) return .5*(2.+x*x*(-5.+x*3.));
X    if (x<2.) return .5*(4.+x*(-8.+x*(5.-x)));
X    return 0.;
X}
X
Xdouble Filt1dGaussian(x)	/* Gaussian (infinite) */
Xdouble x;
X{
X    return exp(-2.*x*x)*sqrt(2./PI);
X}
X
Xdouble Filt1dSinc(x)		/* Sinc, perfect lowpass filter (infinite) */
Xdouble x;
X{
X    return x==0. ? 1. : sin(PI*x)/(PI*x);
X}
X
X/* See Pratt "Digital Image Processing" p. 97 for Bessel functions */
X
Xdouble Filt1dBessel(x)		/* Bessel (for circularly symm. 2-d filt, inf)*/
Xdouble x;
X{
X    return x==0. ? PI/4. : j1(PI*x)/(2.*x);
X}
X
X/*----------------------------------------------------------------------*/
X
X/*
X * normal distribution: has unit area, but it's not for unit spaced samples
X * Normal(x) = Gaussian(x/2)/2
X */
X
Xdouble Filt1dNormal(x)
Xdouble x;
X{
X    return exp(-x*x/2.)/sqrt(2.*PI);
X}
EOF10753
sed 's/^X//' <<'EOF10754' >demo.c
X/* demo.c: simple program to demonstrate use of Filt1d package */
X
X#define DT .1		/* a small increment */
X#include <stdio.h>
X#include "filt1d.h"
X
Xmain(ac, av)
Xint ac;
Xchar **av;
X{
X    char *filtname;
X    double t, x;
X    Filt1d *fp;
X
X    if (ac!=2) {
X	fprintf(stderr, "Usage: demo <filtname>\n");
X	exit(1);
X    }
X    filtname = av[1];
X
X    fp = Filt1dFind(filtname);		/* get filter descriptor ptr */
X    if (!fp) {
X	fprintf(stderr, "filter %s unknown\n", filtname);
X	exit(1);
X    }
X    printf("%s is %sfinite, support diam=%g\n",
X	fp->name, fp->infinite ? "in" : "", fp->support*2.);
X
X    /* evaluate filter function over range of support at some increment */
X    for (t= -fp->support; t<=fp->support; t+=DT) {
X	x = (*fp->func)(t);
X	printf("%s(%6.3f) = %8.5f\n", fp->name, t, x);
X    }
X}
EOF10754
sed 's/^X//' <<'EOF10755' >Makefile
Xdemo: demo.o filt1d.o
X	cc -o demo demo.o filt1d.o -lm
EOF10755
exit