ph@miro.Berkeley.EDU (Paul Heckbert) (08/11/89)
# to unpack, cut here and run the following shell archive through sh
# contents: README Makefile zoom/zoom.1 zoom/zoom.c zoom/Makefile
#
mkdir libsys libpic zoom
echo extracting README
sed 's/^X//' <<'EOF13085' >README
XThis is C source code for a program to do filtered zoom (resize) of 8
Xand 24-bit raster images. The program supports arbitrary floating
Xpoint scale and translation with independent control of x and y scales,
Xupsampling (interpolation) or downsampling (decimation), 1-channel
Xblack and white or 3-channel color pictures, overlapping source and
Xdestination windows within the same frame buffer, and separable
Xfiltering with a large choice of filters. It can also be used for
Xin-place separable filtering of images.
X
XThis directory contains:
X libsys:
X one header file
X libpic:
X Subroutine package for device-independent picture I/O.
X Modify this to add picture I/O capability for your own frame buffers
X and file formats to the pic package and thereby to the zoom program.
X Read the man page pic.3 to learn how.
X zoom:
X Source for the main program zoom, and a set of useful filter routines.
X This code is graphics device independent, so you shouldn't have
X to modify it.
X
XRun "make install" and it should create subdirectories include, lib, and bin,
Xleaving the zoom program in bin/zoom. Note that you shouldn't expect this
Xto compile as is, as it's currently configured to link in the
XSilicon Graphics iris library. Again, see the notes under
X"CREATING A DEVICE LIBRARY" in libpic/pic.3
X
XPlease email me any enhancements, such as new device libraries for libpic;
XI'll collect and redistribute. (I already have libpic device libraries for
Xmemory frame buffers, various MIT, NYIT, Pixar, Xerox PARC, and UC Berkeley
Xpicture file formats, Sunraster format, Utah raster toolkit RLE format,
Xa fairly standard dump format, and hp98721 and iris frame buffers,
Xwhich I may distribute in the future.)
X
XPaul Heckbert 10 Aug 1989
XUC Berkeley
Xph@miro.berkeley.edu
EOF13085
echo extracting Makefile
sed 's/^X//' <<'EOF13086' >Makefile
Xinstall:
X -mkdir include lib bin
X cd libsys; make install
X cd libpic; make install
X cd zoom; make install
X
Xclean:
X -rm -fr include lib bin
X cd libsys; make clean
X cd libpic; make clean
X cd zoom; make clean
EOF13086
echo extracting zoom/zoom.1
sed 's/^X//' <<'EOF13087' >zoom/zoom.1
X.\" $Header$
X.\" a few macros
X.de Cs \" code start
X.DS
X.ft B
X.ta 9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n
X..
X.de Ce \" code end
X.ft R
X.DE
X..
X.de Ss \" subroutine definition start
X.nf
X.ft B
X.ta 9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n,+9n
X..
X.de Sd \" subroutine documentation
X.ft R
X.fi
X.in +4n
X..
X.de Se \" subroutine definition&documentation end
X.in -4n
X..
X.de DS
X.nf
X.in +4n
X.sp .5v
X..
X.de DE
X.sp .5v
X.in -4n
X.fi
X..
X.TH ZOOM 1 "10 August 1989"
X.SH NAME
Xzoom \- filtered image scale and translate
X.SH SYNOPSIS
X.DS
XUsage: zoom [options]
X-src %s source filename
X-dst %s dest filename
X-s %d%d%d%d source box (x0 y0 xsize ysize)
X-d %d%d%d%d dest box
X-sw %d%d%d%d window (x0 y0 x1 y1)
X-dw %d%d%d%d dest window
X-map %f%f%f%f scale and translate: sx sy tx ty (src to dst by default)
X-square use square mapping (don't stretch pixels)
X-intscale use integer scale factors
X-filt %s[%s filter name in x and y (default=triangle)
X-supp %f[%f filter support radius
X-blur %f[%f blur factor: >1 is blurry, <1 is sharp
X-window %s[%s window an IIR filter (default=blackman)
X-mono monochrome mode (1 channel)
X-dev print list of known picture devices/file formats
XWhere %d denotes integer, %f denotes float, %s denotes string,
Xand '[' marks optional following args
X.DE
X.SH DESCRIPTION
X\fIzoom\fP zooms a raster image
Xfrom one picture file or frame buffer to another.
XThe "zoom" operation, also known as "resize",
Xconsists of a scale and translation.
XThe program supports arbitrary floating point scale and translation
Xwith independent control of x and y scales,
Xupsampling (interpolation) or downsampling (decimation),
X1-channel black and white or 3-channel color pictures,
Xoverlapping source and destination windows within the same frame buffer,
Xand separable filtering with a large choice of filters.
XIt can also be used for in-place separable filtering of images.
XThe algorithm is not limited to integer or rational scale factors;
Xit can scale by any floating point number.
XThe program uses the \fIpic\fP package for all picture I/O,
Xso the \fIzoom\fP source code is device-independent,
Xand the program can read and write any of the picture file formats and
Xframe buffers known by \fIpic\fP.
XThe program is optimized to be nearly as fast as special-purpose code
Xfor important special cases such as point sampling.
XThe memory used by the algorithm is modest: it is proportional to the
Xpicture width times the filter width,
Xnot proportional to picture area.
X.PP
XTo run \fIzoom\fP, the user specifies the source and destination file
Xor frame buffer names, specifies the mapping (zoom transformation)
Xand filter.
XMost of these options have reasonable defaults.
XThe program will read and write the full-screen area of the source and
Xdestination pictures by default;
Xthe user can specify subrectangles to read and write if desired.
XRectangles can be defined in either of two ways: as a "box",
Xwhich consists of \fIxmin, ymin, xsize, ysize\fP,
Xor as a "window", which consists of \fIxmin, ymin, xmax, ymax\fP.
XBy default, the source window is zoomed into the destination window.
XAlternatively, a mapping can be specified directly with the \fB-map\fP option.
XA mapping can be constrained to have equal scale in x and y, or to have
Xinteger scale factors, using the \fB-square\fP and \fB-intscale\fP options,
Xrespectively.
X.PP
XFilters are selected with the \fB-filt\fP option.
XIf given one filter name, \fB-filt\fP will use that filter in both x and y;
Xif given two names, different filters can be used in the two dimensions.
XThe command \fBzoom -filt '?'\fP prints
Xthe list of filters currently known to \fIzoom\fP.
XThat list is:
X.DS
XNAME SUPPORT
Xpoint 0
Xbox 0.5
Xtriangle 1
Xquadratic 1.5
Xcubic 2
Xcatrom 2
Xmitchell 2
Xgaussian 1.25
Xsinc 4
Xbessel 3.24
X.DE
XThe option \fB-filt point\fP gives fast pixel replication,
X\fB-filt triangle\fP (the default) gives bilinear interpolation,
Xwhich is suitable for most purposes,
Xand \fB-filt mitchell\fP gives slower, very high quality results.
XThe other filters are provided for experimentation.
XThe time required by a filter is proportional to its support.
X.PP
XMost of the filters known to \fIzoom\fP are FIR (finite impulse response),
Xwith an intrinsic support (width),
Xbut some of them (gaussian, sinc, bessel) are IIR (infinite impulse response),
Xand must be truncated at some arbitrary support.
XThis can be done with the \fB-supp\fP option.
XAgain, the defaults are reasonable.
XThe IIR filters can be windowed (brought down to zero) with various window
Xfunctions listed below:
X.Cs
Xhanning
Xhamming
Xblackman
Xkaiser
X.Ce
XThe sinc and bessel filters are blackman-windowed by default.
XFilters can be scaled artificially to blur or sharpen them with the \fB-blur\fP
Xoption.
X.SH EXAMPLES
X.Ss
Xzoom -src mandrill.dump
X.Sd
XZoom the mandrill from picture file \fBmandrill.dump\fP
Xinto the default destination device (whatever that is),
Xmapping the file's rectangle into the device's rectangle,
Xwith a triangle filter.
X.Se
X
X.Ss
Xzoom -src mandrill.dump -filt point -square -intscale
X.Sd
XZoom the mandrill to the default device's full screen,
Xbut maintain the picture's aspect ratio, and zoom it up by an integer factor
Xwith point sampling (pixel replication)
X.Se
X
X.Ss
Xzoom -src mandrill.dump -dst iris -d 50 75 100 100
X.Sd
XZoom the mandrill into an iris window at position (50,70) with xsize=ysize=100.
X.Se
X
X.Ss
Xzoom -src hp -dst hp -s 100 100 640 512 -d 0 0 1280 1024 -filt mitchell
X.Sd
XZoom from hp to itself with overlapping source and destination windows,
Xusing a mitchell filter.
X.Se
X
X.Ss
Xzoom -src hp -dst hp -s 0 0 500 500 -d 0 0 500 500 -filt cubic -blur 2
X.Sd
XBlur (low pass filter) an image in-place.
X.Se
X.SH SEE ALSO
Xpic(3),
X\fIDiscrete Time Signal Processing\fP, Alan Oppenheim, Ronald Schafer, 1989.
X.SH AUTHOR
XPaul Heckbert, UC Berkeley, August 1989.
Xph@miro.berkeley.edu
EOF13087
echo extracting zoom/zoom.c
sed 's/^X//' <<'EOF13088' >zoom/zoom.c
X/*
X * zoom: subroutines for in-place one-pass filtered image zooming
X *
X * Zooms (resizes) one image into another with independent control of
X * scale, translation, and filtering in x and y. Magnifies or minifies.
X * Filters with an arbitrary separable finite impulse response filter.
X * (A filter h(x,y) is called separable if h(x,y)=f(x)*g(y)).
X * See the filt package regarding filtering.
X * The code is graphics device independent, using only generic scanline
X * i/o operations.
X *
X * The program makes one pass over the image using a moving window of
X * scanlines, so memory usage is proportional to imagewidth*filterheight,
X * not imagewidth*imageheight, as you'd get if the entire image were
X * buffered. The separability of the filter is also exploited, to get
X * a cpu time proportional to npixels*(filterwidth+filterheight),
X * not npixels*filterwidth*filterheight as with direct 2-D convolution.
X *
X * The source and destination images can be identical, with overlapping
X * windows, in which case the algorithm finds the magic split point and
X * uses the appropriate scanning directions to avoid feedback (recursion)
X * in the filtering.
X *
X * Paul Heckbert 12 Sept 1988
X * UC Berkeley
X * ph@miro.berkeley.edu
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/*
X * routines in this source file:
X *
X * PUBLIC:
X * zoom most common entry point; window-to-window zoom
X * zoom_continuous fancier: explicit, continuous scale&tran control
X *
X * SORTA PRIVATE:
X * zoom_unfiltered_mono point-sampled zooming, called by zoom_continuous
X * zoom_unfiltered_rgb point-sampled zooming, called by zoom_continuous
X * zoom_filtered filtered zooming, called by zoom_continuous
X * zoom_filtered_xy zoom, filtering x before y, called by zoom_filtered
X * zoom_filtered_yx zoom, filtering y before x, called by zoom_filtered
X *
X * make_weighttab make lookup table of filter coefficients
X * make_map_table in-place tricks; order scanlines to avoid feedback
X */
X
Xstatic char rcsid[] = "$Header: zoom.c,v 3.0 88/10/10 13:52:07 ph Locked $";
X
X#include <math.h>
X
X#include <simple.h>
X#include <pic.h>
X#include "filt.h"
X#include "scanline.h"
X#include "zoom.h"
X
X#define EPSILON 1e-7 /* error tolerance */
X#define UNDEF PIC_UNDEFINED
X
X#define WEIGHTBITS 14 /* # bits in filter coefficients */
X#define FINALSHIFT (2*WEIGHTBITS-CHANBITS) /* shift after x&y filter passes */
X#define WEIGHTONE (1<<WEIGHTBITS) /* filter weight of one */
X
X#define INTEGER(x) (fabs((x)-floor((x)+.5)) < EPSILON) /* is x an integer? */
X#define FRAC(x) fabs((x)-floor((x)+.5)) /* diff from closest integer */
X
Xtypedef struct { /* ZOOM-SPECIFIC FILTER PARAMETERS */
X double scale; /* filter scale (spacing between centers in a space) */
X double supp; /* scaled filter support radius */
X int wid; /* filter width: max number of nonzero samples */
X} Filtpar;
X
Xint zoom_debug = 0; /* debug level: 0=none, 1=scanline i/o, 2=filters */
Xint zoom_coerce = 1; /* simplify filters if possible? */
Xint zoom_xy = UNDEF; /* filter x before y (1) or vice versa (0)? */
Xint zoom_trimzeros = 1; /* trim zeros in x filter weight tables? */
X
Xstatic int nzero = 0, ntot = 0, nzmax = 0;
X
XFilt *simplify_filter();
X
X/*----------------------------------------------------------------------
X * The mapping from source coordinates to dest coordinates:
X * (Notation: prefix 'a' denotes source coords, 'b' denotes destination coords)
X *
X * bx = sx*ax+tx
X * by = sy*ay+ty
X *
X * where (ax,ay) and (bx,by) are DISCRETE source and dest coords, respectively.
X *
X * An important fine point on terminology: there are two kinds of pixel coords:
X * DISCRETE COORDINATES take on integer values at pixel centers
X * CONTINUOUS COORDINATES take on integer values halfway between pixels
X * For example, an image with discrete coordinate range [0..511] has a
X * continuous coordinate range of [0..512]. The pixel with
X * discrete coords (x,y) has its center at continuous coords (x+.5,y+.5)
X * and its continuous coord domain is [x..x+1] in X and [y..y+1] in Y.
X * Note: discrete coords are not always stored in ints and continuous coords
X * are not always stored in floats.
X *
X * conversion:
X * if c = continuous coord and d = discrete coord, then
X * c = d+.5
X * d = floor(c)
X *
X * To map a discrete src interval [a0..a1] to a discrete dst interval [b0..b1]:
X *
X * b-b0+.5 = bn/an*(a-a0+.5)
X * or b = (bn/an)*a + (b0-.5)-(bn/an)*(a0-.5)
X * = scale*a+tran
X *
X * where a and b are the discrete source and dest coords (either x or y)
X * and an and bn are the interval lengths: an=a1-a0+1, bn=b1-b0+1.
X * a0, an, b0, bn are the mapping parameters used by the zoom routine below.
X * In general, however, for sub-pixel scale and translate control, we allow
X * any real numbers for scale and tran (although there should usually be some
X * correspondence between the mapping function and the src and dst windows).
X *
X * We'll often want the inverse mapping, from discrete dst coord b to
X * continuous src coord ac, relative to the interval origins a0 and b0:
X *
X * b+b0 = s*(ac+a0-.5)+t
X * so ac = (b+b0-s*(a0-.5)-t)/s
X * = (b+offset)/scale = MAP(b, scale, offset)
X *
X * The mapping fields ux and uy in the Mapping structure
X * are these offsets in x & y respectively.
X */
X
X/* the mapping from discrete dest coord b to continuous source coord: */
X#define MAP(b, scale, offset) ((b)+(offset))/(scale)
X
X/*----------------------------------------------------------------------
X * zoom: simplest entry point; window-to-window zoom.
X * Zoom window a of apic to window b of bpic using the specified filters.
X * Window a represents the pixels with x in [a->x0..a->x1]
X * and y in [a->y0..a->y1], (ranges are inclusive). Likewise for b.
X * apic and bpic may be equal and a and b may overlap, no sweat.
X */
X
Xzoom(apic, a, bpic, b, xfilt, yfilt)
XPic *apic, *bpic; /* source and dest pictures */
XWindow_box *a, *b; /* source and dest windows */
XFilt *xfilt, *yfilt; /* filters for x and y */
X{
X Mapping map;
X
X if (b->x0==UNDEF) {
X printf("zoom: undefined window for dest %s: ",
X pic_get_name(bpic));
X window_print("", &b);
X printf("\n");
X exit(1);
X }
X window_box_set_size(a);
X window_box_set_size(b);
X if (a->nx<=0 || a->ny<=0) return;
X map.sx = (double)b->nx/a->nx;
X map.sy = (double)b->ny/a->ny;
X map.tx = b->x0-.5 - map.sx*(a->x0-.5);
X map.ty = b->y0-.5 - map.sy*(a->y0-.5);
X zoom_continuous(apic, a, bpic, b, &map, xfilt, yfilt);
X}
X
X/*
X * zoom_opt: window-to-window zoom with options.
X * Like zoom() but has options to make mapping square (preserve pixel shape)
X * and round to the next lowest integer scale factor.
X */
X
Xzoom_opt(apic, a, bpic, b, xfilt, yfilt, square, intscale)
XPic *apic, *bpic; /* source and dest pictures */
XWindow_box *a, *b; /* source and dest windows */
XFilt *xfilt, *yfilt; /* filters for x and y */
Xint square, intscale; /* square mapping? integer scales? */
X{
X Mapping map;
X
X if (b->x0==UNDEF) {
X printf("zoom_opt: undefined window for dest %s: ",
X pic_get_name(bpic));
X window_print("", &b);
X printf("\n");
X exit(1);
X }
X window_box_set_size(a);
X window_box_set_size(b);
X if (a->nx<=0 || a->ny<=0) return;
X map.sx = (double)b->nx/a->nx;
X map.sy = (double)b->ny/a->ny;
X if (square) {
X if (map.sx>map.sy) { /* use the smaller scale for both sx and sy */
X map.sx = map.sy;
X b->nx = ceil(a->nx*map.sx);
X }
X else {
X map.sy = map.sx;
X b->ny = ceil(a->ny*map.sy);
X }
X }
X if (intscale) {
X if (map.sx>1.) {
X map.sx = floor(map.sx);
X b->nx = ceil(a->nx*map.sx);
X }
X if (map.sy>1.) {
X map.sy = floor(map.sy);
X b->ny = ceil(a->ny*map.sy);
X }
X }
X window_box_set_max(b);
X map.tx = b->x0-.5 - map.sx*(a->x0-.5);
X map.ty = b->y0-.5 - map.sy*(a->y0-.5);
X zoom_continuous(apic, a, bpic, b, &map, xfilt, yfilt);
X}
X
X/*
X * zoom_continuous: zoom with explicit, continuous scale&tran control.
X * Zoom according to the mapping defined in m,
X * reading from window awin of apic and writing to window bwin of bpic.
X * Like zoom but allows continuous, subpixel scales and translates.
X * Meaning of m is explained above. Pixels in the dest window bwin outside
X * of the support of the mapped source window will be left unchanged.
X * Scales in m should be positive.
X */
X
Xzoom_continuous(apic, awin, bpic, bwin, m, xfilt, yfilt)
XPic *apic, *bpic;
XWindow_box *awin, *bwin;
XMapping *m;
XFilt *xfilt, *yfilt;
X{
X int xy;
X Window_box a, b, bc, t;
X Filtpar ax, ay;
X
X /* we consider 3 channel and 4 channel to be the same */
X if ((pic_get_nchan(apic)==1) != (pic_get_nchan(bpic)==1)) {
X fprintf(stderr, "zoom_continuous: source has %d channels, dest has %d channels\n",
X pic_get_nchan(apic), pic_get_nchan(bpic));
X return;
X }
X if (m->sx<=0. || m->sy<=0.) {
X fprintf(stderr, "zoom_continuous: negative scales not allowed\n");
X return;
X }
X
X a = *awin;
X b = *bwin;
X window_box_set_size(&a);
X window_box_set_size(&b);
X
X /*
X * find scale of filter in a space (source space)
X * when minifying, ascale=1/scale, but when magnifying, ascale=1
X */
X ax.scale = xfilt->blur*MAX(1., 1./m->sx);
X ay.scale = yfilt->blur*MAX(1., 1./m->sy);
X /*
X * find support radius of scaled filter
X * if ax.supp and ay.supp are both <=.5 then we've got point sampling.
X * Point sampling is essentially a special filter whose width is fixed
X * at one source pixel.
X */
X ax.supp = MAX(.5, ax.scale*xfilt->supp);
X ay.supp = MAX(.5, ay.scale*yfilt->supp);
X ax.wid = ceil(2.*ax.supp);
X ay.wid = ceil(2.*ay.supp);
X
X /* clip source and dest windows against their respective pictures */
X window_box_intersect(&a, (Window_box *)pic_get_window(apic, &t), &a);
X if (pic_get_window(bpic, &t)->x0 != UNDEF)
X window_box_intersect(&b, &t, &b);
X
X /* find valid dest window (transformed source + support margin) */
X bc.x0 = ceil(m->sx*(a.x0-ax.supp)+m->tx+EPSILON);
X bc.y0 = ceil(m->sy*(a.y0-ay.supp)+m->ty+EPSILON);
X bc.x1 = floor(m->sx*(a.x1+ax.supp)+m->tx-EPSILON);
X bc.y1 = floor(m->sy*(a.y1+ay.supp)+m->ty-EPSILON);
X window_box_set_size(&bc);
X
X if (b.x0<bc.x0 || b.x1>bc.x1 || b.y0<bc.y0 || b.y1>bc.y1) {
X /* requested dest window lies outside the valid dest, so clip dest */
X window_box_print(" clipping odst=", &b);
X window_box_print(" to ", &bc);
X printf("\n");
X window_box_intersect(&b, &bc, &b);
X }
X pic_set_window(bpic, &b);
X
X /* compute offsets for MAP (these will be .5 if zoom() routine was called)*/
X m->ux = b.x0-m->sx*(a.x0-.5)-m->tx;
X m->uy = b.y0-m->sy*(a.y0-.5)-m->ty;
X
X printf("src=%s:%s", pic_get_dev(apic), pic_get_name(apic));
X window_box_print("", &a);
X printf("\ndst=%s:%s", pic_get_dev(bpic), pic_get_name(bpic));
X window_box_print("", &b);
X printf("\n");
X
X if (a.nx<=0 || a.ny<=0 || b.nx<=0 || b.ny<=0) return;
X
X /* check for high-level simplifications of filter */
X xfilt = simplify_filter("x", zoom_coerce, xfilt, &ax, m->sx, m->ux);
X yfilt = simplify_filter("y", zoom_coerce, yfilt, &ay, m->sy, m->uy);
X
X /*
X * decide which filtering order (xy or yx) is faster for this mapping by
X * counting convolution multiplies
X */
X xy = zoom_xy!=UNDEF ? zoom_xy :
X b.nx*(a.ny*ax.wid+b.ny*ay.wid) < b.ny*(b.nx*ax.wid+a.nx*ay.wid);
X
X /* choose the appropriate filtering routine */
X if (str_eq(xfilt->name, "point") && str_eq(yfilt->name, "point"))
X zoom_unfiltered(apic, &a, bpic, &b, m);
X else if (xy)
X zoom_filtered_xy(apic, &a, bpic, &b, m, xfilt, &ax, yfilt, &ay);
X else
X zoom_filtered_yx(apic, &a, bpic, &b, m, xfilt, &ax, yfilt, &ay);
X}
X
X/*
X * filter_simplify: check if our discrete sampling of an arbitrary continuous
X * filter, parameterized by the filter spacing (a->scale), its radius (a->supp),
X * and the scale and offset of the coordinate mapping (s and u), causes the
X * filter to reduce to point sampling.
X *
X * It reduces if support is less than 1 pixel or
X * if integer scale and translation, and filter is cardinal
X */
X
Xstatic Filt *simplify_filter(dim, coerce, filter, a, s, u)
Xchar *dim;
Xint coerce;
XFilt *filter;
XFiltpar *a;
Xdouble s, u; /* scale and offset */
X{
X if (coerce && (a->supp<=.5 ||
X filter->cardinal && INTEGER(1./a->scale) && INTEGER(1./(s*a->scale))
X && INTEGER((u/s-.5)/a->scale))) {
X if (!str_eq(filter->name, "point"))
X printf("coercing %sfilter=%s to point\n", dim, filter->name);
X filter = filt_find("point");
X a->scale = 1.;
X a->supp = .5;
X a->wid = 1;
X }
X return filter;
X}
X
X/*----------------------------------------------------------------------*/
X
X/* zoom_unfiltered: do unfiltered zoom (point sampling) */
X
Xzoom_unfiltered(apic, a, bpic, b, m)
XPic *apic, *bpic;
XWindow_box *a, *b;
XMapping *m;
X{
X int overlap;
X
X /* do source and dest windows overlap? */
X overlap = apic==bpic && window_box_overlap(a, b);
X
X printf("-map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
X /*
X * note: We have only x-resample before y-resample versions of the
X * unfiltered zoom case; we could optimize further by creating
X * y-resample before x-resample versions.
X */
X if (pic_get_nchan(apic)==1)
X zoom_unfiltered_mono(apic, a, bpic, b, m, overlap);
X else
X zoom_unfiltered_rgb (apic, a, bpic, b, m, overlap);
X}
X
X/* zoom_unfiltered_mono: monochrome unfiltered zoom */
X
Xzoom_unfiltered_mono(apic, a, bpic, b, m, overlap)
XPic *apic, *bpic;
XWindow_box *a, *b;
XMapping *m;
Xint overlap;
X{
X int byi, by, ay, ayold, bx, ax;
X Pixel1 *abuf, *bbuf, *bp; /* source and dest scanline buffers */
X Pixel1 **xmap, **xp; /* mapping from abuf to bbuf (optimization) */
X short *ymap; /* order of dst y coords that avoids feedback */
X char *linewritten; /* has scanline y been written? (debugging) */
X
X ALLOC(abuf, Pixel1, a->nx);
X ALLOC(bbuf, Pixel1, b->nx);
X ALLOC(xmap, Pixel1 *, b->nx);
X ALLOC(ymap, short, b->ny);
X ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
X
X /* if overlapping src & dst, find magic y ordering that avoids feedback */
X make_map_table(m->sy, m->ty, .5, a->y0, b->y0, b->ny, overlap, ymap);
X
X for (bx=0; bx<b->nx; bx++) {
X ax = MAP(bx, m->sx, m->ux);
X xmap[bx] = &abuf[ax];
X }
X
X ayold = -1; /* impossible value for ay */
X for (byi=0; byi<b->ny; byi++) {
X by = ymap[byi];
X ay = MAP(by, m->sy, m->uy);
X /* scan line a.y0+ay goes to b.y0+by */
X if (ay!=ayold) { /* new scan line; read it in */
X ayold = ay;
X if (overlap && linewritten[a->y0+ay])
X fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ay);
X /* read scan line into abuf */
X pic_read_row(apic, a->y0+ay, a->x0, a->nx, abuf);
X /* resample in x */
X for (bp=bbuf, xp=xmap, bx=0; bx<b->nx; bx++)
X *bp++ = **xp++;
X }
X pic_write_row(bpic, b->y0+by, b->x0, b->nx, bbuf);
X linewritten[b->y0+by] = 1;
X }
X free(abuf);
X free(bbuf);
X free(xmap);
X free(ymap);
X}
X
X/* zoom_unfiltered_rgb: color unfiltered zoom */
X
Xzoom_unfiltered_rgb(apic, a, bpic, b, m, overlap)
XPic *apic, *bpic;
XWindow_box *a, *b;
XMapping *m;
Xint overlap;
X{
X int byi, by, ay, ayold, bx, ax;
X Pixel1_rgba *abuf, *bbuf, *bp; /* source and dest scanline buffers */
X Pixel1_rgba **xmap, **xp; /* mapping from abuf to bbuf */
X short *ymap; /* order of dst y coords that avoids feedback */
X char *linewritten; /* has scanline y been written? (debugging) */
X
X ALLOC(abuf, Pixel1_rgba, a->nx);
X ALLOC(bbuf, Pixel1_rgba, b->nx);
X ALLOC(xmap, Pixel1_rgba *, b->nx);
X ALLOC(ymap, short, b->ny);
X ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
X
X /* if overlapping src & dst, find magic y ordering that avoids feedback */
X make_map_table(m->sy, m->ty, .5, a->y0, b->y0, b->ny, overlap, ymap);
X
X for (bx=0; bx<b->nx; bx++) {
X ax = MAP(bx, m->sx, m->ux);
X xmap[bx] = &abuf[ax];
X }
X
X ayold = -1; /* impossible value for ay */
X for (byi=0; byi<b->ny; byi++) {
X by = ymap[byi];
X ay = MAP(by, m->sy, m->uy);
X /* scan line a.y0+ay goes to b.y0+by */
X if (ay!=ayold) { /* new scan line; read it in */
X ayold = ay;
X if (overlap && linewritten[a->y0+ay])
X fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ay);
X /* read scan line into abuf */
X pic_read_row_rgba(apic, a->y0+ay, a->x0, a->nx, abuf);
X /* resample in x */
X for (bp=bbuf, xp=xmap, bx=0; bx<b->nx; bx++)
X *bp++ = **xp++;
X }
X pic_write_row_rgba(bpic, b->y0+by, b->x0, b->nx, bbuf);
X linewritten[b->y0+by] = 1;
X }
X free(abuf);
X free(bbuf);
X free(xmap);
X free(ymap);
X}
X
X/*----------------------------------------------------------------------*/
X
X/*
X * zoom_filtered_xy: filtered zoom, xfilt before yfilt
X *
X * note: when calling make_weighttab, we can trim leading and trailing
X * zeros from the x weight buffers as an optimization,
X * but not for y weight buffers since the split formula is anticipating
X * a constant amount of buffering of source scanlines;
X * trimming zeros in yweight could cause feedback.
X */
X
X
Xzoom_filtered_xy(apic, a, bpic, b, m, xfilt, ax, yfilt, ay)
XPic *apic, *bpic;
XWindow_box *a, *b;
XMapping *m;
XFilt *xfilt, *yfilt; /* x and y filters */
XFiltpar *ax, *ay; /* extra x and y filter parameters */
X{
X int ayf, bx, byi, by, overlap, nchan;
X
X /*PIXELTYPE NBITS RAISON D'ETRE */
X Scanline abbuf; /* PIXEL1 8 src&dst scanline bufs */
X Scanline accum; /* PIXEL4 28 accumulator buf for yfilt */
X Scanline *linebuf, *tbuf; /* PIXEL2 14 circular buf of active lines */
X
X short *ymap; /* order of dst y coords that avoids feedback */
X Weighttab *xweights; /* sampled filter at each dst x pos; for xfilt*/
X short *xweightbuf, *xwp; /* big block of memory addressed by xweights */
X Weighttab yweight; /* a single sampled filter for current y pos */
X char *linewritten; /* has scanline y been written? (debugging) */
X
X nchan = pic_get_nchan(apic)==1 ? PIXEL_MONO : PIXEL_RGB;
X scanline_alloc(&abbuf, PIXEL1|nchan, MAX(a->nx, b->nx));
X scanline_alloc(&accum, PIXEL4|nchan, b->nx);
X ALLOC(linebuf, Scanline, ay->wid);
X /* construct circular buffer of ay->wid intermediate scanlines */
X for (ayf=0; ayf<ay->wid; ayf++) {
X scanline_alloc(&linebuf[ayf], PIXEL2|nchan, b->nx);
X linebuf[ayf].y = -1; /* mark scanline as unread */
X }
X
X ALLOC(ymap, short, b->ny);
X ALLOC(xweights, Weighttab, b->nx);
X ALLOC(xweightbuf, short, b->nx*ax->wid);
X ALLOC(yweight.weight, short, ay->wid);
X ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
X
X /* do source and dest windows overlap? */
X overlap = apic==bpic && window_box_overlap(a, b);
X
X printf("-xy -map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
X printf("X: filt=%s support=%g, scale=%g scaledsupport=%g width=%d\n",
X xfilt->name, xfilt->supp, ax->scale, ax->supp, ax->wid);
X printf("Y: filt=%s support=%g, scale=%g scaledsupport=%g width=%d\n",
X yfilt->name, yfilt->supp, ay->scale, ay->supp, ay->wid);
X
X /*
X * prepare a weighttab (a sampled filter for source pixels) for
X * each dest x position
X */
X for (xwp=xweightbuf, bx=0; bx<b->nx; bx++, xwp+=ax->wid) {
X xweights[bx].weight = xwp;
X make_weighttab(b->x0+bx, MAP(bx, m->sx, m->ux),
X xfilt, ax, a->nx, zoom_trimzeros, &xweights[bx]);
X }
X
X /* if overlapping src & dst, find magic y ordering that avoids feedback */
X make_map_table(m->sy, m->ty, ay->supp, a->y0, b->y0, b->ny, overlap,
X ymap);
X
X for (byi=0; byi<b->ny; byi++) { /* loop over dest scanlines */
X by = ymap[byi];
X if (zoom_debug) printf("by=%d: ", b->y0+by);
X /* prepare a weighttab for dest y position by */
X make_weighttab(b->y0+by, MAP(by, m->sy, m->uy),
X yfilt, ay, a->ny, 0, &yweight);
X if (zoom_debug) printf("ay=%d-%d, reading: ",
X a->y0+yweight.i0, a->y0+yweight.i1-1);
X scanline_zero(&accum);
X
X /* loop over source scanlines that influence this dest scanline */
X for (ayf=yweight.i0; ayf<yweight.i1; ayf++) {
X tbuf = &linebuf[ayf % ay->wid];
X if (tbuf->y != ayf) { /* scanline needs to be read in */
X if (zoom_debug) printf("%d ", a->y0+ayf);
X if (overlap && linewritten[a->y0+ayf])
X fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ayf);
X tbuf->y = ayf;
X /* read source scanline into abbuf */
X abbuf.len = a->nx;
X scanline_read(apic, a->x0, a->y0+ayf, &abbuf);
X /* and filter it into the appropriate line of linebuf (xfilt) */
X scanline_filter(CHANBITS, xweights, &abbuf, tbuf);
X }
X /* add weighted tbuf into accum (these do yfilt) */
X scanline_accum(yweight.weight[ayf-yweight.i0], tbuf, &accum);
X }
X
X /* shift accum down into abbuf and write out dest scanline */
X abbuf.len = b->nx;
X scanline_shift(FINALSHIFT, &accum, &abbuf);
X scanline_write(bpic, b->x0, b->y0+by, &abbuf);
X linewritten[b->y0+by] = 1;
X if (zoom_debug) printf("\n");
X }
X
X scanline_free(&abbuf);
X scanline_free(&accum);
X for (ayf=0; ayf<ay->wid; ayf++)
X scanline_free(&linebuf[ayf]);
X free(ymap);
X free(linebuf);
X free(xweights);
X free(xweightbuf);
X free(yweight.weight);
X free(linewritten);
X statistics();
X}
X
X/* zoom_filtered_yx: filtered zoom, yfilt before xfilt */
X
Xzoom_filtered_yx(apic, a, bpic, b, m, xfilt, ax, yfilt, ay)
XPic *apic, *bpic;
XWindow_box *a, *b;
XMapping *m;
XFilt *xfilt, *yfilt; /* x and y filters */
XFiltpar *ax, *ay; /* extra x and y filter parameters */
X{
X int ayf, bx, byi, by, overlap, nchan;
X
X /*PIXELTYPE NBITS RAISON D'ETRE */
X Scanline bbuf; /* PIXEL1 8 dst scanline buf */
X Scanline accum; /* PIXEL4 22 accumulator buf for yfilt */
X Scanline *linebuf, *tbuf; /* PIXEL1 8 circular buf of active lines */
X
X short *ymap; /* order of dst y coords that avoids feedback */
X Weighttab *xweights; /* sampled filter at each dst x pos; for xfilt*/
X short *xweightbuf, *xwp; /* big block of memory addressed by xweights */
X Weighttab yweight; /* a single sampled filter for current y pos */
X char *linewritten; /* has scanline y been written? (debugging) */
X
X nchan = pic_get_nchan(apic)==1 ? PIXEL_MONO : PIXEL_RGB;
X scanline_alloc(&bbuf, PIXEL1|nchan, b->nx);
X scanline_alloc(&accum, PIXEL4|nchan, a->nx);
X ALLOC(linebuf, Scanline, ay->wid);
X /* construct circular buffer of ay->wid intermediate scanlines */
X for (ayf=0; ayf<ay->wid; ayf++) {
X scanline_alloc(&linebuf[ayf], PIXEL1|nchan, a->nx);
X linebuf[ayf].y = -1; /* mark scanline as unread */
X }
X
X ALLOC(ymap, short, b->ny);
X ALLOC(xweights, Weighttab, b->nx);
X ALLOC(xweightbuf, short, b->nx*ax->wid);
X ALLOC(yweight.weight, short, ay->wid);
X ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
X
X /* do source and dest windows overlap? */
X overlap = apic==bpic && window_box_overlap(a, b);
X
X printf("-yx -map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
X printf("X: filt=%s supp=%g scale=%g scaledsupp=%g wid=%d\n",
X xfilt->name, xfilt->supp, ax->scale, ax->supp, ax->wid);
X printf("Y: filt=%s supp=%g scale=%g scaledsupp=%g wid=%d\n",
X yfilt->name, yfilt->supp, ay->scale, ay->supp, ay->wid);
X
X /*
X * prepare a weighttab (a sampled filter for source pixels) for
X * each dest x position
X */
X for (xwp=xweightbuf, bx=0; bx<b->nx; bx++, xwp+=ax->wid) {
X xweights[bx].weight = xwp;
X make_weighttab(b->x0+bx, MAP(bx, m->sx, m->ux),
X xfilt, ax, a->nx, zoom_trimzeros, &xweights[bx]);
X }
X
X /* if overlapping src & dst, find magic y ordering that avoids feedback */
X make_map_table(m->sy, m->ty, ay->supp, a->y0, b->y0, b->ny, overlap,
X ymap);
X
X for (byi=0; byi<b->ny; byi++) { /* loop over dest scanlines */
X by = ymap[byi];
X if (zoom_debug) printf("by=%d: ", b->y0+by);
X /* prepare a weighttab for dest y position by */
X make_weighttab(b->y0+by, MAP(by, m->sy, m->uy),
X yfilt, ay, a->ny, 0, &yweight);
X if (zoom_debug) printf("ay=%d-%d, reading: ",
X a->y0+yweight.i0, a->y0+yweight.i1-1);
X scanline_zero(&accum);
X
X /* loop over source scanlines that influence this dest scanline */
X for (ayf=yweight.i0; ayf<yweight.i1; ayf++) {
X tbuf = &linebuf[ayf % ay->wid];
X if (tbuf->y != ayf) { /* scanline needs to be read in */
X if (zoom_debug) printf("%d ", a->y0+ayf);
X if (overlap && linewritten[a->y0+ayf])
X fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ayf);
X tbuf->y = ayf;
X /* read source scanline into linebuf */
X scanline_read(apic, a->x0, a->y0+ayf, tbuf);
X }
X /* add weighted tbuf into accum (these do yfilt) */
X scanline_accum(yweight.weight[ayf-yweight.i0], tbuf, &accum);
X }
X
X /* and filter it into the appropriate line of linebuf (xfilt) */
X scanline_filter(FINALSHIFT, xweights, &accum, &bbuf);
X /* and write out dest scanline in bbuf */
X scanline_write(bpic, b->x0, b->y0+by, &bbuf);
X linewritten[b->y0+by] = 1;
X if (zoom_debug) printf("\n");
X }
X
X scanline_free(&bbuf);
X scanline_free(&accum);
X for (ayf=0; ayf<ay->wid; ayf++)
X scanline_free(&linebuf[ayf]);
X free(ymap);
X free(linebuf);
X free(xweights);
X free(xweightbuf);
X free(yweight.weight);
X free(linewritten);
X statistics();
X}
X
X/*
X * make_weighttab: sample the continuous filter, scaled by ap->scale and
X * positioned at continuous source coordinate cen, for source coordinates in
X * the range [0..len-1], writing the weights into wtab.
X * Scale the weights so they sum to WEIGHTONE, and trim leading and trailing
X * zeros if trimzeros!=0.
X * b is the dest coordinate (for diagnostics).
X */
X
Xstatic make_weighttab(b, cen, filter, ap, len, trimzeros, wtab)
Xint b, len;
Xdouble cen;
XFilt *filter;
XFiltpar *ap;
Xint trimzeros;
XWeighttab *wtab;
X{
X int i0, i1, i, sum, t, stillzero, lastnonzero, nz;
X short *wp;
X double den, sc, tr;
X
X /* find the source coord range of this positioned filter: [i0..i1-1] */
X i0 = cen-ap->supp+.5;
X i1 = cen+ap->supp+.5;
X if (i0<0) i0 = 0;
X if (i1>len) i1 = len;
X if (i0>=i1) {
X fprintf(stderr, "make_weighttab: null filter at %d\n", b);
X exit(1);
X }
X /* the range of source samples to buffer: */
X wtab->i0 = i0;
X wtab->i1 = i1;
X
X /* find scale factor sc to normalize the filter */
X for (den=0, i=i0; i<i1; i++)
X den += filt_func(filter, (i+.5-cen)/ap->scale);
X /* set sc so that sum of sc*func() is approximately WEIGHTONE */
X sc = den==0. ? WEIGHTONE : WEIGHTONE/den;
X if (zoom_debug>1) printf(" b=%d cen=%g scale=%g [%d..%d) sc=%g: ",
X b, cen, ap->scale, i0, i1, sc);
X
X /* compute the discrete, sampled filter coefficients */
X stillzero = trimzeros;
X for (sum=0, wp=wtab->weight, i=i0; i<i1; i++) {
X
X /* evaluate the filter function: */
X tr = sc*filt_func(filter, (i+.5-cen)/ap->scale);
X
X if (tr<MINSHORT || tr>MAXSHORT) {
X fprintf(stderr, "tr=%g at %d\n", tr, b);
X exit(1);
X }
X t = floor(tr+.5);
X if (stillzero && t==0) i0++; /* find first nonzero */
X else {
X stillzero = 0;
X *wp++ = t; /* add weight to table */
X sum += t;
X if (t!=0) lastnonzero = i; /* find last nonzero */
X }
X }
X ntot += wtab->i1-wtab->i0;
X if (sum==0) {
X nz = wtab->i1-wtab->i0;
X /* fprintf(stderr, "sum=0 at %d\n", b); */
X wtab->i0 = wtab->i0+wtab->i1 >> 1;
X wtab->i1 = wtab->i0+1;
X wtab->weight[0] = WEIGHTONE;
X }
X else {
X if (trimzeros) { /* skip leading and trailing zeros */
X /* set wtab->i0 and ->i1 to the nonzero support of the filter */
X nz = wtab->i1-wtab->i0-(lastnonzero-i0+1);
X wtab->i0 = i0;
X wtab->i1 = i1 = lastnonzero+1;
X }
X else /* keep leading and trailing zeros */
X nz = 0;
X if (sum!=WEIGHTONE) {
X /*
X * Fudge the center slightly to make sum=WEIGHTONE exactly.
X * Is this the best way to normalize a discretely sampled
X * continuous filter?
X */
X i = cen+.5;
X if (i<i0) i = i0; else if (i>=i1) i = i1-1;
X t = WEIGHTONE-sum;
X if (zoom_debug>1) printf("[%d]+=%d ", i, t);
X wtab->weight[i-i0] += t; /* fudge center sample */
X }
X }
X if (nz>nzmax) nzmax = nz;
X nzero += nz;
X if (zoom_debug>1) {
X printf("\t");
X for (wp=wtab->weight, i=i0; i<i1; i++, wp++)
X printf("%5d ", *wp);
X printf("\n");
X }
X}
X
X/*
X * make_map_table:
X * construct a table for the mapping a = (b-t)/s for b-b0 in [0..bn-1]
X * ordered so that the buffer resampling operation
X * for (bi=0; bi<bn; bi++) {
X * b = map[bi];
X * a = MAP(b, scale, offset);
X * buf[b] = buf[a];
X * }
X * can work IN PLACE without feedback
X *
X * a and b here are source and dest coordinates, in either x or y.
X * This is needed only if there is overlap between source and dest windows.
X */
X
Xstatic make_map_table(scale, tran, asupp, a0, b0, bn, overlap, map)
Xdouble scale, tran; /* scale and translate */
Xdouble asupp; /* filter support radius in source space */
Xint a0, b0, bn; /* source and dest positions, dest length */
Xint overlap; /* do source and dest overlap? */
Xshort *map; /* mapping table */
X{
X int split, b, i0;
X double z, u;
X
X /* find fixed point of mapping; let split=b-b0 at the point where a=b */
X
X if (overlap) { /* source and dest windows overlap */
X if (scale==1.) /* no scale change, translation only */
X /* if moving left then scan to right, and vice versa */
X split = a0<b0 ? 0 : bn;
X else { /* magnify or minify */
X
X /* THE MAGIC SPLIT FORMULA: */
X
X if (scale>1.) { /* magnify */
X /*
X * choose integer nearest midpoint of valid interval:
X * x < split <= x+s/(s-1)
X * where x=(tran+scale*asupp)/(1-scale)-b0
X */
X split = floor((tran+scale*asupp-.5)/(1.-scale)-b0+1.);
X }
X else { /* minify */
X /*
X * only one valid split pt in this case:
X * x <= split < x+1
X * so we take extra care (x as above)
X */
X split = ceil((tran+scale*asupp)/(1.-scale)-b0);
X /*
X * The above formula is perfect for exact arithmetic,
X * but hardware roundoff can cause split to be one too big.
X * Check for roundoff by simulating precisely the calculation
X * of i0 in make_weighttab.
X */
X u = b0-scale*(a0-.5)-tran; /* recalculate offset */
X z = MAP(split-1, scale, u); /* cen at b=split-1 */
X z = z-asupp+.5;
X i0 = z; /* i0 at b=split-1 */
X if (a0+i0>=b0+split) /* feedback at b=split-1? */
X split--; /* correct for roundoff */
X }
X if (split<0) split = 0;
X else if (split>bn) split = bn;
X printf("split at y=%d\n", b0+split);
X }
X
X if (scale>=1.) { /* magnify: scan in toward split */
X for (b=0; b<split; b++) *map++ = b;
X for (b=bn-1; b>=split; b--) *map++ = b;
X }
X else { /* minify: scan out away from split */
X for (b=split-1; b>=0; b--) *map++ = b;
X for (b=split; b<bn; b++) *map++ = b;
X }
X }
X else /* no overlap; either order will work */
X for (b=0; b<bn; b++) *map++ = b; /* we opt for forward order */
X}
X
Xstatic statistics()
X{
X# ifdef DEBUG
X printf("%d/%d=%.2f%% of filter coeffs were zero, nzmax=%d\n",
X nzero, ntot, 100.*nzero/ntot, nzmax);
X# endif
X}
EOF13088
echo extracting zoom/Makefile
sed 's/^X//' <<'EOF13089' >zoom/Makefile
X# $Header: Makefile,v 3.0 88/10/10 13:52:15 ph Locked $
X
XDEST = ..
XL = $(DEST)/lib
XBIN = $(DEST)/bin
XCOPTS = -g
XIPATH = -I$(DEST)/include
XCFLAGS = $(IPATH) $(COPTS)
XOBJ = zoom_main.o zoom.o scanline.o filt.o
XSRC = zoom_main.c zoom.c scanline.c filt.c
XEXTRALIB = -lgl_s
X# to pick up iris graphics library, use "make EXTRALIB=-lgl_s ..."
X
Xzoom: $(OBJ) $L/libpic.a
X cc $(COPTS) -o zoom $(OBJ) $L/libpic.a $(EXTRALIB) -lm
X
Xzoom_main.o:
X cc $(CFLAGS) -c zoom_main.c \
X -DDEFAULT_FILE='"iris"'
X# "note: we'll pick up pic_list from libpic.a(pic_all.o)"
X
Xzoom.lint: GHOST
X lint $(IPATH) $(SRC) >zoom.lint
X
Xinstall: zoom
X mv zoom $(BIN)
X
XGHOST:
X
Xzoom_main.o: zoom.h filt.h
Xzoom.o: zoom.h scanline.h filt.h
Xscanline.o: scanline.h
Xfilt.o: filt.h
X
Xclean:
X -rm -f zoom *.o *.lint
EOF13089
exit