#ifndef lint
static char sccsid[] = "@(#)image.c 9.11 88/02/23 Copyright 1985 Sun Micro";
#endif

/*
 * Copyright (c) 1985 by Sun Microsystems, Inc.
 */

/*-
	Image manipulation routines

	image.c, Mon Jul  7 16:54:34 1986

 */

#include <psio.h>
#include <sys/types.h>
#ifdef REF
#include <ref/config.h>
#endif
#include <pixrect/pixrect.h>
#include <pixrect/pr_io.h>
#include <pixrect/pr_util.h>
#include <rasterfile.h>
#include <pixrect/pr_io.h>
#include <pixrect/memvar.h>
#include "shape.h"
#include "canvas.h"
#include "cursor.h"
#include "cscript.h"


extern unsigned char
            red_map[256],
            green_map[256],
            blue_map[256],
            red_inverse[256],
            green_inverse[256],
            blue_inverse[256],
            grey_inverse[256];


#ifdef LITTLEENDIAN
#define PIXBIT(x) (1<<(x&0xF))
#define PIXLEFT(m) ((m)>>1)
#define PIXRIGHT(m) (unsigned short)((m)<<1)
#else
#define PIXBIT(x) (0x8000>>(x&0xF))
#define PIXLEFT(m) (unsigned short)((m)<<1)
#define PIXRIGHT(m) ((m)>>1)
#endif

#define RIGHTSTEP(p,m) ((m = PIXRIGHT(m)) == 0 ? m = PIXBIT(0), p++ : 0)
#define LEFTSTEP(p,m) ((m = PIXLEFT(m)) == 0 ? m = PIXBIT(15), p-- : 0)

struct pixrect *cv_dither24to8();
struct pixrect *cv_dither24to1();
struct pixrect *cv_dither8to1();
struct pixrect *cs_scalepixrect();
struct pixel {
    unsigned char red,
                green,
                blue;
};

static colormap_t colormap = {1, 256, red_map, green_map, blue_map};

struct spoint rotated_upper_left;
struct pixrect *cs_scale_and_rotate();

cs_imagecanvas(gc, cv, imagekind)
    register struct graphics_context *gc;
    struct pixrect *cv;
    int         imagekind;	/* 0=>image; 1 implies mask; -1=>inverted mask */
{
    extern      mem_rop();
#ifdef REF
#ifdef SUNREF
    extern	magic_hack_bgcolor;
#else
    int	        magic_hack_bgcolor;
#endif
#else
    extern      magic_hack_bgcolor;
#endif
    struct pixrect *in;
    struct pixrect *new = 0;
    int         opcode;
    struct shape *outline;
    cs_cursorup();
    if (imagekind == 0)
	opcode = PIX_COLOR(black_color.hint) | PIX_SRC;
    else if (imagekind > 0)
	if (gc->canvas->pr_depth == 1)
	    if (gc->func == 0)
		opcode = PIX_NOT(PIX_SRC) & PIX_DST;
	    else
		opcode = PIX_SRC | PIX_DST;
	else
	    opcode = PIX_COLOR(gc->color.hint) | PIX_SRC | PIX_DST;
    else if (gc->canvas->pr_depth == 1)
	if (gc->func == 0)
	    opcode = PIX_SRC & PIX_DST;
	else
	    opcode = PIX_NOT(PIX_SRC) | PIX_DST;
    else
	opcode = PIX_COLOR(gc->color.hint) | PIX_NOT(PIX_SRC) | PIX_DST;
#ifndef REF 
    magic_hack_bgcolor = white_color.hint;
#endif
    if (gc->clip != ((struct canvas *) gc->canvas)->lastclip
	    && cv_is_canvas(gc->canvas))
	cv_validate_clip(gc->canvas, gc->clip);
    if (cv_is_canvas(cv)) {
	if (((struct canvas *) cv)->lastclip != ((struct canvas *) cv)->outerclip)
	    cv_validate_clip(cv, ((struct canvas *) cv)->outerclip);
	cs_checkcursordims(((struct canvas *) cv),
			   0, 0,
			   cv->pr_width, cv->pr_height);
    }
    if (cv->pr_ops->pro_rop != mem_rop) {
	in = mem_create(cv->pr_width, cv->pr_height, cv->pr_depth);
	pr_rop(in, 0, 0, cv->pr_width, cv->pr_height, PIX_SRC, cv, 0, 0);
    }
    else
	in = cv;
    new = cs_scale_and_rotate(in, gc, &outline);
    if (new) {
	if (outline) {
	    cs_checkcursordims(gc->canvas, outline->pos.x, outline->pos.y,
			       outline->size.x, outline->size.y);
	    pr_shaperop(gc->canvas, outline, opcode, new, outline->pos.x, outline->pos.y);
	    sh_decref(outline);
	}
	else {
	    cs_checkcursordims(gc->canvas, cfloorfr(gc->transform.matrix[2][0]),
		       cfloorfr(gc->transform.matrix[2][1]) - new->pr_height,
			       new->pr_width, new->pr_height);
	    pr_rop(gc->canvas,
		   cfloorfr(gc->transform.matrix[2][0]),
		   cfloorfr(gc->transform.matrix[2][1]) - new->pr_height,
		   new->pr_width, new->pr_height, opcode, new, 0, 0);
	}
    }
    if (in && in != cv)
	pr_destroy(in);
    if (new && new != in && new != cv)
	pr_destroy(new);
    return new != 0;
}

