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