#include "i.h"


// big tables in separate files
#include "rgb.h"
#include "ycbcr.h"

enum {
	CLAMPBOFF = 300,
	NCLAMPB = CLAMPBOFF+256+CLAMPBOFF,
	CLAMPNOFF = 64,
	NCLAMPN = CLAMPNOFF+256+CLAMPNOFF
};

int	dbgimg = 0;

static int		progressive = 0;

static uchar	clampb[NCLAMPB];
static int		clampn_b[NCLAMPN];
static int		clampn_g[NCLAMPN];
static int		clampn_r[NCLAMPN];
static int		cmap_r[256];
static int		cmap_g[256];
static int		cmap_b[256];

static void	imgerror(ImageSource* is, int err);
static int	getc(ImageSource* is);
static void	ungetc(ImageSource* is);
static void	ungetc2(ImageSource* is);
static int	getlew(ImageSource* is);
static int	getbew(ImageSource* is);
static void	readis(ImageSource* is, uchar* buf, int n);
static int	getn(ImageSource* is, int n);
static Image*	newimage(ImageSource* is, int w, int h);
static MaskedImage*	newmi(Image* im);
static void	setdims(ImageSource* is);
static void	tracearray(int* a, int n, char* name);
static void	getxbitmaphdr(ImageSource* is);
static MaskedImage*	getxbitmapmim(ImageSource* is);
static int	getxbitmapdefine(ImageSource* is, int* pval);
static void	get_to_char(ImageSource* is, int cterm);
static int	get_hexbyte(ImageSource* is);
static int	hexdig(int c);
static void	getgifhdr(ImageSource* is);
static int	gifgettoimage(ImageSource* is);
static MaskedImage*	getgifmim(ImageSource* is);
static uchar*	gifreadcmap(ImageSource* is, int bpe, int* pmaplen);
static void	getjpeghdr(ImageSource* is);
static int	jpegmarker(ImageSource* is);
static void	jpegtabmisc(ImageSource* is, int* pm, int* pn);
static void	jpeghuffmantables(ImageSource* is, int n);
static int	jpeghuffmantable(ImageSource* is);
static void	jpegquanttables(ImageSource* is, int n);
static int	jpegquanttable(ImageSource* is);
static MaskedImage*	getjpegmim(ImageSource* is);
static uchar**	jpegbaselinescan(ImageSource* is, int* pnchans);
static void	jrestart(ImageSource* is, int mcu);
static void	jpegprogressivescan(ImageSource* is);
static uchar**	jprogressiveIDCT(ImageSource* is, int* pnchans);
static void	jprogressiveinit(ImageSource* is, Jpegstate* h);
static void	jprogressivedc(ImageSource* is, int comp);
static void	jprogressiveac(ImageSource* is, int comp);
static void	jprogressiveacinc(ImageSource* is, int comp);
static void	jincrement(ImageSource* is, int* acc, int k, int Pt);
static void	colormap1(Jpegstate* h, uchar* pic, int* data, int mcu, int nacross);
static void	colormapall1(Jpegstate* h, uchar** chans, int* data0, int* data1, int* data2, int mcu, int nacross);
static void	colormap(Jpegstate* h, uchar** chans, int** data0, int** data1, int** data2, int mcu, int nacross, int Hmax, int Vmax, int* H, int* V);
static int	jdecode(ImageSource* is, Huffman* t);
static int	jnextbyte(ImageSource* is);
static int	jnextborm(ImageSource* is);
static int	jreceive(ImageSource* is, int s);
static int	jreceiveEOB(ImageSource* is, int s);
static int	jreceivebit(ImageSource* is);
static void	nibbles(int c, int* p1, int* p2);
static void	idct(int* b);
static void	init_tabs(void);
static void	remap1(uchar* pic, int dx, int dy, uchar* cmap, int cmaplen);
static void	remapgrey(uchar* pic, int dx, int dy);
static void	remaprgb(uchar** chans, int dx, int dy);
static uchar*	resample(uchar* src, int sw, int sh, int dw, int dh);

void
imginit(void)
{
	init_tabs();
}

// Return true if mtype is an image type we can handle
int
supported(int mtype)
{
	switch(mtype) {
	case ImageJpeg:
	case ImageGif:
	case ImageXXBitmap:
		return 1;
	}
	return 0;
}

// w,h passed in are specified width and height.
// Result will be resampled if they don't match the dimensions
// in the decoded picture (if only one of w,h is specified, the other
// dimension is scaled by the same factor).
ImageSource*
newimagesource(ByteSource* bs, int w, int h)
{
	int	mtype;
	ImageSource*	is;

	dbgimg = config.dbg['i'];
	mtype = UnknownType;
	if(bs->hdr != nil)
		mtype = bs->hdr->mtype;
	is = (ImageSource*)emallocz(sizeof(ImageSource));
	is->width = w;
	is->height = h;
	is->mtype = mtype;
	is->bs = bs;
	return is;
}

int
getmim(ImageSource* is, MaskedImage** pans)
{
	MaskedImage*	ans;
	int	ret;
	int	v;
	int	ed;

	if(dbgimg)
		trace("img: getmim %U %dx%d\n", is->bs->req->url, is->width, is->height);
	if(dbgev)
		logtime("IMAGE_GETMIM", is->width*is->height);
	ans = nil;
	*pans = nil;
	ret = Mimnone;
	ed = is->bs->edata;

	// To avoid the need to check for error returns on a lot of low level functions,
	// we setjmp here, and call imgerror to return control to this place if
	// there is an error in processing the image.
	if((v = setjmp(is->getmimsave))) {
		is->err = v;
		ret = Mimerror;
		if(dbgimg)
			trace("getmim got err: %S\n", errphrase(is->err));
	}
	else {
		if(is->bs->hdr == nil)
			return Mimnone;
		else if(is->mtype == UnknownType)
			is->mtype = is->bs->hdr->mtype;
		// temporary hack: wait until whole file is here first
		if(is->bs->edata == is->bs->dalloclen) {
			switch(is->mtype) {
			case ImageJpeg:
				ans = getjpegmim(is);
				break;
			case ImageGif:
				ans = getgifmim(is);
				break;
			case ImageXXBitmap:
				ans = getxbitmapmim(is);
				break;
			default:
				is->err = ERRunsupimg;
				ret = Mimerror;
				ans = nil;
				break;
			}
			if(ans != nil)
				ret = Mimdone;
		}
	}
	if(dbgimg)
		trace("img: getmim returns (%d,%p)\n", ret, (void*)ans);
	if(dbgev)
		logtime("IMAGE_GETMIM_END", 0);
	is->bs->lim = ed;
	if(is->i > ed)
		is->bs->lim = is->i;
	*pans = ans;
	return ret;
}

// transfer control back to getmim, with given error code
static void
imgerror(ImageSource* is, int err)
{
	if(dbgimg)
		trace("Image error: %S\n", errphrase(err));
	longjmp(is->getmimsave, err);
}

// Get next char or call imgerror if cannot
static int
getc(ImageSource* is)
{
	if(is->i >= is->bs->dalloclen) {
		imgerror(is, ERReof);
	}
	return is->bs->data[is->i++];
}

// Unget the last character.
// When called before any other getting routines, we
// know the buffer still has that character in it.
static void
ungetc(ImageSource* is)
{
	assert(is->i > 0);
	is->i--;
}

// Like ungetc, but ungets two bytes
static void
ungetc2(ImageSource* is)
{
	assert(is->i > 1);
	is->i -= 2;
}

// Get 2 bytes and return the 16-bit value, little-endian order.
static int
getlew(ImageSource* is)
{
	int	c0;
	int	c1;

	c0 = getc(is);
	c1 = getc(is);
	return c0 + (c1 << 8);
}

// Get 2 bytes and return the 16-bit value, big-endian order.
static int
getbew(ImageSource* is)
{
	int	c0;
	int	c1;

	c0 = getc(is);
	c1 = getc(is);
	return (c0 << 8) + c1;
}

// Copy next n bytes of input into buf
// or call imgerror if cannot.
static void
readis(ImageSource* is, uchar* buf, int n)
{
	if(is->i + n < is->bs->dalloclen) {
		memmove(buf, is->bs->data+is->i, n);
		is->i += n;
	}
	else
		imgerror(is, ERReof);
}

// Caller needs n bytes.
// Return index into is->bs->data where at least
// the next n bytes can be found.
// There might be a "premature eof" exception.
static int
getn(ImageSource* is, int n)
{
	int	i;

	i = is->i;
	if(i + n <= is->bs->dalloclen)
		is->i += n;
	else
		imgerror(is, ERReof);
	return i;
}

// draw's allocimage but with some defaults.
// calls imgerror if out of memory.
static Image*
newimage(ImageSource* is, int w, int h)
{
	Image*	im;

	if(!imcacheneed(w*h))
		imgerror(is, ERRnomem);
	im = allocimage(display, Rect(0, 0, w, h), CMAP8, 1, DWhite);
	if(im == nil)
		imgerror(is, ERRnomem);
	return im;
}

static MaskedImage*
newmi(Image* im)
{
	MaskedImage* mi;

	mi = (MaskedImage*)emallocz(sizeof(MaskedImage));
	mi->im = im;
	return mi;
}

// Call this after origw and origh are set to set the width and height
// to our desired (rescaled) answer dimensions.
// If only one of the dimensions is specified, the other is scaled by
// the same factor.
static void
setdims(ImageSource* is)
{
	int	sw;
	int	sh;
	int	dw;
	int	dh;

	sw = is->origw;
	sh = is->origh;
	dw = is->width;
	dh = is->height;
	if(dw == 0 && dh == 0) {
		dw = sw;
		dh = sh;
	}
	else if(dw == 0 || dh == 0) {
		if(dw == 0) {
			dw = (int)((((double)sw)*((double)dh/(double)sh)));
			if(dw == 0)
				dw = 1;
		}
		else {
			dh = (int)((((double)sh)*((double)dw/(double)sw)));
			if(dh == 0)
				dh = 1;
		}
	}
	is->width = dw;
	is->height = dh;
}

// for debugging
static void
tracearray(int* a, int n, char* name)
{
	int	i;

	trace("%s:", name);
	for(i = 0; i < n; i++) {
		if((i%10) == 0)
			trace("\n%5d: ", i);
		trace("%6d", a[i]);
	}
	trace("\n");
}

// XBitmaps

static void
getxbitmaphdr(ImageSource* is)
{
	int	fnd;

	fnd = getxbitmapdefine(is, &is->origw);
	if(fnd)
		fnd = getxbitmapdefine(is, &is->origh);
	if(!fnd)
		imgerror(is, ERRimgbad);
	fnd = getxbitmapdefine(is, nil);
	if(fnd)
		getxbitmapdefine(is, nil);
	get_to_char(is, '{');
}