struct pixrect *
cs_scale_and_rotate(in, gc, outlinep)
    struct pixrect *in;
    struct graphics_context *gc;
    struct shape **outlinep;
{
    struct pixrect *skew1 = 0;	/* the rotated image */
    *outlinep = 0;
    if (gc->transform.matrix[1][0] == 0
	    && gc->transform.matrix[0][0] > 0
	    && gc->transform.matrix[1][1] < 0
	    && gc->transform.matrix[0][1] == 0) {
	skew1 = cs_scalepixrect(in,
				cfloorfr(gc->transform.matrix[0][0]),
				-cfloorfr(gc->transform.matrix[1][1]),
				gc->canvas->pr_depth);
    }
    else {			/* rotation! */
	short       xv,
	            yv,		/* unit up (vertical) vector */
	            xh,
	            yh;		/* unit right (horizontal) vector */
	short       axv,
	            ayv,
	            axh,
	            ayh;	/* abs versions of prev */
	short       x0,
	            y0,		/* upper right corner of result image */
	            hsize,
	            vsize,	/* horizontal and vertical size of the image
				 * before rotation */
	            verror;	/* bresenham error used for stepping down the
				 * left edge */
	short       hxmajor,	/* true if the horizontal vector uses x-major
				 * bresenham */
	            vxmajor,	/* true if the vertical vector uses x-major
				 * bresenham */
	            iherror;	/* initial horizontal bresenham error */
	short       sline,	/* bytes per raster line in source */
	            dline,	/* bytes per raster line in dest */
	            vcnt;	/* counts down the vertical vector */
	int         area;	/* The area of the rotated image */
	struct pixrect *new;	/* The scaled but unrotated image */
	struct shape *outline;	/* the shape of the rotated image */
	cs_dtransform(gc,
		      0, 1,	/* unit up (vertical) vector */
		      &xv, &yv);
	cs_dtransform(gc,
		      1, 0,	/* unit across (horizontal) vector */
		      &xh, &yh);
	cs_gsave(gc);
	cs_newpath(gc);
	cs_frmoveto(gc, 0, 0);	/* construct the unit square */
	cs_frlineto(gc, 0, fracti(1));
	cs_frlineto(gc, fracti(1), fracti(1));
	cs_frlineto(gc, fracti(1), 0);
	cs_closepath(gc);
	outline = cv_pathtoshape(gc, ~1);
	cs_grestore(gc);
	if (outline == 0)
	    return 0;
	hxmajor = vxmajor = 1;
	axv = xv >= 0 ? xv : -xv;
	ayv = yv >= 0 ? yv : -yv;
	if (axv <= ayv)
	    vxmajor = 0;
	axh = xh >= 0 ? xh : -xh;
	ayh = yh >= 0 ? yh : -yh;
	if (axh > ayh)
	    hsize = axh;
	else
	    hsize = ayh, hxmajor = 0;
	area = xh * yv - xv * yh;
	if (area < 0)
	    area = -area;
	if (hsize==0 || area==0) {
	    sh_decref(outline);
	    return 0;
	}
	vsize = (area + (hsize >> 1)) / hsize;
	new = cs_scalepixrect(in, hsize, vsize, gc->canvas->pr_depth);
	if (new) {
	    skew1 = mem_create(outline->size.x, outline->size.y, new->pr_depth);
	    x0 = cfloorfr(gc->transform.matrix[2][0]) + xv - outline->pos.x;
	    rotated_upper_left.x = x0;
	    y0 = cfloorfr(gc->transform.matrix[2][1]) + yv - outline->pos.y;
	    if (yv > 0)		/* HACK! */
		y0--;
	    rotated_upper_left.y = y0;
	    verror = 0;
	    iherror = 0;
	    switch (new->pr_depth) {
	    case 8:
		{
		    unsigned char *sleft,	/* source left edge (initial
						 * value of sp) */
		               *dleft;	/* destination left edge (initial
					 * value of dp) */
		    register unsigned char
		               *sp,	/* source pointer */
		               *dp;	/* destination pointer */
		    assert(new->pr_depth == 8);
		    sline = mpr_d(new)->md_linebytes;
		    sleft = mprd8_addr(mpr_d(new), 0, 0, 8);
		    dline = mpr_d(skew1)->md_linebytes;
		    dleft = mprd8_addr(mpr_d(skew1), x0, y0, 8);
		    for (vcnt = vsize; --vcnt >= 0;) {
			register    hcnt;
			register    herror = iherror;
			sp = sleft;
			dp = dleft;

#ifdef undef
			for (hcnt = hsize; --hcnt >= 0;) {
			    *dp = *sp++;
			    if (hxmajor) {
				if (xh > 0)
				    dp++;
				else
				    dp--;
				herror += ayh;
				if (herror > axh) {
				    herror -= axh;
				    if (yh > 0)
					dp += dline;
				    else
					dp -= dline;
				}
			    }
			    else {	/* hymajor */
				if (yh > 0)
				    dp += dline;
				else
				    dp -= dline;
				herror += axh;
				if (herror > ayh) {
				    herror -= ayh;
				    if (xh > 0)
					dp++;
				    else
					dp--;
				}
			    }
			}
#endif

			if (hxmajor)
			    if (xh > 0)
				for (hcnt = hsize; --hcnt >= 0;) {
				    *dp = *sp++;
				    dp++;
				    herror += ayh;
				    if (herror > axh) {
					herror -= axh;
					if (yh > 0)
					    dp += dline;
					else
					    dp -= dline;
				    }
				}
			    else/* xh<=0 */
				for (hcnt = hsize; --hcnt >= 0;) {
				    *dp = *sp++;
				    dp--;
				    herror += ayh;
				    if (herror > axh) {
					herror -= axh;
					if (yh > 0)
					    dp += dline;
					else
					    dp -= dline;
				    }
				}
			else
			 /* !hxmajor */ if (yh > 0)
			    for (hcnt = hsize; --hcnt >= 0;) {
				*dp = *sp++;
				dp += dline;
				herror += axh;
				if (herror > ayh) {
				    herror -= ayh;
				    if (xh > 0)
					dp++;
				    else
					dp--;
				}
			    }
			else	/* yh<=0 */
			    for (hcnt = hsize; --hcnt >= 0;) {
				*dp = *sp++;
				dp -= dline;
				herror += axh;
				if (herror > ayh) {
				    herror -= ayh;
				    if (xh > 0)
					dp++;
				    else
					dp--;
				}
			    }
			sleft += sline;
			{
			    int         ihchange = 0;
			    if (vxmajor) {
				if (xv > 0)
				    dleft--;
				else
				    dleft++;
				verror += ayv;
				if (verror > axv) {
				    verror -= axv;
				    if (yv > 0) {
					ihchange = -1;
					dleft -= dline;
				    }
				    else {
					dleft += dline;
					ihchange = 1;
				    }
				}
			    }
			    else {	/* vymajor */
				if (yv > 0)
				    dleft -= dline;
				else
				    dleft += dline;
				verror += axv;
				if (verror > ayv) {
				    verror -= ayv;
				    if (xv > 0) {
					ihchange = -1;
					dleft--;
				    }
				    else {
					dleft++;
					ihchange = 1;
				    }
				}
			    }
			    if (ihchange) {
				if (((vxmajor ? yv : xv) ^ (hxmajor ? xh : yh)) > 0
					 /* ihchange < 0 */ ) {
				    if (hxmajor) {
					iherror -= ayh;
					if (iherror >= 0)
					    continue;
					iherror += axh;
				    }
				    else {
					iherror -= axh;
					if (iherror >= 0)
					    continue;
					iherror += ayh;
				    }
				}
				else {	/* if (ihchange > 0) */
				    if (hxmajor) {
					iherror += ayh;
					if (iherror < axh)
					    continue;
					iherror -= axh;
				    }
				    else {
					iherror += axh;
					if (iherror < ayh)
					    continue;
					iherror -= ayh;
				    }
				}
				if (vxmajor) {
				    verror -= ayv;
				    if (xv <= 0)
					dleft--;
				    else
					dleft++;
				}
				else {	/* vymajor */
				    verror -= axv;
				    if (yv <= 0)
					dleft -= dline;
				    else
					dleft += dline;
				}
			    }
			}
		    }
		}
		break;
	    case 1:
		{
		    unsigned short *sleft;	/* source left edge (initial
						 * value of sp) */
		    register unsigned short
		               *sp,	/* source pointer */
		               *dp;	/* destination pointer */

#ifdef DEBUG
		    unsigned short *slim = (unsigned short *)
		    ((int) mpr_d(new)->md_image
		     + mpr_d(new)->md_linebytes * new->pr_height);
#endif
		    unsigned short *dbot = (unsigned short *) mpr_d(skew1)->md_image;
		    unsigned short *dlim = (unsigned short *)
		    ((int) mpr_d(skew1)->md_image
		     + mpr_d(skew1)->md_linebytes * skew1->pr_height);
		    sline = mpr_d(new)->md_linebytes / sizeof *sp;
		    sleft = (unsigned short *) mprd_addr(mpr_d(new), 0, 0);
		    dline = mpr_d(skew1)->md_linebytes / sizeof *dp;
		    for (vcnt = vsize; --vcnt >= 0;) {
			register    hcnt;
			register    herror = iherror;
			register unsigned short dmask,
			            smask;
			smask = PIXBIT(mprd_skew(mpr_d(new), 0, 0));
			dmask = PIXBIT(mprd_skew(mpr_d(skew1), x0, y0));
			sp = sleft;
			dp = (unsigned short *) mprd_addr(mpr_d(skew1), x0, y0);
			for (hcnt = hsize; --hcnt >= 0;) {
			    assert(sp >= (unsigned short *) mpr_d(new)->md_image && sp < slim);
/*-			    assert(dp>=(unsigned short*)mpr_d(skew1)->md_image
				   && dp<dlim); */
			    /*
			     * The range test on dp in the following line is a
			     * terrible hack.  There is a bug in the
			     * computation of the alignment between the
			     * character and the destination pixrect that
			     * occasionally causes off-by-one errors.  The
			     * off-by-one errors lead to touching pixels
			     * outside of the pixrect. It's a numerical
			     * stability issue.  I don't have the time to fix
			     * it now, so I'll just leave this test in.
			     */
			    if (*sp & smask && dbot <= dp && dp < dlim)
				*dp |= dmask;
			    RIGHTSTEP(sp, smask);
			    if (hxmajor) {
				if (xh > 0)
				    RIGHTSTEP(dp, dmask);
				else
				    LEFTSTEP(dp, dmask);
				herror += ayh;
				if (herror > axh) {
				    herror -= axh;
				    if (yh > 0)
					dp += dline;
				    else
					dp -= dline;
				}
			    }
			    else {	/* hymajor */
				if (yh > 0)
				    dp += dline;
				else
				    dp -= dline;
				herror += axh;
				if (herror > ayh) {
				    herror -= ayh;
				    if (xh > 0)
					RIGHTSTEP(dp, dmask);
				    else
					LEFTSTEP(dp, dmask);
				}
			    }
			}
			sleft += sline;
			{
			    int         ihchange = 0;
			    if (vxmajor) {
				if (xv > 0)
				    x0--;
				else
				    x0++;
				verror += ayv;
				if (verror > axv) {
				    verror -= axv;
				    if (yv > 0) {
					ihchange = -1;
					y0--;
				    }
				    else {
					y0++;
					ihchange = 1;
				    }
				}
			    }
			    else {	/* vymajor */
				if (yv > 0)
				    y0--;
				else
				    y0++;
				verror += axv;
				if (verror > ayv) {
				    verror -= ayv;
				    if (xv > 0) {
					ihchange = -1;
					x0--;
				    }
				    else {
					x0++;
					ihchange = 1;
				    }
				}
			    }
			    if (ihchange) {
				if (((vxmajor ? yv : xv) ^ (hxmajor ? xh : yh)) > 0
					 /* ihchange < 0 */ ) {
				    if (hxmajor) {
					iherror -= ayh;
					if (iherror >= 0)
					    continue;
					iherror += axh;
				    }
				    else {
					iherror -= axh;
					if (iherror >= 0)
					    continue;
					iherror += ayh;
				    }
				}
				else {	/* if (ihchange > 0) */
				    if (hxmajor) {
					iherror += ayh;
					if (iherror < axh)
					    continue;
					iherror -= axh;
				    }
				    else {
					iherror += axh;
					if (iherror < ayh)
					    continue;
					iherror -= ayh;
				    }
				}
				if (vxmajor) {
				    verror -= ayv;
				    if (xv <= 0)
					x0--;
				    else
					x0++;
				}
				else {	/* vymajor */
				    verror -= axv;
				    if (yv <= 0)
					y0--;
				    else
					y0++;
				}
			    }
			}
		    }
		}
		break;
	    }
	}
	*outlinep = outline;
	if (new && new != in)
	    pr_destroy(new);
    }
    return skew1;
}

