Path: utzoo!attcan!uunet!tut.cis.ohio-state.edu!ucbvax!pasteur!miro.Berkeley.EDU!ph From: ph@miro.Berkeley.EDU (Paul Heckbert) Newsgroups: comp.graphics Subject: source code for filtered image zoom, part 2/3 Message-ID: <16204@pasteur.Berkeley.EDU> Date: 11 Aug 89 02:44:04 GMT Sender: news@pasteur.Berkeley.EDU Reply-To: ph@miro.Berkeley.EDU (Paul Heckbert) Organization: University of California at Berkeley Lines: 1594 # to unpack, cut here and run the following shell archive through sh # contents: zoom/scanline.h zoom/scanline.c zoom/zoom.h zoom/zoom_main.c # zoom/filt.c zoom/filt.h libpic/dump.c # echo extracting zoom/scanline.h sed 's/^X//' <<'EOF13099' >zoom/scanline.h X/* scanline.h: definitions for generic scanline data type and routines */ X X#ifndef SCANLINE_HDR X#define SCANLINE_HDR X X/* $Header: scanline.h,v 3.0 88/10/10 13:52:04 ph Locked $ */ X X#include X X/* scanline type is composed of #bytes per channel and #channels */ X X/* alternatives for number of bytes per channel */ X#define PIXEL_BYTESMASK 3 X#define PIXEL1 0 X#define PIXEL2 1 X#define PIXEL4 2 X X/* alternatives for number of channels per pixel */ X#define PIXEL_CHANMASK 4 X#define PIXEL_MONO 0 X#define PIXEL_RGB 4 X X/* X * for 1-byte rgb we use Pixel1_rgba, not Pixel1_rgb, since some compilers X * (e.g. mips) don't pad 3-byte structures to 4 bytes as most compilers do X */ Xtypedef Pixel1_rgba Scanline_rgb1; Xtypedef Pixel2_rgba Scanline_rgb2; Xtypedef Pixel4_rgb Scanline_rgb4; X Xtypedef struct { /* GENERIC SCANLINE */ X int type; /* scanline type: #bytes ored with #channels */ X int len; /* length of row array */ X int y; /* y coordinate of this scanline (if we care) */ X union { /* one of the following is valid dep. on type */ X Pixel1 *row1; X Pixel2 *row2; X Pixel4 *row4; X Scanline_rgb1 *row1_rgb; X Scanline_rgb2 *row2_rgb; X Scanline_rgb4 *row4_rgb; X } u; X} Scanline; X Xtypedef struct { /* SAMPLED FILTER WEIGHT TABLE */ X short i0, i1; /* range of samples is [i0..i1-1] */ X short *weight; /* weight[i] goes with pixel at i0+i */ X} Weighttab; X X#define CHANBITS 8 /* # bits per image channel */ X#define CHANWHITE 255 /* pixel value of white */ X X#endif EOF13099 echo extracting zoom/scanline.c sed 's/^X//' <<'EOF13100' >zoom/scanline.c X/* X * scanline: package of generic scanline operations X * Performs various filtering operations on scanlines of pixels. X * black&white or color, at 1, 2, or 4 bytes per pixel. X * Operations are currently implemented for only a subset of the possible data X * type combinations. X * X * Paul Heckbert 12 Sept 1988 X */ X Xstatic char rcsid[] = "$Header: scanline.c,v 3.0 88/10/10 13:51:11 ph Locked $"; X#include X X#include X#include X#include "scanline.h" X X/* some machines/compilers don't do unsigned chars correctly (PDP 11's) */ X/* #define BAD_UCHAR /* define this if your machine has bad uchars */ X#ifdef BAD_UCHAR X# define UCHAR(x) ((x)&255) X#else X# define UCHAR(x) (x) X#endif X X#define BYTES(buf) ((buf)->type&PIXEL_BYTESMASK) X#define CHAN(buf) ((buf)->type&PIXEL_CHANMASK) X#define TYPECOMB(type1, type2) ((type1)<<3 | (type2)) X#define PEG(val, t) (t = (val), t < 0 ? 0 : t > CHANWHITE ? CHANWHITE : t) X X#define YOU_LOSE(routine) { \ X fprintf(stderr, "scanline_%s: bad scanline type: %d\n", \ X routine, buf->type); \ X exit(1); \ X} X X#define YOU_BLEW_IT(routine) { \ X fprintf(stderr, "scanline_%s: bad type combination: %d & %d\n", \ X routine, abuf->type, bbuf->type); \ X exit(1); \ X} X X/* X * scanline_alloc: allocate memory for scanline buf of the given type and length X */ X Xscanline_alloc(buf, type, len) XScanline *buf; Xint type, len; X{ X buf->type = type; X buf->len = len; X switch (type) { X case PIXEL_MONO|PIXEL1: ALLOC(buf->u.row1, Pixel1, len); break; X case PIXEL_MONO|PIXEL2: ALLOC(buf->u.row2, Pixel2, len); break; X case PIXEL_MONO|PIXEL4: ALLOC(buf->u.row4, Pixel4, len); break; X case PIXEL_RGB |PIXEL1: X ALLOC(buf->u.row1_rgb, Scanline_rgb1, len); break; X case PIXEL_RGB |PIXEL2: X ALLOC(buf->u.row2_rgb, Scanline_rgb2, len); break; X case PIXEL_RGB |PIXEL4: X ALLOC(buf->u.row4_rgb, Scanline_rgb4, len); break; X default: YOU_LOSE("alloc"); X } X} X X/* X * scanline_free: free the memory used by buf (but not buf itself) X */ X Xscanline_free(buf) XScanline *buf; X{ X switch (buf->type) { X case PIXEL_MONO|PIXEL1: free(buf->u.row1); break; X case PIXEL_MONO|PIXEL2: free(buf->u.row2); break; X case PIXEL_MONO|PIXEL4: free(buf->u.row4); break; X case PIXEL_RGB |PIXEL1: free(buf->u.row1_rgb); break; X case PIXEL_RGB |PIXEL2: free(buf->u.row2_rgb); break; X case PIXEL_RGB |PIXEL4: free(buf->u.row4_rgb); break; X default: YOU_LOSE("free"); X } X} X X/* X * scanline_zero: zero a scanline X */ X Xscanline_zero(buf) XScanline *buf; X{ X switch (buf->type) { X case PIXEL_MONO|PIXEL4: X bzero(buf->u.row4, buf->len*sizeof(Pixel4)); X break; X case PIXEL_RGB |PIXEL4: X bzero(buf->u.row4_rgb, buf->len*sizeof(Scanline_rgb4)); X break; X default: X YOU_LOSE("zero"); X } X} X X/* X * scanline_read: read buf->len pixels into buf starting at (x0,y0) in picture p X */ X Xscanline_read(p, x0, y0, buf) XPic *p; Xint x0, y0; XScanline *buf; X{ X switch (buf->type) { X case PIXEL_MONO|PIXEL1: X pic_read_row(p, y0, x0, buf->len, buf->u.row1); X break; X case PIXEL_RGB |PIXEL1: X pic_read_row_rgba(p, y0, x0, buf->len, buf->u.row1_rgb); X break; X default: X YOU_LOSE("read"); X } X} X X/* X * scanline_write: write buf->len pixels from buf into pic p starting at (x0,y0) X */ X Xscanline_write(p, x0, y0, buf) XPic *p; Xint x0, y0; XScanline *buf; X{ X switch (buf->type) { X case PIXEL_MONO|PIXEL1: X pic_write_row(p, y0, x0, buf->len, buf->u.row1); X break; X case PIXEL_RGB |PIXEL1: X pic_write_row_rgba(p, y0, x0, buf->len, buf->u.row1_rgb); X break; X default: X YOU_LOSE("write"); X } X} X X/* X * scanline_filter: convolve abuf with the sampled filter in wtab, X * writing the result to bbuf after shifting down by shift bits X * length of arrays taken to be bbuf->len X * X * Note that the Pixel4->Pixel1 routines shift abuf down by 8 bits before X * multiplying, to avoid overflow. X */ X Xscanline_filter(shift, wtab, abuf, bbuf) Xint shift; XWeighttab *wtab; XScanline *abuf, *bbuf; X{ X switch (TYPECOMB(abuf->type, bbuf->type)) { X case TYPECOMB(PIXEL_MONO|PIXEL1, PIXEL_MONO|PIXEL2): X pixel12_filter(bbuf->len, shift, wtab, abuf->u.row1, bbuf->u.row2); X break; X case TYPECOMB(PIXEL_MONO|PIXEL4, PIXEL_MONO|PIXEL1): X pixel41_filter(bbuf->len, shift, wtab, abuf->u.row4, bbuf->u.row1); X break; X case TYPECOMB(PIXEL_RGB |PIXEL1, PIXEL_RGB |PIXEL2): X pixel12_rgb_filter(bbuf->len, shift, wtab, X abuf->u.row1_rgb, bbuf->u.row2_rgb); X break; X case TYPECOMB(PIXEL_RGB |PIXEL4, PIXEL_RGB |PIXEL1): X pixel41_rgb_filter(bbuf->len, shift, wtab, X abuf->u.row4_rgb, bbuf->u.row1_rgb); X break; X default: X YOU_BLEW_IT("filter"); X } X} X Xpixel12_filter(bn, shift, wtab, abuf, bbuf) Xint bn, shift; XWeighttab *wtab; XPixel1 *abuf; XPixel2 *bbuf; X{ X register int af, sum; X register Pixel1 *ap; X register short *wp; X int b; X X for (b=0; bweight, ap= &abuf[wtab->i0], X af=wtab->i1-wtab->i0; af>0; af--) X sum += *wp++ * (short)UCHAR(*ap++); X *bbuf++ = sum>>shift; X } X} X Xpixel41_filter(bn, shift, wtab, abuf, bbuf) Xint bn, shift; XWeighttab *wtab; XPixel4 *abuf; XPixel1 *bbuf; X{ X register int af, sum; X register Pixel4 *ap; X register short *wp; X int t, b; X X for (b=0; bweight, ap= &abuf[wtab->i0], X af=wtab->i1-wtab->i0; af>0; af--) X sum += *wp++ * (short)(*ap++ >> CHANBITS); X *bbuf++ = PEG(sum>>shift, t); X } X} X Xpixel12_rgb_filter(bn, shift, wtab, abuf, bbuf) Xint bn, shift; XWeighttab *wtab; XScanline_rgb1 *abuf; XScanline_rgb2 *bbuf; X{ X register Scanline_rgb1 *ap; X register short *wp; X register int sumr, sumg, sumb, af; X int b; X X for (b=0; bweight, ap= &abuf[wtab->i0], X af=wtab->i1-wtab->i0; af>0; af--, wp++, ap++) { X sumr += *wp * (short)UCHAR(ap->r); X sumg += *wp * (short)UCHAR(ap->g); X sumb += *wp * (short)UCHAR(ap->b); X } X bbuf->r = sumr>>shift; X bbuf->g = sumg>>shift; X bbuf->b = sumb>>shift; X } X} X Xpixel41_rgb_filter(bn, shift, wtab, abuf, bbuf) Xint bn, shift; XWeighttab *wtab; XScanline_rgb4 *abuf; XScanline_rgb1 *bbuf; X{ X register Scanline_rgb4 *ap; X register short *wp; X register int sumr, sumg, sumb, af; X int t, b; X X for (b=0; bweight, ap= &abuf[wtab->i0], X af=wtab->i1-wtab->i0; af>0; af--, wp++, ap++) { X sumr += *wp * (short)(ap->r >> CHANBITS); X sumg += *wp * (short)(ap->g >> CHANBITS); X sumb += *wp * (short)(ap->b >> CHANBITS); X } X bbuf->r = PEG(sumr>>shift, t); X bbuf->g = PEG(sumg>>shift, t); X bbuf->b = PEG(sumb>>shift, t); X } X} X X/* X * scanline_accum: bbuf += weight*accum X * length of arrays taken to be bbuf->len X */ X Xscanline_accum(weight, abuf, bbuf) Xint weight; XScanline *abuf, *bbuf; X{ X if (weight==0) return; X if (BYTES(bbuf)!=PIXEL4 || CHAN(abuf)!=CHAN(bbuf)) X YOU_BLEW_IT("accum"); X switch (abuf->type) { X case PIXEL_MONO|PIXEL1: X pixel14_accum(bbuf->len, weight, abuf->u.row1, bbuf->u.row4); X break; X case PIXEL_MONO|PIXEL2: X pixel24_accum(bbuf->len, weight, abuf->u.row2, bbuf->u.row4); X break; X case PIXEL_RGB |PIXEL1: X pixel14_rgb_accum(bbuf->len, weight, X abuf->u.row1_rgb, bbuf->u.row4_rgb); X break; X case PIXEL_RGB |PIXEL2: X pixel24_rgb_accum(bbuf->len, weight, X abuf->u.row2_rgb, bbuf->u.row4_rgb); X break; X default: X YOU_BLEW_IT("accum"); X } X} X Xpixel14_accum(n, weight, buf, accum) Xregister int n, weight; Xregister Pixel1 *buf; Xregister Pixel4 *accum; X{ X for (; n>0; n--) X *accum++ += (short)weight * (short)UCHAR(*buf++); X} X Xpixel24_accum(n, weight, buf, accum) Xregister int n, weight; Xregister Pixel2 *buf; Xregister Pixel4 *accum; X{ X for (; n>0; n--) X *accum++ += (short)weight * *buf++; X} X Xpixel14_rgb_accum(n, weight, buf, accum) Xregister int n, weight; Xregister Scanline_rgb1 *buf; Xregister Scanline_rgb4 *accum; X{ X for (; n>0; n--, accum++, buf++) { X accum->r += (short)weight * (short)UCHAR(buf->r); X accum->g += (short)weight * (short)UCHAR(buf->g); X accum->b += (short)weight * (short)UCHAR(buf->b); X } X} X Xpixel24_rgb_accum(n, weight, buf, accum) Xregister int n, weight; Xregister Scanline_rgb2 *buf; Xregister Scanline_rgb4 *accum; X{ X for (; n>0; n--, accum++, buf++) { X accum->r += (short)weight * buf->r; X accum->g += (short)weight * buf->g; X accum->b += (short)weight * buf->b; X } X} X X/* X * scanline_shift: bbuf = abuf>>shift X * length of arrays taken to be bbuf->len X */ X Xscanline_shift(shift, abuf, bbuf) Xint shift; XScanline *abuf, *bbuf; X{ X if (BYTES(abuf)!=PIXEL4 || BYTES(bbuf)!=PIXEL1 || CHAN(abuf)!=CHAN(bbuf)) X YOU_BLEW_IT("shift"); X switch (CHAN(abuf)) { X case PIXEL_MONO: X pixel41_shift(bbuf->len, shift, abuf->u.row4, bbuf->u.row1); X break; X case PIXEL_RGB: X pixel41_rgb_shift(bbuf->len, shift, X abuf->u.row4_rgb, bbuf->u.row1_rgb); X break; X } X} X Xpixel41_shift(n, shift, accum, bbuf) Xregister int n, shift; Xregister Pixel4 *accum; Xregister Pixel1 *bbuf; X{ X register int t, half; X X half = 1 << shift-1; X for (; n>0; n--) X *bbuf++ = PEG(*accum+++half >> shift, t); X} X Xpixel41_rgb_shift(n, shift, accum, bbuf) Xregister int n, shift; Xregister Scanline_rgb4 *accum; Xregister Scanline_rgb1 *bbuf; X{ X register int t, half; X X half = 1 << shift-1; X for (; n>0; n--, accum++, bbuf++) { X bbuf->r = PEG(accum->r+half >> shift, t); X bbuf->g = PEG(accum->g+half >> shift, t); X bbuf->b = PEG(accum->b+half >> shift, t); X } X} EOF13100 echo extracting zoom/zoom.h sed 's/^X//' <<'EOF13101' >zoom/zoom.h X/* $Header: zoom.h,v 3.0 88/10/10 13:52:11 ph Locked $ */ X Xtypedef struct { /* SOURCE TO DEST COORDINATE MAPPING */ X double sx, sy; /* x and y scales */ X double tx, ty; /* x and y translations */ X double ux, uy; /* x and y offset used by MAP, private fields */ X} Mapping; X X/* see explanation in zoom.c */ X Xextern int zoom_debug; Xextern int zoom_coerce; /* simplify filters if possible? */ Xextern int zoom_xy; /* filter x before y (1) or vice versa (0)? */ Xextern int zoom_trimzeros; /* trim zeros from filter weight table? */ EOF13101 echo extracting zoom/zoom_main.c sed 's/^X//' <<'EOF13102' >zoom/zoom_main.c X/* X * zoom_main: main program for image zooming X * X * see comments in zoom.c X * X * note: DEFAULT_FILE must be set in Makefile X */ X X#include X X#include X#include X#include "filt.h" X#include "zoom.h" X X#define UNDEF PIC_UNDEFINED X#define FILTER_DEFAULT "triangle" X#define WINDOW_DEFAULT "blackman" X Xdouble atof_check(); X Xstatic char usage[] = "\ XUsage: zoom [options]\n\ X-src %s source filename\n\ X-dst %s dest filename\n\ X-s %d%d%d%d source box (x0 y0 xsize ysize)\n\ X-d %d%d%d%d dest box\n\ X-sw %d%d%d%d window (x0 y0 x1 y1)\n\ X-dw %d%d%d%d dest window\n\ X-map %f%f%f%f scale and translate: sx sy tx ty (src to dst by default)\n\ X-square use square mapping (don't stretch pixels)\n\ X-intscale use integer scale factors\n\ X-filt %s[%s filter name in x and y (default=triangle)\n\ X \"-filt '?'\" prints a filter catalog\n\ X-supp %f[%f filter support radius\n\ X-blur %f[%f blur factor: >1 is blurry, <1 is sharp\n\ X-window %s[%s window an IIR filter (default=blackman)\n\ X-mono monochrome mode (1 channel)\n\ X-debug %d print filter coefficients\n\ X-xy filter x before y\n\ X-yx filter y before x\n\ X-plain disable filter coercion\n\ X-keep0 keep zeros in xfilter\n\ X-dev print list of known picture devices/file formats\n\ XWhere %d denotes integer, %f denotes float, %s denotes string,\n\ Xand '[' marks optional following args\n\ X"; X Xmain(ac, av) Xint ac; Xchar **av; X{ X char *xfiltname = FILTER_DEFAULT, *yfiltname = 0; X char *xwindowname = 0, *ywindowname = 0; X char *srcfile = DEFAULT_FILE; X char *dstfile = DEFAULT_FILE; X int xyflag, yxflag, mono, nocoerce, nchan, square, intscale, keepzeros, i; X double xsupp = -1., ysupp = -1.; X double xblur = -1., yblur = -1.; X Pic *apic, *bpic; X Window_box a; /* src window */ X Window_box b; /* dst window */ X Mapping m; X Filt *xfilt, *yfilt, xf, yf; X X a.x0 = b.x0 = UNDEF; X a.x1 = b.x1 = UNDEF; X m.sx = 0.; X square = 0; X intscale = 0; X mono = 0; X xyflag = 0; X yxflag = 0; X nocoerce = 0; X keepzeros = 0; X X for (i=1; i=0.) xfilt->supp = xsupp; X if (xsupp>=0. && ysupp<0.) ysupp = xsupp; X if (ysupp>=0.) yfilt->supp = ysupp; X if (xblur>=0.) xfilt->blur = xblur; X if (xblur>=0. && yblur<0.) yblur = xblur; X if (yblur>=0.) yfilt->blur = yblur; X X if (!ywindowname) ywindowname = xwindowname; X if (xwindowname || xfilt->windowme) { X if (!xwindowname) xwindowname = WINDOW_DEFAULT; X xfilt = filt_window(xfilt, xwindowname); X } X if (ywindowname || yfilt->windowme) { X if (!ywindowname) ywindowname = WINDOW_DEFAULT; X yfilt = filt_window(yfilt, ywindowname); X } X X if (xfilt->printproc) { X printf("xfilt: "); X filt_print_client(xfilt); X } X if (yfilt->printproc) { X printf("yfilt: "); X filt_print_client(yfilt); X } X X if (m.sx==0.) zoom_opt(apic, &a, bpic, &b, xfilt, yfilt, square, intscale); X else zoom_continuous(apic, &a, bpic, &b, &m, xfilt, yfilt); X X pic_close(apic); X if (bpic!=apic) pic_close(bpic); X} X Xok(enough, option) Xint enough; Xchar *option; X{ X if (!enough) { X fprintf(stderr, "insufficient args to %s\n", option); X exit(1); X } X return 1; X} X X/* atof_check: ascii to float conversion with checking */ X Xdouble atof_check(str) Xchar *str; X{ X char *s; X X for (s=str; *s; s++) X if (strchr("0123456789.+-eE", *s) == NULL) { X fprintf(stderr, "expected numeric argument, not %s\n", str); X exit(1); X } X return atof(str); X} EOF13102 echo extracting zoom/filt.c sed 's/^X//' <<'EOF13103' >zoom/filt.c X/* X * filt: package of 1-d signal filters, both FIR and IIR X * X * Paul Heckbert, ph@miro.berkeley.edu 23 Oct 1986, 10 Sept 1988 X * X * Copyright (c) 1989 Paul S. Heckbert X * This source may be used for peaceful, nonprofit purposes only, unless X * under licence from the author. This notice should remain in the source. X * X * X * Documentation: X * To use these routines, X * #include X * Then call filt_find to select the desired filter, e.g. X * Filt *f; X * f = filt_find("catrom"); X * This filter function (impulse response) is then evaluated by calling X * filt_func(f, x). Each filter is nonzero between -f->supp and f->supp. X * Typically one will use the filter something like this: X * double phase, x, weight; X * for (x=phase-f->supp; xsupp; x+=deltax) { X * weight = filt_func(f, x); X * # do something with weight X * } X * X * Example of windowing an IIR filter: X * f = filt_find("sinc"); X * # if you wanted to vary sinc width, set f->supp = desired_width here X * # e.g. f->supp = 6; X * f = filt_window(f, "blackman"); X * You can then use f as above. X */ X X/* X * references: X * X * A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975 X * X * R.W. Hamming, Digital Filters, Prentice-Hall, Englewood Cliffs, NJ, 1983 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 Xstatic char rcsid[] = "$Header: filt.c,v 2.2 88/12/30 21:32:17 ph Locked $"; X#include X X#include X#include "filt.h" X X#define EPSILON 1e-7 X Xtypedef struct { /* data for parameterized Mitchell filter */ X double p0, p2, p3; X double q0, q1, q2, q3; X} mitchell_data; X Xtypedef struct { /* data for parameterized Kaiser window */ X double a; /* = w*(N-1)/2 in Oppenheim&Schafer notation */ X double i0a; X /* X * typically 4name)!=0) { X fprintf(stderr, "filt_insert: there's already a filter called %s\n", X f->name); X return; X } X if (nfilt>=NFILTMAX) { X fprintf(stderr, "filt_insert: too many filters: %d\n", nfilt+1); X return; X } X filt[nfilt++] = *f; X} X X/* X * filt_catalog: print a filter catalog to stdout X */ X Xvoid filt_catalog() X{ X int i; X Filt *f; X X if (nfilt==0) filt_init(); X for (i=0; iname, f->supp, f->windowme ? " (windowed by default)" : ""); X if (f->printproc) { X printf("\n "); X filt_print_client(f); X } X printf("\n"); X} X X/*-------------------- windowing a filter --------------------*/ X X/* X * filt_window: given an IIR filter f and the name of a window function, X * create a compound filter that is the product of the two: X * wf->func(x) = f->func(x) * w->func(x/s) X * X * note: allocates memory that is (probably) never freed X */ X XFilt *filt_window(f, windowname) XFilt *f; Xchar *windowname; X{ X Filt *w, *wf; X window_data *d; X X if (str_eq(windowname, "box")) return f; /* window with box is NOP */ X w = filt_find(windowname); X ALLOC(wf, Filt, 1); X *wf = *f; X ALLOC(wf->name, char, 50); X sprintf(wf->name, "%s*%s", f->name, w->name); X wf->func = window_func; X wf->initproc = 0; X if (f->printproc || w->printproc) wf->printproc = window_print; X else wf->printproc = 0; X ALLOC(d, window_data, 1); X d->filter = f; X d->window = w; X wf->clientdata = (char *)d; X return wf; X} X Xstatic double window_func(x, d) Xdouble x; Xchar *d; X{ X register window_data *w; X X w = (window_data *)d; X# ifdef DEBUG X printf("%s*%s(%g) = %g*%g = %g\n", X w->filter->name, w->window->name, x); X filt_func(w->filter, x), filt_func(w->window, x/w->filter->supp), X filt_func(w->filter, x) * filt_func(w->window, x/w->filter->supp)); X# endif X return filt_func(w->filter, x) * filt_func(w->window, x/w->filter->supp); X} X Xstatic void window_print(d) Xchar *d; X{ X window_data *w; X X w = (window_data *)d; X if (w->filter->printproc) filt_print_client(w->filter); X if (w->window->printproc) filt_print_client(w->window); X} X X/*--------------- unit-area filters for unit-spaced samples ---------------*/ X X/* all filters centered on 0 */ X Xdouble filt_box(x, d) /* box, pulse, Fourier window, */ Xdouble x; /* 1st order (constant) b-spline */ Xchar *d; X{ X if (x<-.5) return 0.; X if (x<.5) return 1.; X return 0.; X} X Xdouble filt_triangle(x, d) /* triangle, Bartlett window, */ Xdouble x; /* 2nd order (linear) b-spline */ Xchar *d; 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 filt_quadratic(x, d) /* 3rd order (quadratic) b-spline */ Xdouble x; Xchar *d; 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 filt_cubic(x, d) /* 4th order (cubic) b-spline */ Xdouble x; Xchar *d; 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 filt_catrom(x, d) /* Catmull-Rom spline, Overhauser spline */ Xdouble x; Xchar *d; 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 filt_gaussian(x, d) /* Gaussian (infinite) */ Xdouble x; Xchar *d; X{ X return exp(-2.*x*x)*sqrt(2./PI); X} X Xdouble filt_sinc(x, d) /* Sinc, perfect lowpass filter (infinite) */ Xdouble x; Xchar *d; X{ X return x==0. ? 1. : sin(PI*x)/(PI*x); X} X Xdouble filt_bessel(x, d) /* Bessel (for circularly symm. 2-d filt, inf)*/ Xdouble x; Xchar *d; X{ X /* X * See Pratt "Digital Image Processing" p. 97 for Bessel functions X * zeros are at approx x=1.2197, 2.2331, 3.2383, 4.2411, 5.2428, 6.2439, X * 7.2448, 8.2454 X */ X return x==0. ? PI/4. : j1(PI*x)/(2.*x); X} X X/*-------------------- parameterized filters --------------------*/ X Xdouble filt_mitchell(x, d) /* Mitchell & Netravali's two-param cubic */ Xdouble x; Xchar *d; X{ X register mitchell_data *m; X X /* X * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics", X * SIGGRAPH 88 X */ X m = (mitchell_data *)d; X if (x<-2.) return 0.; X if (x<-1.) return m->q0-x*(m->q1-x*(m->q2-x*m->q3)); X if (x<0.) return m->p0+x*x*(m->p2-x*m->p3); X if (x<1.) return m->p0+x*x*(m->p2+x*m->p3); X if (x<2.) return m->q0+x*(m->q1+x*(m->q2+x*m->q3)); X return 0.; X} X Xstatic void mitchell_init(b, c, d) Xdouble b, c; Xchar *d; X{ X mitchell_data *m; X X m = (mitchell_data *)d; X m->p0 = ( 6. - 2.*b ) / 6.; X m->p2 = (-18. + 12.*b + 6.*c) / 6.; X m->p3 = ( 12. - 9.*b - 6.*c) / 6.; X m->q0 = ( 8.*b + 24.*c) / 6.; X m->q1 = ( - 12.*b - 48.*c) / 6.; X m->q2 = ( 6.*b + 30.*c) / 6.; X m->q3 = ( - b - 6.*c) / 6.; X} X Xstatic void mitchell_print(d) Xchar *d; X{ X mitchell_data *m; X X m = (mitchell_data *)d; X printf("mitchell: p0=%g p2=%g p3=%g q0=%g q1=%g q2=%g q3=%g\n", X m->p0, m->p2, m->p3, m->q0, m->q1, m->q2, m->q3); X} X X/*-------------------- window functions --------------------*/ X Xdouble filt_hanning(x, d) /* Hanning window */ Xdouble x; Xchar *d; X{ X return .5+.5*cos(PI*x); X} X Xdouble filt_hamming(x, d) /* Hamming window */ Xdouble x; Xchar *d; X{ X return .54+.46*cos(PI*x); X} X Xdouble filt_blackman(x, d) /* Blackman window */ Xdouble x; Xchar *d; X{ X return .42+.50*cos(PI*x)+.08*cos(2.*PI*x); X} X X/*-------------------- parameterized windows --------------------*/ X Xdouble filt_kaiser(x, d) /* parameterized Kaiser window */ Xdouble x; Xchar *d; X{ X /* from Oppenheim & Schafer, Hamming */ X kaiser_data *k; X X k = (kaiser_data *)d; X return bessel_i0(k->a*sqrt(1.-x*x))*k->i0a; X} X Xstatic void kaiser_init(a, b, d) Xdouble a, b; Xchar *d; X{ X kaiser_data *k; X X k = (kaiser_data *)d; X k->a = a; X k->i0a = 1./bessel_i0(a); X} X Xstatic void kaiser_print(d) Xchar *d; X{ X kaiser_data *k; X X k = (kaiser_data *)d; X printf("kaiser: a=%g i0a=%g\n", k->a, k->i0a); X} X Xdouble bessel_i0(x) Xdouble x; X{ X /* X * modified zeroth order Bessel function of the first kind. X * Are there better ways to compute this than the power series? X */ X register int i; X double sum, y, t; X X sum = 1.; X y = x*x/4.; X t = y; X for (i=2; t>EPSILON; i++) { X sum += t; X t *= (double)y/(i*i); X } X return sum; X} X X/*--------------- filters for non-unit spaced samples ---------------*/ X Xdouble filt_normal(x, d) /* normal distribution (infinite) */ Xdouble x; Xchar *d; 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 return exp(-x*x/2.)/sqrt(2.*PI); X} EOF13103 echo extracting zoom/filt.h sed 's/^X//' <<'EOF13104' >zoom/filt.h X/* filt.h: definitions for filter data types and routines */ X X#ifndef FILT_HDR X#define FILT_HDR X X/* $Header: filt.h,v 2.2 88/12/30 21:32:20 ph Locked $ */ X Xtypedef struct { /* A 1-D FILTER */ X char *name; /* name of filter */ X double (*func)(); /* filter function */ X double supp; /* radius of nonzero portion */ X double blur; /* blur factor (1=normal) */ X char windowme; /* should filter be windowed? */ X char cardinal; /* is this filter cardinal? X ie, does func(x) = (x==0) for integer x? */ X char unitrange; /* does filter stay within the range [0..1] */ X void (*initproc)(); /* initialize client data, if any */ X void (*printproc)(); /* print client data, if any */ X char *clientdata; /* client info to be passed to func */ X} Filt; X X#define filt_func(f, x) (*(f)->func)((x), (f)->clientdata) X#define filt_print_client(f) (*(f)->printproc)((f)->clientdata) X XFilt *filt_find(); XFilt *filt_window(); Xvoid filt_print(); Xvoid filt_catalog(); X X X/* the filter collection: */ X Xdouble filt_box(); /* box, pulse, Fourier window, */ Xdouble filt_triangle(); /* triangle, Bartlett window, */ Xdouble filt_quadratic(); /* 3rd order (quadratic) b-spline */ Xdouble filt_cubic(); /* 4th order (cubic) b-spline */ Xdouble filt_catrom(); /* Catmull-Rom spline, Overhauser spline */ Xdouble filt_mitchell(); /* Mitchell & Netravali's two-param cubic */ Xdouble filt_gaussian(); /* Gaussian (infinite) */ Xdouble filt_sinc(); /* Sinc, perfect lowpass filter (infinite) */ Xdouble filt_bessel(); /* Bessel (for circularly symm. 2-d filt, inf)*/ X Xdouble filt_hanning(); /* Hanning window */ Xdouble filt_hamming(); /* Hamming window */ Xdouble filt_blackman(); /* Blackman window */ Xdouble filt_kaiser(); /* parameterized Kaiser window */ X Xdouble filt_normal(); /* normal distribution (infinite) */ X X X/* support routines */ Xdouble bessel_i0(); X X#endif EOF13104 echo extracting libpic/dump.c sed 's/^X//' <<'EOF13105' >libpic/dump.c X/* X * dump: subroutine package to read and write my dump picture file format X * X * Paul Heckbert 1 Oct 1988 X */ X Xstatic char rcsid[] = "$Header: dump.c,v 2.1 88/11/01 21:09:43 ph Locked $"; X X#include X#include "dump.h" X X#ifdef vax /* swap bytes in header if little-endian */ X# define BYTESWAP X#endif X X#define UNDEF PIXEL_UNDEFINED X X#define CHECK_HEAD_UNWRITTEN(p, subrname) { \ X if ((p)->headwritten) { \ X fprintf(stderr, "%s: can't change state once writing commences\n", \ X subrname); \ X exit(1); \ X } \ X} X#define CHECK_HEAD_WRITTEN(p, subrname) \ X if (!(p)->headwritten) dump_write_head(p, subrname); else X XDump *dump_open_read(), *dump_open_write(); X XDump *dump_open(file, mode) Xchar *file, *mode; X{ X if (str_eq(mode, "r")) return dump_open_read(file); X if (str_eq(mode, "w")) return dump_open_write(file); X fprintf(stderr, "dump_open: can't do mode %s\n", mode); X exit(1); /*NOTREACHED*/ X} X Xstatic Dump *dump_open_write(file) Xchar *file; X{ X Dump *p; X X ALLOC(p, Dump, 1); X if ((p->fp = fopen(file, "w")) == NULL) { X fprintf(stderr, "dump_open_write: can't write %s\n", file); X free(p); X return 0; X } X p->h.magic = DUMP_MAGIC; X p->h.nchan = 3; X p->h.dx = UNDEF; X strcpy(p->name, file); X p->headsize = sizeof(Dump_head); X p->headwritten = 0; X p->curx = p->cury = 0; X return p; X} X Xstatic void dump_write_head(p, subrname) XDump *p; Xchar *subrname; X{ X if (p->h.dx==UNDEF) { X fprintf(stderr, "%s: size of %s is uninitialized\n", p->name); X exit(1); X } X# ifdef BYTESWAP X headswap(&p->h); X# endif X if (fwrite(&p->h, sizeof(Dump_head), 1, p->fp) != 1) { X fprintf(stderr, "%s: write error on %s\n", subrname, p->name); X exit(1); X } X printf("%s: %dx%d %d-chan\n", p->name, p->h.dx, p->h.dy, p->h.nchan); X p->headwritten = 1; X} X Xstatic Dump *dump_open_read(file) Xchar *file; X{ X int badhead; X Dump *p; X X ALLOC(p, Dump, 1); X if ((p->fp = fopen(file, "r")) == NULL) { X fprintf(stderr, "dump_open_read: can't find %s\n", file); X free(p); X return 0; X } X X badhead = fread(&p->h, sizeof(Dump_head), 1, p->fp) != 1; X# ifdef BYTESWAP X headswap(&p->h); X# endif X if (badhead || p->h.magic!=DUMP_MAGIC) { X fprintf(stderr, "dump_open_read: %s is not a dump file\n", file); X free(p); X return 0; X } X printf("%s: %dx%d %d-chan\n", file, p->h.dx, p->h.dy, p->h.nchan); X strcpy(p->name, file); X p->headsize = sizeof(Dump_head); X p->headwritten = 1; X p->curx = p->cury = 0; X return p; X} X X#ifdef BYTESWAP X Xstatic headswap(h) XDump_head *h; X{ X swap_short(&h->magic); X swap_short(&h->nchan); X swap_short(&h->dx); X swap_short(&h->dy); X} X X#endif X Xvoid dump_close(p) XDump *p; X{ X if (p->fp) fclose(p->fp); X free(p); X} X Xchar *dump_get_name(p) XDump *p; X{ X return p->name; X} X Xvoid dump_clear(p, pv) XDump *p; XPixel1 pv; X{ X fprintf(stderr, "dump_clear: unimplemented\n"); X} X Xvoid dump_clear_rgba(p, r, g, b, a) XDump *p; XPixel1 r, g, b, a; X{ X fprintf(stderr, "dump_clear_rgba: unimplemented\n"); X} X X/*-------------------- file writing routines --------------------*/ X Xvoid dump_set_nchan(p, nchan) XDump *p; Xint nchan; X{ X CHECK_HEAD_UNWRITTEN(p, "dump_set_nchan"); X if (nchan!=1 && nchan!=3) { X fprintf(stderr, "dump_set_nchan: can't handle nchan=%d\n", nchan); X exit(1); X } X p->h.nchan = nchan; X} X Xvoid dump_set_box(p, ox, oy, dx, dy) XDump *p; Xint ox, oy, dx, dy; X{ X CHECK_HEAD_UNWRITTEN(p, "dump_set_box"); X /* ignore ox, oy */ X p->h.dx = dx; X p->h.dy = dy; X} X Xvoid dump_write_pixel(p, x, y, pv) XDump *p; Xint x, y; XPixel1 pv; X{ X fprintf(stderr, "dump_write_pixel: unimplemented\n"); X} X Xvoid dump_write_pixel_rgba(p, x, y, r, g, b, a) XDump *p; Xint x, y; XPixel1 r, g, b, a; X{ X fprintf(stderr, "dump_write_pixel_rgba: unimplemented\n"); X} X Xvoid dump_write_row(p, y, x0, nx, buf) Xregister Dump *p; Xint y, x0, nx; XPixel1 *buf; X{ X CHECK_HEAD_WRITTEN(p, "dump_write_row"); X if (x0!=p->curx || y!=p->cury) X dump_jump_to_pixel(p, x0, y); X if (fwrite(buf, nx*sizeof(Pixel1), 1, p->fp) != 1) { X fprintf(stderr, "dump_write_row: write error on %s\n", p->name); X exit(1); X } X dump_advance(p, nx); X} X Xvoid dump_write_row_rgba(p, y, x0, nx, buf) Xregister Dump *p; Xint y, x0, nx; Xregister Pixel1_rgba *buf; X{ X register int x; X X CHECK_HEAD_WRITTEN(p, "dump_write_row_rgba"); X if (x0!=p->curx || y!=p->cury) X dump_jump_to_pixel(p, x0, y); X for (x=nx; --x>=0; buf++) { X putc(buf->r, p->fp); X putc(buf->g, p->fp); X putc(buf->b, p->fp); X } X dump_advance(p, nx); X} X X/*-------------------- file reading routines --------------------*/ X Xint dump_get_nchan(p) XDump *p; X{ X return p->h.nchan; X} X Xvoid dump_get_box(p, ox, oy, dx, dy) XDump *p; Xint *ox, *oy, *dx, *dy; X{ X if (p->h.dx==UNDEF) { X *ox = UNDEF; /* used by some programs (e.g. zoom) */ X *oy = UNDEF; X } X else { X *ox = 0; X *oy = 0; X } X *dx = p->h.dx; X *dy = p->h.dy; X} X XPixel1 dump_read_pixel(p, x, y) XDump *p; Xint x, y; X{ X fprintf(stderr, "dump_read_pixel: unimplemented\n"); X} X Xvoid dump_read_pixel_rgba(p, x, y, pv) XDump *p; Xint x, y; XPixel1_rgba *pv; X{ X fprintf(stderr, "dump_read_pixel_rgba: unimplemented\n"); X} X Xvoid dump_read_row(p, y, x0, nx, buf) Xregister Dump *p; Xint y, x0, nx; XPixel1 *buf; X{ X if (x0!=p->curx || y!=p->cury) X dump_jump_to_pixel(p, x0, y); X if (fread(buf, nx*sizeof(Pixel1), 1, p->fp) != 1) { X fprintf(stderr, "dump_read_row: read error on %s\n", p->name); X exit(1); X } X dump_advance(p, nx); X} X Xvoid dump_read_row_rgba(p, y, x0, nx, buf) Xregister Dump *p; Xint y, x0, nx; Xregister Pixel1_rgba *buf; X{ X register int x; X X if (x0!=p->curx || y!=p->cury) X dump_jump_to_pixel(p, x0, y); X for (x=nx; --x>=0; buf++) { X buf->r = getc(p->fp); X buf->g = getc(p->fp); X buf->b = getc(p->fp); X buf->a = PIXEL1_MAX; X } X dump_advance(p, nx); X} X Xvoid dump_jump_to_pixel(p, x, y) XDump *p; Xint x, y; X{ X /* fprintf(stderr, "jumping from (%d,%d) to (%d,%d) in %s\n", */ X /* p->curx, p->cury, x, y, p->name); */ X p->curx = x; X p->cury = y; X assert(fseek(p->fp, p->headsize+(y*p->h.dx+x)*p->h.nchan*sizeof(Pixel1), 0) == 0); X} X Xvoid dump_advance(p, nx) Xregister Dump *p; Xint nx; X{ X p->curx += nx; X if (p->curx >= p->h.dx) { X p->curx -= p->h.dx; X p->cury++; X } X} EOF13105 exit