1 #include <u.h>
2 #include <libc.h>
3 #include <bio.h>
4 #include <draw.h>
5 #include "imagefile.h"
6 
7 enum {
8 	/* Constants, all preceded by byte 0xFF */
9 	SOF	=0xC0,	/* Start of Frame */
10 	SOF2=0xC2,	/* Start of Frame; progressive Huffman */
11 	JPG	=0xC8,	/* Reserved for JPEG extensions */
12 	DHT	=0xC4,	/* Define Huffman Tables */
13 	DAC	=0xCC,	/* Arithmetic coding conditioning */
14 	RST	=0xD0,	/* Restart interval termination */
15 	RST7	=0xD7,	/* Restart interval termination (highest value) */
16 	SOI	=0xD8,	/* Start of Image */
17 	EOI	=0xD9,	/* End of Image */
18 	SOS	=0xDA,	/* Start of Scan */
19 	DQT	=0xDB,	/* Define quantization tables */
20 	DNL	=0xDC,	/* Define number of lines */
21 	DRI	=0xDD,	/* Define restart interval */
22 	DHP	=0xDE,	/* Define hierarchical progression */
23 	EXP	=0xDF,	/* Expand reference components */
24 	APPn	=0xE0,	/* Reserved for application segments */
25 	JPGn	=0xF0,	/* Reserved for JPEG extensions */
26 	COM	=0xFE,	/* Comment */
27 
28 	CLAMPOFF	= 300,
29 	NCLAMP		= CLAMPOFF+700
30 };
31 
32 typedef struct Framecomp Framecomp;
33 typedef struct Header Header;
34 typedef struct Huffman Huffman;
35 
36 struct Framecomp	/* Frame component specifier from SOF marker */
37 {
38 	int	C;
39 	int	H;
40 	int	V;
41 	int	Tq;
42 };
43 
44 struct Huffman
45 {
46 	int	*size;	/* malloc'ed */
47 	int	*code;	/* malloc'ed */
48 	int	*val;		/* malloc'ed */
49 	int	mincode[17];
50 	int	maxcode[17];
51 	int	valptr[17];
52 	/* fast lookup */
53 	int	value[256];
54 	int	shift[256];
55 };
56 
57 
58 struct Header
59 {
60 	Biobuf	*fd;
61 	char		err[256];
62 	jmp_buf	errlab;
63 	/* variables in i/o routines */
64 	int		sr;	/* shift register, right aligned */
65 	int		cnt;	/* # bits in right part of sr */
66 	uchar	*buf;
67 	int		nbuf;
68 	int		peek;
69 
70 	int		Nf;
71 
72 	Framecomp	comp[3];
73 	uchar	mode;
74 	int		X;
75 	int		Y;
76 	int		qt[4][64];		/* quantization tables */
77 	Huffman	dcht[4];
78 	Huffman	acht[4];
79 	int		**data[3];
80 	int		ndata[3];
81 
82 	uchar	*sf;	/* start of frame; do better later */
83 	uchar	*ss;	/* start of scan; do better later */
84 	int		ri;	/* restart interval */
85 
86 	/* progressive scan */
87 	Rawimage *image;
88 	Rawimage **array;
89 	int		*dccoeff[3];
90 	int		**accoeff[3];	/* only need 8 bits plus quantization */
91 	int		naccoeff[3];
92 	int		nblock[3];
93 	int		nacross;
94 	int		ndown;
95 	int		Hmax;
96 	int		Vmax;
97 };
98 
99 static	uchar	clamp[NCLAMP];
100 
101 static	Rawimage	*readslave(Header*, int);
102 static	int			readsegment(Header*, int*);
103 static	void			quanttables(Header*, uchar*, int);
104 static	void			huffmantables(Header*, uchar*, int);
105 static	void			soiheader(Header*);
106 static	int			nextbyte(Header*, int);
107 static	int			int2(uchar*, int);
108 static	void			nibbles(int, int*, int*);
109 static	int			receive(Header*, int);
110 static	int			receiveEOB(Header*, int);
111 static	int			receivebit(Header*);
112 static	void			restart(Header*, int);
113 static	int			decode(Header*, Huffman*);
114 static	Rawimage*	baselinescan(Header*, int);
115 static	void			progressivescan(Header*, int);
116 static	Rawimage*	progressiveIDCT(Header*, int);
117 static	void			idct(int*);
118 static	void			colormap1(Header*, int, Rawimage*, int*, int, int);
119 static	void			colormapall1(Header*, int, Rawimage*, int*, int*, int*, int, int);
120 static	void			colormap(Header*, int, Rawimage*, int**, int**, int**, int, int, int, int, int*, int*);
121 static	void			jpgerror(Header*, char*, ...);
122 
123 static	char		readerr[] = "ReadJPG: read error: %r";
124 static	char		memerr[] = "ReadJPG: malloc failed: %r";
125 
126 static	int zig[64] = {
127 	0, 1, 8, 16, 9, 2, 3, 10, 17, /* 0-7 */
128 	24, 32, 25, 18, 11, 4, 5, /* 8-15 */
129 	12, 19, 26, 33, 40, 48, 41, 34, /* 16-23 */
130 	27, 20, 13, 6, 7, 14, 21, 28, /* 24-31 */
131 	35, 42, 49, 56, 57, 50, 43, 36, /* 32-39 */
132 	29, 22, 15, 23, 30, 37, 44, 51, /* 40-47 */
133 	58, 59, 52, 45, 38, 31, 39, 46, /* 48-55 */
134 	53, 60, 61, 54, 47, 55, 62, 63 /* 56-63 */
135 };
136 
137 static
138 void
139 jpginit(void)
140 {
141 	int k;
142 	static int inited;
143 
144 	if(inited)
145 		return;
146 	inited = 1;
147 	for(k=0; k<CLAMPOFF; k++)
148 		clamp[k] = 0;
149 	for(; k<CLAMPOFF+256; k++)
150 		clamp[k] = k-CLAMPOFF;
151 	for(; k<NCLAMP; k++)
152 		clamp[k] = 255;
153 }
154 
155 static
156 void*
157 jpgmalloc(Header *h, int n, int clear)
158 {
159 	void *p;
160 
161 	p = malloc(n);
162 	if(p == nil)
163 		jpgerror(h, memerr);
164 	if(clear)
165 		memset(p, 0, n);
166 	return p;
167 }
168 
169 static
170 void
171 clear(void *pp)
172 {
173 	void **p = (void**)pp;
174 
175 	if(*p){
176 		free(*p);
177 		*p = nil;
178 	}
179 }
180 
181 static
182 void
183 jpgfreeall(Header *h, int freeimage)
184 {
185 	int i, j;
186 
187 	clear(&h->buf);
188 	if(h->dccoeff[0])
189 		for(i=0; i<3; i++)
190 			clear(&h->dccoeff[i]);
191 	if(h->accoeff[0])
192 		for(i=0; i<3; i++){
193 			if(h->accoeff[i])
194 				for(j=0; j<h->naccoeff[i]; j++)
195 					clear(&h->accoeff[i][j]);
196 			clear(&h->accoeff[i]);
197 		}
198 	for(i=0; i<4; i++){
199 		clear(&h->dcht[i].size);
200 		clear(&h->acht[i].size);
201 		clear(&h->dcht[i].code);
202 		clear(&h->acht[i].code);
203 		clear(&h->dcht[i].val);
204 		clear(&h->acht[i].val);
205 	}
206 	if(h->data[0])
207 		for(i=0; i<3; i++){
208 			if(h->data[i])
209 				for(j=0; j<h->ndata[i]; j++)
210 					clear(&h->data[i][j]);
211 			clear(&h->data[i]);
212 		}
213 	if(freeimage && h->image!=nil){
214 		clear(&h->array);
215 		clear(&h->image->cmap);
216 		for(i=0; i<3; i++)
217 			clear(&h->image->chans[i]);
218 		clear(&h->image);
219 	}
220 }
221 
222 static
223 void
224 jpgerror(Header *h, char *fmt, ...)
225 {
226 	va_list arg;
227 
228 	va_start(arg, fmt);
229 	vseprint(h->err, h->err+sizeof h->err, fmt, arg);
230 	va_end(arg);
231 
232 	werrstr(h->err);
233 	jpgfreeall(h, 1);
234 	longjmp(h->errlab, 1);
235 }
236 
237 Rawimage**
238 Breadjpg(Biobuf *b, int colorspace)
239 {
240 	Rawimage *r, **array;
241 	Header *h;
242 	char buf[ERRMAX];
243 
244 	buf[0] = '\0';
245 	if(colorspace!=CYCbCr && colorspace!=CRGB){
246 		errstr(buf, sizeof buf);	/* throw it away */
247 		werrstr("ReadJPG: unknown color space");
248 		return nil;
249 	}
250 	jpginit();
251 	h = malloc(sizeof(Header));
252 	array = malloc(sizeof(Header));
253 	if(h==nil || array==nil){
254 		free(h);
255 		free(array);
256 		return nil;
257 	}
258 	h->array = array;
259 	memset(h, 0, sizeof(Header));
260 	h->fd = b;
261 	errstr(buf, sizeof buf);	/* throw it away */
262 	if(setjmp(h->errlab))
263 		r = nil;
264 	else
265 		r = readslave(h, colorspace);
266 	jpgfreeall(h, 0);
267 	free(h);
268 	array[0] = r;
269 	array[1] = nil;
270 	return array;
271 }
272 
273 Rawimage**
274 readjpg(int fd, int colorspace)
275 {
276 	Rawimage** a;
277 	Biobuf b;
278 
279 	if(Binit(&b, fd, OREAD) < 0)
280 		return nil;
281 	a = Breadjpg(&b, colorspace);
282 	Bterm(&b);
283 	return a;
284 }
285 
286 static
287 Rawimage*
288 readslave(Header *header, int colorspace)
289 {
290 	Rawimage *image;
291 	int nseg, i, H, V, m, n;
292 	uchar *b;
293 
294 	soiheader(header);
295 	nseg = 0;
296 	image = nil;
297 
298 	header->buf = jpgmalloc(header, 4096, 0);
299 	header->nbuf = 4096;
300 	while(header->err[0] == '\0'){
301 		nseg++;
302 		n = readsegment(header, &m);
303 		b = header->buf;
304 		switch(m){
305 		case -1:
306 			return image;
307 
308 		case APPn+0:
309 			if(nseg==1 && strncmp((char*)b, "JFIF", 4)==0)  /* JFIF header; check version */
310 				if(b[5]>1 || b[6]>2)
311 					sprint(header->err, "ReadJPG: can't handle JFIF version %d.%2d", b[5], b[6]);
312 			break;
313 
314 		case APPn+1: case APPn+2: case APPn+3: case APPn+4: case APPn+5:
315 		case APPn+6: case APPn+7: case APPn+8: case APPn+9: case APPn+10:
316 		case APPn+11: case APPn+12: case APPn+13: case APPn+14: case APPn+15:
317 			break;
318 
319 		case DQT:
320 			quanttables(header, b, n);
321 			break;
322 
323 		case SOF:
324 		case SOF2:
325 			header->Y = int2(b, 1);
326 			header->X = int2(b, 3);
327 			header->Nf =b[5];
328 			for(i=0; i<header->Nf; i++){
329 				header->comp[i].C = b[6+3*i+0];
330 				nibbles(b[6+3*i+1], &H, &V);
331 				if(H<=0 || V<=0)
332 					jpgerror(header, "non-positive sampling factor (Hsamp or Vsamp)");
333 				header->comp[i].H = H;
334 				header->comp[i].V = V;
335 				header->comp[i].Tq = b[6+3*i+2];
336 			}
337 			header->mode = m;
338 			header->sf = b;
339 			break;
340 
341 		case  SOS:
342 			header->ss = b;
343 			switch(header->mode){
344 			case SOF:
345 				image = baselinescan(header, colorspace);
346 				break;
347 			case SOF2:
348 				progressivescan(header, colorspace);
349 				break;
350 			default:
351 				sprint(header->err, "unrecognized or unspecified encoding %d", header->mode);
352 				break;
353 			}
354 			break;
355 
356 		case  DHT:
357 			huffmantables(header, b, n);
358 			break;
359 
360 		case  DRI:
361 			header->ri = int2(b, 0);
362 			break;
363 
364 		case  COM:
365 			break;
366 
367 		case EOI:
368 			if(header->mode == SOF2)
369 				image = progressiveIDCT(header, colorspace);
370 			return image;
371 
372 		default:
373 			sprint(header->err, "ReadJPG: unknown marker %.2x", m);
374 			break;
375 		}
376 	}
377 	return image;
378 }
379 
380 /* readsegment is called after reading scan, which can have */
381 /* read ahead a byte.  so we must check peek here */
382 static
383 int
384 readbyte(Header *h)
385 {
386 	uchar x;
387 
388 	if(h->peek >= 0){
389 		x = h->peek;
390 		h->peek = -1;
391 	}else if(Bread(h->fd, &x, 1) != 1)
392 		jpgerror(h, readerr);
393 	return x;
394 }
395 
396 static
397 int
398 marker(Header *h)
399 {
400 	int c;
401 
402 	while((c=readbyte(h)) == 0)
403 		fprint(2, "ReadJPG: skipping zero byte at offset %lld\n", Boffset(h->fd));
404 	if(c != 0xFF)
405 		jpgerror(h, "ReadJPG: expecting marker; found 0x%x at offset %lld\n", c, Boffset(h->fd));
406 	while(c == 0xFF)
407 		c = readbyte(h);
408 	return c;
409 }
410 
411 static
412 int
413 int2(uchar *buf, int n)
414 {
415 	return (buf[n]<<8) + buf[n+1];
416 }
417 
418 static
419 void
420 nibbles(int b, int *p0, int *p1)
421 {
422 	*p0 = (b>>4) & 0xF;
423 	*p1 = b & 0xF;
424 }
425 
426 static
427 void
428 soiheader(Header *h)
429 {
430 	h->peek = -1;
431 	if(marker(h) != SOI)
432 		jpgerror(h, "ReadJPG: unrecognized marker in header");
433 	h->err[0] = '\0';
434 	h->mode = 0;
435 	h->ri = 0;
436 }
437 
438 static
439 int
440 readsegment(Header *h, int *markerp)
441 {
442 	int m, n;
443 	uchar tmp[2];
444 
445 	m = marker(h);
446 	switch(m){
447 	case EOI:
448 		*markerp = m;
449 		return 0;
450 	case 0:
451 		jpgerror(h, "ReadJPG: expecting marker; saw %.2x at offset %lld", m, Boffset(h->fd));
452 	}
453 	if(Bread(h->fd, tmp, 2) != 2)
454     Readerr:
455 		jpgerror(h, readerr);
456 	n = int2(tmp, 0);
457 	if(n < 2)
458 		goto Readerr;
459 	n -= 2;
460 	if(n > h->nbuf){
461 		free(h->buf);
462 		h->buf = jpgmalloc(h, n+1, 0); /* +1 for sentinel */
463 		h->nbuf = n;
464 	}
465 	if(Bread(h->fd, h->buf, n) != n)
466 		goto Readerr;
467 	*markerp = m;
468 	return n;
469 }
470 
471 static
472 int
473 huffmantable(Header *h, uchar *b)
474 {
475 	Huffman *t;
476 	int Tc, th, n, nsize, i, j, k, v, cnt, code, si, sr, m;
477 	int *maxcode;
478 
479 	nibbles(b[0], &Tc, &th);
480 	if(Tc > 1)
481 		jpgerror(h, "ReadJPG: unknown Huffman table class %d", Tc);
482 	if(th>3 || (h->mode==SOF && th>1))
483 		jpgerror(h, "ReadJPG: unknown Huffman table index %d", th);
484 	if(Tc == 0)
485 		t = &h->dcht[th];
486 	else
487 		t = &h->acht[th];
488 
489 	/* flow chart C-2 */
490 	nsize = 0;
491 	for(i=0; i<16; i++)
492 		nsize += b[1+i];
493 	t->size = jpgmalloc(h, (nsize+1)*sizeof(int), 1);
494 	k = 0;
495 	for(i=1; i<=16; i++){
496 		n = b[i];
497 		for(j=0; j<n; j++)
498 			t->size[k++] = i;
499 	}
500 	t->size[k] = 0;
501 
502 	/* initialize HUFFVAL */
503 	t->val = jpgmalloc(h, nsize*sizeof(int), 1);
504 	for(i=0; i<nsize; i++)
505 		t->val[i] = b[17+i];
506 
507 	/* flow chart C-3 */
508 	t->code = jpgmalloc(h, (nsize+1)*sizeof(int), 1);
509 	k = 0;
510 	code = 0;
511 	si = t->size[0];
512 	for(;;){
513 		do
514 			t->code[k++] = code++;
515 		while(t->size[k] == si);
516 		if(t->size[k] == 0)
517 			break;
518 		do{
519 			code <<= 1;
520 			si++;
521 		}while(t->size[k] != si);
522 	}
523 
524 	/* flow chart F-25 */
525 	i = 0;
526 	j = 0;
527 	for(;;){
528 		for(;;){
529 			i++;
530 			if(i > 16)
531 				goto outF25;
532 			if(b[i] != 0)
533 				break;
534 			t->maxcode[i] = -1;
535 		}
536 		t->valptr[i] = j;
537 		t->mincode[i] = t->code[j];
538 		j += b[i]-1;
539 		t->maxcode[i] = t->code[j];
540 		j++;
541 	}
542 outF25:
543 
544 	/* create byte-indexed fast path tables */
545 	maxcode = t->maxcode;
546 	/* stupid startup algorithm: just run machine for each byte value */
547 	for(v=0; v<256; ){
548 		cnt = 7;
549 		m = 1<<7;
550 		code = 0;
551 		sr = v;
552 		i = 1;
553 		for(;;i++){
554 			if(sr & m)
555 				code |= 1;
556 			if(code <= maxcode[i])
557 				break;
558 			code <<= 1;
559 			m >>= 1;
560 			if(m == 0){
561 				t->shift[v] = 0;
562 				t->value[v] = -1;
563 				goto continueBytes;
564 			}
565 			cnt--;
566 		}
567 		t->shift[v] = 8-cnt;
568 		t->value[v] = t->val[t->valptr[i]+(code-t->mincode[i])];
569 
570     continueBytes:
571 		v++;
572 	}
573 
574 	return nsize;
575 }
576 
577 static
578 void
579 huffmantables(Header *h, uchar *b, int n)
580 {
581 	int l, mt;
582 
583 	for(l=0; l<n; l+=17+mt)
584 		mt = huffmantable(h, &b[l]);
585 }
586 
587 static
588 int
589 quanttable(Header *h, uchar *b)
590 {
591 	int i, pq, tq, *q;
592 
593 	nibbles(b[0], &pq, &tq);
594 	if(pq > 1)
595 		jpgerror(h, "ReadJPG: unknown quantization table class %d", pq);
596 	if(tq > 3)
597 		jpgerror(h, "ReadJPG: unknown quantization table index %d", tq);
598 	q = h->qt[tq];
599 	for(i=0; i<64; i++){
600 		if(pq == 0)
601 			q[i] = b[1+i];
602 		else
603 			q[i] = int2(b, 1+2*i);
604 	}
605 	return 64*(1+pq);
606 }
607 
608 static
609 void
610 quanttables(Header *h, uchar *b, int n)
611 {
612 	int l, m;
613 
614 	for(l=0; l<n; l+=1+m)
615 		m = quanttable(h, &b[l]);
616 }
617 
618 static
619 Rawimage*
620 baselinescan(Header *h, int colorspace)
621 {
622 	int Ns, z, k, m, Hmax, Vmax, comp;
623 	int allHV1, nblock, ri, mcu, nacross, nmcu;
624 	Huffman *dcht, *acht;
625 	int block, t, diff, *qt;
626 	uchar *ss;
627 	Rawimage *image;
628 	int Td[3], Ta[3], H[3], V[3], DC[3];
629 	int ***data, *zz;
630 
631 	ss = h->ss;
632 	Ns = ss[0];
633 	if((Ns!=3 && Ns!=1) || Ns!=h->Nf)
634 		jpgerror(h, "ReadJPG: can't handle scan not 3 components");
635 
636 	image = jpgmalloc(h, sizeof(Rawimage), 1);
637 	h->image = image;
638 	image->r = Rect(0, 0, h->X, h->Y);
639 	image->cmap = nil;
640 	image->cmaplen = 0;
641 	image->chanlen = h->X*h->Y;
642 	image->fields = 0;
643 	image->gifflags = 0;
644 	image->gifdelay = 0;
645 	image->giftrindex = 0;
646 	if(Ns == 3)
647 		image->chandesc = colorspace;
648 	else
649 		image->chandesc = CY;
650 	image->nchans = h->Nf;
651 	for(k=0; k<h->Nf; k++)
652 		image->chans[k] = jpgmalloc(h, h->X*h->Y, 0);
653 
654 	/* compute maximum H and V */
655 	Hmax = 0;
656 	Vmax = 0;
657 	for(comp=0; comp<Ns; comp++){
658 		if(h->comp[comp].H > Hmax)
659 			Hmax = h->comp[comp].H;
660 		if(h->comp[comp].V > Vmax)
661 			Vmax = h->comp[comp].V;
662 	}
663 
664 	/* initialize data structures */
665 	allHV1 = 1;
666 	data = h->data;
667 	for(comp=0; comp<Ns; comp++){
668 		/* JPEG requires scan components to be in same order as in frame, */
669 		/* so if both have 3 we know scan is Y Cb Cr and there's no need to */
670 		/* reorder */
671 		nibbles(ss[2+2*comp], &Td[comp], &Ta[comp]);
672 		H[comp] = h->comp[comp].H;
673 		V[comp] = h->comp[comp].V;
674 		nblock = H[comp]*V[comp];
675 		if(nblock != 1)
676 			allHV1 = 0;
677 		data[comp] = jpgmalloc(h, nblock*sizeof(int*), 0);
678 		h->ndata[comp] = nblock;
679 		DC[comp] = 0;
680 		for(m=0; m<nblock; m++)
681 			data[comp][m] = jpgmalloc(h, 8*8*sizeof(int), 0);
682 	}
683 
684 	ri = h->ri;
685 
686 	h->cnt = 0;
687 	h->sr = 0;
688 	h->peek = -1;
689 	nacross = ((h->X+(8*Hmax-1))/(8*Hmax));
690 	nmcu = ((h->Y+(8*Vmax-1))/(8*Vmax))*nacross;
691 	for(mcu=0; mcu<nmcu; ){
692 		for(comp=0; comp<Ns; comp++){
693 			dcht = &h->dcht[Td[comp]];
694 			acht = &h->acht[Ta[comp]];
695 			qt = h->qt[h->comp[comp].Tq];
696 
697 			for(block=0; block<H[comp]*V[comp]; block++){
698 				/* F-22 */
699 				t = decode(h, dcht);
700 				diff = receive(h, t);
701 				DC[comp] += diff;
702 
703 				/* F-23 */
704 				zz = data[comp][block];
705 				memset(zz, 0, 8*8*sizeof(int));
706 				zz[0] = qt[0]*DC[comp];
707 				k = 1;
708 
709 				for(;;){
710 					t = decode(h, acht);
711 					if((t&0x0F) == 0){
712 						if((t&0xF0) != 0xF0)
713 							break;
714 						k += 16;
715 					}else{
716 						k += t>>4;
717 						z = receive(h, t&0xF);
718 						zz[zig[k]] = z*qt[k];
719 						if(k == 63)
720 							break;
721 						k++;
722 					}
723 				}
724 
725 				idct(zz);
726 			}
727 		}
728 
729 		/* rotate colors to RGB and assign to bytes */
730 		if(Ns == 1) /* very easy */
731 			colormap1(h, colorspace, image, data[0][0], mcu, nacross);
732 		else if(allHV1) /* fairly easy */
733 			colormapall1(h, colorspace, image, data[0][0], data[1][0], data[2][0], mcu, nacross);
734 		else /* miserable general case */
735 			colormap(h, colorspace, image, data[0], data[1], data[2], mcu, nacross, Hmax, Vmax, H, V);
736 		/* process restart marker, if present */
737 		mcu++;
738 		if(ri>0 && mcu<nmcu && mcu%ri==0){
739 			restart(h, mcu);
740 			for(comp=0; comp<Ns; comp++)
741 				DC[comp] = 0;
742 		}
743 	}
744 	return image;
745 }
746 
747 static
748 void
749 restart(Header *h, int mcu)
750 {
751 	int rest, rst, nskip;
752 
753 	rest = mcu/h->ri-1;
754 	nskip = 0;
755 	do{
756 		do{
757 			rst = nextbyte(h, 1);
758 			nskip++;
759 		}while(rst>=0 && rst!=0xFF);
760 		if(rst == 0xFF){
761 			rst = nextbyte(h, 1);
762 			nskip++;
763 		}
764 	}while(rst>=0 && (rst&~7)!=RST);
765 	if(nskip != 2)
766 		sprint(h->err, "ReadJPG: skipped %d bytes at restart %d\n", nskip-2, rest);
767 	if(rst < 0)
768 		jpgerror(h, readerr);
769 	if((rst&7) != (rest&7))
770 		jpgerror(h, "ReadJPG: expected RST%d got %d", rest&7, rst&7);
771 	h->cnt = 0;
772 	h->sr = 0;
773 }
774 
775 static
776 Rawimage*
777 progressiveIDCT(Header *h, int colorspace)
778 {
779 	int k, m, comp, block, Nf, bn;
780 	int allHV1, nblock, mcu, nmcu;
781 	int H[3], V[3], blockno[3];
782 	int *dccoeff, **accoeff;
783 	int ***data, *zz;
784 
785 	Nf = h->Nf;
786 	allHV1 = 1;
787 	data = h->data;
788 
789 	for(comp=0; comp<Nf; comp++){
790 		H[comp] = h->comp[comp].H;
791 		V[comp] = h->comp[comp].V;
792 		nblock = h->nblock[comp];
793 		if(nblock != 1)
794 			allHV1 = 0;
795 		h->ndata[comp] = nblock;
796 		data[comp] = jpgmalloc(h, nblock*sizeof(int*), 0);
797 		for(m=0; m<nblock; m++)
798 			data[comp][m] = jpgmalloc(h, 8*8*sizeof(int), 0);
799 	}
800 
801 	memset(blockno, 0, sizeof blockno);
802 	nmcu = h->nacross*h->ndown;
803 	for(mcu=0; mcu<nmcu; mcu++){
804 		for(comp=0; comp<Nf; comp++){
805 			dccoeff = h->dccoeff[comp];
806 			accoeff = h->accoeff[comp];
807 			bn = blockno[comp];
808 			for(block=0; block<h->nblock[comp]; block++){
809 				zz = data[comp][block];
810 				memset(zz, 0, 8*8*sizeof(int));
811 				zz[0] = dccoeff[bn];
812 
813 				for(k=1; k<64; k++)
814 					zz[zig[k]] = accoeff[bn][k];
815 
816 				idct(zz);
817 				bn++;
818 			}
819 			blockno[comp] = bn;
820 		}
821 
822 		/* rotate colors to RGB and assign to bytes */
823 		if(Nf == 1) /* very easy */
824 			colormap1(h, colorspace, h->image, data[0][0], mcu, h->nacross);
825 		else if(allHV1) /* fairly easy */
826 			colormapall1(h, colorspace, h->image, data[0][0], data[1][0], data[2][0], mcu, h->nacross);
827 		else /* miserable general case */
828 			colormap(h, colorspace, h->image, data[0], data[1], data[2], mcu, h->nacross, h->Hmax, h->Vmax, H, V);
829 	}
830 
831 	return h->image;
832 }
833 
834 static
835 void
836 progressiveinit(Header *h, int colorspace)
837 {
838 	int Nf, Ns, j, k, nmcu, comp;
839 	uchar *ss;
840 	Rawimage *image;
841 
842 	ss = h->ss;
843 	Ns = ss[0];
844 	Nf = h->Nf;
845 	if((Ns!=3 && Ns!=1) || Ns!=Nf)
846 		jpgerror(h, "ReadJPG: image must have 1 or 3 components");
847 
848 	image = jpgmalloc(h, sizeof(Rawimage), 1);
849 	h->image = image;
850 	image->r = Rect(0, 0, h->X, h->Y);
851 	image->cmap = nil;
852 	image->cmaplen = 0;
853 	image->chanlen = h->X*h->Y;
854 	image->fields = 0;
855 	image->gifflags = 0;
856 	image->gifdelay = 0;
857 	image->giftrindex = 0;
858 	if(Nf == 3)
859 		image->chandesc = colorspace;
860 	else
861 		image->chandesc = CY;
862 	image->nchans = h->Nf;
863 	for(k=0; k<Nf; k++){
864 		image->chans[k] = jpgmalloc(h, h->X*h->Y, 0);
865 		h->nblock[k] = h->comp[k].H*h->comp[k].V;
866 	}
867 
868 	/* compute maximum H and V */
869 	h->Hmax = 0;
870 	h->Vmax = 0;
871 	for(comp=0; comp<Nf; comp++){
872 		if(h->comp[comp].H > h->Hmax)
873 			h->Hmax = h->comp[comp].H;
874 		if(h->comp[comp].V > h->Vmax)
875 			h->Vmax = h->comp[comp].V;
876 	}
877 	h->nacross = ((h->X+(8*h->Hmax-1))/(8*h->Hmax));
878 	h->ndown = ((h->Y+(8*h->Vmax-1))/(8*h->Vmax));
879 	nmcu = h->nacross*h->ndown;
880 
881 	for(k=0; k<Nf; k++){
882 		h->dccoeff[k] = jpgmalloc(h, h->nblock[k]*nmcu * sizeof(int), 1);
883 		h->accoeff[k] = jpgmalloc(h, h->nblock[k]*nmcu * sizeof(int*), 1);
884 		h->naccoeff[k] = h->nblock[k]*nmcu;
885 		for(j=0; j<h->nblock[k]*nmcu; j++)
886 			h->accoeff[k][j] = jpgmalloc(h, 64*sizeof(int), 1);
887 	}
888 
889 }
890 
891 static
892 void
893 progressivedc(Header *h, int comp, int Ah, int Al)
894 {
895 	int Ns, z, ri, mcu,  nmcu;
896 	int block, t, diff, qt, *dc, bn;
897 	Huffman *dcht;
898 	uchar *ss;
899 	int Td[3], DC[3], blockno[3];
900 
901 	ss= h->ss;
902 	Ns = ss[0];
903 	if(Ns!=h->Nf)
904 		jpgerror(h, "ReadJPG: can't handle progressive with Nf!=Ns in DC scan");
905 
906 	/* initialize data structures */
907 	h->cnt = 0;
908 	h->sr = 0;
909 	h->peek = -1;
910 	for(comp=0; comp<Ns; comp++){
911 		/*
912 		 * JPEG requires scan components to be in same order as in frame,
913 		 * so if both have 3 we know scan is Y Cb Cr and there's no need to
914 		 * reorder
915 		 */
916 		nibbles(ss[2+2*comp], &Td[comp], &z);	/* z is ignored */
917 		DC[comp] = 0;
918 	}
919 
920 	ri = h->ri;
921 
922 	nmcu = h->nacross*h->ndown;
923 	memset(blockno, 0, sizeof blockno);
924 	for(mcu=0; mcu<nmcu; ){
925 		for(comp=0; comp<Ns; comp++){
926 			dcht = &h->dcht[Td[comp]];
927 			qt = h->qt[h->comp[comp].Tq][0];
928 			dc = h->dccoeff[comp];
929 			bn = blockno[comp];
930 
931 			for(block=0; block<h->nblock[comp]; block++){
932 				if(Ah == 0){
933 					t = decode(h, dcht);
934 					diff = receive(h, t);
935 					DC[comp] += diff;
936 					dc[bn] = qt*DC[comp]<<Al;
937 				}else
938 					dc[bn] |= qt*receivebit(h)<<Al;
939 				bn++;
940 			}
941 			blockno[comp] = bn;
942 		}
943 
944 		/* process restart marker, if present */
945 		mcu++;
946 		if(ri>0 && mcu<nmcu && mcu%ri==0){
947 			restart(h, mcu);
948 			for(comp=0; comp<Ns; comp++)
949 				DC[comp] = 0;
950 		}
951 	}
952 }
953 
954 static
955 void
956 progressiveac(Header *h, int comp, int Al)
957 {
958 	int Ns, Ss, Se, z, k, eobrun, x, y, nver, tmcu, blockno, *acc, rs;
959 	int ri, mcu, nacross, ndown, nmcu, nhor;
960 	Huffman *acht;
961 	int *qt, rrrr, ssss, q;
962 	uchar *ss;
963 	int Ta, H, V;
964 
965 	ss = h->ss;
966 	Ns = ss[0];
967 	if(Ns != 1)
968 		jpgerror(h, "ReadJPG: illegal Ns>1 in progressive AC scan");
969 	Ss = ss[1+2];
970 	Se = ss[2+2];
971 	H = h->comp[comp].H;
972 	V = h->comp[comp].V;
973 
974 	nacross = h->nacross*H;
975 	ndown = h->ndown*V;
976 	q = 8*h->Hmax/H;
977 	nhor = (h->X+q-1)/q;
978 	q = 8*h->Vmax/V;
979 	nver = (h->Y+q-1)/q;
980 
981 	/* initialize data structures */
982 	h->cnt = 0;
983 	h->sr = 0;
984 	h->peek = -1;
985 	nibbles(ss[1+1], &z, &Ta);	/* z is thrown away */
986 
987 	ri = h->ri;
988 
989 	eobrun = 0;
990 	acht = &h->acht[Ta];
991 	qt = h->qt[h->comp[comp].Tq];
992 	nmcu = nacross*ndown;
993 	mcu = 0;
994 	for(y=0; y<nver; y++){
995 		for(x=0; x<nhor; x++){
996 			/* Figure G-3  */
997 			if(eobrun > 0){
998 				--eobrun;
999 				continue;
1000 			}
1001 
1002 			/* arrange blockno to be in same sequence as original scan calculation. */
1003 			tmcu = x/H + (nacross/H)*(y/V);
1004 			blockno = tmcu*H*V + H*(y%V) + x%H;
1005 			acc = h->accoeff[comp][blockno];
1006 			k = Ss;
1007 			for(;;){
1008 				rs = decode(h, acht);
1009 				/* XXX remove rrrr ssss as in baselinescan */
1010 				nibbles(rs, &rrrr, &ssss);
1011 				if(ssss == 0){
1012 					if(rrrr < 15){
1013 						eobrun = 0;
1014 						if(rrrr > 0)
1015 							eobrun = receiveEOB(h, rrrr)-1;
1016 						break;
1017 					}
1018 					k += 16;
1019 				}else{
1020 					k += rrrr;
1021 					z = receive(h, ssss);
1022 					acc[k] = z*qt[k]<<Al;
1023 					if(k == Se)
1024 						break;
1025 					k++;
1026 				}
1027 			}
1028 		}
1029 
1030 		/* process restart marker, if present */
1031 		mcu++;
1032 		if(ri>0 && mcu<nmcu && mcu%ri==0){
1033 			restart(h, mcu);
1034 			eobrun = 0;
1035 		}
1036 	}
1037 }
1038 
1039 static
1040 void
1041 increment(Header *h, int acc[], int k, int Pt)
1042 {
1043 	if(acc[k] == 0)
1044 		return;
1045 	if(receivebit(h) != 0)
1046 		if(acc[k] < 0)
1047 			acc[k] -= Pt;
1048 		else
1049 			acc[k] += Pt;
1050 }
1051 
1052 static
1053 void
1054 progressiveacinc(Header *h, int comp, int Al)
1055 {
1056 	int Ns, i, z, k, Ss, Se, Ta, **ac, H, V;
1057 	int ri, mcu, nacross, ndown, nhor, nver, eobrun, nzeros, pending, x, y, tmcu, blockno, q, nmcu;
1058 	Huffman *acht;
1059 	int *qt, rrrr, ssss, *acc, rs;
1060 	uchar *ss;
1061 
1062 	ss = h->ss;
1063 	Ns = ss[0];
1064 	if(Ns != 1)
1065 		jpgerror(h, "ReadJPG: illegal Ns>1 in progressive AC scan");
1066 	Ss = ss[1+2];
1067 	Se = ss[2+2];
1068 	H = h->comp[comp].H;
1069 	V = h->comp[comp].V;
1070 
1071 	nacross = h->nacross*H;
1072 	ndown = h->ndown*V;
1073 	q = 8*h->Hmax/H;
1074 	nhor = (h->X+q-1)/q;
1075 	q = 8*h->Vmax/V;
1076 	nver = (h->Y+q-1)/q;
1077 
1078 	/* initialize data structures */
1079 	h->cnt = 0;
1080 	h->sr = 0;
1081 	h->peek = -1;
1082 	nibbles(ss[1+1], &z, &Ta);	/* z is thrown away */
1083 	ri = h->ri;
1084 
1085 	eobrun = 0;
1086 	ac = h->accoeff[comp];
1087 	acht = &h->acht[Ta];
1088 	qt = h->qt[h->comp[comp].Tq];
1089 	nmcu = nacross*ndown;
1090 	mcu = 0;
1091 	pending = 0;
1092 	nzeros = -1;
1093 	for(y=0; y<nver; y++){
1094 		for(x=0; x<nhor; x++){
1095 			/* Figure G-7 */
1096 
1097 			/*  arrange blockno to be in same sequence as original scan calculation. */
1098 			tmcu = x/H + (nacross/H)*(y/V);
1099 			blockno = tmcu*H*V + H*(y%V) + x%H;
1100 			acc = ac[blockno];
1101 			if(eobrun > 0){
1102 				if(nzeros > 0)
1103 					jpgerror(h, "ReadJPG: zeros pending at block start");
1104 				for(k=Ss; k<=Se; k++)
1105 					increment(h, acc, k, qt[k]<<Al);
1106 				--eobrun;
1107 				continue;
1108 			}
1109 
1110 			for(k=Ss; k<=Se; ){
1111 				if(nzeros >= 0){
1112 					if(acc[k] != 0)
1113 						increment(h, acc, k, qt[k]<<Al);
1114 					else if(nzeros-- == 0)
1115 						acc[k] = pending;
1116 					k++;
1117 					continue;
1118 				}
1119 				rs = decode(h, acht);
1120 				nibbles(rs, &rrrr, &ssss);
1121 				if(ssss == 0){
1122 					if(rrrr < 15){
1123 						eobrun = 0;
1124 						if(rrrr > 0)
1125 							eobrun = receiveEOB(h, rrrr)-1;
1126 						while(k <= Se){
1127 							increment(h, acc, k, qt[k]<<Al);
1128 							k++;
1129 						}
1130 						break;
1131 					}
1132 					for(i=0; i<16; k++){
1133 						increment(h, acc, k, qt[k]<<Al);
1134 						if(acc[k] == 0)
1135 							i++;
1136 					}
1137 					continue;
1138 				}else if(ssss != 1)
1139 					jpgerror(h, "ReadJPG: ssss!=1 in progressive increment");
1140 				nzeros = rrrr;
1141 				pending = receivebit(h);
1142 				if(pending == 0)
1143 					pending = -1;
1144 				pending *= qt[k]<<Al;
1145 			}
1146 		}
1147 
1148 		/* process restart marker, if present */
1149 		mcu++;
1150 		if(ri>0 && mcu<nmcu && mcu%ri==0){
1151 			restart(h, mcu);
1152 			eobrun = 0;
1153 			nzeros = -1;
1154 		}
1155 	}
1156 }
1157 
1158 static
1159 void
1160 progressivescan(Header *h, int colorspace)
1161 {
1162 	uchar *ss;
1163 	int Ns, Ss, Ah, Al, c, comp, i;
1164 
1165 	if(h->dccoeff[0] == nil)
1166 		progressiveinit(h, colorspace);
1167 
1168 	ss = h->ss;
1169 	Ns = ss[0];
1170 	Ss = ss[1+2*Ns];
1171 	nibbles(ss[3+2*Ns], &Ah, &Al);
1172 	c = ss[1];
1173 	comp = -1;
1174 	for(i=0; i<h->Nf; i++)
1175 		if(h->comp[i].C == c)
1176 			comp = i;
1177 	if(comp == -1)
1178 		jpgerror(h, "ReadJPG: bad component index in scan header");
1179 
1180 	if(Ss == 0){
1181 		progressivedc(h, comp, Ah, Al);
1182 		return;
1183 	}
1184 	if(Ah == 0){
1185 		progressiveac(h, comp, Al);
1186 		return;
1187 	}
1188 	progressiveacinc(h, comp, Al);
1189 }
1190 
1191 enum {
1192 	c1 = 2871,	/* 1.402 * 2048 */
1193 	c2 = 705,		/* 0.34414 * 2048 */
1194 	c3 = 1463,	/* 0.71414 * 2048 */
1195 	c4 = 3629	/* 1.772 * 2048 */
1196 };
1197 
1198 static
1199 void
1200 colormap1(Header *h, int colorspace, Rawimage *image, int data[8*8], int mcu, int nacross)
1201 {
1202 	uchar *pic;
1203 	int x, y, dx, dy, minx, miny;
1204 	int r, k, pici;
1205 
1206 	USED(colorspace);
1207 	pic = image->chans[0];
1208 	minx = 8*(mcu%nacross);
1209 	dx = 8;
1210 	if(minx+dx > h->X)
1211 		dx = h->X-minx;
1212 	miny = 8*(mcu/nacross);
1213 	dy = 8;
1214 	if(miny+dy > h->Y)
1215 		dy = h->Y-miny;
1216 	pici = miny*h->X+minx;
1217 	k = 0;
1218 	for(y=0; y<dy; y++){
1219 		for(x=0; x<dx; x++){
1220 			r = clamp[(data[k+x]+128)+CLAMPOFF];
1221 			pic[pici+x] = r;
1222 		}
1223 		pici += h->X;
1224 		k += 8;
1225 	}
1226 }
1227 
1228 static
1229 void
1230 colormapall1(Header *h, int colorspace, Rawimage *image, int data0[8*8], int data1[8*8], int data2[8*8], int mcu, int nacross)
1231 {
1232 	uchar *rpic, *gpic, *bpic, *rp, *gp, *bp;
1233 	int *p0, *p1, *p2;
1234 	int x, y, dx, dy, minx, miny;
1235 	int r, g, b, k, pici;
1236 	int Y, Cr, Cb;
1237 
1238 	rpic = image->chans[0];
1239 	gpic = image->chans[1];
1240 	bpic = image->chans[2];
1241 	minx = 8*(mcu%nacross);
1242 	dx = 8;
1243 	if(minx+dx > h->X)
1244 		dx = h->X-minx;
1245 	miny = 8*(mcu/nacross);
1246 	dy = 8;
1247 	if(miny+dy > h->Y)
1248 		dy = h->Y-miny;
1249 	pici = miny*h->X+minx;
1250 	k = 0;
1251 	for(y=0; y<dy; y++){
1252 		p0 = data0+k;
1253 		p1 = data1+k;
1254 		p2 = data2+k;
1255 		rp = rpic+pici;
1256 		gp = gpic+pici;
1257 		bp = bpic+pici;
1258 		if(colorspace == CYCbCr)
1259 			for(x=0; x<dx; x++){
1260 				*rp++ = clamp[*p0++ + 128 + CLAMPOFF];
1261 				*gp++ = clamp[*p1++ + 128 + CLAMPOFF];
1262 				*bp++ = clamp[*p2++ + 128 + CLAMPOFF];
1263 			}
1264 		else
1265 			for(x=0; x<dx; x++){
1266 				Y = (*p0++ + 128) << 11;
1267 				Cb = *p1++;
1268 				Cr = *p2++;
1269 				r = Y+c1*Cr;
1270 				g = Y-c2*Cb-c3*Cr;
1271 				b = Y+c4*Cb;
1272 				*rp++ = clamp[(r>>11)+CLAMPOFF];
1273 				*gp++ = clamp[(g>>11)+CLAMPOFF];
1274 				*bp++ = clamp[(b>>11)+CLAMPOFF];
1275 			}
1276 		pici += h->X;
1277 		k += 8;
1278 	}
1279 }
1280 
1281 static
1282 void
1283 colormap(Header *h, int colorspace, Rawimage *image, int *data0[8*8], int *data1[8*8], int *data2[8*8], int mcu, int nacross, int Hmax, int Vmax,  int *H, int *V)
1284 {
1285 	uchar *rpic, *gpic, *bpic;
1286 	int x, y, dx, dy, minx, miny;
1287 	int r, g, b, pici, H0, H1, H2;
1288 	int t, b0, b1, b2, y0, y1, y2, x0, x1, x2;
1289 	int Y, Cr, Cb;
1290 
1291 	rpic = image->chans[0];
1292 	gpic = image->chans[1];
1293 	bpic = image->chans[2];
1294 	minx = 8*Hmax*(mcu%nacross);
1295 	dx = 8*Hmax;
1296 	if(minx+dx > h->X)
1297 		dx = h->X-minx;
1298 	miny = 8*Vmax*(mcu/nacross);
1299 	dy = 8*Vmax;
1300 	if(miny+dy > h->Y)
1301 		dy = h->Y-miny;
1302 	pici = miny*h->X+minx;
1303 	H0 = H[0];
1304 	H1 = H[1];
1305 	H2 = H[2];
1306 	for(y=0; y<dy; y++){
1307 		t = y*V[0];
1308 		b0 = H0*(t/(8*Vmax));
1309 		y0 = 8*((t/Vmax)&7);
1310 		t = y*V[1];
1311 		b1 = H1*(t/(8*Vmax));
1312 		y1 = 8*((t/Vmax)&7);
1313 		t = y*V[2];
1314 		b2 = H2*(t/(8*Vmax));
1315 		y2 = 8*((t/Vmax)&7);
1316 		x0 = 0;
1317 		x1 = 0;
1318 		x2 = 0;
1319 		for(x=0; x<dx; x++){
1320 			if(colorspace == CYCbCr){
1321 				rpic[pici+x] = clamp[data0[b0][y0+x0++*H0/Hmax] + 128 + CLAMPOFF];
1322 				gpic[pici+x] = clamp[data1[b1][y1+x1++*H1/Hmax] + 128 + CLAMPOFF];
1323 				bpic[pici+x] = clamp[data2[b2][y2+x2++*H2/Hmax] + 128 + CLAMPOFF];
1324 			}else{
1325 				Y = (data0[b0][y0+x0++*H0/Hmax]+128)<<11;
1326 				Cb = data1[b1][y1+x1++*H1/Hmax];
1327 				Cr = data2[b2][y2+x2++*H2/Hmax];
1328 				r = Y+c1*Cr;
1329 				g = Y-c2*Cb-c3*Cr;
1330 				b = Y+c4*Cb;
1331 				rpic[pici+x] = clamp[(r>>11)+CLAMPOFF];
1332 				gpic[pici+x] = clamp[(g>>11)+CLAMPOFF];
1333 				bpic[pici+x] = clamp[(b>>11)+CLAMPOFF];
1334 			}
1335 			if(x0*H0/Hmax >= 8){
1336 				x0 = 0;
1337 				b0++;
1338 			}
1339 			if(x1*H1/Hmax >= 8){
1340 				x1 = 0;
1341 				b1++;
1342 			}
1343 			if(x2*H2/Hmax >= 8){
1344 				x2 = 0;
1345 				b2++;
1346 			}
1347 		}
1348 		pici += h->X;
1349 	}
1350 }
1351 
1352 /*
1353  * decode next 8-bit value from entropy-coded input.  chart F-26
1354  */
1355 static
1356 int
1357 decode(Header *h, Huffman *t)
1358 {
1359 	int code, v, cnt, m, sr, i;
1360 	int *maxcode;
1361 	static int badcode;
1362 
1363 	maxcode = t->maxcode;
1364 	if(h->cnt < 8)
1365 		nextbyte(h, 0);
1366 	/* fast lookup */
1367 	code = (h->sr>>(h->cnt-8))&0xFF;
1368 	v = t->value[code];
1369 	if(v >= 0){
1370 		h->cnt -= t->shift[code];
1371 		return v;
1372 	}
1373 
1374 	h->cnt -= 8;
1375 	if(h->cnt == 0)
1376 		nextbyte(h, 0);
1377 	h->cnt--;
1378 	cnt = h->cnt;
1379 	m = 1<<cnt;
1380 	sr = h->sr;
1381 	code <<= 1;
1382 	i = 9;
1383 	for(;;i++){
1384 		if(sr & m)
1385 			code |= 1;
1386 		if(code <= maxcode[i])
1387 			break;
1388 		code <<= 1;
1389 		m >>= 1;
1390 		if(m == 0){
1391 			sr = nextbyte(h, 0);
1392 			m = 0x80;
1393 			cnt = 8;
1394 		}
1395 		cnt--;
1396 	}
1397 	if(i >= 17){
1398 		if(badcode == 0)
1399 			fprint(2, "badly encoded %dx%d JPEG file; ignoring bad value\n", h->X, h->Y);
1400 		badcode = 1;
1401 		i = 0;
1402 	}
1403 	h->cnt = cnt;
1404 	return t->val[t->valptr[i]+(code-t->mincode[i])];
1405 }
1406 
1407 /*
1408  * load next byte of input
1409  */
1410 static
1411 int
1412 nextbyte(Header *h, int marker)
1413 {
1414 	int b, b2;
1415 
1416 	if(h->peek >= 0){
1417 		b = h->peek;
1418 		h->peek = -1;
1419 	}else{
1420 		b = Bgetc(h->fd);
1421 		if(b == Beof)
1422 			jpgerror(h, "truncated file");
1423 		b &= 0xFF;
1424 	}
1425 
1426 	if(b == 0xFF){
1427 		if(marker)
1428 			return b;
1429 		b2 = Bgetc(h->fd);
1430 		if(b2 != 0){
1431 			if(b2 == Beof)
1432 				jpgerror(h, "truncated file");
1433 			b2 &= 0xFF;
1434 			if(b2 == DNL)
1435 				jpgerror(h, "ReadJPG: DNL marker unimplemented");
1436 			/* decoder is reading into marker; satisfy it and restore state */
1437 			Bungetc(h->fd);
1438 			h->peek = b;
1439 		}
1440 	}
1441 	h->cnt += 8;
1442 	h->sr = (h->sr<<8) | b;
1443 	return b;
1444 }
1445 
1446 /*
1447  * return next s bits of input, MSB first, and level shift it
1448  */
1449 static
1450 int
1451 receive(Header *h, int s)
1452 {
1453 	int v, m;
1454 
1455 	while(h->cnt < s)
1456 		nextbyte(h, 0);
1457 	h->cnt -= s;
1458 	v = h->sr >> h->cnt;
1459 	m = (1<<s);
1460 	v &= m-1;
1461 	/* level shift */
1462 	if(v < (m>>1))
1463 		v += ~(m-1)+1;
1464 	return v;
1465 }
1466 
1467 /*
1468  * return next s bits of input, decode as EOB
1469  */
1470 static
1471 int
1472 receiveEOB(Header *h, int s)
1473 {
1474 	int v, m;
1475 
1476 	while(h->cnt < s)
1477 		nextbyte(h, 0);
1478 	h->cnt -= s;
1479 	v = h->sr >> h->cnt;
1480 	m = (1<<s);
1481 	v &= m-1;
1482 	/* level shift */
1483 	v += m;
1484 	return v;
1485 }
1486 
1487 /*
1488  * return next bit of input
1489  */
1490 static
1491 int
1492 receivebit(Header *h)
1493 {
1494 	if(h->cnt < 1)
1495 		nextbyte(h, 0);
1496 	h->cnt--;
1497 	return (h->sr >> h->cnt) & 1;
1498 }
1499 
1500 /*
1501  *  Scaled integer implementation.
1502  *  inverse two dimensional DCT, Chen-Wang algorithm
1503  * (IEEE ASSP-32, pp. 803-816, Aug. 1984)
1504  * 32-bit integer arithmetic (8 bit coefficients)
1505  * 11 mults, 29 adds per DCT
1506  *
1507  * coefficients extended to 12 bit for IEEE1180-1990 compliance
1508  */
1509 
1510 enum {
1511 	W1		= 2841,	/* 2048*sqrt(2)*cos(1*pi/16)*/
1512 	W2		= 2676,	/* 2048*sqrt(2)*cos(2*pi/16)*/
1513 	W3		= 2408,	/* 2048*sqrt(2)*cos(3*pi/16)*/
1514 	W5		= 1609,	/* 2048*sqrt(2)*cos(5*pi/16)*/
1515 	W6		= 1108,	/* 2048*sqrt(2)*cos(6*pi/16)*/
1516 	W7		= 565,	/* 2048*sqrt(2)*cos(7*pi/16)*/
1517 
1518 	W1pW7	= 3406,	/* W1+W7*/
1519 	W1mW7	= 2276,	/* W1-W7*/
1520 	W3pW5	= 4017,	/* W3+W5*/
1521 	W3mW5	= 799,	/* W3-W5*/
1522 	W2pW6	= 3784,	/* W2+W6*/
1523 	W2mW6	= 1567,	/* W2-W6*/
1524 
1525 	R2		= 181	/* 256/sqrt(2)*/
1526 };
1527 
1528 static
1529 void
1530 idct(int b[8*8])
1531 {
1532 	int x, y, eighty, v;
1533 	int x0, x1, x2, x3, x4, x5, x6, x7, x8;
1534 	int *p;
1535 
1536 	/* transform horizontally*/
1537 	for(y=0; y<8; y++){
1538 		eighty = y<<3;
1539 		/* if all non-DC components are zero, just propagate the DC term*/
1540 		p = b+eighty;
1541 		if(p[1]==0)
1542 		if(p[2]==0 && p[3]==0)
1543 		if(p[4]==0 && p[5]==0)
1544 		if(p[6]==0 && p[7]==0){
1545 			v = p[0]<<3;
1546 			p[0] = v;
1547 			p[1] = v;
1548 			p[2] = v;
1549 			p[3] = v;
1550 			p[4] = v;
1551 			p[5] = v;
1552 			p[6] = v;
1553 			p[7] = v;
1554 			continue;
1555 		}
1556 		/* prescale*/
1557 		x0 = (p[0]<<11)+128;
1558 		x1 = p[4]<<11;
1559 		x2 = p[6];
1560 		x3 = p[2];
1561 		x4 = p[1];
1562 		x5 = p[7];
1563 		x6 = p[5];
1564 		x7 = p[3];
1565 		/* first stage*/
1566 		x8 = W7*(x4+x5);
1567 		x4 = x8 + W1mW7*x4;
1568 		x5 = x8 - W1pW7*x5;
1569 		x8 = W3*(x6+x7);
1570 		x6 = x8 - W3mW5*x6;
1571 		x7 = x8 - W3pW5*x7;
1572 		/* second stage*/
1573 		x8 = x0 + x1;
1574 		x0 -= x1;
1575 		x1 = W6*(x3+x2);
1576 		x2 = x1 - W2pW6*x2;
1577 		x3 = x1 + W2mW6*x3;
1578 		x1 = x4 + x6;
1579 		x4 -= x6;
1580 		x6 = x5 + x7;
1581 		x5 -= x7;
1582 		/* third stage*/
1583 		x7 = x8 + x3;
1584 		x8 -= x3;
1585 		x3 = x0 + x2;
1586 		x0 -= x2;
1587 		x2 = (R2*(x4+x5)+128)>>8;
1588 		x4 = (R2*(x4-x5)+128)>>8;
1589 		/* fourth stage*/
1590 		p[0] = (x7+x1)>>8;
1591 		p[1] = (x3+x2)>>8;
1592 		p[2] = (x0+x4)>>8;
1593 		p[3] = (x8+x6)>>8;
1594 		p[4] = (x8-x6)>>8;
1595 		p[5] = (x0-x4)>>8;
1596 		p[6] = (x3-x2)>>8;
1597 		p[7] = (x7-x1)>>8;
1598 	}
1599 	/* transform vertically*/
1600 	for(x=0; x<8; x++){
1601 		/* if all non-DC components are zero, just propagate the DC term*/
1602 		p = b+x;
1603 		if(p[8*1]==0)
1604 		if(p[8*2]==0 && p[8*3]==0)
1605 		if(p[8*4]==0 && p[8*5]==0)
1606 		if(p[8*6]==0 && p[8*7]==0){
1607 			v = (p[8*0]+32)>>6;
1608 			p[8*0] = v;
1609 			p[8*1] = v;
1610 			p[8*2] = v;
1611 			p[8*3] = v;
1612 			p[8*4] = v;
1613 			p[8*5] = v;
1614 			p[8*6] = v;
1615 			p[8*7] = v;
1616 			continue;
1617 		}
1618 		/* prescale*/
1619 		x0 = (p[8*0]<<8)+8192;
1620 		x1 = p[8*4]<<8;
1621 		x2 = p[8*6];
1622 		x3 = p[8*2];
1623 		x4 = p[8*1];
1624 		x5 = p[8*7];
1625 		x6 = p[8*5];
1626 		x7 = p[8*3];
1627 		/* first stage*/
1628 		x8 = W7*(x4+x5) + 4;
1629 		x4 = (x8+W1mW7*x4)>>3;
1630 		x5 = (x8-W1pW7*x5)>>3;
1631 		x8 = W3*(x6+x7) + 4;
1632 		x6 = (x8-W3mW5*x6)>>3;
1633 		x7 = (x8-W3pW5*x7)>>3;
1634 		/* second stage*/
1635 		x8 = x0 + x1;
1636 		x0 -= x1;
1637 		x1 = W6*(x3+x2) + 4;
1638 		x2 = (x1-W2pW6*x2)>>3;
1639 		x3 = (x1+W2mW6*x3)>>3;
1640 		x1 = x4 + x6;
1641 		x4 -= x6;
1642 		x6 = x5 + x7;
1643 		x5 -= x7;
1644 		/* third stage*/
1645 		x7 = x8 + x3;
1646 		x8 -= x3;
1647 		x3 = x0 + x2;
1648 		x0 -= x2;
1649 		x2 = (R2*(x4+x5)+128)>>8;
1650 		x4 = (R2*(x4-x5)+128)>>8;
1651 		/* fourth stage*/
1652 		p[8*0] = (x7+x1)>>14;
1653 		p[8*1] = (x3+x2)>>14;
1654 		p[8*2] = (x0+x4)>>14;
1655 		p[8*3] = (x8+x6)>>14;
1656 		p[8*4] = (x8-x6)>>14;
1657 		p[8*5] = (x0-x4)>>14;
1658 		p[8*6] = (x3-x2)>>14;
1659 		p[8*7] = (x7-x1)>>14;
1660 	}
1661 }