struct pixrect *
cs_scalepixrect(pr, dwidth, dheight, ddepth)
    register struct pixrect *pr;
{
    struct pixrect *prn = 0;
    struct pixrect *temp = 0;
    if (dwidth <= 0)
	dwidth = pr->pr_width;
    if (dheight <= 0)
	dheight = pr->pr_height;
    if (ddepth <= 0)
	ddepth = pr->pr_depth;
    if (pr->pr_depth == 4) {
	temp = mem_create(dwidth, dheight, 8);
	cv_quickscale4to8(pr, temp);
	pr = temp;
    }
    else if (dwidth != pr->pr_width || dheight != pr->pr_height)
	if (ddepth <= 8 && pr->pr_depth == 24) {
	    temp = mem_create(dwidth, dheight, ddepth);
	    cv_slowscale24to8(pr, temp, &colormap);
	    pr = temp;
	}
	else {
	    temp = mem_create(dwidth, dheight, pr->pr_depth);
	    if (temp == 0)
		return 0;
	    switch (pr->pr_depth) {
	    case 24:
		cv_slowscale24(pr, temp);
		break;
	    case 8:
		cv_quickscale8(pr, temp);
		break;
	    case 1:
		cv_quickscale1(pr, temp);
		break;
	    default:
		pr_destroy(temp);
		return 0;
	    }
	    pr = temp;
	}
    if (pr->pr_depth != ddepth) {
	switch (pr->pr_depth * 100 + ddepth) {
	case 2408:
	    prn = cv_dither24to8(pr, &colormap);
	    break;
	case 801:
	    prn = cv_dither8to1(pr, &colormap);
	    break;
	case 2401:
	    prn = cv_dither24to1(pr, &colormap);
	    break;
	case 108:
	    prn = pr;
	    break;
	default:
	    prn = 0;
	}
    }
    else
	prn = pr;
    if (temp && prn != temp)
	pr_destroy(temp);
    return prn;
}