static MaskedImage*
getxbitmapmim(ImageSource* is)
{
	int	bytesperline;
	uchar*	pixels;
	int	npixels;
	int	pixi;
	int	i;
	Image*	im;
	int	v;
	int	kend;
	int	k;
	int	j;

	getxbitmaphdr(is);
	setdims(is);
	bytesperline = (is->origw + 7)/8;
	npixels = is->origw*is->origh;
	pixels = (uchar*)emalloc(npixels);
	pixi = 0;
	for(i = 0; i < is->origh; i++) {
		for(j = 0; j < bytesperline; j++) {
			v = get_hexbyte(is);
			kend = 7;
			if(j == bytesperline - 1)
				kend = (is->origw - 1)%8;
			for(k = 0; k <= kend; k++) {
				if(v&(1 << k))
					pixels[pixi] = Black;
				else
					pixels[pixi] = White;
				pixi++;
			}
		}
	}
	if(is->width != is->origw || is->height != is->origh)
		pixels = resample(pixels, is->origw, is->origh, is->width, is->height);
	im = newimage(is, is->width, is->height);
	loadimage(im, im->r, pixels, npixels);
	return newmi(im);
}

// get a line, which should be of form
//	'//define fieldname val'
// and return 1 if found.
// If pval != nil, put val in *pval.
static int
getxbitmapdefine(ImageSource* is, int* pval)
{
	int	fnd;
	int	n;
	int	c;

	fnd = 0;
	n = 0;
	c = getc(is);
	if(c == '#') {
		get_to_char(is, ' ');
		get_to_char(is, ' ');
		c = getc(is);
		while(c >= '0' && c <= '9') {
			fnd = 1;
			n = n*10 + c - '0';
			c = getc(is);
		}
	}
	else
		ungetc(is);
	get_to_char(is, '\n');
	if(pval != nil)
		*pval = n;
	return fnd;
}

// read fd until get char cterm
// (getc will call imgerror if eof hit first)
static void
get_to_char(ImageSource* is, int cterm)
{
	while(1) {
		if(getc(is) == cterm)
			return;
	}
}

// read fd until get xDD, were DD are hex digits.
// (getc will call imgerror if eof hit first)
static int
get_hexbyte(ImageSource* is)
{
	int	n1;
	int	n2;

	get_to_char(is, 'x');
	n1 = hexdig(getc(is));
	n2 = hexdig(getc(is));
	if(n1 < 0 || n2 < 0)
		imgerror(is, ERRimgbad);
	return (n1 << 4) + n2;
}

static int
hexdig(int c)
{
	if('0' <= c && c <= '9')
		c -= 48;
	else if('a' <= c && c <= 'f')
		c += 10 - 'a';
	else if('A' <= c && c <= 'F')
		c += 10 - 'A';
	else
		c = -1;
	return c;
}

// GIF

// GIF flags
enum {
	TRANSP = 1,
	INPUT = 2,
	DISPMASK = (7<<2),
	HASCMAP = 0x80,
	INTERLACED = 0x40
};

static void
getgifhdr(ImageSource* is)
{
	Gifstate*	h;
	int	i;
	char*	vers;

	if(dbgimg)
		trace("img: getgifhdr\n");
	h = (Gifstate*)emallocz(sizeof(Gifstate));
	i = getn(is, 6);
	vers = (char*)(is->bs->data+i);
	if(strncmp(vers, "GIF87a", 6) && strncmp(vers, "GIF89a", 6))
		imgerror(is, ERRimgbad);
	is->origw = getlew(is);
	is->origh = getlew(is);
	h->fields = getc(is);
	h->bgrnd = getc(is);
	h->aspect = getc(is);
	setdims(is);
	if(dbgimg)
		trace("img: getgifhdr has origw=%d, origh=%d, w=%d, h=%d, fields=16r%x, bgrnd=%d, aspect=%d\n",
			is->origw, is->origh, is->width, is->height, h->fields, h->bgrnd, h->aspect);
	h->tbl = (GifEntry*)emalloc(4096 * sizeof(GifEntry));
	for(i = 0; i < 258; i++) {
		h->tbl[i].prefix = -1;
		h->tbl[i].exten = i;
	}
	if(h->fields&HASCMAP)
		h->globalcmap = gifreadcmap(is, (h->fields&7) + 1, &h->globalcmaplen);
	is->gstate = h;
	if(warn && h->aspect != 0)
		trace("warning: non-standard aspect ratio in GIF image ignored\n");
	if(!gifgettoimage(is))
		imgerror(is, ERRimgbad);
}

static int
gifgettoimage(ImageSource* is)
{
	Gifstate*	h;
	int	nbytes;
	int	i;
	int	hsize;
	int	hasdata;
	int	c;

	if(dbgimg)
		trace("img: gifgettoimage\n");
	h = is->gstate;
	while(1) {
		// some GIFs omit Trailer
		if(is->i >= is->bs->dalloclen)
			break;

		switch(c = getc(is)) {
		case 0x2C:			// Image Descriptor
			return 1;
			break;

		case 0x21:			// Extension
			hsize = 0;
			hasdata = 0;

			switch(getc(is)) {
			case 0x01:			// Plain Text Extension
				hsize = 14;
				hasdata = 1;
				if(dbgimg)
					trace("gifgettoimage: text extension\n");
				break;
			case 0xF9:			// Graphic Control Extension
				getc(is);			// blocksize (should be 4)
				h->flags = getc(is);
				h->delay = getlew(is);
				h->trindex = getc(is);
				getc(is);			// block terminator (should be 0)
				if(dbgimg)
					trace("gifgettoimage: graphic control flags=16r%x, delay=%d, trindex=%d\n",
						h->flags, h->delay, h->trindex);
				break;
			case 0xFE:			// Comment Extension
				if(dbgimg)
					trace("gifgettoimage: comment extension\n");
				hasdata = 1;
				break;
			case 0xFF:			// Application Extension
				if(dbgimg)
					trace("gifgettoimage: application extension\n");
				hsize = getc(is);
				// standard says this must be 11, but Adobe likes to put out 10-byte ones,
				// so we pay attention to the field.
				hasdata = 1;
				break;
			default:
				imgerror(is, ERRimgbad);
				break;
			}
			if(hsize > 0)
				getn(is, hsize);
			if(hasdata) {
				while(1) {
					if((nbytes = getc(is)) == 0)
						break;
					i = getn(is, nbytes);
					if(dbgimg)
						trace("extension data: '%S'\n", toStr(is->bs->data+i, nbytes, UTF_8));
				}
			}
			break;

		case 0x3B:		// Trailer
			goto loop_done;
			break;

		default:
			if(c == 0)
				continue;				// Fix for some buggy gifs
			imgerror(is, ERRimgbad);
			break;
		}
	}
loop_done:
	return 0;
}

