treadjpg.c - plan9port - [fork] Plan 9 from user space
HTML git clone git://src.adamsgaard.dk/plan9port
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
treadjpg.c (33031B)
---
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 }