struct pixrect *
cv_dither24to8(ipr, map)
    struct pixrect *ipr;
    colormap_t *map;
{
    struct pixrect *pr;
    register    x,
                y;
    struct error {
	short       red,
	            green,
	            blue;
    };
    struct error *errormap;
    register struct pixel *pixel;
    register struct error *error;
    register unsigned char *dst;
    errormap = (struct error *) snoopcalloc("image.c", sizeof(struct error), ipr->pr_width + 2);
    pr = mem_create(ipr->pr_width, ipr->pr_height, 8);
    dst = (unsigned char *) mpr_d(pr)->md_image;
    for (y = 0; y < ipr->pr_height; y++) {
	error = errormap + 1;
	pixel = (struct pixel *) mprd_addr(mpr_d(ipr), 0, y);
	for (x = 0; x < ipr->pr_width; x++) {
	    register    red = pixel->red + error->red;
	    register    green = pixel->green + error->green;
	    register    blue = pixel->blue + error->blue;
	    register    index;
	    register    T;

	    if (red > 255)
		error->red = red - 255, red = 255;
	    else if (red < 0)
		error->red = red, red = 0;
	    else
		error->red = 0;

	    if (green > 255)
		error->green = green - 255, green = 255;
	    else if (green < 0)
		error->green = green, green = 0;
	    else
		error->green = 0;

	    if (blue > 255)
		error->blue = blue - 255, blue = 255;
	    else if (blue < 0)
		error->blue = blue, blue = 0;
	    else
		error->blue = 0;
	    index = red_inverse[red] + green_inverse[green] + blue_inverse[blue];
	    error->red -= red_map[index] - red;
	    error->green -= green_map[index] - green;
	    error->blue -= blue_map[index] - blue;
	    *dst++ = index;
	    T = (error->red * 3 + 4) >> 3;
	    error[1].red += T;
	    error[0].red -= T + T;
	    error[-1].red += T;
	    T = (error->green * 3 + 4) >> 3;
	    error[1].green += T;
	    error[0].green -= T + T;
	    error[-1].green += T;
	    T = (error->blue * 3 + 4) >> 3;
	    error[1].blue += T;
	    error[0].blue -= T + T;
	    error[-1].blue += T;
	    error++;
	    /*
 	     * add 3 to value of pixel
	     */
	    pixel = (struct pixel *) ((char *)pixel + 3);
	}
	if (ipr->pr_width & 1)
	    *dst++ = 0;
    }
    free(errormap);
    return pr;
}

