[comp.graphics] source code for filtered image zoom, part 2/3

ph@miro.Berkeley.EDU (Paul Heckbert) (08/11/89)

# 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 <pixel.h>
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 <stdio.h>
X
X#include <simple.h>
X#include <pic.h>
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; b<bn; b++, wtab++) {
X	/* start sum at 1<<shift-1 for rounding */
X	for (sum=1 << shift-1, wp=wtab->weight, 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; b<bn; b++, wtab++) {
X	/* start sum at 1<<shift-1 for rounding */
X	for (sum=1 << shift-1, wp=wtab->weight, 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; b<bn; b++, wtab++, bbuf++) {
X	/* start sum at 1<<shift-1 for rounding */
X	sumr = sumg = sumb = 1 << shift-1;
X	for (wp=wtab->weight, 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; b<bn; b++, wtab++, bbuf++) {
X	/* start sum at 1<<shift-1 for rounding */
X	sumr = sumg = sumb = 1 << shift-1;
X	for (wp=wtab->weight, 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 <math.h>
X
X#include <simple.h>
X#include <pic.h>
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<ac; i++)
X	if (av[i][0]=='-')
X	    if (str_eq(av[i], "-src") && ok(i+1<ac, "-src"))
X		srcfile = av[++i];
X	    else if (str_eq(av[i], "-dst") && ok(i+1<ac, "-dst"))
X		dstfile = av[++i];
X	    else if (str_eq(av[i], "-s") && ok(i+4<ac, "-s")) {
X		a.x0 = atof_check(av[++i]);
X		a.y0 = atof_check(av[++i]);
X		a.nx = atof_check(av[++i]);
X		a.ny = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-d") && ok(i+4<ac, "-d")) {
X		b.x0 = atof_check(av[++i]);
X		b.y0 = atof_check(av[++i]);
X		b.nx = atof_check(av[++i]);
X		b.ny = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-sw") && ok(i+4<ac, "-sw")) {
X		a.x0 = atof_check(av[++i]);
X		a.y0 = atof_check(av[++i]);
X		a.x1 = atof_check(av[++i]);
X		a.y1 = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-dw") && ok(i+4<ac, "-dw")) {
X		b.x0 = atof_check(av[++i]);
X		b.y0 = atof_check(av[++i]);
X		b.x1 = atof_check(av[++i]);
X		b.y1 = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-map") && ok(i+4<ac, "-map")) {
X		m.sx = atof_check(av[++i]);
X		m.sy = atof_check(av[++i]);
X		m.tx = atof_check(av[++i]);
X		m.ty = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-square"))
X		square = 1;
X	    else if (str_eq(av[i], "-intscale"))
X		intscale = 1;
X	    else if (str_eq(av[i], "-filt") && ok(i+1<ac, "-filt")) {
X		xfiltname = av[++i];
X		if (i+1<ac && av[i+1][0]!='-') yfiltname = av[++i];
X	    }
X	    else if (str_eq(av[i], "-supp") && ok(i+1<ac, "-supp")) {
X		xsupp = atof_check(av[++i]);
X		if (i+1<ac && av[i+1][0]!='-') ysupp = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-blur") && ok(i+1<ac, "-blur")) {
X		xblur = atof_check(av[++i]);
X		if (i+1<ac && av[i+1][0]!='-') yblur = atof_check(av[++i]);
X	    }
X	    else if (str_eq(av[i], "-window") && ok(i+1<ac, "-window")) {
X		xwindowname = av[++i];
X		if (i+1<ac && av[i+1][0]!='-') ywindowname = av[++i];
X	    }
X	    else if (str_eq(av[i], "-mono"))
X		mono = 1;
X	    else if (str_eq(av[i], "-debug") && ok(i+1<ac, "-debug"))
X		zoom_debug = atof_check(av[++i]);
X	    else if (str_eq(av[i], "-xy"))
X		xyflag = 1;
X	    else if (str_eq(av[i], "-yx"))
X		yxflag = 1;
X	    else if (str_eq(av[i], "-plain"))
X		nocoerce = 1;
X	    else if (str_eq(av[i], "-keep0"))
X		keepzeros = 1;
X	    else if (str_eq(av[i], "-dev")) {
X		pic_catalog();
X		exit(0);
X	    }
X	    else {
X		if (!str_eq(av[i], "-"))
X		    fprintf(stderr, "unrecognized argument: %s\n", av[i]);
X		fputs(usage, stderr);
X		exit(1);
X	    }
X
X    if (str_eq(xfiltname, "?")) {
X	filt_catalog();
X	exit(0);
X    }
X    if (xyflag) zoom_xy = 1;
X    if (yxflag) zoom_xy = 0;
X    zoom_coerce = !nocoerce;
X    zoom_trimzeros = !keepzeros;
X
X    bpic = pic_open(dstfile, "w");
X    if (!bpic) {
X	fprintf(stderr, "can't get %s\n", dstfile);
X	exit(1);
X    }
X    apic = str_eq(srcfile, dstfile) ? bpic : pic_open(srcfile, "r");
X    if (!apic) {
X	fprintf(stderr, "can't get %s\n", srcfile);
X	exit(1);
X    }
X
X    /* how many channels are src and dst?  maybe src knows... */
X    nchan = pic_get_nchan(apic);
X    if (nchan==UNDEF) {
X	nchan = mono ? 1 : 3;
X	pic_set_nchan(apic, nchan);
X    }
X    pic_set_nchan(bpic, nchan);
X
X    /*
X     * set source and dest subwindows a and b of apic and bpic
X     * note: pic_get_window doesn't set nx, ny fields of Window_box struct
X     */
X    if (a.x0==UNDEF) pic_get_window(apic, (Window_box *)&a);
X    if (a.x1==UNDEF) window_box_set_max(&a);
X    if (b.x0==UNDEF) pic_get_window(bpic, (Window_box *)&b);
X    if (b.x1==UNDEF) window_box_set_max(&b);
X    /*
X     * nx and ny uninitialized at this point
X     * pic_get_window might return x0=UNDEF for bpic
X     */
X
X    if (!yfiltname) yfiltname = xfiltname;
X    xfilt = filt_find(xfiltname);
X    yfilt = filt_find(yfiltname);
X    if (!xfilt || !yfilt) {
X	fprintf(stderr, "can't find filters %s and %s\n",
X	    xfiltname, yfiltname);
X	exit(1);
X    }
X    /* copy the filters before modifying them */
X    xf = *xfilt; xfilt = &xf;
X    yf = *yfilt; yfilt = &yf;
X    if (xsupp>=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 <filt.h>
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; x<f->supp; 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 <math.h>
X
X#include <simple.h>
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 4<a<9
X     * param a trades off main lobe width (sharpness)
X     * for side lobe amplitude (ringing)
X     */
X} kaiser_data;
X
Xtypedef struct {	/* data for windowed function compound filter */
X    Filt *filter;
X    Filt *window;
X} window_data;
X
Xvoid mitchell_init(), mitchell_print();
Xvoid kaiser_init(), kaiser_print();
Xvoid window_print();
Xdouble window_func();
Xstatic mitchell_data md;
Xstatic kaiser_data kd;
X
X#define NFILTMAX 30
Xstatic int nfilt = 0;
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
Xstatic Filt filt[NFILTMAX] = {
X/*  NAME	FUNC		SUPP	BLUR WIN CARD UNIT  OPT... */
X   {"point",	filt_box,	0.0,	1.0,  0,  1,  1},
X   {"box",	filt_box,       0.5,	1.0,  0,  1,  1},
X   {"triangle",	filt_triangle,  1.0,	1.0,  0,  1,  1},
X   {"quadratic",filt_quadratic, 1.5,	1.0,  0,  0,  1},
X   {"cubic",	filt_cubic,     2.0,	1.0,  0,  0,  1},
X
X   {"catrom",	filt_catrom,    2.0,	1.0,  0,  1,  0},
X   {"mitchell",	filt_mitchell,	2.0,	1.0,  0,  0,  0,
X	mitchell_init, mitchell_print, (char *)&md},
X
X   {"gaussian",	filt_gaussian,  1.25,	1.0,  0,  0,  1},
X   {"sinc",	filt_sinc,      4.0,	1.0,  1,  1,  0},
X   {"bessel",	filt_bessel,    3.2383,	1.0,  1,  0,  0},
X
X   {"hanning",	filt_hanning,	1.0,	1.0,  0,  1,  1},
X   {"hamming",	filt_hamming,	1.0,	1.0,  0,  1,  1},
X   {"blackman",	filt_blackman,	1.0,	1.0,  0,  1,  1},
X   {"kaiser",	filt_kaiser,	1.0,	1.0,  0,  1,  1,
X	kaiser_init, kaiser_print, (char *)&kd},
X   {0}
X};
X
X/*-------------------- general filter routines --------------------*/
X
Xstatic filt_init()
X{
X    /* count the filters to set nfilt */
X    for (nfilt=0; nfilt<NFILTMAX && filt[nfilt].name; nfilt++);
X    /* better way?? */
X    mitchell_init(1./3., 1./3., (char *)&md);
X    kaiser_init(6.5, 0., (char *)&kd);
X}
X
X/*
X * filt_find: return ptr to filter descriptor given filter name
X */
X
XFilt *filt_find(name)
Xchar *name;
X{
X    int i;
X
X    if (nfilt==0) filt_init();
X    for (i=0; i<nfilt; i++)
X	if (str_eq(name, filt[i].name))
X	    return &filt[i];
X    return 0;
X}
X
X/*
X * filt_insert: insert filter f in filter collection
X */
X
Xvoid filt_insert(f)
XFilt *f;
X{
X    if (nfilt==0) filt_init();
X    if (filt_find(f->name)!=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; i<nfilt; i++)
X	filt_print(&filt[i]);
X}
X
X/*
X * filt_print: print info about a filter to stdout
X */
X
Xvoid filt_print(f)
XFilt *f;
X{
X    printf("%-9s\t%4.2f%s",
X	f->name, 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 <simple.h>
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