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         if(h->image != nil)
 230                 fprint(2, "jpg: partial image: %s\n", h->err);
 231         werrstr("%s", h->err);
 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 = h->image;
 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 1 or 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)
 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 i, Td[3], DC[3], blockno[3];
 901
 902         ss= h->ss;
 903         Ns = ss[0];
 904         if(Ns!=1 && Ns!=h->Nf)
 905                 jpgerror(h, "ReadJPG: can't handle progressive with Ns!=1 and Nf!=Ns in DC scan");
 906
 907         /* initialize data structures */
 908         h->cnt = 0;
 909         h->sr = 0;
 910         h->peek = -1;
 911
 912         for(i=0; i<Ns; i++){
 913                 /*
 914                  * JPEG requires scan components to be in same order as in frame,
 915                  * so if both have 3 we know scan is Y Cb Cr and there's no need to
 916                  * reorder
 917                  */
 918                 nibbles(ss[2+2*i], &Td[i], &z); /* z is ignored */
 919                 DC[i] = 0;
 920         }
 921
 922         ri = h->ri;
 923
 924         nmcu = h->nacross*h->ndown;
 925         memset(blockno, 0, sizeof blockno);
 926         for(mcu=0; mcu<nmcu; ){
 927                 for(i=0; i<Ns; i++){
 928                         if(Ns != 1) comp = i;
 929
 930                         dcht = &h->dcht[Td[i]];
 931                         qt = h->qt[h->comp[comp].Tq][0];
 932                         dc = h->dccoeff[comp];
 933                         bn = blockno[i];
 934
 935                         for(block=0; block<h->nblock[comp]; block++){
 936                                 if(Ah == 0){
 937                                         t = decode(h, dcht);
 938                                         diff = receive(h, t);
 939                                         DC[i] += diff;
 940                                         dc[bn] = qt*DC[i]<<Al;
 941                                 }else
 942                                         dc[bn] |= qt*receivebit(h)<<Al;
 943                                 bn++;
 944                         }
 945                         blockno[i] = bn;
 946                 }
 947
 948                 /* process restart marker, if present */
 949                 mcu++;
 950                 if(ri>0 && mcu<nmcu && mcu%ri==0){
 951                         restart(h, mcu);
 952                         for(i=0; i<Ns; i++)
 953                                 DC[i] = 0;
 954                 }
 955         }
 956 }
 957
 958 static
 959 void
 960 progressiveac(Header *h, int comp, int Al)
 961 {
 962         int Ns, Ss, Se, z, k, eobrun, x, y, nver, tmcu, blockno, *acc, rs;
 963         int ri, mcu, nacross, ndown, nmcu, nhor;
 964         Huffman *acht;
 965         int *qt, rrrr, ssss, q;
 966         uchar *ss;
 967         int Ta, H, V;
 968
 969         ss = h->ss;
 970         Ns = ss[0];
 971         if(Ns != 1)
 972                 jpgerror(h, "ReadJPG: illegal Ns>1 in progressive AC scan");
 973         Ss = ss[1+2];
 974         Se = ss[2+2];
 975         H = h->comp[comp].H;
 976         V = h->comp[comp].V;
 977
 978         nacross = h->nacross*H;
 979         ndown = h->ndown*V;
 980         q = 8*h->Hmax/H;
 981         nhor = (h->X+q-1)/q;
 982         q = 8*h->Vmax/V;
 983         nver = (h->Y+q-1)/q;
 984
 985         /* initialize data structures */
 986         h->cnt = 0;
 987         h->sr = 0;
 988         h->peek = -1;
 989         nibbles(ss[1+1], &z, &Ta);      /* z is thrown away */
 990
 991         ri = h->ri;
 992
 993         eobrun = 0;
 994         acht = &h->acht[Ta];
 995         qt = h->qt[h->comp[comp].Tq];
 996         nmcu = nacross*ndown;
 997         mcu = 0;
 998         for(y=0; y<nver; y++){
 999                 for(x=0; x<nhor; x++){
1000                         /* Figure G-3  */
1001                         if(eobrun > 0){
1002                                 --eobrun;
1003                                 continue;
1004                         }
1005
1006                         /* arrange blockno to be in same sequence as original scan calculation. */
1007                         tmcu = x/H + (nacross/H)*(y/V);
1008                         blockno = tmcu*H*V + H*(y%V) + x%H;
1009                         acc = h->accoeff[comp][blockno];
1010                         k = Ss;
1011                         do {
1012                                 rs = decode(h, acht);
1013                                 /* XXX remove rrrr ssss as in baselinescan */
1014                                 nibbles(rs, &rrrr, &ssss);
1015                                 if(ssss == 0){
1016                                         if(rrrr < 15){
1017                                                 eobrun = 0;
1018                                                 if(rrrr > 0)
1019                                                         eobrun = receiveEOB(h, rrrr)-1;
1020                                                 break;
1021                                         }
1022                                         k += 16;
1023                                 }else{
1024                                         z = receive(h, ssss);
1025                                         k += rrrr;
1026                                         if(k > Se)
1027                                                 break;
1028                                         acc[k] = z*qt[k]<<Al;
1029                                         k++;
1030                                 }
1031                         } while(k <= Se);
1032                 }
1033
1034                 /* process restart marker, if present */
1035                 mcu++;
1036                 if(ri>0 && mcu<nmcu && mcu%ri==0){
1037                         restart(h, mcu);
1038                         eobrun = 0;
1039                 }
1040         }
1041 }
1042
1043 static
1044 void
1045 increment(Header *h, int acc[], int k, int Pt)
1046 {
1047         if(acc[k] == 0)
1048                 return;
1049         if(receivebit(h) != 0)
1050                 if(acc[k] < 0)
1051                         acc[k] -= Pt;
1052                 else
1053                         acc[k] += Pt;
1054 }
1055
1056 static
1057 void
1058 progressiveacinc(Header *h, int comp, int Al)
1059 {
1060         int Ns, i, z, k, Ss, Se, Ta, **ac, H, V;
1061         int ri, mcu, nacross, ndown, nhor, nver, eobrun, nzeros, pending, x, y, tmcu, blockno, q, nmcu;
1062         Huffman *acht;
1063         int *qt, rrrr, ssss, *acc, rs;
1064         uchar *ss;
1065
1066         ss = h->ss;
1067         Ns = ss[0];
1068         if(Ns != 1)
1069                 jpgerror(h, "ReadJPG: illegal Ns>1 in progressive AC scan");
1070         Ss = ss[1+2];
1071         Se = ss[2+2];
1072         H = h->comp[comp].H;
1073         V = h->comp[comp].V;
1074
1075         nacross = h->nacross*H;
1076         ndown = h->ndown*V;
1077         q = 8*h->Hmax/H;
1078         nhor = (h->X+q-1)/q;
1079         q = 8*h->Vmax/V;
1080         nver = (h->Y+q-1)/q;
1081
1082         /* initialize data structures */
1083         h->cnt = 0;
1084         h->sr = 0;
1085         h->peek = -1;
1086         nibbles(ss[1+1], &z, &Ta);      /* z is thrown away */
1087         ri = h->ri;
1088
1089         eobrun = 0;
1090         ac = h->accoeff[comp];
1091         acht = &h->acht[Ta];
1092         qt = h->qt[h->comp[comp].Tq];
1093         nmcu = nacross*ndown;
1094         mcu = 0;
1095         pending = 0;
1096         nzeros = -1;
1097         for(y=0; y<nver; y++){
1098                 for(x=0; x<nhor; x++){
1099                         /* Figure G-7 */
1100
1101                         /*  arrange blockno to be in same sequence as original scan calculation. */
1102                         tmcu = x/H + (nacross/H)*(y/V);
1103                         blockno = tmcu*H*V + H*(y%V) + x%H;
1104                         acc = ac[blockno];
1105                         if(eobrun > 0){
1106                                 if(nzeros > 0)
1107                                         jpgerror(h, "ReadJPG: zeros pending at block start");
1108                                 for(k=Ss; k<=Se; k++)
1109                                         increment(h, acc, k, qt[k]<<Al);
1110                                 --eobrun;
1111                                 continue;
1112                         }
1113
1114                         for(k=Ss; k<=Se; ){
1115                                 if(nzeros >= 0){
1116                                         if(acc[k] != 0)
1117                                                 increment(h, acc, k, qt[k]<<Al);
1118                                         else if(nzeros-- == 0)
1119                                                 acc[k] = pending;
1120                                         k++;
1121                                         continue;
1122                                 }
1123                                 rs = decode(h, acht);
1124                                 nibbles(rs, &rrrr, &ssss);
1125                                 if(ssss == 0){
1126                                         if(rrrr < 15){
1127                                                 eobrun = 0;
1128                                                 if(rrrr > 0)
1129                                                         eobrun = receiveEOB(h, rrrr)-1;
1130                                                 while(k <= Se){
1131                                                         increment(h, acc, k, qt[k]<<Al);
1132                                                         k++;
1133                                                 }
1134                                                 break;
1135                                         }
1136                                         for(i=0; i<16; k++){
1137                                                 increment(h, acc, k, qt[k]<<Al);
1138                                                 if(acc[k] == 0)
1139                                                         i++;
1140                                         }
1141                                         continue;
1142                                 }else if(ssss != 1)
1143                                         jpgerror(h, "ReadJPG: ssss!=1 in progressive increment");
1144                                 nzeros = rrrr;
1145                                 pending = receivebit(h);
1146                                 if(pending == 0)
1147                                         pending = -1;
1148                                 pending *= qt[k]<<Al;
1149                         }
1150                 }
1151
1152                 /* process restart marker, if present */
1153                 mcu++;
1154                 if(ri>0 && mcu<nmcu && mcu%ri==0){
1155                         restart(h, mcu);
1156                         eobrun = 0;
1157                         nzeros = -1;
1158                 }
1159         }
1160 }
1161
1162 static
1163 void
1164 progressivescan(Header *h, int colorspace)
1165 {
1166         uchar *ss;
1167         int Ns, Ss, Ah, Al, c, comp, i;
1168
1169         if(h->dccoeff[0] == nil)
1170                 progressiveinit(h, colorspace);
1171
1172         ss = h->ss;
1173         Ns = ss[0];
1174         Ss = ss[1+2*Ns];
1175         nibbles(ss[3+2*Ns], &Ah, &Al);
1176         c = ss[1];
1177         comp = -1;
1178         for(i=0; i<h->Nf; i++)
1179                 if(h->comp[i].C == c)
1180                         comp = i;
1181         if(comp == -1)
1182                 jpgerror(h, "ReadJPG: bad component index in scan header");
1183
1184         if(Ss == 0){
1185                 progressivedc(h, comp, Ah, Al);
1186                 return;
1187         }
1188         if(Ah == 0){
1189                 progressiveac(h, comp, Al);
1190                 return;
1191         }
1192         progressiveacinc(h, comp, Al);
1193 }
1194
1195 enum {
1196         c1 = 2871,      /* 1.402 * 2048 */
1197         c2 = 705,               /* 0.34414 * 2048 */
1198         c3 = 1463,      /* 0.71414 * 2048 */
1199         c4 = 3629,      /* 1.772 * 2048 */
1200 };
1201
1202 static
1203 void
1204 colormap1(Header *h, int colorspace, Rawimage *image, int data[8*8], int mcu, int nacross)
1205 {
1206         uchar *pic;
1207         int x, y, dx, dy, minx, miny;
1208         int r, k, pici;
1209
1210         USED(colorspace);
1211         pic = image->chans[0];
1212         minx = 8*(mcu%nacross);
1213         dx = 8;
1214         if(minx+dx > h->X)
1215                 dx = h->X-minx;
1216         miny = 8*(mcu/nacross);
1217         dy = 8;
1218         if(miny+dy > h->Y)
1219                 dy = h->Y-miny;
1220         pici = miny*h->X+minx;
1221         k = 0;
1222         for(y=0; y<dy; y++){
1223                 for(x=0; x<dx; x++){
1224                         r = clamp[(data[k+x]+128)+CLAMPOFF];
1225                         pic[pici+x] = r;
1226                 }
1227                 pici += h->X;
1228                 k += 8;
1229         }
1230 }
1231
1232 static
1233 void
1234 colormapall1(Header *h, int colorspace, Rawimage *image, int data0[8*8], int data1[8*8], int data2[8*8], int mcu, int nacross)
1235 {
1236         uchar *rpic, *gpic, *bpic, *rp, *gp, *bp;
1237         int *p0, *p1, *p2;
1238         int x, y, dx, dy, minx, miny;
1239         int r, g, b, k, pici;
1240         int Y, Cr, Cb;
1241
1242         rpic = image->chans[0];
1243         gpic = image->chans[1];
1244         bpic = image->chans[2];
1245         minx = 8*(mcu%nacross);
1246         dx = 8;
1247         if(minx+dx > h->X)
1248                 dx = h->X-minx;
1249         miny = 8*(mcu/nacross);
1250         dy = 8;
1251         if(miny+dy > h->Y)
1252                 dy = h->Y-miny;
1253         pici = miny*h->X+minx;
1254         k = 0;
1255         for(y=0; y<dy; y++){
1256                 p0 = data0+k;
1257                 p1 = data1+k;
1258                 p2 = data2+k;
1259                 rp = rpic+pici;
1260                 gp = gpic+pici;
1261                 bp = bpic+pici;
1262                 if(colorspace == CYCbCr)
1263                         for(x=0; x<dx; x++){
1264                                 *rp++ = clamp[*p0++ + 128 + CLAMPOFF];
1265                                 *gp++ = clamp[*p1++ + 128 + CLAMPOFF];
1266                                 *bp++ = clamp[*p2++ + 128 + CLAMPOFF];
1267                         }
1268                 else
1269                         for(x=0; x<dx; x++){
1270                                 Y = (*p0++ + 128) << 11;
1271                                 Cb = *p1++;
1272                                 Cr = *p2++;
1273                                 r = Y+c1*Cr;
1274                                 g = Y-c2*Cb-c3*Cr;
1275                                 b = Y+c4*Cb;
1276                                 *rp++ = clamp[(r>>11)+CLAMPOFF];
1277                                 *gp++ = clamp[(g>>11)+CLAMPOFF];
1278                                 *bp++ = clamp[(b>>11)+CLAMPOFF];
1279                         }
1280                 pici += h->X;
1281                 k += 8;
1282         }
1283 }
1284
1285 static
1286 void
1287 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)
1288 {
1289         uchar *rpic, *gpic, *bpic;
1290         int x, y, dx, dy, minx, miny;
1291         int r, g, b, pici, H0, H1, H2;
1292         int t, b0, b1, b2, y0, y1, y2, x0, x1, x2;
1293         int Y, Cr, Cb;
1294
1295         rpic = image->chans[0];
1296         gpic = image->chans[1];
1297         bpic = image->chans[2];
1298         minx = 8*Hmax*(mcu%nacross);
1299         dx = 8*Hmax;
1300         if(minx+dx > h->X)
1301                 dx = h->X-minx;
1302         miny = 8*Vmax*(mcu/nacross);
1303         dy = 8*Vmax;
1304         if(miny+dy > h->Y)
1305                 dy = h->Y-miny;
1306         pici = miny*h->X+minx;
1307         H0 = H[0];
1308         H1 = H[1];
1309         H2 = H[2];
1310         for(y=0; y<dy; y++){
1311                 t = y*V[0];
1312                 b0 = H0*(t/(8*Vmax));
1313                 y0 = 8*((t/Vmax)&7);
1314                 t = y*V[1];
1315                 b1 = H1*(t/(8*Vmax));
1316                 y1 = 8*((t/Vmax)&7);
1317                 t = y*V[2];
1318                 b2 = H2*(t/(8*Vmax));
1319                 y2 = 8*((t/Vmax)&7);
1320                 x0 = 0;
1321                 x1 = 0;
1322                 x2 = 0;
1323                 for(x=0; x<dx; x++){
1324                         if(colorspace == CYCbCr){
1325                                 rpic[pici+x] = clamp[data0[b0][y0+x0++*H0/Hmax] + 128 + CLAMPOFF];
1326                                 gpic[pici+x] = clamp[data1[b1][y1+x1++*H1/Hmax] + 128 + CLAMPOFF];
1327                                 bpic[pici+x] = clamp[data2[b2][y2+x2++*H2/Hmax] + 128 + CLAMPOFF];
1328                         }else{
1329                                 Y = (data0[b0][y0+x0++*H0/Hmax]+128)<<11;
1330                                 Cb = data1[b1][y1+x1++*H1/Hmax];
1331                                 Cr = data2[b2][y2+x2++*H2/Hmax];
1332                                 r = Y+c1*Cr;
1333                                 g = Y-c2*Cb-c3*Cr;
1334                                 b = Y+c4*Cb;
1335                                 rpic[pici+x] = clamp[(r>>11)+CLAMPOFF];
1336                                 gpic[pici+x] = clamp[(g>>11)+CLAMPOFF];
1337                                 bpic[pici+x] = clamp[(b>>11)+CLAMPOFF];
1338                         }
1339                         if(x0*H0/Hmax >= 8){
1340                                 x0 = 0;
1341                                 b0++;
1342                         }
1343                         if(x1*H1/Hmax >= 8){
1344                                 x1 = 0;
1345                                 b1++;
1346                         }
1347                         if(x2*H2/Hmax >= 8){
1348                                 x2 = 0;
1349                                 b2++;
1350                         }
1351                 }
1352                 pici += h->X;
1353         }
1354 }
1355
1356 /*
1357  * decode next 8-bit value from entropy-coded input.  chart F-26
1358  */
1359 static
1360 int
1361 decode(Header *h, Huffman *t)
1362 {
1363         int code, v, cnt, sr, i;
1364         int *maxcode;
1365
1366         maxcode = t->maxcode;
1367         if(h->cnt < 8)
1368                 nextbyte(h, 0);
1369         /* fast lookup */
1370         code = (h->sr>>(h->cnt-8))&0xFF;
1371         v = t->value[code];
1372         if(v >= 0){
1373                 h->cnt -= t->shift[code];
1374                 return v;
1375         }
1376
1377         h->cnt -= 8;
1378         if(h->cnt == 0)
1379                 nextbyte(h, 0);
1380         cnt = h->cnt;
1381         sr = h->sr;
1382         code <<= 1;
1383         for(i = 9; i<17; i++){
1384                 cnt--;
1385                 if(sr & (1<<cnt))
1386                         code |= 1;
1387                 if(code <= maxcode[i])
1388                         break;
1389                 code <<= 1;
1390                 if(cnt == 0){
1391                         sr = nextbyte(h, 0);
1392                         cnt = 8;
1393                 }
1394         }
1395         if(i >= 17){
1396                 /* bad code */
1397                 code = 0;
1398                 i = 0;
1399         }
1400         h->cnt = cnt;
1401         return t->val[t->valptr[i]+(code-t->mincode[i])];
1402 }
1403
1404 /*
1405  * load next byte of input
1406  */
1407 static
1408 int
1409 nextbyte(Header *h, int marker)
1410 {
1411         int b, b2;
1412
1413         if(h->peek >= 0){
1414                 b = h->peek;
1415                 h->peek = -1;
1416         }else{
1417                 b = Bgetc(h->fd);
1418                 if(b == Beof)
1419                         jpgerror(h, "truncated file");
1420                 b &= 0xFF;
1421         }
1422
1423         if(b == 0xFF){
1424                 if(marker)
1425                         return b;
1426                 b2 = Bgetc(h->fd);
1427                 if(b2 != 0){
1428                         if(b2 == Beof)
1429                                 jpgerror(h, "truncated file");
1430                         b2 &= 0xFF;
1431                         if(b2 == DNL)
1432                                 jpgerror(h, "ReadJPG: DNL marker unimplemented");
1433                         /* decoder is reading into marker; satisfy it and restore state */
1434                         Bungetc(h->fd);
1435                         h->peek = b;
1436                 }
1437         }
1438         h->cnt += 8;
1439         h->sr = (h->sr<<8) | b;
1440         return b;
1441 }
1442
1443 /*
1444  * return next s bits of input, MSB first, and level shift it
1445  */
1446 static
1447 int
1448 receive(Header *h, int s)
1449 {
1450         int v, m;
1451
1452         while(h->cnt < s)
1453                 nextbyte(h, 0);
1454         h->cnt -= s;
1455         v = h->sr >> h->cnt;
1456         m = (1<<s);
1457         v &= m-1;
1458         /* level shift */
1459         if(v < (m>>1))
1460                 v += ~(m-1)+1;
1461         return v;
1462 }
1463
1464 /*
1465  * return next s bits of input, decode as EOB
1466  */
1467 static
1468 int
1469 receiveEOB(Header *h, int s)
1470 {
1471         int v, m;
1472
1473         while(h->cnt < s)
1474                 nextbyte(h, 0);
1475         h->cnt -= s;
1476         v = h->sr >> h->cnt;
1477         m = (1<<s);
1478         v &= m-1;
1479         /* level shift */
1480         v += m;
1481         return v;
1482 }
1483
1484 /*
1485  * return next bit of input
1486  */
1487 static
1488 int
1489 receivebit(Header *h)
1490 {
1491         if(h->cnt < 1)
1492                 nextbyte(h, 0);
1493         h->cnt--;
1494         return (h->sr >> h->cnt) & 1;
1495 }
1496
1497 /*
1498  *  Scaled integer implementation.
1499  *  inverse two dimensional DCT, Chen-Wang algorithm
1500  * (IEEE ASSP-32, pp. 803-816, Aug. 1984)
1501  * 32-bit integer arithmetic (8 bit coefficients)
1502  * 11 mults, 29 adds per DCT
1503  *
1504  * coefficients extended to 12 bit for IEEE1180-1990 compliance
1505  */
1506
1507 enum {
1508         W1              = 2841, /* 2048*sqrt(2)*cos(1*pi/16)*/
1509         W2              = 2676, /* 2048*sqrt(2)*cos(2*pi/16)*/
1510         W3              = 2408, /* 2048*sqrt(2)*cos(3*pi/16)*/
1511         W5              = 1609, /* 2048*sqrt(2)*cos(5*pi/16)*/
1512         W6              = 1108, /* 2048*sqrt(2)*cos(6*pi/16)*/
1513         W7              = 565,  /* 2048*sqrt(2)*cos(7*pi/16)*/
1514
1515         W1pW7   = 3406, /* W1+W7*/
1516         W1mW7   = 2276, /* W1-W7*/
1517         W3pW5   = 4017, /* W3+W5*/
1518         W3mW5   = 799,  /* W3-W5*/
1519         W2pW6   = 3784, /* W2+W6*/
1520         W2mW6   = 1567, /* W2-W6*/
1521
1522         R2              = 181   /* 256/sqrt(2)*/
1523 };
1524
1525 static
1526 void
1527 idct(int b[8*8])
1528 {
1529         int x, y, eighty, v;
1530         int x0, x1, x2, x3, x4, x5, x6, x7, x8;
1531         int *p;
1532
1533         /* transform horizontally*/
1534         for(y=0; y<8; y++){
1535                 eighty = y<<3;
1536                 /* if all non-DC components are zero, just propagate the DC term*/
1537                 p = b+eighty;
1538                 if(p[1]==0)
1539                 if(p[2]==0 && p[3]==0)
1540                 if(p[4]==0 && p[5]==0)
1541                 if(p[6]==0 && p[7]==0){
1542                         v = p[0]<<3;
1543                         p[0] = v;
1544                         p[1] = v;
1545                         p[2] = v;
1546                         p[3] = v;
1547                         p[4] = v;
1548                         p[5] = v;
1549                         p[6] = v;
1550                         p[7] = v;
1551                         continue;
1552                 }
1553                 /* prescale*/
1554                 x0 = (p[0]<<11)+128;
1555                 x1 = p[4]<<11;
1556                 x2 = p[6];
1557                 x3 = p[2];
1558                 x4 = p[1];
1559                 x5 = p[7];
1560                 x6 = p[5];
1561                 x7 = p[3];
1562                 /* first stage*/
1563                 x8 = W7*(x4+x5);
1564                 x4 = x8 + W1mW7*x4;
1565                 x5 = x8 - W1pW7*x5;
1566                 x8 = W3*(x6+x7);
1567                 x6 = x8 - W3mW5*x6;
1568                 x7 = x8 - W3pW5*x7;
1569                 /* second stage*/
1570                 x8 = x0 + x1;
1571                 x0 -= x1;
1572                 x1 = W6*(x3+x2);
1573                 x2 = x1 - W2pW6*x2;
1574                 x3 = x1 + W2mW6*x3;
1575                 x1 = x4 + x6;
1576                 x4 -= x6;
1577                 x6 = x5 + x7;
1578                 x5 -= x7;
1579                 /* third stage*/
1580                 x7 = x8 + x3;
1581                 x8 -= x3;
1582                 x3 = x0 + x2;
1583                 x0 -= x2;
1584                 x2 = (R2*(x4+x5)+128)>>8;
1585                 x4 = (R2*(x4-x5)+128)>>8;
1586                 /* fourth stage*/
1587                 p[0] = (x7+x1)>>8;
1588                 p[1] = (x3+x2)>>8;
1589                 p[2] = (x0+x4)>>8;
1590                 p[3] = (x8+x6)>>8;
1591                 p[4] = (x8-x6)>>8;
1592                 p[5] = (x0-x4)>>8;
1593                 p[6] = (x3-x2)>>8;
1594                 p[7] = (x7-x1)>>8;
1595         }
1596         /* transform vertically*/
1597         for(x=0; x<8; x++){
1598                 /* if all non-DC components are zero, just propagate the DC term*/
1599                 p = b+x;
1600                 if(p[8*1]==0)
1601                 if(p[8*2]==0 && p[8*3]==0)
1602                 if(p[8*4]==0 && p[8*5]==0)
1603                 if(p[8*6]==0 && p[8*7]==0){
1604                         v = (p[8*0]+32)>>6;
1605                         p[8*0] = v;
1606                         p[8*1] = v;
1607                         p[8*2] = v;
1608                         p[8*3] = v;
1609                         p[8*4] = v;
1610                         p[8*5] = v;
1611                         p[8*6] = v;
1612                         p[8*7] = v;
1613                         continue;
1614                 }
1615                 /* prescale*/
1616                 x0 = (p[8*0]<<8)+8192;
1617                 x1 = p[8*4]<<8;
1618                 x2 = p[8*6];
1619                 x3 = p[8*2];
1620                 x4 = p[8*1];
1621                 x5 = p[8*7];
1622                 x6 = p[8*5];
1623                 x7 = p[8*3];
1624                 /* first stage*/
1625                 x8 = W7*(x4+x5) + 4;
1626                 x4 = (x8+W1mW7*x4)>>3;
1627                 x5 = (x8-W1pW7*x5)>>3;
1628                 x8 = W3*(x6+x7) + 4;
1629                 x6 = (x8-W3mW5*x6)>>3;
1630                 x7 = (x8-W3pW5*x7)>>3;
1631                 /* second stage*/
1632                 x8 = x0 + x1;
1633                 x0 -= x1;
1634                 x1 = W6*(x3+x2) + 4;
1635                 x2 = (x1-W2pW6*x2)>>3;
1636                 x3 = (x1+W2mW6*x3)>>3;
1637                 x1 = x4 + x6;
1638                 x4 -= x6;
1639                 x6 = x5 + x7;
1640                 x5 -= x7;
1641                 /* third stage*/
1642                 x7 = x8 + x3;
1643                 x8 -= x3;
1644                 x3 = x0 + x2;
1645                 x0 -= x2;
1646                 x2 = (R2*(x4+x5)+128)>>8;
1647                 x4 = (R2*(x4-x5)+128)>>8;
1648                 /* fourth stage*/
1649                 p[8*0] = (x7+x1)>>14;
1650                 p[8*1] = (x3+x2)>>14;
1651                 p[8*2] = (x0+x4)>>14;
1652                 p[8*3] = (x8+x6)>>14;
1653                 p[8*4] = (x8-x6)>>14;
1654                 p[8*5] = (x0-x4)>>14;
1655                 p[8*6] = (x3-x2)>>14;
1656                 p[8*7] = (x7-x1)>>14;
1657         }
1658 }