struct pixrect *
cv_dither24to1(ipr, map)
    struct pixrect *ipr;
    colormap_t *map;
{
    struct pixrect *pr;
    register    x,
                y;
    short      *errormap;
    register struct pixel *pixel;
    register short *error;
    register unsigned char *dst;
    errormap = (short *) snoopcalloc("image.c", sizeof(short), ipr->pr_width + 2);
    pr = mem_create(ipr->pr_width, ipr->pr_height, 1);
    for (y = 0; y < ipr->pr_height; y++) {
	error = errormap + 1;
	pixel = (struct pixel *) mprd_addr(mpr_d(ipr), 0, y);
	dst = (unsigned char *) mprd_addr(mpr_d(pr), 0, y);
	for (x = 0; x < ipr->pr_width; x++) {
	    register    Tn;
	    register    T;
	    T = pixel->red + pixel->green + pixel->blue;
	    T += *error;
	    if (T > 0x80 * 3)
		T -= 0x80 * 3;
	    else
		dst[x >> 3] |= 0x80 >> (x & 7);
	    Tn = (T * 3 + 4) >> 3;
	    error[1] += Tn;
	    error[-1] += Tn;
	    *error++ = T - 2 * Tn;
	    /*
	     * add 3 to the value of pixel
	     */
	    pixel = (struct pixel *) ((char *) pixel + 3);
	}
    }
    free(errormap);
    return pr;
}

struct pixrect *
cv_dither8to1(ipr, map)
    struct pixrect *ipr;
    colormap_t *map;
{
    struct pixrect *pr;
    register    x,
                y;
    short      *errormap;
    register unsigned char *pixel;
    register short *error;
    register unsigned short *dst;
    unsigned char intensity_map[256];
    errormap = (short *) snoopcalloc("image.c", sizeof(short), ipr->pr_width + 2);
    for (y = map->length; --y >= 0;) {

#ifdef undef
	register    v = map->map[0][y];
	x = map->map[1][y];
	if (x > v)
	    v = x;
	x = map->map[2][y];
	if (x > v)
	    v = x;
	intensity_map[y] = v;
#endif

	intensity_map[y] = (map->map[0][y] + map->map[1][y] + map->map[2][y]) / 3;
    }
    pr = mem_create(ipr->pr_width, ipr->pr_height, 1);
    for (y = 0; y < ipr->pr_height; y++) {
	register    mask;
	error = errormap + 1;
	pixel = (unsigned char *) mprd_addr(mpr_d(ipr), 0, y);
	dst = (unsigned short *) mprd_addr(mpr_d(pr), 0, y);
	mask = PIXBIT(mprd_skew(mpr_d(pr), 0, y));
	for (x = 0; x < ipr->pr_width; x++) {
	    register    Tn;
	    register    T;
	    T = intensity_map[*pixel++] + *error;
	    if (T > 0x80)
		T -= 0xFF;
	    else
		*dst |= mask;
	    RIGHTSTEP(dst, mask);
	    Tn = (T * 3 + 4) >> 3;
	    error[1] += Tn;
	    error[-1] += Tn;
	    *error++ = T - 2 * Tn;
	}
    }
    free(errormap);
    return pr;
}

static unsigned char dither_matrix[16][16];

static
make_dither_matrix(n)
    register    n;
{
    register    i,
                j;
    register    halfn;
    if (n <= 0)
	return;
    halfn = n >> 1;
    make_dither_matrix(halfn);
    for (i = 0; i < halfn; i++)
	for (j = 0; j < halfn; j++) {
	    register T;
	    register unsigned char *dmp = &dither_matrix[i][j];
	    T = *dmp << 2;
	    *dmp = T;
	    dither_matrix[i + halfn][j + halfn] = T + 1;
	    dither_matrix[i + halfn][j] = T + 2;
	    dither_matrix[i][j + halfn] = T + 3;
	}
}

/* matrix dither */
struct pixrect *
cv_mdither8to1(ipr, map)
    struct pixrect *ipr;
    colormap_t *map;
{
    struct pixrect *pr;
    register    x,
                y;
    unsigned char intensity_map[256];
    for (y = map->length; --y >= 0;) {
	register    v = map->map[0][y];
	x = map->map[1][y];
	if (x > v)
	    v = x;
	x = map->map[2][y];
	if (x > v)
	    v = x;
	intensity_map[y] = v;
    }
    pr = mem_create(ipr->pr_width, ipr->pr_height, 1);
    if (dither_matrix[1][1] == 0)
	make_dither_matrix(16);
    for (y = 0; y < ipr->pr_height; y++) {
	register unsigned char *src,
	           *dst;
	register unsigned char *dither = dither_matrix[y & 0xF];
	src = (unsigned char *) mprd_addr(mpr_d(ipr), 0, y);
	dst = (unsigned char *) mprd_addr(mpr_d(pr), 0, y);
	for (x = 0; x < ipr->pr_width; x++)
	    if (intensity_map[*src++] <= dither[x & 0xF])
		dst[x >> 3] |= 0x80 >> (x & 7);
    }
    return pr;
}


