sys/src/cmd/jpg/readjpg.c

   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("%s", 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(2*sizeof(Rawimage*));
 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                                 /* hack: colormap1() doesnt handle resampling */
 332                                 if(header->Nf == 1)
 333                                         H = V = 1;
 334                                 header->comp[i].H = H;
 335                                 header->comp[i].V = V;
 336                                 header->comp[i].Tq = b[6+3*i+2];
 337                         }
 338                         header->mode = m;
 339                         header->sf = b;
 340                         break;
 341
 342                 case  SOS:
 343                         header->ss = b;
 344                         switch(header->mode){
 345                         case SOF:
 346                                 image = baselinescan(header, colorspace);
 347                                 break;
 348                         case SOF2:
 349                                 progressivescan(header, colorspace);
 350                                 break;
 351                         default:
 352                                 sprint(header->err, "unrecognized or unspecified encoding %d", header->mode);
 353                                 break;
 354                         }
 355                         break;
 356
 357                 case  DHT:
 358                         huffmantables(header, b, n);
 359                         break;
 360
 361                 case  DRI:
 362                         header->ri = int2(b, 0);
 363                         break;
 364
 365                 case  COM:
 366                         break;
 367
 368                 case EOI:
 369                         if(header->mode == SOF2)
 370                                 image = progressiveIDCT(header, colorspace);
 371                         return image;
 372
 373                 default:
 374                         sprint(header->err, "ReadJPG: unknown marker %.2x", m);
 375                         break;
 376                 }
 377         }
 378         return image;
 379 }
 380
 381 /* readsegment is called after reading scan, which can have */
 382 /* read ahead a byte.  so we must check peek here */
 383 static
 384 int
 385 readbyte(Header *h)
 386 {
 387         uchar x;
 388
 389         if(h->peek >= 0){
 390                 x = h->peek;
 391                 h->peek = -1;
 392         }else if(Bread(h->fd, &x, 1) != 1)
 393                 jpgerror(h, readerr);
 394         return x;
 395 }
 396
 397 static
 398 int
 399 marker(Header *h)
 400 {
 401         int c;
 402
 403 Again:
 404         while((c=readbyte(h)) == 0)
 405                 ;
 406         if(c != 0xFF)
 407                 goto Again;
 408         while(c == 0xFF)
 409                 c = readbyte(h);
 410         if(c == 0)
 411                 goto Again;
 412         return c;
 413 }
 414
 415 static
 416 int
 417 int2(uchar *buf, int n)
 418 {
 419         return (buf[n]<<8) + buf[n+1];
 420 }
 421
 422 static
 423 void
 424 nibbles(int b, int *p0, int *p1)
 425 {
 426         *p0 = (b>>4) & 0xF;
 427         *p1 = b & 0xF;
 428 }
 429
 430 static
 431 void
 432 soiheader(Header *h)
 433 {
 434         h->peek = -1;
 435         if(marker(h) != SOI)
 436                 jpgerror(h, "ReadJPG: unrecognized marker in header");
 437         h->err[0] = '\0';
 438         h->mode = 0;
 439         h->ri = 0;
 440 }
 441
 442 static
 443 int
 444 readsegment(Header *h, int *markerp)
 445 {
 446         int m, n;
 447         uchar tmp[2];
 448
 449         if((m = marker(h)) == EOI){
 450                 *markerp = m;
 451                 return 0;
 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                 /* zero in case of short read later */
 463                 h->buf = jpgmalloc(h, n+1, 1); /* +1 for sentinel */
 464                 h->nbuf = n;
 465         }
 466         /* accept short reads to cope with some real-world jpegs */
 467         if(Bread(h->fd, h->buf, n) < 0)
 468                 goto Readerr;
 469         *markerp = m;
 470         return n;
 471 }
 472
 473 static
 474 int
 475 huffmantable(Header *h, uchar *b)
 476 {
 477         Huffman *t;
 478         int Tc, th, n, nsize, i, j, k, v, cnt, code, si, sr;
 479         int *maxcode;
 480
 481         nibbles(b[0], &Tc, &th);
 482         if(Tc > 1)
 483                 jpgerror(h, "ReadJPG: unknown Huffman table class %d", Tc);
 484         if(th>3 || (h->mode==SOF && th>1))
 485                 jpgerror(h, "ReadJPG: unknown Huffman table index %d", th);
 486         if(Tc == 0)
 487                 t = &h->dcht[th];
 488         else
 489                 t = &h->acht[th];
 490
 491         /* flow chart C-2 */
 492         nsize = 0;
 493         for(i=1; i<=16; i++)
 494                 nsize += b[i];
 495         if(nsize == 0)
 496                 return 0;
 497         t->size = jpgmalloc(h, (nsize+1)*sizeof(int), 1);
 498         k = 0;
 499         for(i=1; i<=16; i++){
 500                 n = b[i];
 501                 for(j=0; j<n; j++)
 502                         t->size[k++] = i;
 503         }
 504         t->size[k] = 0;
 505
 506         /* initialize HUFFVAL */
 507         t->val = jpgmalloc(h, nsize*sizeof(int), 1);
 508         for(i=0; i<nsize; i++)
 509                 t->val[i] = b[17+i];
 510
 511         /* flow chart C-3 */
 512         t->code = jpgmalloc(h, (nsize+1)*sizeof(int), 1);
 513         k = 0;
 514         code = 0;
 515         si = t->size[0];
 516         for(;;){
 517                 do
 518                         t->code[k++] = code++;
 519                 while(t->size[k] == si);
 520                 if(t->size[k] == 0)
 521                         break;
 522                 do{
 523                         code <<= 1;
 524                         si++;
 525                 }while(t->size[k] != si);
 526         }
 527
 528         /* flow chart F-25 */
 529         i = 0;
 530         j = 0;
 531         for(;;){
 532                 for(;;){
 533                         i++;
 534                         if(i > 16)
 535                                 goto outF25;
 536                         if(b[i] != 0)
 537                                 break;
 538                         t->maxcode[i] = -1;
 539                 }
 540                 t->valptr[i] = j;
 541                 t->mincode[i] = t->code[j];
 542                 j += b[i]-1;
 543                 t->maxcode[i] = t->code[j];
 544                 j++;
 545         }
 546 outF25:
 547
 548         /* create byte-indexed fast path tables */
 549         maxcode = t->maxcode;
 550         /* stupid startup algorithm: just run machine for each byte value */
 551         for(v=0; v<256; ){
 552                 code = 0;
 553                 cnt = 8;
 554                 sr = v;
 555                 for(i=1;;i++){
 556                         cnt--;
 557                         if(sr & (1<<cnt))
 558                                 code |= 1;
 559                         if(code <= maxcode[i])
 560                                 break;
 561                         code <<= 1;
 562                         if(cnt == 0){
 563                                 t->shift[v] = 0;
 564                                 t->value[v] = -1;
 565                                 goto continueBytes;
 566                         }
 567                 }
 568                 t->shift[v] = 8-cnt;
 569                 t->value[v] = t->val[t->valptr[i]+(code-t->mincode[i])];
 570
 571     continueBytes:
 572                 v++;
 573         }
 574
 575         return nsize;
 576 }
 577
 578 static
 579 void
 580 huffmantables(Header *h, uchar *b, int n)
 581 {
 582         int l, mt;
 583
 584         for(l=0; l<n; l+=17+mt)
 585                 mt = huffmantable(h, &b[l]);
 586 }
 587
 588 static
 589 int
 590 quanttable(Header *h, uchar *b)
 591 {
 592         int i, pq, tq, *q;
 593
 594         nibbles(b[0], &pq, &tq);
 595         if(pq > 1)
 596                 jpgerror(h, "ReadJPG: unknown quantization table class %d", pq);
 597         if(tq > 3)
 598                 jpgerror(h, "ReadJPG: unknown quantization table index %d", tq);
 599         q = h->qt[tq];
 600         for(i=0; i<64; i++){
 601                 if(pq == 0)
 602                         q[i] = b[1+i];
 603                 else
 604                         q[i] = int2(b, 1+2*i);
 605         }
 606         return 64*(1+pq);
 607 }
 608
 609 static
 610 void
 611 quanttables(Header *h, uchar *b, int n)
 612 {
 613         int l, m;
 614
 615         for(l=0; l<n; l+=1+m)
 616                 m = quanttable(h, &b[l]);
 617 }
 618
 619 static
 620 Rawimage*
 621 baselinescan(Header *h, int colorspace)
 622 {
 623         int Ns, z, k, m, Hmax, Vmax, comp;
 624         int allHV1, nblock, ri, mcu, nacross, nmcu;
 625         Huffman *dcht, *acht;
 626         int block, t, diff, *qt;
 627         uchar *ss;
 628         Rawimage *image;
 629         int Td[3], Ta[3], H[3], V[3], DC[3];
 630         int ***data, *zz;
 631
 632         ss = h->ss;
 633         Ns = ss[0];
 634         if((Ns!=3 && Ns!=1) || Ns!=h->Nf)
 635                 jpgerror(h, "ReadJPG: can't handle scan not 3 components");
 636
 637         image = jpgmalloc(h, sizeof(Rawimage), 1);
 638         h->image = image;
 639         image->r = Rect(0, 0, h->X, h->Y);
 640         image->cmap = nil;
 641         image->cmaplen = 0;
 642         image->chanlen = h->X*h->Y;
 643         image->fields = 0;
 644         image->gifflags = 0;
 645         image->gifdelay = 0;
 646         image->giftrindex = 0;
 647         if(Ns == 3)
 648                 image->chandesc = colorspace;
 649         else
 650                 image->chandesc = CY;
 651         image->nchans = h->Nf;
 652         for(k=0; k<h->Nf; k++)
 653                 image->chans[k] = jpgmalloc(h, h->X*h->Y, 0);
 654
 655         /* compute maximum H and V */
 656         Hmax = 0;
 657         Vmax = 0;
 658         for(comp=0; comp<Ns; comp++){
 659                 if(h->comp[comp].H > Hmax)
 660                         Hmax = h->comp[comp].H;
 661                 if(h->comp[comp].V > Vmax)
 662                         Vmax = h->comp[comp].V;
 663         }
 664
 665         /* initialize data structures */
 666         allHV1 = 1;
 667         data = h->data;
 668         for(comp=0; comp<Ns; comp++){
 669                 /* JPEG requires scan components to be in same order as in frame, */
 670                 /* so if both have 3 we know scan is Y Cb Cr and there's no need to */
 671                 /* reorder */
 672                 nibbles(ss[2+2*comp], &Td[comp], &Ta[comp]);
 673                 H[comp] = h->comp[comp].H;
 674                 V[comp] = h->comp[comp].V;
 675                 nblock = H[comp]*V[comp];
 676                 if(nblock != 1)
 677                         allHV1 = 0;
 678                 data[comp] = jpgmalloc(h, nblock*sizeof(int*), 0);
 679                 h->ndata[comp] = nblock;
 680                 DC[comp] = 0;
 681                 for(m=0; m<nblock; m++)
 682                         data[comp][m] = jpgmalloc(h, 8*8*sizeof(int), 0);
 683         }
 684
 685         ri = h->ri;
 686
 687         h->cnt = 0;
 688         h->sr = 0;
 689         h->peek = -1;
 690         nacross = ((h->X+(8*Hmax-1))/(8*Hmax));
 691         nmcu = ((h->Y+(8*Vmax-1))/(8*Vmax))*nacross;
 692         for(mcu=0; mcu<nmcu; ){
 693                 for(comp=0; comp<Ns; comp++){
 694                         dcht = &h->dcht[Td[comp]];
 695                         acht = &h->acht[Ta[comp]];
 696                         qt = h->qt[h->comp[comp].Tq];
 697
 698                         for(block=0; block<H[comp]*V[comp]; block++){
 699                                 /* F-22 */
 700                                 t = decode(h, dcht);
 701                                 diff = receive(h, t);
 702                                 DC[comp] += diff;
 703
 704                                 /* F-23 */
 705                                 zz = data[comp][block];
 706                                 memset(zz, 0, 8*8*sizeof(int));
 707                                 zz[0] = qt[0]*DC[comp];
 708                                 k = 1;
 709                                 do{
 710                                         t = decode(h, acht);
 711                                         assert(t >= 0);
 712                                         if((t&0x0F) == 0){
 713                                                 if((t&0xF0) != 0xF0)
 714                                                         break;
 715                                                 k += 16;
 716                                         }else{
 717                                                 z = receive(h, t&0xF);
 718                                                 k += t>>4;
 719                                                 if(k >= 64)
 720                                                         break;
 721                                                 zz[zig[k]] = z*qt[k];
 722                                                 k++;
 723                                         }
 724                                 } while(k < 64);
 725
 726                                 idct(zz);
 727                         }
 728                 }
 729
 730                 /* rotate colors to RGB and assign to bytes */
 731                 if(Ns == 1) /* very easy */
 732                         colormap1(h, colorspace, image, data[0][0], mcu, nacross);
 733                 else if(allHV1) /* fairly easy */
 734                         colormapall1(h, colorspace, image, data[0][0], data[1][0], data[2][0], mcu, nacross);
 735                 else /* miserable general case */
 736                         colormap(h, colorspace, image, data[0], data[1], data[2], mcu, nacross, Hmax, Vmax, H, V);
 737                 /* process restart marker, if present */
 738                 mcu++;
 739                 if(ri>0 && mcu<nmcu && mcu%ri==0){
 740                         restart(h, mcu);
 741                         for(comp=0; comp<Ns; comp++)
 742                                 DC[comp] = 0;
 743                 }
 744         }
 745         return image;
 746 }
 747
 748 static
 749 void
 750 restart(Header *h, int mcu)
 751 {
 752         int rest, rst, nskip;
 753
 754         rest = mcu/h->ri-1;
 755         nskip = 0;
 756         do{
 757                 do{
 758                         rst = nextbyte(h, 1);
 759                         nskip++;
 760                 }while(rst>=0 && rst!=0xFF);
 761                 if(rst == 0xFF){
 762                         rst = nextbyte(h, 1);
 763                         nskip++;
 764                 }
 765         }while(rst>=0 && (rst&~7)!=RST);
 766         if(nskip != 2)
 767                 sprint(h->err, "ReadJPG: skipped %d bytes at restart %d\n", nskip-2, rest);
 768         if(rst < 0)
 769                 jpgerror(h, readerr);
 770         if((rst&7) != (rest&7))
 771                 jpgerror(h, "ReadJPG: expected RST%d got %d", rest&7, rst&7);
 772         h->cnt = 0;
 773         h->sr = 0;
 774 }
 775
 776 static
 777 Rawimage*
 778 progressiveIDCT(Header *h, int colorspace)
 779 {
 780         int k, m, comp, block, Nf, bn;
 781         int allHV1, nblock, mcu, nmcu;
 782         int H[3], V[3], blockno[3];
 783         int *dccoeff, **accoeff;
 784         int ***data, *zz;
 785
 786         Nf = h->Nf;
 787         allHV1 = 1;
 788         data = h->data;
 789
 790         for(comp=0; comp<Nf; comp++){
 791                 H[comp] = h->comp[comp].H;
 792                 V[comp] = h->comp[comp].V;
 793                 nblock = h->nblock[comp];
 794                 if(nblock != 1)
 795                         allHV1 = 0;
 796                 h->ndata[comp] = nblock;
 797                 data[comp] = jpgmalloc(h, nblock*sizeof(int*), 0);
 798                 for(m=0; m<nblock; m++)
 799                         data[comp][m] = jpgmalloc(h, 8*8*sizeof(int), 0);
 800         }
 801
 802         memset(blockno, 0, sizeof blockno);
 803         nmcu = h->nacross*h->ndown;
 804         for(mcu=0; mcu<nmcu; mcu++){
 805                 for(comp=0; comp<Nf; comp++){
 806                         dccoeff = h->dccoeff[comp];
 807                         accoeff = h->accoeff[comp];
 808                         bn = blockno[comp];
 809                         for(block=0; block<h->nblock[comp]; block++){
 810                                 zz = data[comp][block];
 811                                 memset(zz, 0, 8*8*sizeof(int));
 812                                 zz[0] = dccoeff[bn];
 813
 814                                 for(k=1; k<64; k++)
 815                                         zz[zig[k]] = accoeff[bn][k];
 816
 817                                 idct(zz);
 818                                 bn++;
 819                         }
 820                         blockno[comp] = bn;
 821                 }
 822
 823                 /* rotate colors to RGB and assign to bytes */
 824                 if(Nf == 1) /* very easy */
 825                         colormap1(h, colorspace, h->image, data[0][0], mcu, h->nacross);
 826                 else if(allHV1) /* fairly easy */
 827                         colormapall1(h, colorspace, h->image, data[0][0], data[1][0], data[2][0], mcu, h->nacross);
 828                 else /* miserable general case */
 829                         colormap(h, colorspace, h->image, data[0], data[1], data[2], mcu, h->nacross, h->Hmax, h->Vmax, H, V);
 830         }
 831
 832         return h->image;
 833 }
 834
 835 static
 836 void
 837 progressiveinit(Header *h, int colorspace)
 838 {
 839         int Nf, Ns, j, k, nmcu, comp;
 840         uchar *ss;
 841         Rawimage *image;
 842
 843         ss = h->ss;
 844         Ns = ss[0];
 845         Nf = h->Nf;
 846         if((Ns!=3 && Ns!=1) || Ns!=Nf)
 847                 jpgerror(h, "ReadJPG: image must have 1 or 3 components");
 848
 849         image = jpgmalloc(h, sizeof(Rawimage), 1);
 850         h->image = image;
 851         image->r = Rect(0, 0, h->X, h->Y);
 852         image->cmap = nil;
 853         image->cmaplen = 0;
 854         image->chanlen = h->X*h->Y;
 855         image->fields = 0;
 856         image->gifflags = 0;
 857         image->gifdelay = 0;
 858         image->giftrindex = 0;
 859         if(Nf == 3)
 860                 image->chandesc = colorspace;
 861         else
 862                 image->chandesc = CY;
 863         image->nchans = h->Nf;
 864         for(k=0; k<Nf; k++){
 865                 image->chans[k] = jpgmalloc(h, h->X*h->Y, 0);
 866                 h->nblock[k] = h->comp[k].H*h->comp[k].V;
 867         }
 868
 869         /* compute maximum H and V */
 870         h->Hmax = 0;
 871         h->Vmax = 0;
 872         for(comp=0; comp<Nf; comp++){
 873                 if(h->comp[comp].H > h->Hmax)
 874                         h->Hmax = h->comp[comp].H;
 875                 if(h->comp[comp].V > h->Vmax)
 876                         h->Vmax = h->comp[comp].V;
 877         }
 878         h->nacross = ((h->X+(8*h->Hmax-1))/(8*h->Hmax));
 879         h->ndown = ((h->Y+(8*h->Vmax-1))/(8*h->Vmax));
 880         nmcu = h->nacross*h->ndown;
 881
 882         for(k=0; k<Nf; k++){
 883                 h->dccoeff[k] = jpgmalloc(h, h->nblock[k]*nmcu * sizeof(int), 1);
 884                 h->accoeff[k] = jpgmalloc(h, h->nblock[k]*nmcu * sizeof(int*), 1);
 885                 h->naccoeff[k] = h->nblock[k]*nmcu;
 886                 for(j=0; j<h->nblock[k]*nmcu; j++)
 887                         h->accoeff[k][j] = jpgmalloc(h, 64*sizeof(int), 1);
 888         }
 889
 890 }
 891
 892 static
 893 void
 894 progressivedc(Header *h, int comp, int Ah, int Al)
 895 {
 896         int Ns, z, ri, mcu,  nmcu;
 897         int block, t, diff, qt, *dc, bn;
 898         Huffman *dcht;
 899         uchar *ss;
 900         int Td[3], DC[3], blockno[3];
 901
 902         ss= h->ss;
 903         Ns = ss[0];
 904         if(Ns!=h->Nf)
 905                 jpgerror(h, "ReadJPG: can't handle progressive with Nf!=Ns in DC scan");
 906
 907         /* initialize data structures */
 908         h->cnt = 0;
 909         h->sr = 0;
 910         h->peek = -1;
 911         for(comp=0; comp<Ns; comp++){
 912                 /*
 913                  * JPEG requires scan components to be in same order as in frame,
 914                  * so if both have 3 we know scan is Y Cb Cr and there's no need to
 915                  * reorder
 916                  */
 917                 nibbles(ss[2+2*comp], &Td[comp], &z);   /* z is ignored */
 918                 DC[comp] = 0;
 919         }
 920
 921         ri = h->ri;
 922
 923         nmcu = h->nacross*h->ndown;
 924         memset(blockno, 0, sizeof blockno);
 925         for(mcu=0; mcu<nmcu; ){
 926                 for(comp=0; comp<Ns; comp++){
 927                         dcht = &h->dcht[Td[comp]];
 928                         qt = h->qt[h->comp[comp].Tq][0];
 929                         dc = h->dccoeff[comp];
 930                         bn = blockno[comp];
 931
 932                         for(block=0; block<h->nblock[comp]; block++){
 933                                 if(Ah == 0){
 934                                         t = decode(h, dcht);
 935                                         diff = receive(h, t);
 936                                         DC[comp] += diff;
 937                                         dc[bn] = qt*DC[comp]<<Al;
 938                                 }else
 939                                         dc[bn] |= qt*receivebit(h)<<Al;
 940                                 bn++;
 941                         }
 942                         blockno[comp] = bn;
 943                 }
 944
 945                 /* process restart marker, if present */
 946                 mcu++;
 947                 if(ri>0 && mcu<nmcu && mcu%ri==0){
 948                         restart(h, mcu);
 949                         for(comp=0; comp<Ns; comp++)
 950                                 DC[comp] = 0;
 951                 }
 952         }
 953 }
 954
 955 static
 956 void
 957 progressiveac(Header *h, int comp, int Al)
 958 {
 959         int Ns, Ss, Se, z, k, eobrun, x, y, nver, tmcu, blockno, *acc, rs;
 960         int ri, mcu, nacross, ndown, nmcu, nhor;
 961         Huffman *acht;
 962         int *qt, rrrr, ssss, q;
 963         uchar *ss;
 964         int Ta, H, V;
 965
 966         ss = h->ss;
 967         Ns = ss[0];
 968         if(Ns != 1)
 969                 jpgerror(h, "ReadJPG: illegal Ns>1 in progressive AC scan");
 970         Ss = ss[1+2];
 971         Se = ss[2+2];
 972         H = h->comp[comp].H;
 973         V = h->comp[comp].V;
 974
 975         nacross = h->nacross*H;
 976         ndown = h->ndown*V;
 977         q = 8*h->Hmax/H;
 978         nhor = (h->X+q-1)/q;
 979         q = 8*h->Vmax/V;
 980         nver = (h->Y+q-1)/q;
 981
 982         /* initialize data structures */
 983         h->cnt = 0;
 984         h->sr = 0;
 985         h->peek = -1;
 986         nibbles(ss[1+1], &z, &Ta);      /* z is thrown away */
 987
 988         ri = h->ri;
 989
 990         eobrun = 0;
 991         acht = &h->acht[Ta];
 992         qt = h->qt[h->comp[comp].Tq];
 993         nmcu = nacross*ndown;
 994         mcu = 0;
 995         for(y=0; y<nver; y++){
 996                 for(x=0; x<nhor; x++){
 997                         /* Figure G-3  */
 998                         if(eobrun > 0){
 999                                 --eobrun;
1000                                 continue;
1001                         }
1002
1003                         /* arrange blockno to be in same sequence as original scan calculation. */
1004                         tmcu = x/H + (nacross/H)*(y/V);
1005                         blockno = tmcu*H*V + H*(y%V) + x%H;
1006                         acc = h->accoeff[comp][blockno];
1007                         k = Ss;
1008                         do {
1009                                 rs = decode(h, acht);
1010                                 /* XXX remove rrrr ssss as in baselinescan */
1011                                 nibbles(rs, &rrrr, &ssss);
1012                                 if(ssss == 0){
1013                                         if(rrrr < 15){
1014                                                 eobrun = 0;
1015                                                 if(rrrr > 0)
1016                                                         eobrun = receiveEOB(h, rrrr)-1;
1017                                                 break;
1018                                         }
1019                                         k += 16;
1020                                 }else{
1021                                         z = receive(h, ssss);
1022                                         k += rrrr;
1023                                         if(k > Se)
1024                                                 break;
1025                                         acc[k] = z*qt[k]<<Al;
1026                                         k++;
1027                                 }
1028                         } while(k <= Se);
1029                 }
1030
1031                 /* process restart marker, if present */
1032                 mcu++;
1033                 if(ri>0 && mcu<nmcu && mcu%ri==0){
1034                         restart(h, mcu);
1035                         eobrun = 0;
1036                 }
1037         }
1038 }
1039
1040 static
1041 void
1042 increment(Header *h, int acc[], int k, int Pt)
1043 {
1044         if(acc[k] == 0)
1045                 return;
1046         if(receivebit(h) != 0)
1047                 if(acc[k] < 0)
1048                         acc[k] -= Pt;
1049                 else
1050                         acc[k] += Pt;
1051 }
1052
1053 static
1054 void
1055 progressiveacinc(Header *h, int comp, int Al)
1056 {
1057         int Ns, i, z, k, Ss, Se, Ta, **ac, H, V;
1058         int ri, mcu, nacross, ndown, nhor, nver, eobrun, nzeros, pending, x, y, tmcu, blockno, q, nmcu;
1059         Huffman *acht;
1060         int *qt, rrrr, ssss, *acc, rs;
1061         uchar *ss;
1062
1063         ss = h->ss;
1064         Ns = ss[0];
1065         if(Ns != 1)
1066                 jpgerror(h, "ReadJPG: illegal Ns>1 in progressive AC scan");
1067         Ss = ss[1+2];
1068         Se = ss[2+2];
1069         H = h->comp[comp].H;
1070         V = h->comp[comp].V;
1071
1072         nacross = h->nacross*H;
1073         ndown = h->ndown*V;
1074         q = 8*h->Hmax/H;
1075         nhor = (h->X+q-1)/q;
1076         q = 8*h->Vmax/V;
1077         nver = (h->Y+q-1)/q;
1078
1079         /* initialize data structures */
1080         h->cnt = 0;
1081         h->sr = 0;
1082         h->peek = -1;
1083         nibbles(ss[1+1], &z, &Ta);      /* z is thrown away */
1084         ri = h->ri;
1085
1086         eobrun = 0;
1087         ac = h->accoeff[comp];
1088         acht = &h->acht[Ta];
1089         qt = h->qt[h->comp[comp].Tq];
1090         nmcu = nacross*ndown;
1091         mcu = 0;
1092         pending = 0;
1093         nzeros = -1;
1094         for(y=0; y<nver; y++){
1095                 for(x=0; x<nhor; x++){
1096                         /* Figure G-7 */
1097
1098                         /*  arrange blockno to be in same sequence as original scan calculation. */
1099                         tmcu = x/H + (nacross/H)*(y/V);
1100                         blockno = tmcu*H*V + H*(y%V) + x%H;
1101                         acc = ac[blockno];
1102                         if(eobrun > 0){
1103                                 if(nzeros > 0)
1104                                         jpgerror(h, "ReadJPG: zeros pending at block start");
1105                                 for(k=Ss; k<=Se; k++)
1106                                         increment(h, acc, k, qt[k]<<Al);
1107                                 --eobrun;
1108                                 continue;
1109                         }
1110
1111                         for(k=Ss; k<=Se; ){
1112                                 if(nzeros >= 0){
1113                                         if(acc[k] != 0)
1114                                                 increment(h, acc, k, qt[k]<<Al);
1115                                         else if(nzeros-- == 0)
1116                                                 acc[k] = pending;
1117                                         k++;
1118                                         continue;
1119                                 }
1120                                 rs = decode(h, acht);
1121                                 nibbles(rs, &rrrr, &ssss);
1122                                 if(ssss == 0){
1123                                         if(rrrr < 15){
1124                                                 eobrun = 0;
1125                                                 if(rrrr > 0)
1126                                                         eobrun = receiveEOB(h, rrrr)-1;
1127                                                 while(k <= Se){
1128                                                         increment(h, acc, k, qt[k]<<Al);
1129                                                         k++;
1130                                                 }
1131                                                 break;
1132                                         }
1133                                         for(i=0; i<16; k++){
1134                                                 increment(h, acc, k, qt[k]<<Al);
1135                                                 if(acc[k] == 0)
1136                                                         i++;
1137                                         }
1138                                         continue;
1139                                 }else if(ssss != 1)
1140                                         jpgerror(h, "ReadJPG: ssss!=1 in progressive increment");
1141                                 nzeros = rrrr;
1142                                 pending = receivebit(h);
1143                                 if(pending == 0)
1144                                         pending = -1;
1145                                 pending *= qt[k]<<Al;
1146                         }
1147                 }
1148
1149                 /* process restart marker, if present */
1150                 mcu++;
1151                 if(ri>0 && mcu<nmcu && mcu%ri==0){
1152                         restart(h, mcu);
1153                         eobrun = 0;
1154                         nzeros = -1;
1155                 }
1156         }
1157 }
1158
1159 static
1160 void
1161 progressivescan(Header *h, int colorspace)
1162 {
1163         uchar *ss;
1164         int Ns, Ss, Ah, Al, c, comp, i;
1165
1166         if(h->dccoeff[0] == nil)
1167                 progressiveinit(h, colorspace);
1168
1169         ss = h->ss;
1170         Ns = ss[0];
1171         Ss = ss[1+2*Ns];
1172         nibbles(ss[3+2*Ns], &Ah, &Al);
1173         c = ss[1];
1174         comp = -1;
1175         for(i=0; i<h->Nf; i++)
1176                 if(h->comp[i].C == c)
1177                         comp = i;
1178         if(comp == -1)
1179                 jpgerror(h, "ReadJPG: bad component index in scan header");
1180
1181         if(Ss == 0){
1182                 progressivedc(h, comp, Ah, Al);
1183                 return;
1184         }
1185         if(Ah == 0){
1186                 progressiveac(h, comp, Al);
1187                 return;
1188         }
1189         progressiveacinc(h, comp, Al);
1190 }
1191
1192 enum {
1193         c1 = 2871,      /* 1.402 * 2048 */
1194         c2 = 705,               /* 0.34414 * 2048 */
1195         c3 = 1463,      /* 0.71414 * 2048 */
1196         c4 = 3629,      /* 1.772 * 2048 */
1197 };
1198
1199 static
1200 void
1201 colormap1(Header *h, int colorspace, Rawimage *image, int data[8*8], int mcu, int nacross)
1202 {
1203         uchar *pic;
1204         int x, y, dx, dy, minx, miny;
1205         int r, k, pici;
1206
1207         USED(colorspace);
1208         pic = image->chans[0];
1209         minx = 8*(mcu%nacross);
1210         dx = 8;
1211         if(minx+dx > h->X)
1212                 dx = h->X-minx;
1213         miny = 8*(mcu/nacross);
1214         dy = 8;
1215         if(miny+dy > h->Y)
1216                 dy = h->Y-miny;
1217         pici = miny*h->X+minx;
1218         k = 0;
1219         for(y=0; y<dy; y++){
1220                 for(x=0; x<dx; x++){
1221                         r = clamp[(data[k+x]+128)+CLAMPOFF];
1222                         pic[pici+x] = r;
1223                 }
1224                 pici += h->X;
1225                 k += 8;
1226         }
1227 }
1228
1229 static
1230 void
1231 colormapall1(Header *h, int colorspace, Rawimage *image, int data0[8*8], int data1[8*8], int data2[8*8], int mcu, int nacross)
1232 {
1233         uchar *rpic, *gpic, *bpic, *rp, *gp, *bp;
1234         int *p0, *p1, *p2;
1235         int x, y, dx, dy, minx, miny;
1236         int r, g, b, k, pici;
1237         int Y, Cr, Cb;
1238
1239         rpic = image->chans[0];
1240         gpic = image->chans[1];
1241         bpic = image->chans[2];
1242         minx = 8*(mcu%nacross);
1243         dx = 8;
1244         if(minx+dx > h->X)
1245                 dx = h->X-minx;
1246         miny = 8*(mcu/nacross);
1247         dy = 8;
1248         if(miny+dy > h->Y)
1249                 dy = h->Y-miny;
1250         pici = miny*h->X+minx;
1251         k = 0;
1252         for(y=0; y<dy; y++){
1253                 p0 = data0+k;
1254                 p1 = data1+k;
1255                 p2 = data2+k;
1256                 rp = rpic+pici;
1257                 gp = gpic+pici;
1258                 bp = bpic+pici;
1259                 if(colorspace == CYCbCr)
1260                         for(x=0; x<dx; x++){
1261                                 *rp++ = clamp[*p0++ + 128 + CLAMPOFF];
1262                                 *gp++ = clamp[*p1++ + 128 + CLAMPOFF];
1263                                 *bp++ = clamp[*p2++ + 128 + CLAMPOFF];
1264                         }
1265                 else
1266                         for(x=0; x<dx; x++){
1267                                 Y = (*p0++ + 128) << 11;
1268                                 Cb = *p1++;
1269                                 Cr = *p2++;
1270                                 r = Y+c1*Cr;
1271                                 g = Y-c2*Cb-c3*Cr;
1272                                 b = Y+c4*Cb;
1273                                 *rp++ = clamp[(r>>11)+CLAMPOFF];
1274                                 *gp++ = clamp[(g>>11)+CLAMPOFF];
1275                                 *bp++ = clamp[(b>>11)+CLAMPOFF];
1276                         }
1277                 pici += h->X;
1278                 k += 8;
1279         }
1280 }
1281
1282 static
1283 void
1284 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)
1285 {
1286         uchar *rpic, *gpic, *bpic;
1287         int x, y, dx, dy, minx, miny;
1288         int r, g, b, pici, H0, H1, H2;
1289         int t, b0, b1, b2, y0, y1, y2, x0, x1, x2;
1290         int Y, Cr, Cb;
1291
1292         rpic = image->chans[0];
1293         gpic = image->chans[1];
1294         bpic = image->chans[2];
1295         minx = 8*Hmax*(mcu%nacross);
1296         dx = 8*Hmax;
1297         if(minx+dx > h->X)
1298                 dx = h->X-minx;
1299         miny = 8*Vmax*(mcu/nacross);
1300         dy = 8*Vmax;
1301         if(miny+dy > h->Y)
1302                 dy = h->Y-miny;
1303         pici = miny*h->X+minx;
1304         H0 = H[0];
1305         H1 = H[1];
1306         H2 = H[2];
1307         for(y=0; y<dy; y++){
1308                 t = y*V[0];
1309                 b0 = H0*(t/(8*Vmax));
1310                 y0 = 8*((t/Vmax)&7);
1311                 t = y*V[1];
1312                 b1 = H1*(t/(8*Vmax));
1313                 y1 = 8*((t/Vmax)&7);
1314                 t = y*V[2];
1315                 b2 = H2*(t/(8*Vmax));
1316                 y2 = 8*((t/Vmax)&7);
1317                 x0 = 0;
1318                 x1 = 0;
1319                 x2 = 0;
1320                 for(x=0; x<dx; x++){
1321                         if(colorspace == CYCbCr){
1322                                 rpic[pici+x] = clamp[data0[b0][y0+x0++*H0/Hmax] + 128 + CLAMPOFF];
1323                                 gpic[pici+x] = clamp[data1[b1][y1+x1++*H1/Hmax] + 128 + CLAMPOFF];
1324                                 bpic[pici+x] = clamp[data2[b2][y2+x2++*H2/Hmax] + 128 + CLAMPOFF];
1325                         }else{
1326                                 Y = (data0[b0][y0+x0++*H0/Hmax]+128)<<11;
1327                                 Cb = data1[b1][y1+x1++*H1/Hmax];
1328                                 Cr = data2[b2][y2+x2++*H2/Hmax];
1329                                 r = Y+c1*Cr;
1330                                 g = Y-c2*Cb-c3*Cr;
1331                                 b = Y+c4*Cb;
1332                                 rpic[pici+x] = clamp[(r>>11)+CLAMPOFF];
1333                                 gpic[pici+x] = clamp[(g>>11)+CLAMPOFF];
1334                                 bpic[pici+x] = clamp[(b>>11)+CLAMPOFF];
1335                         }
1336                         if(x0*H0/Hmax >= 8){
1337                                 x0 = 0;
1338                                 b0++;
1339                         }
1340                         if(x1*H1/Hmax >= 8){
1341                                 x1 = 0;
1342                                 b1++;
1343                         }
1344                         if(x2*H2/Hmax >= 8){
1345                                 x2 = 0;
1346                                 b2++;
1347                         }
1348                 }
1349                 pici += h->X;
1350         }
1351 }
1352
1353 /*
1354  * decode next 8-bit value from entropy-coded input.  chart F-26
1355  */
1356 static
1357 int
1358 decode(Header *h, Huffman *t)
1359 {
1360         int code, v, cnt, sr, i;
1361         int *maxcode;
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         cnt = h->cnt;
1378         sr = h->sr;
1379         code <<= 1;
1380         for(i = 9; i<17; i++){
1381                 cnt--;
1382                 if(sr & (1<<cnt))
1383                         code |= 1;
1384                 if(code <= maxcode[i])
1385                         break;
1386                 code <<= 1;
1387                 if(cnt == 0){
1388                         sr = nextbyte(h, 0);
1389                         cnt = 8;
1390                 }
1391         }
1392         if(i >= 17){
1393                 /* bad code */
1394                 code = 0;
1395                 i = 0;
1396         }
1397         h->cnt = cnt;
1398         return t->val[t->valptr[i]+(code-t->mincode[i])];
1399 }
1400
1401 /*
1402  * load next byte of input
1403  */
1404 static
1405 int
1406 nextbyte(Header *h, int marker)
1407 {
1408         int b, b2;
1409
1410         if(h->peek >= 0){
1411                 b = h->peek;
1412                 h->peek = -1;
1413         }else{
1414                 b = Bgetc(h->fd);
1415                 if(b == Beof)
1416                         jpgerror(h, "truncated file");
1417                 b &= 0xFF;
1418         }
1419
1420         if(b == 0xFF){
1421                 if(marker)
1422                         return b;
1423                 b2 = Bgetc(h->fd);
1424                 if(b2 != 0){
1425                         if(b2 == Beof)
1426                                 jpgerror(h, "truncated file");
1427                         b2 &= 0xFF;
1428                         if(b2 == DNL)
1429                                 jpgerror(h, "ReadJPG: DNL marker unimplemented");
1430                         /* decoder is reading into marker; satisfy it and restore state */
1431                         Bungetc(h->fd);
1432                         h->peek = b;
1433                 }
1434         }
1435         h->cnt += 8;
1436         h->sr = (h->sr<<8) | b;
1437         return b;
1438 }
1439
1440 /*
1441  * return next s bits of input, MSB first, and level shift it
1442  */
1443 static
1444 int
1445 receive(Header *h, int s)
1446 {
1447         int v, m;
1448
1449         while(h->cnt < s)
1450                 nextbyte(h, 0);
1451         h->cnt -= s;
1452         v = h->sr >> h->cnt;
1453         m = (1<<s);
1454         v &= m-1;
1455         /* level shift */
1456         if(v < (m>>1))
1457                 v += ~(m-1)+1;
1458         return v;
1459 }
1460
1461 /*
1462  * return next s bits of input, decode as EOB
1463  */
1464 static
1465 int
1466 receiveEOB(Header *h, int s)
1467 {
1468         int v, m;
1469
1470         while(h->cnt < s)
1471                 nextbyte(h, 0);
1472         h->cnt -= s;
1473         v = h->sr >> h->cnt;
1474         m = (1<<s);
1475         v &= m-1;
1476         /* level shift */
1477         v += m;
1478         return v;
1479 }
1480
1481 /*
1482  * return next bit of input
1483  */
1484 static
1485 int
1486 receivebit(Header *h)
1487 {
1488         if(h->cnt < 1)
1489                 nextbyte(h, 0);
1490         h->cnt--;
1491         return (h->sr >> h->cnt) & 1;
1492 }
1493
1494 /*
1495  *  Scaled integer implementation.
1496  *  inverse two dimensional DCT, Chen-Wang algorithm
1497  * (IEEE ASSP-32, pp. 803-816, Aug. 1984)
1498  * 32-bit integer arithmetic (8 bit coefficients)
1499  * 11 mults, 29 adds per DCT
1500  *
1501  * coefficients extended to 12 bit for IEEE1180-1990 compliance
1502  */
1503
1504 enum {
1505         W1              = 2841, /* 2048*sqrt(2)*cos(1*pi/16)*/
1506         W2              = 2676, /* 2048*sqrt(2)*cos(2*pi/16)*/
1507         W3              = 2408, /* 2048*sqrt(2)*cos(3*pi/16)*/
1508         W5              = 1609, /* 2048*sqrt(2)*cos(5*pi/16)*/
1509         W6              = 1108, /* 2048*sqrt(2)*cos(6*pi/16)*/
1510         W7              = 565,  /* 2048*sqrt(2)*cos(7*pi/16)*/
1511
1512         W1pW7   = 3406, /* W1+W7*/
1513         W1mW7   = 2276, /* W1-W7*/
1514         W3pW5   = 4017, /* W3+W5*/
1515         W3mW5   = 799,  /* W3-W5*/
1516         W2pW6   = 3784, /* W2+W6*/
1517         W2mW6   = 1567, /* W2-W6*/
1518
1519         R2              = 181   /* 256/sqrt(2)*/
1520 };
1521
1522 static
1523 void
1524 idct(int b[8*8])
1525 {
1526         int x, y, eighty, v;
1527         int x0, x1, x2, x3, x4, x5, x6, x7, x8;
1528         int *p;
1529
1530         /* transform horizontally*/
1531         for(y=0; y<8; y++){
1532                 eighty = y<<3;
1533                 /* if all non-DC components are zero, just propagate the DC term*/
1534                 p = b+eighty;
1535                 if(p[1]==0)
1536                 if(p[2]==0 && p[3]==0)
1537                 if(p[4]==0 && p[5]==0)
1538                 if(p[6]==0 && p[7]==0){
1539                         v = p[0]<<3;
1540                         p[0] = v;
1541                         p[1] = v;
1542                         p[2] = v;
1543                         p[3] = v;
1544                         p[4] = v;
1545                         p[5] = v;
1546                         p[6] = v;
1547                         p[7] = v;
1548                         continue;
1549                 }
1550                 /* prescale*/
1551                 x0 = (p[0]<<11)+128;
1552                 x1 = p[4]<<11;
1553                 x2 = p[6];
1554                 x3 = p[2];
1555                 x4 = p[1];
1556                 x5 = p[7];
1557                 x6 = p[5];
1558                 x7 = p[3];
1559                 /* first stage*/
1560                 x8 = W7*(x4+x5);
1561                 x4 = x8 + W1mW7*x4;
1562                 x5 = x8 - W1pW7*x5;
1563                 x8 = W3*(x6+x7);
1564                 x6 = x8 - W3mW5*x6;
1565                 x7 = x8 - W3pW5*x7;
1566                 /* second stage*/
1567                 x8 = x0 + x1;
1568                 x0 -= x1;
1569                 x1 = W6*(x3+x2);
1570                 x2 = x1 - W2pW6*x2;
1571                 x3 = x1 + W2mW6*x3;
1572                 x1 = x4 + x6;
1573                 x4 -= x6;
1574                 x6 = x5 + x7;
1575                 x5 -= x7;
1576                 /* third stage*/
1577                 x7 = x8 + x3;
1578                 x8 -= x3;
1579                 x3 = x0 + x2;
1580                 x0 -= x2;
1581                 x2 = (R2*(x4+x5)+128)>>8;
1582                 x4 = (R2*(x4-x5)+128)>>8;
1583                 /* fourth stage*/
1584                 p[0] = (x7+x1)>>8;
1585                 p[1] = (x3+x2)>>8;
1586                 p[2] = (x0+x4)>>8;
1587                 p[3] = (x8+x6)>>8;
1588                 p[4] = (x8-x6)>>8;
1589                 p[5] = (x0-x4)>>8;
1590                 p[6] = (x3-x2)>>8;
1591                 p[7] = (x7-x1)>>8;
1592         }
1593         /* transform vertically*/
1594         for(x=0; x<8; x++){
1595                 /* if all non-DC components are zero, just propagate the DC term*/
1596                 p = b+x;
1597                 if(p[8*1]==0)
1598                 if(p[8*2]==0 && p[8*3]==0)
1599                 if(p[8*4]==0 && p[8*5]==0)
1600                 if(p[8*6]==0 && p[8*7]==0){
1601                         v = (p[8*0]+32)>>6;
1602                         p[8*0] = v;
1603                         p[8*1] = v;
1604                         p[8*2] = v;
1605                         p[8*3] = v;
1606                         p[8*4] = v;
1607                         p[8*5] = v;
1608                         p[8*6] = v;
1609                         p[8*7] = v;
1610                         continue;
1611                 }
1612                 /* prescale*/
1613                 x0 = (p[8*0]<<8)+8192;
1614                 x1 = p[8*4]<<8;
1615                 x2 = p[8*6];
1616                 x3 = p[8*2];
1617                 x4 = p[8*1];
1618                 x5 = p[8*7];
1619                 x6 = p[8*5];
1620                 x7 = p[8*3];
1621                 /* first stage*/
1622                 x8 = W7*(x4+x5) + 4;
1623                 x4 = (x8+W1mW7*x4)>>3;
1624                 x5 = (x8-W1pW7*x5)>>3;
1625                 x8 = W3*(x6+x7) + 4;
1626                 x6 = (x8-W3mW5*x6)>>3;
1627                 x7 = (x8-W3pW5*x7)>>3;
1628                 /* second stage*/
1629                 x8 = x0 + x1;
1630                 x0 -= x1;
1631                 x1 = W6*(x3+x2) + 4;
1632                 x2 = (x1-W2pW6*x2)>>3;
1633                 x3 = (x1+W2mW6*x3)>>3;
1634                 x1 = x4 + x6;
1635                 x4 -= x6;
1636                 x6 = x5 + x7;
1637                 x5 -= x7;
1638                 /* third stage*/
1639                 x7 = x8 + x3;
1640                 x8 -= x3;
1641                 x3 = x0 + x2;
1642                 x0 -= x2;
1643                 x2 = (R2*(x4+x5)+128)>>8;
1644                 x4 = (R2*(x4-x5)+128)>>8;
1645                 /* fourth stage*/
1646                 p[8*0] = (x7+x1)>>14;
1647                 p[8*1] = (x3+x2)>>14;
1648                 p[8*2] = (x0+x4)>>14;
1649                 p[8*3] = (x8+x6)>>14;
1650                 p[8*4] = (x8-x6)>>14;
1651                 p[8*5] = (x0-x4)>>14;
1652                 p[8*6] = (x3-x2)>>14;
1653                 p[8*7] = (x7-x1)>>14;
1654         }
1655 }