static MaskedImage*
getgifmim(ImageSource* is)
{
	Gifstate*	h;
	int	left;
	int	top;
	int	width;
	int	height;
	int	totpix;
	int	c;
	int	incode;
	int	codesize;
	int	CTM;
	int	EOD;
	uchar*	pic;
	int	pici;
	uchar*	data;
	int	datai;
	int	dataend;
	int	nbits;
	int	sreg;
	uchar*	stack;
	int	stacki;
	int	fc;
	GifEntry*	tbl;
	Image*	mask;
	int	bgcolor;
	int	i;
	Image*	im;
	MaskedImage*	mi;
	int	dispmeth;
	int	nbytes;
	int	code;
	int	nb;
	int	early;
	int	csize;
	int	csmask;
	int	nentry;
	int	maxentry;
	int	first;
	int	ocode;
	int	start;
	int	by;
	uchar*	ipic;
	int	ipici;
	int	ipiclim;
	int	w2;
	int	w4;
	int	w8;
	int	k;
	double	wscale;
	double	hscale;
	int	owidth;
	int	oheight;
	uchar	v;
	int	x;
	int	bytesperrow;
	uchar	trpix;
	uchar*	mpic;
	int	mpici;
	int	y;
	Image*	newmask;
	Rectangle	r;
	int	col;
	Image*	newim;
	int	ilstart[4];
	int	ilby[4];

	if(is->gstate == nil)
		getgifhdr(is);

	// At this point, should just have read Image Descriptor marker byte
	h = is->gstate;
	left = getlew(is);
	top = getlew(is);
	width = getlew(is);
	height = getlew(is);
	h->fields = getc(is);
	totpix = width*height;
	h->cmap = h->globalcmap;
	h->cmaplen = h->globalcmaplen;
	if(dbgimg)
		trace("getgifmim, left=%d, top=%d, width=%d, height=%d, fields=16r%x\n",
			left, top, width, height, h->fields);
	if(dbgev)
		logtime("IMAGE_GETGIFMIM", 0);
	if(h->fields&HASCMAP)
		h->cmap = gifreadcmap(is, (h->fields&7) + 1, &h->cmaplen);
	if(h->cmap == nil)
		imgerror(is, ERRimgbad);

	// now decode the image
	codesize = getc(is);
	if(codesize > 8)
		imgerror(is, ERRimgbad);
	if(h->cmaplen != 3*(1 << codesize) && (codesize != 2 || h->cmaplen != 3*2))
		imgerror(is, ERRimgbad);

	CTM = 1 << codesize;
	EOD = CTM + 1;

	pic = (uchar*)emalloc(totpix);
	pici = 0;
	data = is->bs->data;
	datai = 0;
	dataend = 0;

	nbits = 0;
	sreg = 0;
	stack = (uchar*)emalloc(4096);
	fc = 0;
	tbl = h->tbl;
	while(1) {
Decode_continue:
		csize = codesize + 1;
		csmask = ((1 << csize) - 1);
		nentry = EOD + 1;
		maxentry = csmask;
		first = 1;

		for(ocode = -1; ; ocode = incode) {
			while(nbits < csize) {
				if(datai == dataend) {
					nbytes = getc(is);
					if(nbytes == 0)
						goto Decode_done;
					datai = getn(is, nbytes);
					dataend = datai + nbytes;
				}
				c = data[datai++];
				sreg |= c << nbits;
				nbits += 8;
			}
			code = sreg&csmask;
			sreg >>= csize;
			nbits -= csize;

			if(code == EOD) {
				nbytes = getc(is);
				if(nbytes != 0 && warn)
					trace("warning: unexpected data past EOD\n");
				goto Decode_done;
			}
			if(code == CTM)
				goto Decode_continue;

			stacki = 4096 - 1;

			incode = code;

			// special case for KwKwK 
			if(code == nentry) {
				stack[stacki--] = fc;
				code = ocode;
			}

			if(code > nentry)
				imgerror(is, ERRimgbad);

			for(c = code; c >= 0; c = tbl[c].prefix)
				stack[stacki--] = tbl[c].exten;

			nb = 4096 - (stacki + 1);
			if(pici + nb > totpix) {
				// this common error is harmless
				// we have to keep reading to keep the blocks in sync
				;
			}
			else {
				memmove(pic+pici, stack+stacki+1, nb);
				pici += nb;
			}

			fc = stack[stacki + 1];

			if(first) {
				first = 0;
				continue;
			}

			early = 0; // peculiar tiff feature here for reference
			if(nentry == maxentry - early) {
				if(csize >= 12)
					continue;
				csize++;
				maxentry = (1 << csize);
				csmask = maxentry - 1;
				if(csize < 12)
					maxentry--;
			}
			tbl[nentry].prefix = ocode;
			tbl[nentry].exten = fc;
			nentry++;
		}
	}
Decode_done:
	while(pici < totpix) {
		// shouldn't happen, but sometimes get buggy gifs
		pic[pici++] = 0;
	}

	if(h->fields&INTERLACED) {
		if(dbgimg)
			trace("getgifmim uninterlacing\n");
		if(dbgev)
			logtime("IMAGE_GETGIFMIM_INTERLACE_START", 0);
		// (TODO: Could un-interlace in place.
		// Decompose permutation into cycles,
		// then need one double-copy of a line
		// per cycle).
		ipic = (uchar*)emalloc(totpix);
		pici = 0;
		ipiclim = totpix - width;
		w2 = width + width;
		w4 = w2 + w2;
		w8 = w4 + w4;
		ilstart[0] = 0;
		ilstart[1] = w4;
		ilstart[2] = w2;
		ilstart[3] = width;
		ilby[0] = w8;
		ilby[1] = w8;
		ilby[2] = w4;
		ilby[3] = w2;
		for(k = 0; k < 4; k++) {
			start = ilstart[k];
			by = ilby[k];
			for(ipici = start; ipici <= ipiclim; ipici += by) {
				memmove(ipic+ipici, pic+pici, width);
				pici += width;
			}
		}
		pic = ipic;
		if(dbgev)
			logtime("IMAGE_GETGIFMIM_INTERLACE_END", 0);
	}
	if(is->width != is->origw || is->height != is->origh) {
		// need to resample, using same factors as original image
		wscale = (double)is->width/(double)is->origw;
		hscale = (double)is->height/(double)is->origh;
		owidth = width;
		oheight = height;
		width = (int)((wscale*(double)width)+.5);
		if(width == 0)
			width = 1;
		height = (int)((hscale*(double)height)+.5);
		if(height == 0)
			height = 1;
		left = (int)((wscale*(double)left)+.5);
		top = (int)((hscale*(double)top)+.5);
		pic = resample(pic, owidth, oheight, width, height);
	}
	mask = nil;
	if(h->flags&TRANSP) {
		if(dbgimg)
			trace("getgifmim making mask, trindex=%d\n", h->trindex);
		if(dbgev)
			logtime("IMAGE_GETGIFMIM_MASK_START", 0);
		// make a 1-bit deep bitmap for mask
		// expect most mask bits will be 1
		bytesperrow = (width + 7)/8;
		trpix = h->trindex;
		nb = bytesperrow*height;
		mpic = (uchar*)emalloc(nb);
		mpici = 0;
		pici = 0;
		for(y = 0; y < height; y++) {
			v = 0xFF;
			k = 0;
			for(x = 0; x < width; x++) {
				if(pic[pici++] == trpix)
					v &= ~(0x80 >> k);
				if(++k == 8) {
					k = 0;
					mpic[mpici++] = v;
					v = 0xFF;
				}
			}
			if(k != 0)
				mpic[mpici++] = v;
		}
		if(!imcacheneed(nb))
			imgerror(is, ERRnomem);
		mask = allocimage(display, Rect(0, 0, width, height), GREY1, 0, DOpaque);
		if(mask == nil)
			imgerror(is, ERRnomem);
		loadimage(mask, mask->r, mpic, nb);
		if(dbgev)
			logtime("IMAGE_GETGIFMIM_MASK_END", 0);
	}
	if(dbgev)
		logtime("IMAGE_GETGIFMIM_REMAP_START", 0);
	remap1(pic, width, height, h->cmap, h->cmaplen);
	if(dbgev)
		logtime("IMAGE_GETGIFMIM_REMAP_END", 0);
	bgcolor = -1;
	i = h->bgrnd;
	if(i >= 0 && 3*i + 2 < h->cmaplen) {
		bgcolor = ((h->cmap[3*i]) << 16)
			| ((h->cmap[3*i + 1]) << 8)
			| (h->cmap[3*i + 2]);
	}
	im = newimage(is, width, height);
	loadimage(im, im->r, pic, totpix);
	if(is->curframe == 0) {
		// make sure first frame fills up whole rectangle
		if(is->width != width || is->height != height || left != 0 || top != 0) {
			r = Rect(left, top, left + width, top + height);
			col = DWhite;
			if(bgcolor != -1)
				col = (bgcolor<<8)|0xFF;
			newim = allocimage(display, Rect(0, 0, is->width, is->height), im->chan, 0, col);
			if(newim == nil)
				imgerror(is, ERRnomem);
			draw(newim, r, im, mask, ZP);
			im = newim;
			if(mask != nil) {
				newmask = allocimage(display, Rect(0, 0, is->width, is->height), GREY1, 0, DTransparent);
				if(newmask == nil)
					imgerror(is, ERRnomem);
				draw(newmask, r, mask, nil, ZP);
				mask = newmask;
			}
			left = 0;
			top = 0;
		}
	}
	mi = newmi(im);
	mi->mask = mask;
	mi->delay = h->delay*10;	// convert centiseconds to milliseconds
	mi->origin = Pt(left, top);
	dispmeth = (h->flags >> 2)&7;
	if(dispmeth == 2) {
		// reset to background color after displaying this frame
		mi->bgcolor = bgcolor;
	}
	else if(dispmeth == 3) {
		// Supposed to "reset to previous", which appears to
		// mean the previous frame that didn't have a "reset to previous".
		// Signal this special case to layout by setting bgcolor to -2.
		mi->bgcolor = -2;
	}
	if(gifgettoimage(is)) {
		mi->more = 1;
		is->curframe++;
		// have to reinitialize table for next time
		for(i = 0; i < 258; i++) {
			h->tbl[i].prefix = -1;
			h->tbl[i].exten = i;
		}
	}
	if(dbgev)
		logtime("IMAGE_GETGIFMIM_END", 0);
	return mi;
}

// Read a GIF colormap, where bpe is number of bits in an entry.
// Causes a 'premature eof' imgerror if can't get the whole map.
static uchar*
gifreadcmap(ImageSource* is, int bpe, int* pmaplen)
{
	int	size;
	uchar*	map;

	size = 3*(1 << bpe);
	map = (uchar*)emalloc(size);
	if(dbgimg > 1)
		trace("gifreadcmap wants %d bytes\n", size);
	readis(is, map, size);
	*pmaplen = size;
	return map;
}

// JPG

// Constants
enum {
	// All of these are preceded by 0xFF byte
	SOF = 0xC0,	// Start of Frame
	SOF2 = 0xC2,	// Start of Frame; progressive Huffman
	JPG = 0xC8,	// Reserved for JPEG extensions
	DHT = 0xC4,	// Define Huffman Tables
	DAC = 0xCC,	// Arithmetic coding conditioning
	RST = 0xD0,	// Restart interval termination
	RST7 = 0xD7,	// Restart interval termination (highest value)
	SOI = 0xD8,	// Start of Image
	EOI = 0xD9,	// End of Image
	SOS = 0xDA,	// Start of Scan
	DQT = 0xDB,	// Define quantization tables
	DNL = 0xDC,	// Define number of lines
	DRI = 0xDD,	// Define restart interval
	DHP = 0xDE,	// Define hierarchical progression
	EXP = 0xDF,	// Expand reference components
	APPn = 0xE0,	// Reserved for application segments
	JPGn = 0xF0,	// Reserved for JPEG extensions
	COM = 0xFE,	// Comment

	NBUF = 16*1024,

	// channel descriptions
	CRGB = 0,		// three channels, R, G, B, no map
	CY = 1,		// one channel, luminance
	CRGB1 = 2,	// one channel, map present
	CYCbCr = 3,	// three channels, Y, Cb, Cr, no map

	jpegcolorspace = CYCbCr,	// default space to convert jpeg to

	// Some floats converted to fixed point (scaling by 2048)
	jc1 = 2871,		// 1.402 * 2048
	jc2 = 705,			// 0.34414 * 2048
	jc3 = 1463,		// 0.71414 * 2048
	jc4 = 3629,		// 1.772 * 2048

	W1 = 2841,		// 2048*sqrt(2)*cos(1*pi/16)
	W2 = 2676,		// 2048*sqrt(2)*cos(2*pi/16)
	W3 = 2408,		// 2048*sqrt(2)*cos(3*pi/16)
	W5 = 1609,		// 2048*sqrt(2)*cos(5*pi/16)
	W6 = 1108,		// 2048*sqrt(2)*cos(6*pi/16)
	W7 = 565,			// 2048*sqrt(2)*cos(7*pi/16)

	W1pW7 = 3406,	// W1+W7
	W1mW7 = 2276,	// W1-W7
	W3pW5 = 4017,	// W3+W5
	W3mW5 = 799,		// W3-W5
	W2pW6 = 3784,	// W2+W6
	W2mW6 = 1567,	// W2-W6

	R2 = 181,			// 256/sqrt(2)
};

static void
getjpeghdr(ImageSource* is)
{
	Jpegstate*	h;
	int	m;
	int	n;
	int	i;
	int	H;
	int	V;

	if(dbgimg)
		trace("getjpeghdr\n");
	h = (Jpegstate*)emallocz(sizeof(Jpegstate));
	is->jstate = h;
	if(jpegmarker(is) != SOI)
		imgerror(is, ERRimgbad);
	jpegtabmisc(is, &m, &n);
	if(!(m == SOF || m == SOF2))
		imgerror(is, ERRimgbad);
	getc(is);			// sample precision
	h->Y = getbew(is);
	h->X = getbew(is);
	h->Nf = getc(is);
	if(dbgimg)
		trace("start of frame, Y=%d, X=%d, Nf=%d\n", h->Y, h->X, h->Nf);
	h->comp = (Framecomp*)emalloc(h->Nf * sizeof(Framecomp));
	h->nblock = (int*)emalloc(h->Nf * sizeof(int));
	for(i = 0; i < h->Nf; i++) {
		h->comp[i].C = getc(is);
		nibbles(getc(is), &H, &V);
		h->comp[i].H = H;
		h->comp[i].V = V;
		h->comp[i].Tq = getc(is);
		h->nblock[i] = H*V;
		if(dbgimg)
			trace("comp[%d]: C=%d, H=%d, V=%d, Tq=%d\n",
				i, h->comp[i].C, H, V, h->comp[i].Tq);
	}
	h->mode = m;
	is->origw = h->X;
	is->origh = h->Y;
	setdims(is);
	if(n != 6 + 3*h->Nf)
		imgerror(is, ERRimgbad);
}

static int
jpegmarker(ImageSource* is)
{
	if(getc(is) != 0xFF)
		imgerror(is, ERRimgbad);
	return getc(is);
}

// Consume tables and miscellaneous marker segments,
// returning the marker id and length of the first non-such-segment
// (after having consumed the marker) in *pm and *pn.
static void
jpegtabmisc(ImageSource* is, int *pm, int *pn)
{
	Jpegstate*	h;
	int	m;
	int	n;
	int	vers0;
	int	vers1;
	uchar*	buf;
	int	i;

	h = is->jstate;
	while(1) {
		h->nseg++;
		m = jpegmarker(is);
		n = 0;
		if(m != EOI)
			n = getbew(is) - 2;
		if(dbgimg > 1)
			trace("jpegtabmisc reading segment, got m=%x, n=%d\n", m, n);

		switch(m) {
		case SOF:
		case SOF2:
		case SOS:
		case EOI:
			goto Loop_done;
			break;

		case APPn+0:
			if(h->nseg == 1 && n >= 6) {
				i = getn(is, 6);
				buf = is->bs->data;
				n -= 6;
				if(!strncmp((char*)(buf+i), "JFIF", 4)) {
					vers0 = buf[i + 5];
					vers1 = buf[i + 6];
					if(vers0 > 1 || vers1 > 2)
						imgerror(is, ERRimgbad);		// unimplemented version
				}
			}
			break;

		case DQT:
			jpegquanttables(is, n);
			n = 0;
			break;

		case DHT:
			jpeghuffmantables(is, n);
			n = 0;
			break;
		case DRI:
			h->ri = getbew(is);
			n -= 2;
			break;

		case COM:
			break;

		default:
			if(!(m >= APPn+1 && m <= APPn+15))
				imgerror(is, ERRimgbad);		// unexpected marker
			break;
		}
		if(n > 0)
			getn(is, n);
	}
Loop_done:
	*pm = m;
	*pn = n;
}

