home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / 3dscale_geom / scale_g.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-30  |  28.4 KB  |  855 lines
  1. /*
  2.  * zoom: subroutines for in-place one-pass filtered image zooming
  3.  *
  4.  * Zooms (resizes) one image into another with independent control of
  5.  * scale, translation, and filtering in x and y.  Magnifies or minifies.
  6.  * Filters with an arbitrary separable finite impulse response filter.
  7.  * (A filter h(x,y) is called separable if h(x,y)=f(x)*g(y)).
  8.  * See the filt package regarding filtering.
  9.  * The code is graphics device independent, using only generic scanline
  10.  * i/o operations.
  11.  *
  12.  * The program makes one pass over the image using a moving window of
  13.  * scanlines, so memory usage is proportional to imagewidth*filterheight,
  14.  * not imagewidth*imageheight, as you'd get if the entire image were
  15.  * buffered.  The separability of the filter is also exploited, to get
  16.  * a cpu time proportional to npixels*(filterwidth+filterheight),
  17.  * not npixels*filterwidth*filterheight as with direct 2-D convolution.
  18.  *
  19.  * The source and destination images can be identical, with overlapping
  20.  * windows, in which case the algorithm finds the magic split point and
  21.  * uses the appropriate scanning directions to avoid feedback (recursion)
  22.  * in the filtering.
  23.  *
  24.  * Paul Heckbert    12 Sept 1988
  25.  * UC Berkeley
  26.  * ph@miro.berkeley.edu
  27.  *
  28.  * Copyright (c) 1989  Paul S. Heckbert
  29.  * This source may be used for peaceful, nonprofit purposes only, unless
  30.  * under licence from the author. This notice should remain in the source.
  31.  */
  32.  
  33. /*
  34.  * routines in this source file:
  35.  *
  36.  * PUBLIC:
  37.  * zoom            most common entry point; window-to-window zoom
  38.  * zoom_continuous    fancier: explicit, continuous scale&tran control
  39.  *
  40.  * SORTA PRIVATE:
  41.  * zoom_unfiltered_mono    point-sampled zooming, called by zoom_continuous
  42.  * zoom_unfiltered_rgb    point-sampled zooming, called by zoom_continuous
  43.  * zoom_filtered    filtered zooming, called by zoom_continuous
  44.  * zoom_filtered_xy    zoom, filtering x before y, called by zoom_filtered
  45.  * zoom_filtered_yx    zoom, filtering y before x, called by zoom_filtered
  46.  *
  47.  * make_weighttab    make lookup table of filter coefficients
  48.  * make_map_table    in-place tricks; order scanlines to avoid feedback
  49.  */
  50.  
  51. static char rcsid[] = "$Header: zoom.c,v 3.0 88/10/10 13:52:07 ph Locked $";
  52.  
  53.  
  54. #include <math.h>
  55.  
  56. #include "simple.h"
  57. #include "pic.h"
  58. #include "filt.h"
  59. #include "scanline.h"
  60. #include "scaleg.h"
  61.  
  62. #define EPSILON 1e-7            /* error tolerance */
  63.  
  64. #define WEIGHTBITS  14            /* # bits in filter coefficients */
  65. #define FINALSHIFT  (2*WEIGHTBITS-CHANBITS) /* shift after x&y filter passes */
  66. #define WEIGHTONE  (1<<WEIGHTBITS)    /* filter weight of one */
  67.  
  68. #define INTEGER(x) (fabs((x)-floor((x)+.5)) < EPSILON)    /* is x an integer? */
  69. #define FRAC(x) fabs((x)-floor((x)+.5))        /* diff from closest integer */
  70.  
  71. int zoom_debug = 0;    /* debug level: 0=none, 1=scanline i/o, 2=filters */
  72. int zoom_coerce = 1;    /* simplify filters if possible? */
  73. int zoom_xy = UNDEF;    /* filter x before y (1) or vice versa (0)? */
  74. int zoom_trimzeros = 1;    /* trim zeros in x filter weight tables? */
  75.  
  76. static int nzero = 0, ntot = 0, nzmax = 0;
  77.  
  78. Filt *simplify_filter();
  79.  
  80. /*----------------------------------------------------------------------
  81.  * The mapping from source coordinates to dest coordinates:
  82.  * (Notation: prefix 'a' denotes source coords, 'b' denotes destination coords)
  83.  *
  84.  *    bx = sx*ax+tx
  85.  *    by = sy*ay+ty
  86.  *
  87.  * where (ax,ay) and (bx,by) are DISCRETE source and dest coords, respectively.
  88.  *
  89.  * An important fine point on terminology: there are two kinds of pixel coords:
  90.  *    DISCRETE COORDINATES take on integer values at pixel centers
  91.  *    CONTINUOUS COORDINATES take on integer values halfway between pixels
  92.  * For example, an image with discrete coordinate range [0..511] has a
  93.  * continuous coordinate range of [0..512].  The pixel with
  94.  * discrete coords (x,y) has its center at continuous coords (x+.5,y+.5)
  95.  * and its continuous coord domain is [x..x+1] in X and [y..y+1] in Y.
  96.  * Note: discrete coords are not always stored in ints and continuous coords
  97.  * are not always stored in floats.
  98.  *
  99.  * conversion:
  100.  * if c = continuous coord and d = discrete coord, then
  101.  *    c = d+.5
  102.  *    d = floor(c)
  103.  *
  104.  * To map a discrete src interval [a0..a1] to a discrete dst interval [b0..b1]:
  105.  *
  106.  *    b-b0+.5 = bn/an*(a-a0+.5)
  107.  *  or    b = (bn/an)*a + (b0-.5)-(bn/an)*(a0-.5)
  108.  *      = scale*a+tran
  109.  *
  110.  * where a and b are the discrete source and dest coords (either x or y)
  111.  * and an and bn are the interval lengths: an=a1-a0+1, bn=b1-b0+1.
  112.  * a0, an, b0, bn are the mapping parameters used by the zoom routine below.
  113.  * In general, however, for sub-pixel scale and translate control, we allow
  114.  * any real numbers for scale and tran (although there should usually be some
  115.  * correspondence between the mapping function and the src and dst windows).
  116.  *
  117.  * We'll often want the inverse mapping, from discrete dst coord b to
  118.  * continuous src coord ac, relative to the interval origins a0 and b0:
  119.  *
  120.  *    b+b0 = s*(ac+a0-.5)+t
  121.  *  so    ac = (b+b0-s*(a0-.5)-t)/s
  122.  *       = (b+offset)/scale = MAP(b, scale, offset)
  123.  *
  124.  * The mapping fields ux and uy in the Mapping structure
  125.  * are these offsets in x & y respectively.
  126.  */
  127.  
  128. /*----------------------------------------------------------------------
  129.  * zoom: simplest entry point; window-to-window zoom.
  130.  * Zoom window a of apic to window b of bpic using the specified filters.
  131.  * Window a represents the pixels with x in [a->x0..a->x1]
  132.  * and y in [a->y0..a->y1], (ranges are inclusive).  Likewise for b.
  133.  * apic and bpic may be equal and a and b may overlap, no sweat.
  134.  */
  135.  
  136. zoom(apic, a, bpic, b, xfilt, yfilt)
  137. Pic *apic, *bpic;    /* source and dest pictures */
  138. Window_box *a, *b;    /* source and dest windows */
  139. Filt *xfilt, *yfilt;    /* filters for x and y */
  140. {
  141.     Mapping map;
  142.  
  143.     if (b->x0==UNDEF) {
  144.     fprintf(stderr,"zoom: undefined window for dest %s: ",
  145.         pic_get_name(bpic));
  146.     window_print("", &b);
  147.     fprintf(stderr,"\n");
  148.     exit(1);
  149.     }
  150.     window_box_set_size(a);
  151.     window_box_set_size(b);
  152.     if (a->nx<=0 || a->ny<=0) return;
  153.     map.sx = (double)b->nx/a->nx;
  154.     map.sy = (double)b->ny/a->ny;
  155.     map.tx = b->x0-.5 - map.sx*(a->x0-.5);
  156.     map.ty = b->y0-.5 - map.sy*(a->y0-.5);
  157.     zoom_continuous(apic, a, bpic, b, &map, xfilt, yfilt);
  158. }
  159.  
  160. /*
  161.  * zoom_opt: window-to-window zoom with options.
  162.  * Like zoom() but has options to make mapping square (preserve pixel shape)
  163.  * and round to the next lowest integer scale factor.
  164.  */
  165.  
  166. zoom_opt(apic, a, bpic, b, xfilt, yfilt, square, intscale)
  167. Pic *apic, *bpic;    /* source and dest pictures */
  168. Window_box *a, *b;    /* source and dest windows */
  169. Filt *xfilt, *yfilt;    /* filters for x and y */
  170. int square, intscale;    /* square mapping? integer scales? */
  171. {
  172.     Mapping map;
  173.  
  174.     if (b->x0==UNDEF) {
  175.     fprintf(stderr,"zoom_opt: undefined window for dest %s: ",
  176.         pic_get_name(bpic));
  177.     window_print("", &b);
  178.     fprintf(stderr,"\n");
  179.     exit(1);
  180.     }
  181.     window_box_set_size(a);
  182.     window_box_set_size(b);
  183.     if (a->nx<=0 || a->ny<=0) return;
  184.     map.sx = (double)b->nx/a->nx;
  185.     map.sy = (double)b->ny/a->ny;
  186.     if (square) {
  187.     if (map.sx>map.sy) {    /* use the smaller scale for both sx and sy */
  188.         map.sx = map.sy;
  189.         b->nx = ceil(a->nx*map.sx);
  190.     }
  191.     else {
  192.         map.sy = map.sx;
  193.         b->ny = ceil(a->ny*map.sy);
  194.     }
  195.     }
  196.     if (intscale) {
  197.     if (map.sx>1.) {
  198.         map.sx = floor(map.sx);
  199.         b->nx = ceil(a->nx*map.sx);
  200.     }
  201.     if (map.sy>1.) {
  202.         map.sy = floor(map.sy);
  203.         b->ny = ceil(a->ny*map.sy);
  204.     }
  205.     }
  206.     window_box_set_max(b);
  207.     map.tx = b->x0-.5 - map.sx*(a->x0-.5);
  208.     map.ty = b->y0-.5 - map.sy*(a->y0-.5);
  209.     zoom_continuous(apic, a, bpic, b, &map, xfilt, yfilt);
  210. }
  211.  
  212. /*
  213.  * zoom_continuous: zoom with explicit, continuous scale&tran control.
  214.  * Zoom according to the mapping defined in m,
  215.  * reading from window awin of apic and writing to window bwin of bpic.
  216.  * Like zoom but allows continuous, subpixel scales and translates.
  217.  * Meaning of m is explained above.  Pixels in the dest window bwin outside
  218.  * of the support of the mapped source window will be left unchanged.
  219.  * Scales in m should be positive.
  220.  */
  221.  
  222. zoom_continuous(apic, awin, bpic, bwin, m, xfilt, yfilt)
  223. Pic *apic, *bpic;
  224. Window_box *awin, *bwin;
  225. Mapping *m;
  226. Filt *xfilt, *yfilt;
  227. {
  228.     int xy;
  229.     Window_box a, b, bc, t;
  230.     Filtpar ax, ay;
  231.  
  232.     if (m->sx<=0. || m->sy<=0.) {
  233.     fprintf(stderr, "zoom_continuous: negative scales not allowed\n");
  234.     return;
  235.     }
  236.  
  237.     a = *awin;
  238.     b = *bwin;
  239.     window_box_set_size(&a);
  240.     window_box_set_size(&b);
  241.  
  242.     /*
  243.      * find scale of filter in a space (source space)
  244.      * when minifying, ascale=1/scale, but when magnifying, ascale=1
  245.      */
  246.     ax.scale = xfilt->blur*MAX(1., 1./m->sx);
  247.     ay.scale = yfilt->blur*MAX(1., 1./m->sy);
  248.     /*
  249.      * find support radius of scaled filter
  250.      * if ax.supp and ay.supp are both <=.5 then we've got point sampling.
  251.      * Point sampling is essentially a special filter whose width is fixed
  252.      * at one source pixel.
  253.      */
  254.     ax.supp = MAX(.5, ax.scale*xfilt->supp);
  255.     ay.supp = MAX(.5, ay.scale*yfilt->supp);
  256.     ax.wid = ceil(2.*ax.supp);
  257.     ay.wid = ceil(2.*ay.supp);
  258.  
  259.     /* clip source and dest windows against their respective pictures */
  260.     window_box_intersect(&a, (Window_box *)pic_get_window(apic, &t), &a);
  261.     if (pic_get_window(bpic, &t)->x0 != UNDEF)
  262.     window_box_intersect(&b, &t, &b);
  263.  
  264.     /* find valid dest window (transformed source + support margin) */
  265.     bc.x0 =  ceil(m->sx*(a.x0-ax.supp)+m->tx+EPSILON);
  266.     bc.y0 =  ceil(m->sy*(a.y0-ay.supp)+m->ty+EPSILON);
  267.     bc.x1 = floor(m->sx*(a.x1+ax.supp)+m->tx-EPSILON);
  268.     bc.y1 = floor(m->sy*(a.y1+ay.supp)+m->ty-EPSILON);
  269.     window_box_set_size(&bc);
  270.  
  271.     if (b.x0<bc.x0 || b.x1>bc.x1 || b.y0<bc.y0 || b.y1>bc.y1) {
  272.     /* requested dest window lies outside the valid dest, so clip dest */
  273.         fprintf(stderr,"\n Error: destination window not valid. \n\n");
  274.         exit(1);
  275.  
  276.       /* Note: can't clip the dest window because then the hips
  277.           header will be wrong ( Probably could re-write the header
  278.            here, but I'm not going to worry about it until a need arises -BT
  279.  
  280.     window_box_print("    clipping odst=", &b);
  281.     window_box_print(" to ", &bc);
  282.     fprintf(stderr,"\n");
  283.     window_box_intersect(&b, &bc, &b);
  284.       */
  285.     }
  286.     pic_set_window(bpic, &b);
  287.  
  288.     /* compute offsets for MAP (these will be .5 if zoom() routine was called)*/
  289.     m->ux = b.x0-m->sx*(a.x0-.5)-m->tx;
  290.     m->uy = b.y0-m->sy*(a.y0-.5)-m->ty;
  291.  
  292. #ifdef DEBUG
  293.     fprintf(stderr,"src=%s:%s", pic_get_dev(apic), pic_get_name(apic));
  294.     window_box_print("", &a);
  295.     fprintf(stderr,"\ndst=%s:%s", pic_get_dev(bpic), pic_get_name(bpic));
  296.     window_box_print("", &b);
  297.     fprintf(stderr,"\n");
  298. #endif
  299.  
  300.     if (a.nx<=0 || a.ny<=0 || b.nx<=0 || b.ny<=0) {
  301.       fprintf(stderr,"\n Error: negative window size. \n" );
  302.       return;
  303.     }
  304.  
  305.     /* check for high-level simplifications of filter */
  306.     xfilt = simplify_filter("x", zoom_coerce, xfilt, &ax, m->sx, m->ux);
  307.     yfilt = simplify_filter("y", zoom_coerce, yfilt, &ay, m->sy, m->uy);
  308.  
  309.     /*
  310.      * decide which filtering order (xy or yx) is faster for this mapping by
  311.      * counting convolution multiplies
  312.      */
  313.     xy = zoom_xy!=UNDEF ? zoom_xy :
  314.     b.nx*(a.ny*ax.wid+b.ny*ay.wid) < b.ny*(b.nx*ax.wid+a.nx*ay.wid);
  315.  
  316.     /* choose the appropriate filtering routine */
  317.     if (str_eq(xfilt->name, "point") && str_eq(yfilt->name, "point"))
  318.     zoom_unfiltered(apic, &a, bpic, &b, m);
  319.     else if (xy)
  320.     zoom_filtered_xy(apic, &a, bpic, &b, m, xfilt, &ax, yfilt, &ay);
  321.     else
  322.     zoom_filtered_yx(apic, &a, bpic, &b, m, xfilt, &ax, yfilt, &ay);
  323. }
  324.  
  325. /*
  326.  * filter_simplify: check if our discrete sampling of an arbitrary continuous
  327.  * filter, parameterized by the filter spacing (a->scale), its radius (a->supp),
  328.  * and the scale and offset of the coordinate mapping (s and u), causes the
  329.  * filter to reduce to point sampling.
  330.  *
  331.  * It reduces if support is less than 1 pixel or
  332.  * if integer scale and translation, and filter is cardinal
  333.  */
  334.  
  335. static Filt *simplify_filter(dim, coerce, filter, a, s, u)
  336. char *dim;
  337. int coerce;
  338. Filt *filter;
  339. Filtpar *a;
  340. double s, u;    /* scale and offset */
  341. {
  342.     if (coerce && (a->supp<=.5 ||
  343.     filter->cardinal && INTEGER(1./a->scale) && INTEGER(1./(s*a->scale))
  344.     && INTEGER((u/s-.5)/a->scale))) {
  345.         if (!str_eq(filter->name, "point"))
  346.         fprintf(stderr,"coercing %sfilter=%s to point\n", dim, filter->name);
  347.         filter = filt_find("point");
  348.         a->scale = 1.;
  349.         a->supp = .5;
  350.         a->wid = 1;
  351.     }
  352.     return filter;
  353. }
  354.  
  355. /*----------------------------------------------------------------------*/
  356.  
  357. /* zoom_unfiltered: do unfiltered zoom (point sampling) */
  358.  
  359. zoom_unfiltered(apic, a, bpic, b, m)
  360. Pic *apic, *bpic;
  361. Window_box *a, *b;
  362. Mapping *m;
  363. {
  364.     int overlap;
  365.  
  366.     /* do source and dest windows overlap? */
  367.     overlap = apic==bpic && window_box_overlap(a, b);
  368.  
  369. #ifdef DEBUG
  370.     fprintf(stderr,"-map (sx,sy,tx,ty)= %.2g %.2g %.2g %.2g \n",
  371.      m->sx, m->sy, m->tx, m->ty);
  372. #endif
  373.     /*
  374.      * note: We have only x-resample before y-resample versions of the
  375.      * unfiltered zoom case; we could optimize further by creating
  376.      * y-resample before x-resample versions.
  377.      */
  378.  
  379.     zoom_unfiltered_mono(apic, a, bpic, b, m, overlap);
  380. }
  381.  
  382. /* zoom_unfiltered_mono: monochrome unfiltered zoom */
  383.  
  384. zoom_unfiltered_mono(apic, a, bpic, b, m, overlap)
  385. Pic *apic, *bpic;
  386. Window_box *a, *b;
  387. Mapping *m;
  388. int overlap;
  389. {
  390.     int byi, by, ay, ayold, bx, ax;
  391.     Pixel1 *abuf, *bbuf, *bp;    /* source and dest scanline buffers */
  392.     Pixel1 **xmap, **xp;    /* mapping from abuf to bbuf (optimization) */
  393.     short *ymap;        /* order of dst y coords that avoids feedback */
  394.     char *linewritten;        /* has scanline y been written? (debugging) */
  395.  
  396.     ALLOC(abuf, Pixel1, a->nx);
  397.     ALLOC(bbuf, Pixel1, b->nx);
  398.     ALLOC(xmap, Pixel1 *, b->nx);
  399.     ALLOC(ymap, short, b->ny);
  400.     ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
  401.  
  402.     /* if overlapping src & dst, find magic y ordering that avoids feedback */
  403.     make_map_table(m->sy, m->ty, .5, a->y0, b->y0, b->ny, overlap, ymap);
  404.  
  405.     for (bx=0; bx<b->nx; bx++) {
  406.     ax = MAP(bx, m->sx, m->ux);
  407.     xmap[bx] = &abuf[ax];
  408.     }
  409.  
  410.     ayold = -1;        /* impossible value for ay */
  411.     for (byi=0; byi<b->ny; byi++) {
  412.     by = ymap[byi];
  413.     ay = MAP(by, m->sy, m->uy);
  414.     /* scan line a.y0+ay goes to b.y0+by */
  415.     if (ay!=ayold) {        /* new scan line; read it in */
  416.         ayold = ay;
  417.         if (overlap && linewritten[a->y0+ay])
  418.         fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ay);
  419.         /* read scan line into abuf */
  420.         pic_read_row(apic, a->y0+ay, a->x0, a->nx, abuf);
  421.         /* resample in x */
  422.         for (bp=bbuf, xp=xmap, bx=0; bx<b->nx; bx++)
  423.         *bp++ = **xp++;
  424.     }
  425.     pic_write_row(bpic, b->y0+by, b->x0, b->nx, bbuf);
  426.     linewritten[b->y0+by] = 1;
  427.     }
  428.     free(abuf);
  429.     free(bbuf);
  430.     free(xmap);
  431.     free(ymap);
  432. }
  433. /*----------------------------------------------------------------------*/
  434.  
  435. /*
  436.  * zoom_filtered_xy: filtered zoom, xfilt before yfilt
  437.  *
  438.  * note: when calling make_weighttab, we can trim leading and trailing
  439.  * zeros from the x weight buffers as an optimization,
  440.  * but not for y weight buffers since the split formula is anticipating
  441.  * a constant amount of buffering of source scanlines;
  442.  * trimming zeros in yweight could cause feedback.
  443.  */
  444.  
  445.  
  446. zoom_filtered_xy(apic, a, bpic, b, m, xfilt, ax, yfilt, ay)
  447. Pic *apic, *bpic;
  448. Window_box *a, *b;
  449. Mapping *m;
  450. Filt *xfilt, *yfilt;    /* x and y filters */
  451. Filtpar *ax, *ay;    /* extra x and y filter parameters */
  452. {
  453.     int ayf, bx, byi, by, overlap, nchan;
  454.  
  455.                 /*PIXELTYPE NBITS  RAISON D'ETRE */
  456.     Scanline abbuf;        /*   PIXEL1   8  src&dst scanline bufs */
  457.     Scanline accum;        /*   PIXEL4  28  accumulator buf for yfilt */
  458.     Scanline *linebuf, *tbuf;    /*   PIXEL2  14  circular buf of active lines */
  459.  
  460.     short *ymap;        /* order of dst y coords that avoids feedback */
  461.     Weighttab *xweights;    /* sampled filter at each dst x pos; for xfilt*/
  462.     short *xweightbuf, *xwp;    /* big block of memory addressed by xweights */
  463.     Weighttab yweight;        /* a single sampled filter for current y pos */
  464.     char *linewritten;        /* has scanline y been written? (debugging) */
  465.  
  466.     nchan = pic_get_nchan(apic)==1 ? PIXEL_MONO : PIXEL_RGB;
  467.     scanline_alloc(&abbuf, PIXEL1|nchan, MAX(a->nx, b->nx));
  468.     scanline_alloc(&accum, PIXEL4|nchan, b->nx);
  469.     ALLOC(linebuf, Scanline, ay->wid);
  470.     /* construct circular buffer of ay->wid intermediate scanlines */
  471.     for (ayf=0; ayf<ay->wid; ayf++) {
  472.     scanline_alloc(&linebuf[ayf], PIXEL2|nchan, b->nx);
  473.     linebuf[ayf].y = -1;        /* mark scanline as unread */
  474.     }
  475.  
  476.     ALLOC(ymap, short, b->ny);
  477.     ALLOC(xweights, Weighttab, b->nx);
  478.     ALLOC(xweightbuf, short, b->nx*ax->wid);
  479.     ALLOC(yweight.weight, short, ay->wid);
  480.     ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
  481.  
  482.     /* do source and dest windows overlap? */
  483.     overlap = apic==bpic && window_box_overlap(a, b);
  484.  
  485. #ifdef DEBUG
  486.     fprintf(stderr,"-xy -map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
  487.     fprintf(stderr,"X: filt=%s support=%g,  scale=%g scaledsupport=%g width=%d\n",
  488.     xfilt->name, xfilt->supp, ax->scale, ax->supp, ax->wid);
  489.     fprintf(stderr,"Y: filt=%s support=%g,  scale=%g scaledsupport=%g width=%d\n",
  490.     yfilt->name, yfilt->supp, ay->scale, ay->supp, ay->wid);
  491. #endif
  492.  
  493.     /*
  494.      * prepare a weighttab (a sampled filter for source pixels) for
  495.      * each dest x position
  496.      */
  497.     for (xwp=xweightbuf, bx=0; bx<b->nx; bx++, xwp+=ax->wid) {
  498.     xweights[bx].weight = xwp;
  499.     make_weighttab(b->x0+bx, MAP(bx, m->sx, m->ux),
  500.         xfilt, ax, a->nx, zoom_trimzeros, &xweights[bx]);
  501.     }
  502.  
  503.     /* if overlapping src & dst, find magic y ordering that avoids feedback */
  504.     make_map_table(m->sy, m->ty, ay->supp, a->y0, b->y0, b->ny, overlap,
  505.     ymap);
  506.  
  507.     for (byi=0; byi<b->ny; byi++) {    /* loop over dest scanlines */
  508.     by = ymap[byi];
  509.     if (zoom_debug) fprintf(stderr,"by=%d: ", b->y0+by);
  510.     /* prepare a weighttab for dest y position by */
  511.     make_weighttab(b->y0+by, MAP(by, m->sy, m->uy),
  512.         yfilt, ay, a->ny, 0, &yweight);
  513.     if (zoom_debug) fprintf(stderr,"ay=%d-%d, reading: ",
  514.         a->y0+yweight.i0, a->y0+yweight.i1-1);
  515.     scanline_zero(&accum);
  516.  
  517.     /* loop over source scanlines that influence this dest scanline */
  518.     for (ayf=yweight.i0; ayf<yweight.i1; ayf++) {
  519.         tbuf = &linebuf[ayf % ay->wid];
  520.         if (tbuf->y != ayf) {    /* scanline needs to be read in */
  521.         if (zoom_debug) fprintf(stderr,"%d ", a->y0+ayf);
  522.         if (overlap && linewritten[a->y0+ayf])
  523.             fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ayf);
  524.         tbuf->y = ayf;
  525.         /* read source scanline into abbuf */
  526.         abbuf.len = a->nx;
  527.         scanline_read(apic, a->x0, a->y0+ayf, &abbuf);
  528.         /* and filter it into the appropriate line of linebuf (xfilt) */
  529.         scanline_filter(CHANBITS, xweights, &abbuf, tbuf);
  530.         }
  531.         /* add weighted tbuf into accum (these do yfilt) */
  532.         scanline_accum(yweight.weight[ayf-yweight.i0], tbuf, &accum);
  533.     }
  534.  
  535.     /* shift accum down into abbuf and write out dest scanline */
  536.     abbuf.len = b->nx;
  537.     scanline_shift(FINALSHIFT, &accum, &abbuf);
  538.     scanline_write(bpic, b->x0, b->y0+by, &abbuf);
  539.     linewritten[b->y0+by] = 1;
  540.     if (zoom_debug) fprintf(stderr,"\n");
  541.     }
  542.  
  543.     scanline_free(&abbuf);
  544.     scanline_free(&accum);
  545.     for (ayf=0; ayf<ay->wid; ayf++)
  546.     scanline_free(&linebuf[ayf]);
  547.     free(ymap);
  548.     free(linebuf);
  549.     free(xweights);
  550.     free(xweightbuf);
  551.     free(yweight.weight);
  552.     free(linewritten);
  553.     statistics();
  554. }
  555.  
  556. /* zoom_filtered_yx: filtered zoom, yfilt before xfilt */
  557.  
  558. zoom_filtered_yx(apic, a, bpic, b, m, xfilt, ax, yfilt, ay)
  559. Pic *apic, *bpic;
  560. Window_box *a, *b;
  561. Mapping *m;
  562. Filt *xfilt, *yfilt;    /* x and y filters */
  563. Filtpar *ax, *ay;    /* extra x and y filter parameters */
  564. {
  565.     int ayf, bx, byi, by, overlap, nchan;
  566.  
  567.                 /*PIXELTYPE NBITS  RAISON D'ETRE */
  568.     Scanline bbuf;        /*   PIXEL1   8  dst scanline buf */
  569.     Scanline accum;        /*   PIXEL4  22  accumulator buf for yfilt */
  570.     Scanline *linebuf, *tbuf;    /*   PIXEL1   8  circular buf of active lines */
  571.  
  572.     short *ymap;        /* order of dst y coords that avoids feedback */
  573.     Weighttab *xweights;    /* sampled filter at each dst x pos; for xfilt*/
  574.     short *xweightbuf, *xwp;    /* big block of memory addressed by xweights */
  575.     Weighttab yweight;        /* a single sampled filter for current y pos */
  576.     char *linewritten;        /* has scanline y been written? (debugging) */
  577.  
  578.     nchan = pic_get_nchan(apic)==1 ? PIXEL_MONO : PIXEL_RGB;
  579.     scanline_alloc(&bbuf, PIXEL1|nchan, b->nx);
  580.     scanline_alloc(&accum, PIXEL4|nchan, a->nx);
  581.     ALLOC(linebuf, Scanline, ay->wid);
  582.     /* construct circular buffer of ay->wid intermediate scanlines */
  583.     for (ayf=0; ayf<ay->wid; ayf++) {
  584.     scanline_alloc(&linebuf[ayf], PIXEL1|nchan, a->nx);
  585.     linebuf[ayf].y = -1;        /* mark scanline as unread */
  586.     }
  587.  
  588.     ALLOC(ymap, short, b->ny);
  589.     ALLOC(xweights, Weighttab, b->nx);
  590.     ALLOC(xweightbuf, short, b->nx*ax->wid);
  591.     ALLOC(yweight.weight, short, ay->wid);
  592.     ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
  593.  
  594.     /* do source and dest windows overlap? */
  595.     overlap = apic==bpic && window_box_overlap(a, b);
  596.  
  597. #ifdef DEBUG
  598.     fprintf(stderr,"-yx -map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
  599.     fprintf(stderr,"X: filt=%s supp=%g  scale=%g scaledsupp=%g wid=%d\n",
  600.     xfilt->name, xfilt->supp, ax->scale, ax->supp, ax->wid);
  601.     fprintf(stderr,"Y: filt=%s supp=%g  scale=%g scaledsupp=%g wid=%d\n",
  602.     yfilt->name, yfilt->supp, ay->scale, ay->supp, ay->wid);
  603. #endif
  604.  
  605.     /*
  606.      * prepare a weighttab (a sampled filter for source pixels) for
  607.      * each dest x position
  608.      */
  609.     for (xwp=xweightbuf, bx=0; bx<b->nx; bx++, xwp+=ax->wid) {
  610.     xweights[bx].weight = xwp;
  611.     make_weighttab(b->x0+bx, MAP(bx, m->sx, m->ux),
  612.         xfilt, ax, a->nx, zoom_trimzeros, &xweights[bx]);
  613.     }
  614.  
  615.     /* if overlapping src & dst, find magic y ordering that avoids feedback */
  616.     make_map_table(m->sy, m->ty, ay->supp, a->y0, b->y0, b->ny, overlap,
  617.     ymap);
  618.  
  619.     for (byi=0; byi<b->ny; byi++) {     /* loop over dest scanlines */
  620.     by = ymap[byi];
  621.     if (zoom_debug) fprintf(stderr,"by=%d: ", b->y0+by);
  622.     /* prepare a weighttab for dest y position by */
  623.     make_weighttab(b->y0+by, MAP(by, m->sy, m->uy),
  624.         yfilt, ay, a->ny, 0, &yweight);
  625.     if (zoom_debug) fprintf(stderr,"ay=%d-%d, reading: ",
  626.         a->y0+yweight.i0, a->y0+yweight.i1-1);
  627.     scanline_zero(&accum);
  628.  
  629.     /* loop over source scanlines that influence this dest scanline */
  630.     for (ayf=yweight.i0; ayf<yweight.i1; ayf++) {
  631.         tbuf = &linebuf[ayf % ay->wid];
  632.         if (tbuf->y != ayf) {    /* scanline needs to be read in */
  633.         if (zoom_debug) fprintf(stderr,"%d ", a->y0+ayf);
  634.         if (overlap && linewritten[a->y0+ayf])
  635.             fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ayf);
  636.         tbuf->y = ayf;
  637.         /* read source scanline into linebuf */
  638.         scanline_read(apic, a->x0, a->y0+ayf, tbuf);
  639.         }
  640.         /* add weighted tbuf into accum (these do yfilt) */
  641.         scanline_accum(yweight.weight[ayf-yweight.i0], tbuf, &accum);
  642.     }
  643.  
  644.     /* and filter it into the appropriate line of linebuf (xfilt) */
  645.     scanline_filter(FINALSHIFT, xweights, &accum, &bbuf);
  646.     /* and write out dest scanline in bbuf */
  647.     scanline_write(bpic, b->x0, b->y0+by, &bbuf);
  648.     linewritten[b->y0+by] = 1;
  649.     if (zoom_debug) fprintf(stderr,"\n");
  650.     }
  651.  
  652.     scanline_free(&bbuf);
  653.     scanline_free(&accum);
  654.     for (ayf=0; ayf<ay->wid; ayf++)
  655.     scanline_free(&linebuf[ayf]);
  656.     free(ymap);
  657.     free(linebuf);
  658.     free(xweights);
  659.     free(xweightbuf);
  660.     free(yweight.weight);
  661.     free(linewritten);
  662.     statistics();
  663. }
  664.  
  665. /*
  666.  * make_weighttab: sample the continuous filter, scaled by ap->scale and
  667.  * positioned at continuous source coordinate cen, for source coordinates in
  668.  * the range [0..len-1], writing the weights into wtab.
  669.  * Scale the weights so they sum to WEIGHTONE, and trim leading and trailing
  670.  * zeros if trimzeros!=0.
  671.  * b is the dest coordinate (for diagnostics).
  672.  */
  673.  
  674. make_weighttab(b, cen, filter, ap, len, trimzeros, wtab)
  675. int b, len;
  676. double cen;
  677. Filt *filter;
  678. Filtpar *ap;
  679. int trimzeros;
  680. Weighttab *wtab;
  681. {
  682.     int i0, i1, i, sum, t, stillzero, lastnonzero, nz;
  683.     short *wp;
  684.     double den, sc, tr;
  685.  
  686.     /* find the source coord range of this positioned filter: [i0..i1-1] */
  687.     i0 = cen-ap->supp+.5;
  688.     i1 = cen+ap->supp+.5;
  689.     if (i0<0) i0 = 0;
  690.     if (i1>len) i1 = len;
  691.     if (i0>=i1) {
  692.     fprintf(stderr, "make_weighttab: null filter at %d\n", b);
  693.     exit(1);
  694.     }
  695.     /* the range of source samples to buffer: */
  696.     wtab->i0 = i0;
  697.     wtab->i1 = i1;
  698.  
  699.     /* find scale factor sc to normalize the filter */
  700.     for (den=0, i=i0; i<i1; i++)
  701.     den += filt_func(filter, (i+.5-cen)/ap->scale);
  702.     /* set sc so that sum of sc*func() is approximately WEIGHTONE */
  703.     sc = den==0. ? WEIGHTONE : WEIGHTONE/den;
  704.     if (zoom_debug>1) fprintf(stderr,"    b=%d cen=%g scale=%g [%d..%d) sc=%g:  ",
  705.     b, cen, ap->scale, i0, i1, sc);
  706.  
  707.     /* compute the discrete, sampled filter coefficients */
  708.     stillzero = trimzeros;
  709.     for (sum=0, wp=wtab->weight, i=i0; i<i1; i++) {
  710.  
  711.     /* evaluate the filter function: */
  712.     tr = sc*filt_func(filter, (i+.5-cen)/ap->scale);
  713.  
  714.     if (tr<MINSHORT || tr>MAXSHORT) {
  715.         fprintf(stderr, "tr=%g at %d\n", tr, b);
  716.         exit(1);
  717.     }
  718.     t = floor(tr+.5);
  719.     if (stillzero && t==0) i0++;    /* find first nonzero */
  720.     else {
  721.         stillzero = 0;
  722.         *wp++ = t;            /* add weight to table */
  723.         sum += t;
  724.         if (t!=0) lastnonzero = i;    /* find last nonzero */
  725.     }
  726.     }
  727.     ntot += wtab->i1-wtab->i0;
  728.     if (sum==0) {
  729.     nz = wtab->i1-wtab->i0;
  730.     /* fprintf(stderr, "sum=0 at %d\n", b); */
  731.     wtab->i0 = wtab->i0+wtab->i1 >> 1;
  732.     wtab->i1 = wtab->i0+1;
  733.     wtab->weight[0] = WEIGHTONE;
  734.     }
  735.     else {
  736.     if (trimzeros) {        /* skip leading and trailing zeros */
  737.         /* set wtab->i0 and ->i1 to the nonzero support of the filter */
  738.         nz = wtab->i1-wtab->i0-(lastnonzero-i0+1);
  739.         wtab->i0 = i0;
  740.         wtab->i1 = i1 = lastnonzero+1;
  741.     }
  742.     else                /* keep leading and trailing zeros */
  743.         nz = 0;
  744.     if (sum!=WEIGHTONE) {
  745.         /*
  746.          * Fudge the center slightly to make sum=WEIGHTONE exactly.
  747.          * Is this the best way to normalize a discretely sampled
  748.          * continuous filter?
  749.          */
  750.         i = cen+.5;
  751.         if (i<i0) i = i0; else if (i>=i1) i = i1-1;
  752.         t = WEIGHTONE-sum;
  753.         if (zoom_debug>1) fprintf(stderr,"[%d]+=%d ", i, t);
  754.         wtab->weight[i-i0] += t;    /* fudge center sample */
  755.     }
  756.     }
  757.     if (nz>nzmax) nzmax = nz;
  758.     nzero += nz;
  759.     if (zoom_debug>1) {
  760.     fprintf(stderr,"\t");
  761.     for (wp=wtab->weight, i=i0; i<i1; i++, wp++)
  762.         fprintf(stderr,"%5d ", *wp);
  763.     fprintf(stderr,"\n");
  764.     }
  765. }
  766.  
  767. /*
  768.  * make_map_table:
  769.  * construct a table for the mapping a = (b-t)/s for b-b0 in [0..bn-1]
  770.  * ordered so that the buffer resampling operation
  771.  *        for (bi=0; bi<bn; bi++) {
  772.  *        b = map[bi];
  773.  *        a = MAP(b, scale, offset);
  774.  *        buf[b] = buf[a];
  775.  *        }
  776.  * can work IN PLACE without feedback
  777.  *
  778.  * a and b here are source and dest coordinates, in either x or y.
  779.  * This is needed only if there is overlap between source and dest windows.
  780.  */
  781.  
  782. static make_map_table(scale, tran, asupp, a0, b0, bn, overlap, map)
  783. double scale, tran;    /* scale and translate */
  784. double asupp;        /* filter support radius in source space */
  785. int a0, b0, bn;        /* source and dest positions, dest length */
  786. int overlap;        /* do source and dest overlap? */
  787. short *map;        /* mapping table */
  788. {
  789.     int split, b, i0;
  790.     double z, u;
  791.  
  792.     /* find fixed point of mapping; let split=b-b0 at the point where a=b */
  793.  
  794.     if (overlap) {            /* source and dest windows overlap */
  795.     if (scale==1.)            /* no scale change, translation only */
  796.         /* if moving left then scan to right, and vice versa */
  797.         split = a0<b0 ? 0 : bn;
  798.     else {                /* magnify or minify */
  799.  
  800.         /* THE MAGIC SPLIT FORMULA: */
  801.  
  802.         if (scale>1.) {        /* magnify */
  803.         /*
  804.          * choose integer nearest midpoint of valid interval:
  805.          *    x < split <= x+s/(s-1)
  806.          * where x=(tran+scale*asupp)/(1-scale)-b0
  807.          */
  808.         split = floor((tran+scale*asupp-.5)/(1.-scale)-b0+1.);
  809.         }
  810.         else {            /* minify */
  811.         /*
  812.          * only one valid split pt in this case:
  813.          *    x <= split < x+1
  814.          * so we take extra care (x as above)
  815.          */
  816.         split = ceil((tran+scale*asupp)/(1.-scale)-b0);
  817.         /*
  818.          * The above formula is perfect for exact arithmetic,
  819.          * but hardware roundoff can cause split to be one too big.
  820.          * Check for roundoff by simulating precisely the calculation
  821.          * of i0 in make_weighttab.
  822.          */
  823.         u = b0-scale*(a0-.5)-tran;    /* recalculate offset */
  824.         z = MAP(split-1, scale, u);    /* cen at b=split-1 */
  825.         z = z-asupp+.5;
  826.         i0 = z;                /* i0 at b=split-1 */
  827.         if (a0+i0>=b0+split)        /* feedback at b=split-1? */
  828.             split--;            /* correct for roundoff */
  829.         }
  830.         if (split<0) split = 0;
  831.         else if (split>bn) split = bn;
  832.         fprintf(stderr,"split at y=%d\n", b0+split);
  833.     }
  834.  
  835.     if (scale>=1.) {        /* magnify: scan in toward split */
  836.         for (b=0;    b<split;  b++) *map++ = b;
  837.         for (b=bn-1; b>=split; b--) *map++ = b;
  838.     }
  839.     else {                /* minify: scan out away from split */
  840.         for (b=split-1; b>=0;  b--) *map++ = b;
  841.         for (b=split;   b<bn;  b++) *map++ = b;
  842.     }
  843.     }
  844.     else                /* no overlap; either order will work */
  845.     for (b=0; b<bn; b++) *map++ = b;    /* we opt for forward order */
  846. }
  847.  
  848. static statistics()
  849. {
  850. #   ifdef DEBUG
  851.     fprintf(stderr,"%d/%d=%.2f%% of filter coeffs were zero, nzmax=%d\n",
  852.         nzero, ntot, 100.*nzero/ntot, nzmax);
  853. #   endif
  854. }
  855.