/* Scaling of an 8 bit image by pixel replication */
cv_quickscale8(ipr, opr)
    register struct pixrect *ipr,
               *opr;
{
    register unsigned char *src,
               *dst;
    int         y;
    if (opr->pr_depth != 8 || opr->pr_depth != 8)
	return -1;
    if (opr->pr_width == ipr->pr_width
	    && opr->pr_height == ipr->pr_height) {
	mem_rop(opr, 0, 0, opr->pr_width, opr->pr_height,
		PIX_SRC, ipr, 0, 0);
	return 0;
    }
    for (y = opr->pr_height; --y >= 0;) {
	register    err = 0,
	            limit = opr->pr_width,
	            step = ipr->pr_width,
	            cnt = opr->pr_width;
	dst = (unsigned char *) mprd_addr(mpr_d(opr), 0, y);
	src = (unsigned char *) mprd_addr(mpr_d(ipr), 0,
					y * ipr->pr_height / opr->pr_height);
	while (--cnt >= 0) {
	    *dst++ = *src;
	    err += step;
	    while (err >= limit) {
		err -= limit;
		src++;
	    }
	}
    }
    return 0;
}


/* Scaling of a 4 bit image into an 8 bit image by pixel replication */
cv_quickscale4to8(ipr, opr)
    register struct pixrect *ipr,
               *opr;
{
    register unsigned char *src,
               *dst;
    int         y;
    static char grey4_inverse[1 << 4];
    if (ipr->pr_depth != 4 || opr->pr_depth != 8)
	return -1;
    if (grey4_inverse[(1 << 4) - 1] == 0) {
	register unsigned char i;
	for (i = 0; i < 1 << 4; i++)
	    grey4_inverse[i] = grey_inverse[i * 255 / ((1 << 4) - 1)];
    }
    for (y = opr->pr_height; --y >= 0;) {
	register    err = 0,
	            limit = opr->pr_width,
	            cnt = opr->pr_width;
	register    nibble = 0;
	register    step = ipr->pr_width;
	dst = (unsigned char *) mprd_addr(mpr_d(opr), 0, y);
	src = (unsigned char *) mprd_addr(mpr_d(ipr), 0,
					y * ipr->pr_height / opr->pr_height);
	while (--cnt >= 0) {
	    *dst++ = grey4_inverse[nibble & 1 ? *src & 0xF : *src >> 4];
	    err += step;
	    while (err >= limit) {
		err -= limit;
		nibble++;
		if ((nibble & 1) == 0)
		    src++;
	    }
	}
    }
    return 0;
}

/* Scaling of a 1 bit image by pixel replication */
cv_quickscale1(ipr, opr)
    register struct pixrect *ipr,
               *opr;
{
    int         y;
    if (ipr->pr_depth != 1 || opr->pr_depth != 1)
	return -1;
    if (opr->pr_width == ipr->pr_width
	    && opr->pr_height == ipr->pr_height) {
	mem_rop(opr, 0, 0, opr->pr_width, opr->pr_height,
		PIX_SRC, ipr, 0, 0);
	return 0;
    }
    for (y = opr->pr_height; --y >= 0;) {
	register unsigned short *src,
	           *dst;
#ifdef LITTLEENDIAN
	register unsigned short dwrd,		/* These MUST be unsigned */
	            		swrd;		/* for the little endian path */
#else
	register    short	dwrd,		/* and signed for the big */
				swrd;		/* endian path */
#endif
	register    err = 0,
	            limit = opr->pr_width,
	            step = ipr->pr_width,
	            cnt = (opr->pr_width + 15) >> 4,
	            srcleft = 16;
	dst = (unsigned short *) mprd_addr(mpr_d(opr), 0, y);
	src = (unsigned short *) mprd_addr(mpr_d(ipr), 0,
					y * ipr->pr_height / opr->pr_height);
	swrd = *src++;
#ifdef LITTLEENDIAN
	/* On a little endian machine we set bits in the short starting
	 * at the left and moving right, always setting the leftmost 
	 * bit; on a big endian machine we do the reverse and set bits
	 * starting at the right and moving left, always setting the
	 * rightmost bit.
	 */
	while (--cnt >= 0) {
	    int	n;
	    dwrd = 0;

	    for (n = 16; --n >= 0;) {
		dwrd >>= 1;
		if ((swrd&1) != 0)
		    dwrd |= (1<<15);

		err += step;
		while (err >= limit) {
		    err -= limit;
		    swrd >>= 1;
		    if (--srcleft <= 0) {
			swrd = *src++;
			srcleft = 16;
		    }
		}
	    }
	    *dst++ = dwrd;
	}
#else
	while (--cnt >= 0) {
	    int         n;
	    for (n = 16; --n >= 0;) {
		dwrd <<= 1;
		if (swrd < 0)
		    dwrd |= 1;
		err += step;
		while (err >= limit) {
		    err -= limit;
		    swrd <<= 1;
		    if (--srcleft <= 0) {
			swrd = *src++;
			srcleft = 16;
		    }
		}
	    }
	    *dst++ = dwrd;
	}
#endif
    }
    return 0;
}

/*
 * scaling of a 24 bit image by linear interpolation of a pixels four
 * neighbours
 */