// Consume huffman tables
static void
jpeghuffmantables(ImageSource* is, int n)
{
	Jpegstate*	h;
	int	l;

	if(dbgimg)
		trace("jpeghuffmantables\n");
	h = is->jstate;
	if(h->dcht == nil) {
		h->dcht = (Huffman**)emallocz(4 * sizeof(Huffman*));
		h->acht = (Huffman**)emallocz(4 * sizeof(Huffman*));
	}
	l = 0;
	while(l < n)
		l += jpeghuffmantable(is);
	if(l != n)
		imgerror(is, ERRimgbad);		// huffman table bad length
}

static int
jpeghuffmantable(ImageSource* is)
{
	Huffman*	t;
	Jpegstate*	h;
	int	Tc;
	int	th;
	uchar*	b;
	int	bi;
	int	nsize;
	int	i;
	int	k;
	int	code;
	int	si;
	int	j;
	int*	maxcode;
	int	v;
	int	n;
	int	cnt;
	int	m;
	int	sr;
	int	numcodes[16];

	t = (Huffman*)emalloc(sizeof(Huffman));
	h = is->jstate;
	nibbles(getc(is), &Tc, &th);
	if(dbgimg > 1)
		trace("jpeghuffmantable, Tc=%d, th=%d\n", Tc, th);
	if(Tc > 1)
		imgerror(is, ERRimgbad);		// unknown Huffman table class
	if(th > 3 || (h->mode == (uchar)SOF && th > 1))
		imgerror(is, ERRimgbad);		// unknown Huffman table index
	if(Tc == 0)
		h->dcht[th] = t;
	else
		h->acht[th] = t;

	// flow chart C-2
	bi = getn(is, 16);
	b = is->bs->data;
	nsize = 0;
	for(i = 0; i < 16; i++)
		nsize += (numcodes[i] = b[bi + i]);
	t->size = (int*)emalloc((nsize + 1) * sizeof(int));
	k = 0;
	for(i = 1; i <= 16; i++) {
		n = numcodes[i - 1];
		for(j = 0; j < n; j++)
			t->size[k++] = i;
	}
	t->size[k] = 0;

	// initialize HUFFVAL
	t->val = (int*)emalloc(nsize * sizeof(int));
	bi = getn(is, nsize);
	for(i = 0; i < nsize; i++)
		t->val[i] = b[bi++];

	// flow chart C-3
	t->code = (int*)emalloc((nsize + 1) * sizeof(int));
	k = 0;
	code = 0;
	si = t->size[0];
	while(1) {
		do
			t->code[k++] = code++;
		while(t->size[k] == si);
		if(t->size[k] == 0)
			break;
		do {
			code <<= 1;
			si++;
		} while(t->size[k] != si);
	}

	// flow chart F-25
	i = 0;
	j = 0;
	while(1) {
		while(1) {
			i++;
			if(i > 16)
				goto F25_done;
			if(numcodes[i - 1] != 0)
				break;
			t->maxcode[i] = -1;
		}
		t->valptr[i] = j;
		t->mincode[i] = t->code[j];
		j += numcodes[i - 1] - 1;
		t->maxcode[i] = t->code[j];
		j++;
	}
F25_done:

	// stupid startup algorithm: just run machine for each byte value
	maxcode = t->maxcode;
	for(v = 0; v < 256; v++) {
		cnt = 7;
		m = 1 << 7;
		code = 0;
		sr = v;
		i = 1;
		for(; ; i++) {
			if(sr&m)
				code |= 1;
			if(code <= maxcode[i])
				break;
			code <<= 1;
			m >>= 1;
			if(m == 0) {
				t->shift[v] = 0;
				t->value[v] = -1;
				goto Bytes_continue;
			}
			cnt--;
		}
		t->shift[v] = 8 - cnt;
		t->value[v] = t->val[t->valptr[i] + (code - t->mincode[i])];
Bytes_continue:
		;
	}
	if(dbgimg > 2) {
		trace("Huffman table %d:\n", th);
		tracearray(t->size, nsize+1, "size");
		tracearray(t->code, nsize+1, "code");
		tracearray(t->val, nsize, "val");
		tracearray(t->mincode, 17, "mincode");
		tracearray(t->maxcode, 17, "maxcode");
		tracearray(t->value, 256, "value");
		tracearray(t->shift, 256, "shift");
	}
	return nsize + 17;
}

static void
jpegquanttables(ImageSource* is, int n)
{
	Jpegstate*	h;
	int	l;

	if(dbgimg)
		trace("jpegquanttables\n");
	h = is->jstate;
	if(h->qt == nil)
		h->qt = (int**)emalloc(4 * sizeof(int*));
	l = 0;
	while(l < n)
		l += jpegquanttable(is);
	if(l != n)
		imgerror(is, ERRimgbad);
}

static int
jpegquanttable(ImageSource* is)
{
	int	pq;
	int	tq;
	int*	q;
	int	i;

	nibbles(getc(is), &pq, &tq);
	if(dbgimg)
		trace("jpegquanttable pq=%d tq=%d\n", pq, tq);
	if(pq > 1)
		imgerror(is, ERRimgbad);		// unknown quantization table class
	if(tq > 3)
		imgerror(is, ERRimgbad);		// bad quantization table index
	q = (int*)emalloc(64 * sizeof(int));
	is->jstate->qt[tq] = q;
	for(i = 0; i < 64; i++) {
		if(pq == 0)
			q[i] = getc(is);
		else
			q[i] = getbew(is);
	}
	if(dbgimg > 2)
		tracearray(q, 64, "quant table");
	return 1 + (64*(1 + pq));
}

// Have just read Frame header.
// Now expect:
//	((tabl/misc segment(s))* (scan header) (entropy coded segment)+)+ EOI
static MaskedImage*
getjpegmim(ImageSource* is)
{
	Jpegstate*	h;
	uchar**	chans;
	int	width;
	int	height;
	Image*	im;
	int	m;
	int	n;
	Scancomp*	scomp;
	int	i;
	int	k;
	int	nchans;

	if(dbgimg)
		trace("getjpegmim\n");
	if(dbgev)
		logtime("IMAGE_GETJPGMIM", is->width*is->height);
	getjpeghdr(is);
	h = is->jstate;
	chans = nil;
	while(1) {
		jpegtabmisc(is, &m, &n);
		if(m == EOI)
			break;
		if(m != SOS)
			imgerror(is, ERRimgbad);		// expected start of scan
		h->Ns = getc(is);
		if(dbgimg)
			trace("start of scan, Ns=%d\n", h->Ns);
		scomp = (Scancomp*)emalloc(h->Ns * sizeof(Scancomp));
		for(i = 0; i < h->Ns; i++) {
			scomp[i].C = getc(is);
			nibbles(getc(is), &scomp[i].tdc, &scomp[i].tac);
		}
		h->scomp = scomp;
		h->Ss = getc(is);
		h->Se = getc(is);
		nibbles(getc(is), &h->Ah, &h->Al);
		if(n != 4 + h->Ns*2)
			imgerror(is, ERRimgbad);		// SOS header wrong length
		if(h->mode == SOF) {
			if(chans != nil)
				imgerror(is, ERRimgbad);	// baseline has > 1 scan
			chans = jpegbaselinescan(is, &nchans);
		}
		else
			jpegprogressivescan(is);
	}
	if(h->mode == SOF2)
		chans = jprogressiveIDCT(is, &nchans);
	if(chans == nil)
		imgerror(is, ERRimgbad);		// jpeg has no image
	width = is->width;
	height = is->height;
	if(width != h->X || height != h->Y) {
		for(k = 0; k < nchans; k++)
			chans[k] = resample(chans[k], h->X, h->Y, width, height);
	}
	if(dbgev)
		logtime("IMAGE_JPG_REMAP", 0);
	if(nchans == 1)
		remapgrey(chans[0], width, height);
	else
		remaprgb(chans, width, height);
	if(dbgev)
		logtime("IMAGE_JPG_REMAP_END", 0);
	im = newimage(is, width, height);
	loadimage(im, im->r, chans[0], width*height);
	if(dbgev)
		logtime("IMAGE_GETJPGMIM_END", 0);
	return newmi(im);
}

static int zig[64] = {
	0, 1, 8, 16, 9, 2, 3, 10, 17, // 0-7
	24, 32, 25, 18, 11, 4, 5, //# 8-15
	12, 19, 26, 33, 40, 48, 41, 34, // 16-23
	27, 20, 13, 6, 7, 14, 21, 28, // 24-31
	35, 42, 49, 56, 57, 50, 43, 36, // 32-39
	29, 22, 15, 23, 30, 37, 44, 51, // 40-47
	58, 59, 52, 45, 38, 31, 39, 46, // 48-55
	53, 60, 61, 54, 47, 55, 62, 63 // 56-63
};

