Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 **p)
172 
{
173 
	if(*p){
174 
		free(*p);
175 
		*p = nil;
176 
	}
177 
}
178 
 
179 
static
180 
void
181 
jpgfreeall(Header *h, int freeimage)
182 
{
183 
	int i, j;
184 
 
185 
	clear(&h->buf);
186 
	if(h->dccoeff[0])
187 
		for(i=0; i<3; i++)
188 
			clear(&h->dccoeff[i]);
189 
	if(h->accoeff[0])
190 
		for(i=0; i<3; i++){
191 
			if(h->accoeff[i])
192 
				for(j=0; j<h->naccoeff[i]; j++)
193 
					clear(&h->accoeff[i][j]);
194 
			clear(&h->accoeff[i]);
195 
		}
196 
	for(i=0; i<4; i++){
197 
		clear(&h->dcht[i].size);
198 
		clear(&h->acht[i].size);
199 
		clear(&h->dcht[i].code);
200 
		clear(&h->acht[i].code);
201 
		clear(&h->dcht[i].val);
202 
		clear(&h->acht[i].val);
203 
	}
204 
	if(h->data[0])
205 
		for(i=0; i<3; i++){
206 
			if(h->data[i])
207 
				for(j=0; j<h->ndata[i]; j++)
208 
					clear(&h->data[i][j]);
209 
			clear(&h->data[i]);
210 
		}
211 
	if(freeimage && h->image!=nil){
212 
		clear(&h->array);
213 
		clear(&h->image->cmap);
214 
		for(i=0; i<3; i++)
215 
			clear(&h->image->chans[i]);
216 
		clear(&h->image);
217 
	}
218 
}
219 
 
220 
static
221 
void
222 
jpgerror(Header *h, char *fmt, ...)
223 
{
224 
	va_list arg;
225 
 
226 
	va_start(arg, fmt);
227 
	vseprint(h->err, h->err+sizeof h->err, fmt, arg);
228 
	va_end(arg);
229 
 
230 
	werrstr(h->err);
231 
	jpgfreeall(h, 1);
232 
	longjmp(h->errlab, 1);
233 
}
234 
 
235 
Rawimage**
236 
Breadjpg(Biobuf *b, int colorspace)
237 
{
238 
	Rawimage *r, **array;
239 
	Header *h;
240 
	char buf[ERRMAX];
241 
 
242 
	buf[0] = '\0';
243 
	if(colorspace!=CYCbCr && colorspace!=CRGB){
244 
		errstr(buf, sizeof buf);	/* throw it away */
245 
		werrstr("ReadJPG: unknown color space");
246 
		return nil;
247 
	}
248 
	jpginit();
249 
	h = malloc(sizeof(Header));
250 
	array = malloc(sizeof(Header));
251 
	if(h==nil || array==nil){
252 
		free(h);
253 
		free(array);
254 
		return nil;
255 
	}
256 
	h->array = array;
257 
	memset(h, 0, sizeof(Header));
258 
	h->fd = b;
259 
	errstr(buf, sizeof buf);	/* throw it away */
260 
	if(setjmp(h->errlab))
261 
		r = nil;
262 
	else
263 
		r = readslave(h, colorspace);
264 
	jpgfreeall(h, 0);
265 
	free(h);
266 
	array[0] = r;
267 
	array[1] = nil;
268 
	return array;
269 
}
270 
 
271 
Rawimage**
272 
readjpg(int fd, int colorspace)
273 
{
274 
	Rawimage** a;
275 
	Biobuf b;
276 
 
277 
	if(Binit(&b, fd, OREAD) < 0)
278 
		return nil;
279 
	a = Breadjpg(&b, colorspace);
280 
	Bterm(&b);
281 
	return a;
282 
}
283 
 