cv_slowscale24(ipr, opr)
    register struct pixrect *ipr,
               *opr;
{
    int         dy,
                dx;
    int         iy;
    int         y;
    register    lineinc;
    if (opr->pr_depth != 24 || opr->pr_depth != 24)
	return -1;
    if (opr->pr_width == ipr->pr_width
	    && opr->pr_height == ipr->pr_height) {
	mem_rop(opr, 0, 0, opr->pr_width, opr->pr_height,
		PIX_SRC, ipr, 0, 0);
	return 0;
    }
    dy = (ipr->pr_height << 16) / opr->pr_height;
    dx = (ipr->pr_width << 16) / opr->pr_width;
    iy = 0;
    lineinc = mpr_d(ipr)->md_linebytes;
    for (y = 0; y < opr->pr_height; y++) {
	register struct pixel *dst;
	struct pixel *src0;
	register    ix = 0;
	int         cnt = opr->pr_width;
	dst = (struct pixel *) mprd_addr(mpr_d(opr), 0, y);
	src0 = (struct pixel *) mprd_addr(mpr_d(ipr), 0, iy >> 16);
	while (--cnt >= 0) {
	    register struct pixel *src;
	    register unsigned short t1,
		        t2,
		        t3;
	    src = (struct pixel *) ((char *) src0 + (ix >> 16) * 3);
#define m(a,b) ((short) (((short)(a))*((short)(b))))
#define s(o) (*(struct pixel *)(((char *)src) + (o)*3))
#define sn(o) (*(struct pixel *)(((char *)src) + (o)*3+lineinc))
	    t1 = (ix >> 8) & 0xFF;
	    t2 = m(t1, s(1).red);
	    t3 = m(t1, sn(1).red);
	    t1 = 0x100 - t1;
	    t2 += m(t1, s(0).red);
	    t3 += m(t1, sn(0).red);
	    t1 = (iy >> 8) & 0xFF;
	    dst->red = (m(0x100 - t1, t2 >> 8) + m(t1, t3 >> 8)) >> 8;

	    t1 = (ix >> 8) & 0xFF;
	    t2 = m(t1, s(1).green);
	    t3 = m(t1, sn(1).green);
	    t1 = 0x100 - t1;
	    t2 += m(t1, s(0).green);
	    t3 += m(t1, sn(0).green);
	    t1 = (iy >> 8) & 0xFF;
	    dst->green = (m(0x100 - t1, t2 >> 8) + m(t1, t3 >> 8)) >> 8;

	    t1 = (ix >> 8) & 0xFF;
	    t2 = m(t1, s(1).blue);
	    t3 = m(t1, sn(1).blue);
	    t1 = 0x100 - t1;
	    t2 += m(t1, s(0).blue);
	    t3 += m(t1, sn(0).blue);
	    t1 = (iy >> 8) & 0xFF;
	    dst->blue = (m(0x100 - t1, t2 >> 8) + m(t1, t3 >> 8)) >> 8;
	    ix += dx;
	    /*
	     * add 3 to the value of dst
	     */
	    dst = (struct pixel *) ((unsigned char *) dst + 3);
	}
	iy += dy;
    }
    return 0;
}

/*
 * scaling of a 24 bit image into an 8 bit image by linear interpolation of a
 * pixels four neighbours
 */
cv_slowscale24to8(ipr, opr, map)
    struct pixrect *ipr,
               *opr;
    colormap_t *map;
{
    int         dy,
                dx;
    int         iy;
    int         y;
    register    lineinc;
    struct error {
	short       red,
	            green,
	            blue;
    };
    struct error *errormap;
    if (ipr->pr_depth != 24 || opr->pr_depth != 8)
	return -1;
    errormap = (struct error *) snoopcalloc("image.c", sizeof(struct error), opr->pr_width + 2);
    dy = (ipr->pr_height << 16) / opr->pr_height;
    dx = (ipr->pr_width << 16) / opr->pr_width;
    iy = 0;
    lineinc = mpr_d(ipr)->md_linebytes;
    for (y = 0; y < opr->pr_height; y++) {
	register char *dst;
	struct pixel *src0;
	register struct error *error = errormap + 1;
	register    ix = 0;
	int         cnt = opr->pr_width;
	dst = (char *) mprd_addr(mpr_d(opr), 0, y);
	src0 = (struct pixel *) mprd_addr(mpr_d(ipr), 0, iy >> 16);
	while (--cnt >= 0) {
	    register struct pixel *src;
	    register /* unsigned */ short t1,
	                t2,
	                t3;
	    register short red,
	                green;
#define BLUE t3
	    src = (struct pixel *) (((int) src0) + (ix >> 16) * 3);
#undef m
#define m(a,b) ((unsigned short) (((unsigned short)(a))*((unsigned short)(b))))
#define ush(n) (((unsigned short)(n))>>8)
#define s(o) (*(struct pixel *)(((char *)src) + (o)*3))
#define sn(o) (*(struct pixel *)(((char *)src) + (o)*3+lineinc))
	    t1 = (ix >> 8) & 0xFF;
	    t2 = m(t1, s(1).red);
	    t3 = m(t1, sn(1).red);
	    t1 = 0x100 - t1;
	    t2 += m(t1, s(0).red);
	    t3 += m(t1, sn(0).red);
	    t1 = (iy >> 8) & 0xFF;
	    red = (m(0x100 - t1, ush(t2)) + m(t1, ush(t3))) >> 8;

	    t1 = (ix >> 8) & 0xFF;
	    t2 = m(t1, s(1).green);
	    t3 = m(t1, sn(1).green);
	    t1 = 0x100 - t1;
	    t2 += m(t1, s(0).green);
	    t3 += m(t1, sn(0).green);
	    t1 = (iy >> 8) & 0xFF;
	    green = (m(0x100 - t1, ush(t2)) + m(t1, ush(t3))) >> 8;

	    t1 = (ix >> 8) & 0xFF;
	    t2 = m(t1, s(1).blue);
	    t3 = m(t1, sn(1).blue);
	    t1 = 0x100 - t1;
	    t2 += m(t1, s(0).blue);
	    t3 += m(t1, sn(0).blue);
	    t1 = (iy >> 8) & 0xFF;
	    BLUE = (m(0x100 - t1, ush(t2)) + m(t1, ush(t3))) >> 8;
	    ix += dx;
	    red += error->red;
	    green += error->green;
	    BLUE += error->blue;
	    if (red > 255)
		error->red = red - 255, red = 255;
	    else if (red < 0)
		error->red = red, red = 0;
	    else
		error->red = 0;

	    if (green > 255)
		error->green = green - 255, green = 255;
	    else if (green < 0)
		error->green = green, green = 0;
	    else
		error->green = 0;

	    if (BLUE > 255)
		error->blue = BLUE - 255, BLUE = 255;
	    else if (BLUE < 0)
		error->blue = BLUE, BLUE = 0;
	    else
		error->blue = 0;
	    t1 = red_inverse[red] + green_inverse[green] + blue_inverse[BLUE];
	    error->red -= red_map[t1] - red;
	    error->green -= green_map[t1] - green;
	    error->blue -= blue_map[t1] - BLUE;
	    *dst++ = t1;
	    t2 = (error->red * 3 + 4) >> 3;
	    error[1].red += t2;
	    error[0].red -= t2 + t2;
	    error[-1].red += t2;
	    t2 = (error->green * 3 + 4) >> 3;
	    error[1].green += t2;
	    error[0].green -= t2 + t2;
	    error[-1].green += t2;
	    t2 = (error->blue * 3 + 4) >> 3;
	    error[1].blue += t2;
	    error[0].blue -= t2 + t2;
	    error[-1].blue += t2;
	    error++;
#undef BLUE
	}
	iy += dy;
    }
    free(errormap);
    return 0;
}