static uchar**
jpegbaselinescan(ImageSource* is, int *pnchans)
{
	Jpegstate*	h;
	int	Ns;
	uchar**	chans;
	int	k;
	int*	Td;
	int*	Ta;
	int***	data;
	int*	H;
	int*	V;
	int*	DC;
	int	Hmax;
	int	Vmax;
	int	comp;
	int	allHV1;
	int	ri;
	int	nacross;
	int	nmcu;
	int	mcu;
	int	nblock;
	int	m;
	int	z;
	int	rs;
	int	rrrr;
	int	ssss;
	int	t;
	int	diff;
	int*	zz;
	Huffman*	dcht;
	Huffman*	acht;
	int*	qt;
	int	block;

	if(dbgimg)
		trace("jpegbaselinescan\n");
	if(dbgev)
		logtime("IMAGE_JPGBASELINESCAN", 0);
	h = is->jstate;
	Ns = h->Ns;
	if(Ns != h->Nf)
		imgerror(is, ERRimgbad);		// baseline needs Ns==Nf
	if(!(Ns == 3 || Ns == 1))
		imgerror(is, ERRimgbad);		// baseline needs Ns==1 or 3
	chans = (uchar**)emalloc(h->Nf * sizeof(uchar*));
	*pnchans = h->Nf;
	for(k = 0; k < h->Nf; k++)
		chans[k] = (uchar*)emalloc(h->X*h->Y);

	// build per-component arrays
	Td = (int*)emalloc(Ns * sizeof(int));
	Ta = (int*)emalloc(Ns * sizeof(int));
	data = (int***)emalloc(Ns * sizeof(int**));
	H = (int*)emalloc(Ns * sizeof(int));
	V = (int*)emalloc(Ns * sizeof(int));
	DC = (int*)emalloc(Ns * sizeof(int));

	// compute maximum H and V
	Hmax = 0;
	Vmax = 0;
	for(comp = 0; comp < Ns; comp++) {
		if(h->comp[comp].H > Hmax)
			Hmax = h->comp[comp].H;
		if(h->comp[comp].V > Vmax)
			Vmax = h->comp[comp].V;
	}
	if(dbgimg > 1)
		trace("Hmax=%d, Vmax=%d\n", Hmax, Vmax);

	// initialize data structures
	allHV1 = 1;
	for(comp = 0; comp < Ns; comp++) {
		// JPEG requires scan components to be in same order as in frame,
		// so if both have 3 we know scan is Y Cb Cr and there's no need to
		// reorder
		Td[comp] = h->scomp[comp].tdc;
		Ta[comp] = h->scomp[comp].tac;
		H[comp] = h->comp[comp].H;
		V[comp] = h->comp[comp].V;
		nblock = H[comp]*V[comp];
		if(nblock != 1)
			allHV1 = 0;
		data[comp] = (int**)emalloc(nblock * sizeof(int*));
		DC[comp] = 0;
		for(m = 0; m < nblock; m++)
			data[comp][m] = (int*)emalloc(8*8 * sizeof(int));
		if(dbgimg > 2)
			trace("scan comp %d: H=%d, V=%d, nblock=%d, Td=%d, Ta=%d\n",
					comp, H[comp], V[comp], nblock, Td[comp], Ta[comp]);
	}

	ri = h->ri;

	h->cnt = 0;
	h->sr = 0;
	nacross = ((h->X + (8*Hmax - 1))/(8*Hmax));
	nmcu = ((h->Y + (8*Vmax - 1))/(8*Vmax))*nacross;
	if(dbgimg)
		trace("nacross=%d, nmcu=%d\n", nacross, nmcu);
	mcu = 0;
	while(mcu < nmcu) {
		if(dbgimg > 2)
			trace("mcu %d\n", mcu);
		for(comp = 0; comp < Ns; comp++) {
			if(dbgimg > 2)
				trace("comp %d\n", comp);
			dcht = h->dcht[Td[comp]];
			acht = h->acht[Ta[comp]];
			qt = h->qt[h->comp[comp].Tq];

			for(block = 0; block < H[comp]*V[comp]; block++) {
				if(dbgimg > 2)
					trace("block %d\n", block);
				// F-22
				t = jdecode(is, dcht);
				diff = jreceive(is, t);
				DC[comp] += diff;
				if(dbgimg > 2)
					trace("t=%d, diff=%d, DC=%d\n", t, diff, DC[comp]);

				// F=23
				zz = data[comp][block];
				memset(zz, 0, 64*sizeof(int));
				zz[0] = qt[0]*DC[comp];
				k = 1;

				while(1) {
					rs = jdecode(is, acht);
					nibbles(rs, &rrrr, &ssss);
					if(ssss == 0) {
						if(rrrr != 15)
							break;
						k += 16;
					}
					else {
						k += rrrr;
						z = jreceive(is, ssss);
						zz[zig[k]] = z*qt[k];
						if(k == 63)
							break;
						k++;
					}
				}
				idct(zz);
			}
		}

		// rotate colors to RGB and assign to bytes
		if(Ns == 1)
			colormap1(h, chans[0], data[0][0], mcu, nacross);
		else if(allHV1)
			colormapall1(h, chans, data[0][0], data[1][0], data[2][0], mcu, nacross);
		else
			colormap(h, chans, data[0], data[1], data[2], mcu, nacross, Hmax, Vmax, H, V);

		// process restart marker, if present
		mcu++;
		if(ri > 0 && mcu < nmcu && mcu%ri == 0) {
			jrestart(is, mcu);
			for(comp = 0; comp < Ns; comp++)
				DC[comp] = 0;
		}
	}
	if(dbgev)
		logtime("IMAGE_JPGBASELINESCAN_END", 0);
	return chans;
}

static void
jrestart(ImageSource* is, int mcu)
{
	Jpegstate*	h;
	int	ri;
	int	restart;
	int	rst;
	int	nskip;

	h = is->jstate;
	ri = h->ri;
	restart = mcu/ri - 1;
	nskip = 0;
	do {
		do {
			rst = jnextborm(is);
			nskip++;
		} while(rst >= 0 && rst != 0xFF);
		if(rst == 0xFF) {
			rst = jnextborm(is);
			nskip++;
		}
	} while(rst >= 0 && (rst&~7) != RST);
	if(nskip != 2 || rst < 0 || ((rst&7) != (restart&7)))
		imgerror(is, ERRimgbad);		// restart problem
	h->cnt = 0;
	h->sr = 0;
}

static void
jpegprogressivescan(ImageSource* is)
{
	Jpegstate*	h;
	int	c;
	int	comp;
	int	i;

	if(dbgev)
		logtime("IMAGE_JPGPROGSCAN", 0);
	h = is->jstate;
	if(h->dccoeff == nil)
		jprogressiveinit(is, h);

	c = h->scomp[0].C;
	comp = -1;
	for(i = 0; i < h->Nf; i++)
		if(h->comp[i].C == c)
			comp = i;
	if(comp == -1)
		imgerror(is, ERRimgbad);		// bad component index in scan header

	if(h->Ss == 0)
		jprogressivedc(is, comp);
	else if(h->Ah == 0)
		jprogressiveac(is, comp);
	else
		jprogressiveacinc(is, comp);
	if(dbgev)
		logtime("IMAGE_JPGPROGSCAN_END", 0);
}

static uchar**
jprogressiveIDCT(ImageSource* is, int* pnchans)
{
	Jpegstate*	h;
	int	Nf;
	int*	H;
	int*	V;
	int	allHV1;
	int***	data;
	int	comp;
	uchar**	chans;
	int	k;
	int*	blockno;
	int	nmcu;
	int	mcu;
	int	nblock;
	int	m;
	int*	zz;
	int*	dccoeff;
	int**	accoeff;
	int	bn;
	int	block;

	if(dbgev)
		logtime("IMAGE_JPGPROGIDCT", 0);
	h = is->jstate;
	Nf = h->Nf;
	H = (int*)emalloc(Nf * sizeof(int));
	V = (int*)emalloc(Nf * sizeof(int));
	allHV1 = 1;
	data = (int***)emalloc(Nf * sizeof(int**));
	for(comp = 0; comp < Nf; comp++) {
		H[comp] = h->comp[comp].H;
		V[comp] = h->comp[comp].V;
		nblock = h->nblock[comp];
		if(nblock != 1)
			allHV1 = 0;
		data[comp] = (int**)emalloc(nblock * sizeof(int*));
		for(m = 0; m < nblock; m++)
			data[comp][m] = (int*)emalloc(8*8 * sizeof(int));
	}
	chans = (uchar**)emalloc(h->Nf*sizeof(uchar*));
	*pnchans = h->Nf;
	for(k = 0; k < h->Nf; k++)
		chans[k] = (uchar*)emalloc(h->X*h->Y);
	blockno = (int*)emallocz(Nf * sizeof(int));
	nmcu = h->nacross*h->ndown;
	for(mcu = 0; mcu < nmcu; mcu++) {
		for(comp = 0; comp < Nf; comp++) {
			dccoeff = h->dccoeff[comp];
			accoeff = h->accoeff[comp];
			bn = blockno[comp];
			for(block = 0; block < h->nblock[comp]; block++) {
				zz = data[comp][block];
				memset(zz, 0, 8*8 * sizeof(int));
				zz[0] = dccoeff[bn];
				for(k = 1; k < 64; k++)
					zz[zig[k]] = accoeff[bn][k];
				idct(zz);
				bn++;
			}
			blockno[comp] = bn;
		}

		// rotate colors to RGB and assign to bytes
		if(Nf == 1)
			colormap1(h, chans[0], data[0][0], mcu, h->nacross);
		else if(allHV1)
			colormapall1(h, chans, data[0][0], data[1][0], data[2][0], mcu, h->nacross);
		else
			colormap(h, chans, data[0], data[1], data[2], mcu, h->nacross, h->Hmax, h->Vmax, H, V);
	}
	return chans;
}

static void
jprogressiveinit(ImageSource* is, Jpegstate* h)
{
	int	Ns;
	int	Nf;
	int	comp;
	int	nmcu;
	int	k;
	int	j;

	Ns = h->Ns;
	Nf = h->Nf;
	if((Ns != 3 && Ns != 1) || Ns != Nf)
		imgerror(is, ERRimgbad);
	h->Hmax = 0;
	h->Vmax = 0;
	for(comp = 0; comp < Nf; comp++) {
		if(h->comp[comp].H > h->Hmax)
			h->Hmax = h->comp[comp].H;
		if(h->comp[comp].V > h->Vmax)
			h->Vmax = h->comp[comp].V;
	}
	h->nacross = ((h->X + (8*h->Hmax - 1))/(8*h->Hmax));
	h->ndown = ((h->Y + (8*h->Vmax - 1))/(8*h->Vmax));
	nmcu = h->nacross*h->ndown;
	h->dccoeff = (int**)emalloc(Nf * sizeof(int*));
	h->accoeff = (int***)emalloc(Nf * sizeof(int**));
	for(k = 0; k < Nf; k++) {
		h->dccoeff[k] = (int*)emallocz(h->nblock[k]*nmcu * sizeof(int));
		h->accoeff[k] = (int**)emalloc(h->nblock[k]*nmcu * sizeof(int*));
		for(j = 0; j < h->nblock[k]*nmcu; j++)
			h->accoeff[k][j] = (int*)emallocz(64 * sizeof(int));
	}
}

static void
jprogressivedc(ImageSource* is, int comp)
{
	Jpegstate*	h;
	int	Ns;
	int	Ah;
	int	Al;
	int*	Td;
	int*	DC;
	int	ri;
	int	nmcu;
	int*	blockno;
	int	mcu;
	int	t;
	int	diff;
	Huffman*	dcht;
	int	qt;
	int*	dc;
	int	bn;
	int	block;

	h = is->jstate;
	Ns = h->Ns;
	Ah = h->Ah;
	Al = h->Al;
	if(Ns != h->Nf)
		imgerror(is, ERRimgbad);		// progressive with Nf!=Ns in DC scan

	// build per-component arrays
	Td = (int*)emalloc(Ns * sizeof(int));
	DC = (int*)emalloc(Ns * sizeof(int));
	h->cnt = 0;
	h->sr = 0;
	for(comp = 0; comp < Ns; comp++) {
		// JPEG requires scan components to be in same order as in frame,
		// so if both have 3 we know scan is Y Cb Cr and there's no need to
		// reorder
		Td[comp] = h->scomp[comp].tdc;
		DC[comp] = 0;
	}

	ri = h->ri;

	nmcu = h->nacross*h->ndown;
	blockno = (int*)emallocz(Ns * sizeof(int));
	mcu = 0;
	while(mcu < nmcu) {
		for(comp = 0; comp < Ns; comp++) {
			dcht = h->dcht[Td[comp]];
			qt = h->qt[h->comp[comp].Tq][0];
			dc = h->dccoeff[comp];
			bn = blockno[comp];

			for(block = 0; block < h->nblock[comp]; block++) {
				if(Ah == 0) {
					t = jdecode(is, dcht);
					diff = jreceive(is, t);
					DC[comp] += diff;
					dc[bn] = qt*DC[comp] << Al;
				}
				else
					dc[bn] |= qt*jreceivebit(is) << Al;
				bn++;
			}
			blockno[comp] = bn;
		}

		// process restart marker, if present
		mcu++;
		if(ri > 0 && mcu < nmcu && mcu%ri == 0) {
			jrestart(is, mcu);
			for(comp = 0; comp < Ns; comp++)
				DC[comp] = 0;
		}
	}
}

