muffy@violet.berkeley.edu (12/14/87)
I am looking for a program that can generate various types of convolution kernels. If anyone has written such a program or knows where I could get a hold of one, I would appreciate the help. Thanks Kathy Murphy muffy@violet.berkeley.edu
cc4b+@andrew.cmu.edu (Christopher Brian Cox) (12/14/87)
Since there are several requests (and I suspect even more people who are interested) for image processing programs, maybe someone to post them to the net. (comp.sources.?) -Chris cc4b@andrew.cmu.edu
ph@degas.Berkeley.EDU (Paul Heckbert) (12/17/87)
# 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
exitawpaeth@watcgl.waterloo.edu (Alan W. Paeth) (12/18/87)
Subject: Convolution Kernels
/*
* Kaiser window weights/coefficients for rectangular convolution kernels
* Note: this is a 2-D radial distribution, though the reference mentions that
* a rectangular window (Cartesian product of two 1-D filters) is not bad
*
* Kernel coefficients might be negative -- simple-minded filtering software
* take note!
*
* Alan Paeth, University of Waterloo, 1985
* Taken TK5102.5 T93 (book, US Library of Congress)
*/
#include <stdio.h>
#include <math.h>
main(argc, argv)
char **argv;
{
int w, h;
float a, t, x, y, r, z;
if (argc != 5) usage();
w = atoi(argv[1]);
h = atoi(argv[2]);
a = atof(argv[3]);
t = atof(argv[4]);
for(x = -((w-1.0)/2.0); x <= ((w-1.0)/2.0); x += 1.0)
{
for(y = -((h-1.0)/2.0); y <= ((h-1.0)/2.0); y += 1.0)
{
r = sqrt(x*x + y*y);
z = (r>t) ? 0 : j0(a*sqrt(1.0-(r/t)*(r/t)))/j0(a); /* Bessel fn. */
printf("%f ", z);
}
printf("\n");
}
exit(0);
}
usage()
{
fprintf(stderr,"usage: kaiser [width] [height] [slope] [radius]\n");
exit(1);
}shep@datacube.UUCP (12/21/87)
I wrote a program a few years ago for Datacube's VFIR-MKII hardware convolver. Called "vm2", it allows 2-D convolution kernels to be manipulated and viewed in real time. Its novelty is a kernel stack that allows kernels to operate on kernels. For example, the effect of two filters in series can be reduced to a single kernel that is the result of the two kernels convolved with each other. Our applications department is now releasing the program as an application note for that product. Try fernando@datacube.com for more information. Shep Siegel Datacube Inc. DSP Products Group 4 Dearborn Rd. Peabody, Ma. 01960 UUCP: shep@datacube.COM VOICE: (617) 535-6644; FAX: (617) 535-5643; TWX: (710) 347-0125