284 
static
285 
Rawimage*
286 
readslave(Header *header, int colorspace)
287 
{
288 
	Rawimage *image;
289 
	int nseg, i, H, V, m, n;
290 
	uchar *b;
291 
 
292 
	soiheader(header);
293 
	nseg = 0;
294 
	image = nil;
295 
 
296 
	header->buf = jpgmalloc(header, 4096, 0);
297 
	header->nbuf = 4096;
298 
	while(header->err[0] == '\0'){
299 
		nseg++;
300 
		n = readsegment(header, &m);
301 
		b = header->buf;
302 
		switch(m){
303 
		case -1:
304 
			return image;
305 
 
306 
		case APPn+0:
307 
			if(nseg==1 && strncmp((char*)b, "JFIF", 4)==0)  /* JFIF header; check version */
308 
				if(b[5]>1 || b[6]>2)
309 
					sprint(header->err, "ReadJPG: can't handle JFIF version %d.%2d", b[5], b[6]);
310 
			break;
311 
 
312 
		case APPn+1: case APPn+2: case APPn+3: case APPn+4: case APPn+5:
313 
		case APPn+6: case APPn+7: case APPn+8: case APPn+9: case APPn+10:
314 
		case APPn+11: case APPn+12: case APPn+13: case APPn+14: case APPn+15:
315 
			break;
316 
 
317 
		case DQT:
318 
			quanttables(header, b, n);
319 
			break;
320 
 
321 
		case SOF:
322 
		case SOF2:
323 
			header->Y = int2(b, 1);
324 
			header->X = int2(b, 3);
325 
			header->Nf =b[5];
326 
			for(i=0; i<header->Nf; i++){
327 
				header->comp[i].C = b[6+3*i+0];
328 
				nibbles(b[6+3*i+1], &H, &V);
329 
				if(H<=0 || V<=0)
330 
					jpgerror(header, "non-positive sampling factor (Hsamp or Vsamp)");
331 
				header->comp[i].H = H;
332 
				header->comp[i].V = V;
333 
				header->comp[i].Tq = b[6+3*i+2];
334 
			}
335 
			header->mode = m;
336 
			header->sf = b;
337 
			break;
338 
 
339 
		case  SOS:
340 
			header->ss = b;
341 
			switch(header->mode){
342 
			case SOF:
343 
				image = baselinescan(header, colorspace);
344 
				break;
345 
			case SOF2:
346 
				progressivescan(header, colorspace);
347 
				break;
348 
			default:
349 
				sprint(header->err, "unrecognized or unspecified encoding %d", header->mode);
350 
				break;
351 
			}
352 
			break;
353 
 
354 
		case  DHT:
355 
			huffmantables(header, b, n);
356 
			break;
357 
 
358 
		case  DRI:
359 
			header->ri = int2(b, 0);
360 
			break;
361 
 
362 
		case  COM:
363 
			break;
364 
 
365 
		case EOI:
366 
			if(header->mode == SOF2)
367 
				image = progressiveIDCT(header, colorspace);
368 
			return image;
369 
 
370 
		default:
371 
			sprint(header->err, "ReadJPG: unknown marker %.2x", m);
372 
			break;
373 
		}
374 
	}
375 
	return image;
376 
}
377 
 
378 
/* readsegment is called after reading scan, which can have */
379 
/* read ahead a byte.  so we must check peek here */
380 
static
381 
int
382 
readbyte(Header *h)
383 
{
384 
	uchar x;
385 
 
386 
	if(h->peek >= 0){
387 
		x = h->peek;
388 
		h->peek = -1;
389 
	}else if(Bread(h->fd, &x, 1) != 1)
390 
		jpgerror(h, readerr);
391 
	return x;
392 
}
393 
 
394 
static
395 
int
396 
marker(Header *h)
397 
{
398 
	int c;
399 
 
400 
	while((c=readbyte(h)) == 0)
401 
		fprint(2, "ReadJPG: skipping zero byte at offset %lld\n", Boffset(h->fd));
402 
	if(c != 0xFF)
403 
		jpgerror(h, "ReadJPG: expecting marker; found 0x%x at offset %lld\n", c, Boffset(h->fd));
404 
	while(c == 0xFF)
405 
		c = readbyte(h);
406 
	return c;
407 
}
408 
 
409 
static
410 
int
411 
int2(uchar *buf, int n)
412 
{
413 
	return (buf[n]<<8) + buf[n+1];
414 
}
415 
 
416 
static
417 
void
418 
nibbles(int b, int *p0, int *p1)
419 
{
420 
	*p0 = (b>>4) & 0xF;
421 
	*p1 = b & 0xF;
422 
}
423 
 
424 
static
425 
void
426 
soiheader(Header *h)
427 
{
428 
	h->peek = -1;
429 
	if(marker(h) != SOI)
430 
		jpgerror(h, "ReadJPG: unrecognized marker in header");
431 
	h->err[0] = '\0';
432 
	h->mode = 0;
433 
	h->ri = 0;
434 
}
435 
 
436 
static
437 
int
438 
readsegment(Header *h, int *markerp)
439 
{
440 
	int m, n;
441 
	uchar tmp[2];
442 
 
443 
	m = marker(h);
444 
	switch(m){
445 
	case EOI:
446 
		*markerp = m;
447 
		return 0;
448 
	case 0:
449 
		jpgerror(h, "ReadJPG: expecting marker; saw %.2x at offset %lld", m, Boffset(h->fd));
450 
	}
451 
	if(Bread(h->fd, tmp, 2) != 2)
452 
    Readerr:
453 
		jpgerror(h, readerr);
454 
	n = int2(tmp, 0);
455 
	if(n < 2)
456 
		goto Readerr;
457 
	n -= 2;
458 
	if(n > h->nbuf){
459 
		free(h->buf);
460 
		/* zero in case of short read later */
461 
		h->buf = jpgmalloc(h, n+1, 1); /* +1 for sentinel */
462 
		h->nbuf = n;
463 
	}
464 
	/* accept short reads to cope with some real-world jpegs */
465 
	if(Bread(h->fd, h->buf, n) < 0)
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 
}