static void
jprogressiveac(ImageSource* is, int comp)
{
	Jpegstate*	h;
	int	Ns;
	int	Al;
	int	Ss;
	int	Se;
	int	H;
	int	V;
	int	nacross;
	int	ndown;
	int	q;
	int	nhor;
	int	nver;
	int	Ta;
	int	ri;
	int	eobrun;
	Huffman*	acht;
	int*	qt;
	int	nmcu;
	int	mcu;
	int	y;
	int	z;
	int	rs;
	int	rrrr;
	int	ssss;
	int	tmcu;
	int	blockno;
	int*	acc;
	int	k;
	int	x;

	h = is->jstate;
	Ns = h->Ns;
	Al = h->Al;
	if(Ns != 1)
		imgerror(is, ERRimgbad);		//  illegal Ns>1 in progressive AC scan
	Ss = h->Ss;
	Se = h->Se;
	H = h->comp[comp].H;
	V = h->comp[comp].V;

	nacross = h->nacross*H;
	ndown = h->ndown*V;
	q = 8*h->Hmax/H;
	nhor = (h->X + q - 1)/q;
	q = 8*h->Vmax/V;
	nver = (h->Y + q - 1)/q;

	// initialize data structures
	h->cnt = 0;
	h->sr = 0;
	Ta = h->scomp[0].tac;

	ri = h->ri;

	eobrun = 0;
	acht = h->acht[Ta];
	qt = h->qt[h->comp[comp].Tq];
	nmcu = nacross*ndown;
	mcu = 0;
	for(y = 0; y < nver; y++) {
		for(x = 0; x < nhor; x++) {
			// Figure G-3
			if(eobrun > 0) {
				--eobrun;
				continue;
			}

			// arrange blockno to be in same sequence as
			// original scan calculation.
			tmcu = x/H + (nacross/H)*(y/V);
			blockno = tmcu*H*V + H*(y%V) + x%H;
			acc = h->accoeff[comp][blockno];
			k = Ss;
			while(1) {
				rs = jdecode(is, acht);
				nibbles(rs, &rrrr, &ssss);
				if(ssss == 0) {
					if(rrrr < 15) {
						eobrun = 0;
						if(rrrr > 0)
							eobrun = jreceiveEOB(is, rrrr) - 1;
						break;
					}
					k += 16;
				}
				else {
					k += rrrr;
					z = jreceive(is, ssss);
					acc[k] = z*qt[k] << Al;
					if(k == Se)
						break;
					k++;
				}
			}
		}

		// process restart marker, if present
		mcu++;
		if(ri > 0 && mcu < nmcu && mcu%ri == 0) {
			jrestart(is, mcu);
			eobrun = 0;
		}
	}
}

static void
jprogressiveacinc(ImageSource* is, int comp)
{
	Jpegstate*	h;
	int	Ns;
	int	Ss;
	int	Se;
	int	H;
	int	V;
	int	Al;
	int	nacross;
	int	ndown;
	int	q;
	int	nhor;
	int	nver;
	int	Ta;
	int	ri;
	int	eobrun;
	int**	ac;
	Huffman*	acht;
	int*	qt;
	int	nmcu;
	int	mcu;
	int	pending;
	int	nzeros;
	int	y;
	int	k;
	int	i;
	int	rs;
	int	rrrr;
	int	ssss;
	int	tmcu;
	int	blockno;
	int*	acc;
	int	x;

	h = is->jstate;
	Ns = h->Ns;
	if(Ns != 1)
		imgerror(is, ERRimgbad);		// illegal Ns>1 in progressive AC scan
	Ss = h->Ss;
	Se = h->Se;
	H = h->comp[comp].H;
	V = h->comp[comp].V;
	Al = h->Al;

	nacross = h->nacross*H;
	ndown = h->ndown*V;
	q = 8*h->Hmax/H;
	nhor = (h->X + q - 1)/q;
	q = 8*h->Vmax/V;
	nver = (h->Y + q - 1)/q;

	// initialize data structures
	h->cnt = 0;
	h->sr = 0;
	Ta = h->scomp[0].tac;
	ri = h->ri;

	eobrun = 0;
	ac = h->accoeff[comp];
	acht = h->acht[Ta];
	qt = h->qt[h->comp[comp].Tq];
	nmcu = nacross*ndown;
	mcu = 0;
	pending = 0;
	nzeros = -1;
	for(y = 0; y < nver; y++) {
		for(x = 0; x < nhor; x++) {
			// Figure G-7

			// arrange blockno to be in same sequence as
			// original scan calculation.
			tmcu = x/H + (nacross/H)*(y/V);
			blockno = tmcu*H*V + H*(y%V) + x%H;
			acc = ac[blockno];
			if(eobrun > 0) {
				if(nzeros > 0)
					imgerror(is, ERRimgbad);		// zeros pending at block start
				for(k = Ss; k <= Se; k++)
					jincrement(is, acc, k, qt[k] << Al);
				--eobrun;
				continue;
			}

			k = Ss;
			while(k <= Se) {
				if(nzeros >= 0) {
					if(acc[k] != 0)
						jincrement(is, acc, k, qt[k] << Al);
					else if(nzeros-- == 0)
						acc[k] = pending;
					k++;
					continue;
				}
				rs = jdecode(is, acht);
				nibbles(rs, &rrrr, &ssss);
				if(ssss == 0) {
					if(rrrr < 15) {
						eobrun = 0;
						if(rrrr > 0)
							eobrun = jreceiveEOB(is, rrrr) - 1;
						while(k <= Se) {
							jincrement(is, acc, k, qt[k] << Al);
							k++;
						}
						break;
					}
					for(i = 0; i < 16; k++) {
						jincrement(is, acc, k, qt[k] << Al);
						if(acc[k] == 0)
							i++;
					}
					continue;
				}
				else if(ssss != 1)
					imgerror(is, ERRimgbad);		// Jpeg ssss!=1 in progressive increment
				nzeros = rrrr;
				pending = jreceivebit(is);
				if(pending == 0)
					pending = -1;
				pending *= qt[k] << Al;
			}
		}

		// process restart marker, if present
		mcu++;
		if(ri > 0 && mcu < nmcu && mcu%ri == 0) {
			jrestart(is, mcu);
			eobrun = 0;
			nzeros = -1;
		}
	}
}

static void
jincrement(ImageSource* is, int* acc, int k, int Pt)
{
	int	b;

	if(acc[k] == 0)
		return ;
	b = jreceivebit(is);
	if(b != 0)
		if(acc[k] < 0)
			acc[k] -= Pt;
		else
			acc[k] += Pt;
}

// Fills in pixels (x,y) for x = minx=8*(mcu%nacross), minx+1, ..., minx+7 (or h.X-1, if less)
// and for y = miny=8*(mcu/nacross), miny+1, ..., miny+7 (or h.Y-1, if less)
static void
colormap1(Jpegstate* h, uchar* pic, int* data, int mcu, int nacross)
{
	int	minx;
	int	dx;
	int	miny;
	int	dy;
	int	pici;
	int	k;
	int	y;
	int	x;

	minx = 8*(mcu%nacross);
	dx = 8;
	if(minx + dx > h->X)
		dx = h->X - minx;
	miny = 8*(mcu/nacross);
	dy = 8;
	if(miny + dy > h->Y)
		dy = h->Y - miny;
	pici = miny*h->X + minx;
	k = 0;
	for(y = 0; y < dy; y++) {
		for(x = 0; x < dx; x++)
			pic[pici + x] = clampb[(data[k + x] + 128) + CLAMPBOFF];
		pici += h->X;
		k += 8;
	}
}

// Fills in same pixels as colormap1
static void
colormapall1(Jpegstate* h, uchar** chans, int* data0, int* data1, int* data2, int mcu, int nacross)
{
	uchar*	rpic;
	uchar*	gpic;
	uchar*	bpic;
	int	minx;
	int	dx;
	int	miny;
	int	dy;
	int	pici;
	int	k;
	int	y;
	int	Y;
	int	Cb;
	int	Cr;
	int	r;
	int	g;
	int	b;
	int	x;

	rpic = chans[0];
	gpic = chans[1];
	bpic = chans[2];
	minx = 8*(mcu%nacross);
	dx = 8;
	if(minx + dx > h->X)
		dx = h->X - minx;
	miny = 8*(mcu/nacross);
	dy = 8;
	if(miny + dy > h->Y)
		dy = h->Y - miny;
	pici = miny*h->X + minx;
	k = 0;
	for(y = 0; y < dy; y++) {
		for(x = 0; x < dx; x++) {
			if(jpegcolorspace == CYCbCr) {
				rpic[pici + x] = clampb[data0[k + x] + 128 + CLAMPBOFF];
				gpic[pici + x] = clampb[data1[k + x] + 128 + CLAMPBOFF];
				bpic[pici + x] = clampb[data2[k + x] + 128 + CLAMPBOFF];
			}
			else {
				Y = (data0[k + x] + 128) << 11;
				Cb = data1[k + x];
				Cr = data2[k + x];
				r = Y + jc1*Cr;
				g = Y - jc2*Cb - jc3*Cr;
				b = Y + jc4*Cb;
				rpic[pici + x] = clampb[(r >> 11) + CLAMPBOFF];
				gpic[pici + x] = clampb[(g >> 11) + CLAMPBOFF];
				bpic[pici + x] = clampb[(b >> 11) + CLAMPBOFF];
			}
		}
		pici += h->X;
		k += 8;
	}
}