int
cs_writeimage(gc, f, mask, canvasonly)
    register struct graphics_context *gc;
    PSFILE *f;
{
    register struct pixrect *pr;
    register struct shape *outline = cv_pathtoshape(gc, mask);
    register struct canvas *cv = (struct canvas *) gc->canvas;
    struct spoint origin;
    if (f == 0)
	return 0;
    if (outline == 0)
	if (cv_is_canvas(cv)) {
	    outline = cv->outerclip;
	    sh_incref(outline);
	}
	else {
	    qtalloc(outline, (struct shape *));
	    bzero(outline, sizeof *outline);
	    outline->size.x = cv->pixrect.pr_size.x;
	    outline->size.y = cv->pixrect.pr_size.y;
	    outline->refcnt = 1;
	    outline->is_rect = 1;
	}
    if (!cv->validpix)
	cv_validate_pixrects(cv);
    pr = mem_create(outline->size.x, outline->size.y, cv->pixrect.pr_depth);
    origin = outline->pos;
    outline = sh_translate(outline, -outline->pos.x, -outline->pos.y);
    if (canvasonly) {
	cs_cursordown();
    } else {
	cs_cursorup();
    }
    pr_shaperop(pr, outline, PIX_SRC,
		canvasonly || !cv_is_canvas(cv) ? (struct pixrect *) cv
						: cv->screen_part.pixrect,
		-origin.x, -origin.y);
    cs_cursorup();
    sh_decref(outline);
#ifdef LITTLEENDIAN
    pr_dump(pr, f, &colormap, RT_STANDARD, 0);
#else
    pr_dump(pr, f, &colormap, RT_BYTE_ENCODED, 0);
#endif
    pr_destroy(pr);
}

struct pixrect *
cs_readimage(f)
    PSFILE     *f;
{
    colormap_t  map;
    register struct pixrect *pr;
    extern struct pixrect *pr_load();
    register    i;
    register unsigned char *p;
    unsigned char r[256],
                g[256],
                b[256],
                tr[256];
    struct rasterfile rh;
    int         stdmap;
    extern struct pixrect *pr_load_image();
    if (f == 0)
	return 0;
    map.map[0] = r;
    map.map[1] = g;
    map.map[2] = b;
    map.length = 256;
    map.type = 1;
    if (pr_load_header(f, &rh) == PIX_ERR) 
	return 0;
    if (rh.ras_maplength > 0 && pr_load_colormap(f, &rh, &map) == PIX_ERR) 
	return 0;
    pr = pr_load_image(f, &rh, &map);
    if (pr == 0 || pr->pr_depth != 8 || rh.ras_maplength == 0)
	return pr;
    stdmap = 1;
    for (i = 0; i < map.length; i++)
	if (r[i] != red_map[i]
		|| g[i] != green_map[i]
		|| b[i] != blue_map[i]) {
	    stdmap = 0;
	    if (r[i] == g[i] && g[i] == b[i])
		tr[i] = grey_inverse[r[i]];
	    else
		tr[i] = red_inverse[r[i]]
		    + green_inverse[g[i]]
		    + blue_inverse[b[i]];
	}
	else
	    tr[i] = i;
    if (!stdmap) {
	p = (unsigned char *) mpr_d(pr)->md_image;
	i = pr->pr_height * ((pr->pr_width + 1) & ~1);
	while (--i >= 0) {
	    *p = tr[*p];
	    p++;
	}
    }
    return pr;
}

static
isqrt(x)			/* sqrt(x) */
    register    x;
{
    register    root,
                temp;
    if (x <= 0)
	return 0;
    root = temp = x;		/* Compute initial approximation */
    while (temp > 1)
	root >>= 1, temp >>= 2;
    /* Then three Newton/Raphson iterations */
    root = (root + x / root) >> 1;
    root = (root + x / root) >> 1;
    root = (root + x / root) >> 1;
    return root;
}