// Fills in pixels (x,y) for x = minx=8*Hmax*(mcu%nacross), minx+1, ..., minx+8*Hmax-1 (or h.X-1, if less)
// and for y = miny=8*Vmax*(mcu/nacross), miny+1, ..., miny+8*Vmax-1 (or h.Y-1, if less)
static void
colormap(Jpegstate* h, uchar** chans, int** data0, int** data1, int** data2, int mcu, int nacross, int Hmax, int Vmax, int* H, int* V)
{
	uchar*	rpic;
	uchar*	gpic;
	uchar*	bpic;
	int	minx;
	int	dx;
	int	miny;
	int	dy;
	int	pici;
	int	H0;
	int	H1;
	int	H2;
	int	y;
	int	Y;
	int	Cb;
	int	Cr;
	int	r;
	int	g;
	int	b;
	int	t;
	int	b0;
	int	y0;
	int	b1;
	int	y1;
	int	b2;
	int	y2;
	int	x0;
	int	x1;
	int	x2;
	int	x;

	rpic = chans[0];
	gpic = chans[1];
	bpic = chans[2];
	minx = 8*Hmax*(mcu%nacross);
	dx = 8*Hmax;
	if(minx + dx > h->X)
		dx = h->X - minx;
	miny = 8*Vmax*(mcu/nacross);
	dy = 8*Vmax;
	if(miny + dy > h->Y)
		dy = h->Y - miny;
	pici = miny*h->X + minx;
	H0 = H[0];
	H1 = H[1];
	H2 = H[2];
	if(dbgimg > 2)
		trace("colormap, minx=%d, miny=%d, dx=%d, dy=%d, pici=%d, H0=%d, H1=%d, H2=%d\n",
			minx, miny, dx, dy, pici, H0, H1, H2);
	for(y = 0; y < dy; y++) {
		t = y*V[0];
		b0 = H0*(t/(8*Vmax));
		y0 = 8*((t/Vmax)&7);
		t = y*V[1];
		b1 = H1*(t/(8*Vmax));
		y1 = 8*((t/Vmax)&7);
		t = y*V[2];
		b2 = H2*(t/(8*Vmax));
		y2 = 8*((t/Vmax)&7);
		x0 = 0;
		x1 = 0;
		x2 = 0;
		for(x = 0; x < dx; x++) {
			if(jpegcolorspace == CYCbCr) {
				rpic[pici + x] = clampb[data0[b0][y0 + x0++*H0/Hmax] + 128 + CLAMPBOFF];
				gpic[pici + x] = clampb[data1[b1][y1 + x1++*H1/Hmax] + 128 + CLAMPBOFF];
				bpic[pici + x] = clampb[data2[b2][y2 + x2++*H2/Hmax] + 128 + CLAMPBOFF];
			}
			else { // RGB
				Y = (data0[b0][y0 + x0++*H0/Hmax] + 128) << 11;
				Cb = data1[b1][y1 + x1++*H1/Hmax];
				Cr = data2[b2][y2 + x2++*H2/Hmax];
				r = Y + jc1*Cr;
				g = Y - jc2*Cb - jc3*Cr;
				b = Y + jc4*Cb;
				rpic[pici + x] = clampb[(r >> 11) + CLAMPBOFF];
				gpic[pici + x] = clampb[(g >> 11) + CLAMPBOFF];
				bpic[pici + x] = clampb[(b >> 11) + CLAMPBOFF];
			}
			if(x0*H0/Hmax >= 8) {
				x0 = 0;
				b0++;
			}
			if(x1*H1/Hmax >= 8) {
				x1 = 0;
				b1++;
			}
			if(x2*H2/Hmax >= 8) {
				x2 = 0;
				b2++;
			}
		}
		pici += h->X;
	}
}

// decode next 8-bit value from entropy-coded input.  chart F-26
static int
jdecode(ImageSource* is, Huffman* t)
{
	Jpegstate*	h;
	int*	maxcode;
	int	code;
	int	v;
	int	cnt;
	int	m;
	int	sr;
	int	i;

	h = is->jstate;
	maxcode = t->maxcode;
	if(h->cnt < 8)
		jnextbyte(is);
	// fast lookup
	code = (h->sr >> (h->cnt - 8))&255;
	v = t->value[code];
	if(v >= 0) {
		h->cnt -= t->shift[code];
		return v;
	}

	h->cnt -= 8;
	if(h->cnt == 0)
		jnextbyte(is);
	h->cnt--;
	cnt = h->cnt;
	m = 1 << cnt;
	sr = h->sr;
	code <<= 1;
	for(i = 9; ; i++) {
		if(sr&m)
			code |= 1;
		if(code <= maxcode[i])
			break;
		code <<= 1;
		m >>= 1;
		if(m == 0) {
			sr = jnextbyte(is);
			m = 0x80;
			cnt = 8;
		}
		cnt--;
	}
	h->cnt = cnt;
	return t->val[t->valptr[i] + (code - t->mincode[i])];
}

// load next byte of input
static int
jnextbyte(ImageSource* is)
{
	int	b;
	Jpegstate*	h;
	int	b2;

	b = getc(is);
	if(b == 0xFF) {
		b2 = getc(is);
		if(b2 != 0) {
			if(b2 == DNL)
				imgerror(is, ERRimgbad);		// DNL marker unimplemented
			ungetc2(is);
		}
	}
	h = is->jstate;
	h->cnt += 8;
	h->sr = (h->sr << 8)|b;
	return b;
}

// like jnextbyte, but look for marker too
static int
jnextborm(ImageSource* is)
{
	int	b;
	Jpegstate*	h;

	b = getc(is);

	if(b == 0xFF)
		return b;
	h = is->jstate;
	h->cnt += 8;
	h->sr = (h->sr << 8)|b;
	return b;
}

// return next s bits of input, MSB first, and level shift it
static int
jreceive(ImageSource* is, int s)
{
	Jpegstate*	h;
	int	v;
	int	m;

	h = is->jstate;
	while(h->cnt < s)
		jnextbyte(is);
	h->cnt -= s;
	v = h->sr >> h->cnt;
	m = (1 << s);
	v &= m - 1;
	// level shift
	if(v < (m >> 1))
		v += ~(m - 1) + 1;
	return v;
}

// return next s bits of input, decode as EOB
static int
jreceiveEOB(ImageSource* is, int s)
{
	Jpegstate*	h;
	int	v;
	int	m;

	h = is->jstate;
	while(h->cnt < s)
		jnextbyte(is);
	h->cnt -= s;
	v = h->sr >> h->cnt;
	m = (1 << s);
	v &= m - 1;
	// level shift
	v += m;
	return v;
}

// return next bit of input
static int
jreceivebit(ImageSource* is)
{
	Jpegstate*	h;

	h = is->jstate;
	if(h->cnt < 1)
		jnextbyte(is);
	h->cnt--;
	return (h->sr >> h->cnt)&1;
}

static void
nibbles(int c, int* p1, int* p2)
{
	*p1 = c >> 4;
	*p2 = c & 15;
}

// Scaled integer implementation.
// inverse two dimensional DCT, Chen-Wang algorithm
// (IEEE ASSP-32, pp. 803-816, Aug. 1984)
// 32-bit integer arithmetic (8 bit coefficients)
// 11 mults, 29 adds per DCT
//
// coefficients extended to 12 bit for IEEE1180-1990
// compliance

static void
idct(int* b)
{
	int	y;
	int	x;
	int	v;
	int	eighty;
	int	x0;
	int	x1;
	int	x2;
	int	x3;
	int	x4;
	int	x5;
	int	x6;
	int	x7;
	int	x8;

	// transform horizontally
	for(y = 0; y < 8; y++) {
		eighty = y << 3;
		// if all non-DC components are zero, just propagate the DC term
		if(b[eighty + 1] == 0)
		if(b[eighty + 2] == 0 && b[eighty + 3] == 0)
		if(b[eighty + 4] == 0 && b[eighty + 5] == 0)
		if(b[eighty + 6] == 0 && b[eighty + 7] == 0) {
			v = b[eighty] << 3;
			b[eighty + 0] = v;
			b[eighty + 1] = v;
			b[eighty + 2] = v;
			b[eighty + 3] = v;
			b[eighty + 4] = v;
			b[eighty + 5] = v;
			b[eighty + 6] = v;
			b[eighty + 7] = v;
			continue;
		}
		// prescale
		x0 = (b[eighty + 0] << 11) + 128;
		x1 = b[eighty + 4] << 11;
		x2 = b[eighty + 6];
		x3 = b[eighty + 2];
		x4 = b[eighty + 1];
		x5 = b[eighty + 7];
		x6 = b[eighty + 5];
		x7 = b[eighty + 3];
		// first stage
		x8 = W7*(x4 + x5);
		x4 = x8 + W1mW7*x4;
		x5 = x8 - W1pW7*x5;
		x8 = W3*(x6 + x7);
		x6 = x8 - W3mW5*x6;
		x7 = x8 - W3pW5*x7;
		// second stage
		x8 = x0 + x1;
		x0 -= x1;
		x1 = W6*(x3 + x2);
		x2 = x1 - W2pW6*x2;
		x3 = x1 + W2mW6*x3;
		x1 = x4 + x6;
		x4 -= x6;
		x6 = x5 + x7;
		x5 -= x7;
		// third stage
		x7 = x8 + x3;
		x8 -= x3;
		x3 = x0 + x2;
		x0 -= x2;
		x2 = (R2*(x4 + x5) + 128) >> 8;
		x4 = (R2*(x4 - x5) + 128) >> 8;
		// fourth stage
		b[eighty + 0] = (x7 + x1) >> 8;
		b[eighty + 1] = (x3 + x2) >> 8;
		b[eighty + 2] = (x0 + x4) >> 8;
		b[eighty + 3] = (x8 + x6) >> 8;
		b[eighty + 4] = (x8 - x6) >> 8;
		b[eighty + 5] = (x0 - x4) >> 8;
		b[eighty + 6] = (x3 - x2) >> 8;
		b[eighty + 7] = (x7 - x1) >> 8;
	}
	// transform vertically
	for(x = 0; x < 8; x++) {
		// if all non-DC components are zero, just propagate the DC term
		if(b[x + 8*1] == 0)
		if(b[x + 8*2] == 0 && b[x + 8*3] == 0)
		if(b[x + 8*4] == 0 && b[x + 8*5] == 0)
		if(b[x + 8*6] == 0 && b[x + 8*7] == 0) {
			v = (b[x + 8*0] + 32) >> 6;
			b[x + 8*0] = v;
			b[x + 8*1] = v;
			b[x + 8*2] = v;
			b[x + 8*3] = v;
			b[x + 8*4] = v;
			b[x + 8*5] = v;
			b[x + 8*6] = v;
			b[x + 8*7] = v;
			continue;
		}
		// prescale
		x0 = (b[x + 8*0] << 8) + 8192;
		x1 = b[x + 8*4] << 8;
		x2 = b[x + 8*6];
		x3 = b[x + 8*2];
		x4 = b[x + 8*1];
		x5 = b[x + 8*7];
		x6 = b[x + 8*5];
		x7 = b[x + 8*3];
		// first stage
		x8 = W7*(x4 + x5) + 4;
		x4 = (x8 + W1mW7*x4) >> 3;
		x5 = (x8 - W1pW7*x5) >> 3;
		x8 = W3*(x6 + x7) + 4;
		x6 = (x8 - W3mW5*x6) >> 3;
		x7 = (x8 - W3pW5*x7) >> 3;
		// second stage
		x8 = x0 + x1;
		x0 -= x1;
		x1 = W6*(x3 + x2) + 4;
		x2 = (x1 - W2pW6*x2) >> 3;
		x3 = (x1 + W2mW6*x3) >> 3;
		x1 = x4 + x6;
		x4 -= x6;
		x6 = x5 + x7;
		x5 -= x7;
		// third stage
		x7 = x8 + x3;
		x8 -= x3;
		x3 = x0 + x2;
		x0 -= x2;
		x2 = (R2*(x4 + x5) + 128) >> 8;
		x4 = (R2*(x4 - x5) + 128) >> 8;
		// fourth stage
		b[x + 8*0] = (x7 + x1) >> 14;
		b[x + 8*1] = (x3 + x2) >> 14;
		b[x + 8*2] = (x0 + x4) >> 14;
		b[x + 8*3] = (x8 + x6) >> 14;
		b[x + 8*4] = (x8 - x6) >> 14;
		b[x + 8*5] = (x0 - x4) >> 14;
		b[x + 8*6] = (x3 - x2) >> 14;
		b[x + 8*7] = (x7 - x1) >> 14;
	}
}

// Remap colors and dither

int
closest_rgbpix(int r, int g, int b)
{
	int	pix;
	int	dr;
	int	dg;
	int	db;
	int	d;
	int	bestdist;
	int	i;

	pix = closestrgb[((r >> 4) << 8) + ((g >> 4) << 4) + (b >> 4)];
	// If white is the closest but original r,g,b wasn't white,
	// look for another color, because web page designer probably
	// cares more about contrast than actual color
	if(pix == 255 && !(r == 255 && g == 255 && b == 255)) {
		bestdist = 1000000;
		for(i = 1; i < 256; i++) {
			dr = r - rgbvmap_r[i];
			dg = g - rgbvmap_g[i];
			db = b - rgbvmap_b[i];
			d = dr*dr + dg*dg + db*db;
			if(d < bestdist) {
				bestdist = d;
				pix = i;
			}
		}
	}
	return pix;
}

static void
init_tabs(void)
{
	int	j;
	int	t;

	for(j = 0; j < CLAMPNOFF; j++) {
		clampn_b[j] = 0;
		clampn_g[j] = 0;
		clampn_r[j] = 0;
	}
	for(j = 0; j < 256; j++) {
		t = j >> 4;
		clampn_b[CLAMPNOFF + j] = t;
		clampn_g[CLAMPNOFF + j] = t << 4;
		clampn_r[CLAMPNOFF + j] = t << 8;
	}
	for(j = 0; j < CLAMPNOFF; j++) {
		clampn_b[CLAMPNOFF + 256 + j] = 15;
		clampn_g[CLAMPNOFF + 256 + j] = 240;
		clampn_r[CLAMPNOFF + 256 + j] = 3840;
	}
	for(j = 0; j < CLAMPBOFF; j++)
		clampb[j] = (uchar)0;
	for(j = 0; j < 256; j++)
		clampb[CLAMPBOFF + j] = j;
	for(j = 0; j < CLAMPBOFF; j++)
		clampb[CLAMPBOFF + 256 + j] = 255;
}

// Remap pixels of pic[] into the closest colors in the rgbv map,
// and do error diffusion of the result.
// pic is a one-channel image whose rgb values are given by looking
// up values in cmap.
static void
remap1(uchar* pic, int dx, int dy, uchar* cmap, int cmaplen)
{
	int	i;
	int	j;
	int*	ered;
	int*	egrn;
	int*	eblu;
	int	p;
	int	y;
	int	x1;
	int	in;
	int	r;
	int	g;
	int	b;
	int	col;
	int	t;
	int	er;
	int	eg;
	int	eb;
	int	x;

	if(dbgimg)
		trace("remap1, pic len %d, dx=%d, dy=%d\n", dx*dy, dx, dy);
	i = 0;
	for(j = 0; j < 256 && i < cmaplen; j++) {
		cmap_r[j] = cmap[i++];
		cmap_g[j] = cmap[i++];
		cmap_b[j] = cmap[i++];
	}
	// in case input has bad indices
	for(; j < 256; j++) {
		cmap_r[j] = 0;
		cmap_g[j] = 0;
		cmap_b[j] = 0;
	}
	// modified floyd steinberg, coefficients (1 0) 3/16, (0, 1) 3/16, (1, 1) 7/16
	ered = (int*)emallocz((dx + 1) * sizeof(int));
	egrn = (int*)emallocz((dx + 1) * sizeof(int));
	eblu = (int*)emallocz((dx + 1) * sizeof(int));
	p = 0;
	for(y = 0; y < dy; y++) {
		er = 0;
		eg = 0;
		eb = 0;
		x = 0;
		while(x < dx) {
			x1 = x + 1;
			in = pic[p];
			r = cmap_r[in] + ered[x];
			g = cmap_g[in] + egrn[x];
			b = cmap_b[in] + eblu[x];
			col = (pic[p++] = closestrgb[clampn_r[r + CLAMPNOFF]
						+ clampn_g[g + CLAMPNOFF]
						+ clampn_b[b + CLAMPNOFF]]);

			r -= rgbvmap_r[col];
			t = (3*r) >> 4;
			ered[x] = t + er;
			ered[x1] += t;
			er = r - 3*t;

			g -= rgbvmap_g[col];
			t = (3*g) >> 4;
			egrn[x] = t + eg;
			egrn[x1] += t;
			eg = g - 3*t;

			b -= rgbvmap_b[col];
			t = (3*b) >> 4;
			eblu[x] = t + eb;
			eblu[x1] += t;
			eb = b - 3*t;

			x = x1;
		}
	}
}

// Remap pixels of pic[] into the closest greyscale colors in the rgbv map,
// and do error diffusion of the result.
// pic is a one-channel greyscale image.
static void
remapgrey(uchar* pic, int dx, int dy)
{
	int*	e;
	int	p;
	int	y;
	int	x1;
	int	b;
	int	b1;
	int	col;
	int	t;
	int	eb;
	int	x;

	if(dbgimg)
		trace("remapgrey, pic len %d, dx=%d, dy=%d\n", dx*dy, dx, dy);
	// modified floyd steinberg, coefficients (1 0) 3/16, (0, 1) 3/16, (1, 1) 7/16
	e = (int*)emallocz((dx + 1) * sizeof(int));
	p = 0;
	for(y = 0; y < dy; y++) {
		eb = 0;
		x = 0;
		while(x < dx) {
			x1 = x + 1;
			b = pic[p] + e[x];
			b1 = clampn_b[b + CLAMPNOFF];
			col = 17*b1;
			pic[p++] = col;

			b -= rgbvmap_b[col];
			t = (3*b) >> 4;
			e[x] = t + eb;
			e[x1] += t;
			eb = b - 3*t;
			x = x1;
		}
	}
}

// Remap pixels of chans into the closest colors in the rgbv map,
// and do error diffusion of the result.
// chans is a 3-channel image whose channels are either (y,cb,cr) or
// (r,g,b), depending on whether jpegcolorspace is CYCbCr or CRGB.
// Variable names use r,g,b (historical).
// (right now we only call this routine from jpeg code).
static void
remaprgb(uchar** chans, int dx, int dy)
{
	uchar*	rpic;
	uchar*	gpic;
	uchar*	bpic;
	uchar*	pic;
	int*	ered;
	int*	egrn;
	int*	eblu;
	uchar*	closest;
	int*	map0;
	int*	map1;
	int*	map2;
	int	p;
	int	y;
	int	x1;
	int	r;
	int	g;
	int	b;
	int	col;
	int	t;
	int	er;
	int	eg;
	int	eb;
	int	x;

	if(dbgimg)
		trace("remaprgb, pic len %d, dx=%d, dy=%d\n", dx*dy, dx, dy);
	rpic = chans[0];
	gpic = chans[1];
	bpic = chans[2];
	pic = chans[0];
	// modified floyd steinberg, coefficients (1 0) 3/16, (0, 1) 3/16, (1, 1) 7/16
	ered = (int*)emallocz((dx + 1) * sizeof(int));
	egrn = (int*)emallocz((dx + 1) * sizeof(int));
	eblu = (int*)emallocz((dx + 1) * sizeof(int));
	if(jpegcolorspace == CRGB) {
		closest = closestrgb;
		map0 = rgbvmap_r;
		map1 = rgbvmap_g;
		map2 = rgbvmap_b;
	}
	else {
		closest = closestycbcr;
		map0 = rgbvmap_y;
		map1 = rgbvmap_cb;
		map2 = rgbvmap_cr;
	}
	p = 0;
	for(y = 0; y < dy; y++) {
		er = 0;
		eg = 0;
		eb = 0;
		x = 0;
		while(x < dx) {
			x1 = x + 1;
			r = rpic[p] + ered[x];
			g = gpic[p] + egrn[x];
			b = bpic[p] + eblu[x];
			// Errors can be uncorrectable if converting from YCbCr,
			// since we can't guarantee that an extremal value of one of
			// the components selects a color with an extremal value.
			// If we don't, the errors accumulate without bound.  This
			// doesn't happen in RGB because the closest table can guarantee
			// a color on the edge of the gamut, producing a zero error in
			// that component.  For the rotation YCbCr space, there may be
			// no color that can guarantee zero error at the edge.
			// Therefore we must clamp explicitly rather than by assuming
			// an upper error bound of CLAMPOFF.  The performance difference
			// is miniscule anyway.
			if(r < 0)
				r = 0;
			else if(r > 255)
				r = 255;
			if(g < 0)
				g = 0;
			else if(g > 255)
				g = 255;
			if(b < 0)
				b = 0;
			else if(b > 255)
				b = 255;
			col = (pic[p++] = closest[(b >> 4) + 16*((g >> 4) + 16*(r >> 4))]);

			r -= map0[col];
			t = (3*r) >> 4;
			ered[x] = t + er;
			ered[x1] += t;
			er = r - 3*t;

			g -= map1[col];
			t = (3*g) >> 4;
			egrn[x] = t + eg;
			egrn[x1] += t;
			eg = g - 3*t;

			b -= map2[col];
			t = (3*b) >> 4;
			eblu[x] = t + eb;
			eblu[x1] += t;
			eb = b - 3*t;

			x = x1;
		}
	}
}

// Given src array, representing sw*sh pixel values, resample them into
// the returned array, with dimensions dw*dh.
//
// Quick and dirty resampling: just interpolate.
// This lets us resample arrays of pixels indices (e.g., result of gif decoding).
// The filter-based resampling methods need conversion to rgb or grayscale.
// Also, although the results won't look good, people really shouldn't be
// asking the browser to resample except for special purposes (like the common
// case of resizing a 1x1 image to make a spacer).
static uchar*
resample(uchar* src, int sw, int sh, int dw, int dh)
{
	double	xfac;
	double	yfac;
	int	totpix;
	uchar*	dst;
	int	dindex;
	int*	sindices;
	int	dx;
	int	dy;
	int	sx;
	int	sy;
	int	soffset;

	if(dbgev)
		logtime("IMAGE_RESAMPLE_START", 0);
	if(src == nil || sw == 0 || sh == 0 || dw == 0 || dh == 0)
		return src;
	xfac = (double)sw/(double)dw;
	yfac = (double)sh/(double)dh;
	totpix = dw*dh;
	dst = (uchar*)emalloc(totpix);
	dindex = 0;

	// precompute index in src row corresponding to each index in dst row
	sindices = (int*)emalloc(dw * sizeof(int));
	for(dx = 0; dx < dw; dx++) {
		sx = (int)((xfac*(double)dx)+.5);
		if(sx >= sw)
			sx = sw - 1;
		sindices[dx] = sx;
	}
	for(dy = 0; dy < dh; dy++) {
		sy = (int)((yfac*(double)dy)+.5);
		if(sy >= sh)
			sy = sh - 1;
		soffset = sy*sw;
		for(dx = 0; dx < dw; dx++)
			dst[dindex++] = src[soffset + sindices[dx]];
	}
	if(dbgev)
		logtime("IMAGE_RESAMPLE_END", 0);
	return dst;
}