summaryrefslogtreecommitdiff
path: root/appl/lib/readjpg.b
diff options
context:
space:
mode:
Diffstat (limited to 'appl/lib/readjpg.b')
-rw-r--r--appl/lib/readjpg.b973
1 files changed, 973 insertions, 0 deletions
diff --git a/appl/lib/readjpg.b b/appl/lib/readjpg.b
new file mode 100644
index 00000000..d47c15d7
--- /dev/null
+++ b/appl/lib/readjpg.b
@@ -0,0 +1,973 @@
+implement RImagefile;
+
+include "sys.m";
+ sys: Sys;
+
+include "draw.m";
+
+include "bufio.m";
+ bufio: Bufio;
+ Iobuf: import bufio;
+
+include "imagefile.m";
+
+# Constants, all preceded by byte 16rFF
+SOF: con byte 16rC0; # Start of Frame
+SOF2: con byte 16rC2; # Start of Frame; progressive Huffman
+JPG: con byte 16rC8; # Reserved for JPEG extensions
+DHT: con byte 16rC4; # Define Huffman Tables
+DAC: con byte 16rCC; # Arithmetic coding conditioning
+RST: con byte 16rD0; # Restart interval termination
+RST7: con byte 16rD7; # Restart interval termination (highest value)
+SOI: con byte 16rD8; # Start of Image
+EOI: con byte 16rD9; # End of Image
+SOS: con byte 16rDA; # Start of Scan
+DQT: con byte 16rDB; # Define quantization tables
+DNL: con byte 16rDC; # Define number of lines
+DRI: con byte 16rDD; # Define restart interval
+DHP: con byte 16rDE; # Define hierarchical progression
+EXP: con byte 16rDF; # Expand reference components
+APPn: con byte 16rE0; # Reserved for application segments
+JPGn: con byte 16rF0; # Reserved for JPEG extensions
+COM: con byte 16rFE; # Comment
+
+Header: adt
+{
+ fd: ref Iobuf;
+ ch: chan of (ref Rawimage, string);
+ # variables in i/o routines
+ sr: int; # shift register, right aligned
+ cnt: int; # # bits in right part of sr
+ buf: array of byte;
+ bufi: int;
+ nbuf: int;
+
+ Nf: int;
+ comp: array of Framecomp;
+ mode: byte;
+ X,Y: int;
+ qt: array of array of int; # quantization tables
+ dcht: array of ref Huffman;
+ acht: array of ref Huffman;
+ sf: array of byte; # start of frame; do better later
+ ss: array of byte; # start of scan; do better later
+ ri: int;
+};
+
+NBUF: con 16*1024;
+
+Huffman: adt
+{
+ bits: array of int;
+ size: array of int;
+ code: array of int;
+ val: array of int;
+ mincode: array of int;
+ maxcode: array of int;
+ valptr: array of int;
+ # fast lookup
+ value: array of int;
+ shift: array of int;
+};
+
+Framecomp: adt # Frame component specifier from SOF marker
+{
+ C: int;
+ H: int;
+ V: int;
+ Tq: int;
+};
+
+zerobytes: array of byte;
+zeroints: array of int;
+zeroreals: array of real;
+clamp: array of byte;
+NCLAMP: con 1000;
+CLAMPOFF: con 300;
+
+init(iomod: Bufio)
+{
+ if(sys == nil)
+ sys = load Sys Sys->PATH;
+ bufio = iomod;
+ zerobytes = array[8*8] of byte;
+ zeroints = array[8*8] of int;
+ zeroreals = array[8*8] of real;
+ for(k:=0; k<8*8; k++){
+ zerobytes[k] = byte 0;
+ zeroints[k] = 0;
+ zeroreals[k] = 0.0;
+ }
+ clamp = array[NCLAMP] of byte;
+ for(k=0; k<CLAMPOFF; k++)
+ clamp[k] = byte 0;
+ for(; k<CLAMPOFF+256; k++)
+ clamp[k] = byte(k-CLAMPOFF);
+ for(; k<NCLAMP; k++)
+ clamp[k] = byte 255;
+}
+
+read(fd: ref Iobuf): (ref Rawimage, string)
+{
+ # spawn a subprocess so I/O errors can clean up easily
+
+ ch := chan of (ref Rawimage, string);
+ spawn readslave(fd, ch);
+
+ return <-ch;
+}
+
+readmulti(fd: ref Iobuf): (array of ref Rawimage, string)
+{
+ (i, err) := read(fd);
+ if(i != nil){
+ a := array[1] of { i };
+ return (a, err);
+ }
+ return (nil, err);
+}
+
+readslave(fd: ref Iobuf, ch: chan of (ref Rawimage, string))
+{
+ image: ref Rawimage;
+
+ (header, err) := soiheader(fd, ch);
+ if(header == nil){
+ ch <-= (nil, err);
+ exit;
+ }
+ buf := header.buf;
+ nseg := 0;
+
+ Loop:
+ while(err == ""){
+ m: int;
+ b: array of byte;
+ nseg++;
+ (m, b, err) = readsegment(header);
+ case m{
+ -1 =>
+ break Loop;
+
+ int APPn+0 =>
+ if(nseg==1 && string b[0:4]=="JFIF"){ # JFIF header; check version
+ vers0 := int b[5];
+ vers1 := int b[6];
+ if(vers0>1 || vers1>2)
+ err = sys->sprint("ReadJPG: can't handle JFIF version %d.%2d", vers0, vers1);
+ }
+
+ int APPn+1 to int APPn+15 =>
+ ;
+
+ int DQT =>
+ err = quanttables(header, b);
+
+ int SOF =>
+ header.Y = int2(b, 1);
+ header.X = int2(b, 3);
+ header.Nf = int b[5];
+ header.comp = array[header.Nf] of Framecomp;
+ for(i:=0; i<header.Nf; i++){
+ header.comp[i].C = int b[6+3*i+0];
+ (H, V) := nibbles(b[6+3*i+1]);
+ header.comp[i].H = H;
+ header.comp[i].V = V;
+ header.comp[i].Tq = int b[6+3*i+2];
+ }
+ header.mode = SOF;
+ header.sf = b;
+
+ int SOF2 =>
+ err = sys->sprint("ReadJPG: can't handle progressive Huffman mode");
+ break Loop;
+
+ int SOS =>
+ header.ss = b;
+ (image, err) = decodescan(header);
+ if(err != "")
+ break Loop;
+
+ # BUG: THIS SHOULD USE THE LOOP TO FINISH UP
+ x := nextbyte(header, 1);
+ if(x != 16rFF)
+ err = sys->sprint("ReadJPG: didn't see marker at end of scan; saw %x", x);
+ else{
+ x = nextbyte(header, 1);
+ if(x != int EOI)
+ err = sys->sprint("ReadJPG: expected EOI saw %x", x);
+ }
+ break Loop;
+
+ int DHT =>
+ err = huffmantables(header, b);
+
+ int DRI =>
+ header.ri = int2(b, 0);
+
+ int COM =>
+ ;
+
+ int EOI =>
+ break Loop;
+
+ * =>
+ err = sys->sprint("ReadJPG: unknown marker %.2x", m);
+ }
+ }
+ ch <-= (image, err);
+}
+
+readerror(): string
+{
+ return sys->sprint("ReadJPG: read error: %r");
+}
+
+marker(buf: array of byte, n: int): byte
+{
+ if(buf[n] != byte 16rFF)
+ return byte 0;
+ return buf[n+1];
+}
+
+int2(buf: array of byte, n: int): int
+{
+ return (int buf[n]<<8)+(int buf[n+1]);
+}
+
+nibbles(b: byte): (int, int)
+{
+ i := int b;
+ return (i>>4, i&15);
+}
+
+soiheader(fd: ref Iobuf, ch: chan of (ref Rawimage, string)): (ref Header, string)
+{
+ # 1+ for restart preamble (see nextbyte), +1 for sentinel
+ buf := array[1+NBUF+1] of byte;
+ if(fd.read(buf, 2) != 2)
+ return (nil, sys->sprint("ReadJPG: can't read header: %r"));
+ if(marker(buf, 0) != SOI)
+ return (nil, sys->sprint("ReadJPG: unrecognized marker in header"));
+ h := ref Header;
+ h.buf = buf;
+ h.bufi = 0;
+ h.nbuf = 0;
+ h.fd = fd;
+ h.ri = 0;
+ h.ch = ch;
+ return (h, nil);
+}
+
+readsegment(h: ref Header): (int, array of byte, string)
+{
+ if(h.fd.read(h.buf, 2) != 2)
+ return (-1, nil, readerror());
+ m := int marker(h.buf, 0);
+ case m{
+ int EOI =>
+ return (m, nil, nil);
+ 0 =>
+ err := sys->sprint("ReadJPG: expecting marker; saw %.2x%.2x)",
+ int h.buf[0], int h.buf[1]);
+ return (-1, nil, err);
+ }
+ if(h.fd.read(h.buf, 2) != 2)
+ return (-1, nil, readerror());
+ n := int2(h.buf, 0);
+ if(n < 2)
+ return (-1, nil, readerror());
+ n -= 2;
+# if(n > len h.buf){
+# h.buf = array[n+1] of byte; # +1 for sentinel
+# #h.nbuf = n;
+# }
+ b := array[n] of byte;
+ if(h.fd.read(b, n) != n)
+ return (-1, nil, readerror());
+ return (m, b, nil);
+}
+
+huffmantables(h: ref Header, b: array of byte): string
+{
+ if(h.dcht == nil){
+ h.dcht = array[4] of ref Huffman;
+ h.acht = array[4] of ref Huffman;
+ }
+ err: string;
+ mt: int;
+ for(l:=0; l<len b; l+=17+mt){
+ (mt, err) = huffmantable(h, b[l:]);
+ if(err != nil)
+ return err;
+ }
+ return nil;
+}
+
+huffmantable(h: ref Header, b: array of byte): (int, string)
+{
+ t := ref Huffman;
+ (Tc, th) := nibbles(b[0]);
+ if(Tc > 1)
+ return (0, sys->sprint("ReadJPG: unknown Huffman table class %d", Tc));
+ if(th>3 || (h.mode==SOF && th>1))
+ return (0, sys->sprint("ReadJPG: unknown Huffman table index %d", th));
+ if(Tc == 0)
+ h.dcht[th] = t;
+ else
+ h.acht[th] = t;
+
+ # flow chart C-2
+ nsize := 0;
+ for(i:=0; i<16; i++)
+ nsize += int b[1+i];
+ t.size = array[nsize+1] of int;
+ k := 0;
+ for(i=1; i<=16; i++){
+ n := int b[i];
+ for(j:=0; j<n; j++)
+ t.size[k++] = i;
+ }
+ t.size[k] = 0;
+
+ # initialize HUFFVAL
+ t.val = array[nsize] of int;
+ for(i=0; i<nsize; i++){
+ t.val[i] = int b[17+i];
+ }
+
+ # flow chart C-3
+ t.code = array[nsize+1] of int;
+ k = 0;
+ code := 0;
+ si := t.size[0];
+ for(;;){
+ do
+ t.code[k++] = code++;
+ while(t.size[k] == si);
+ if(t.size[k] == 0)
+ break;
+ do{
+ code <<= 1;
+ si++;
+ }while(t.size[k] != si);
+ }
+
+ # flow chart F-25
+ t.mincode = array[17] of int;
+ t.maxcode = array[17] of int;
+ t.valptr = array[17] of int;
+ i = 0;
+ j := 0;
+ F25:
+ for(;;){
+ for(;;){
+ i++;
+ if(i > 16)
+ break F25;
+ if(int b[i] != 0)
+ break;
+ t.maxcode[i] = -1;
+ }
+ t.valptr[i] = j;
+ t.mincode[i] = t.code[j];
+ j += int b[i]-1;
+ t.maxcode[i] = t.code[j];
+ j++;
+ }
+
+ # create byte-indexed fast path tables
+ t.value = array[256] of int;
+ t.shift = array[256] of int;
+ maxcode := t.maxcode;
+ # stupid startup algorithm: just run machine for each byte value
+ Bytes:
+ for(v:=0; v<256; v++){
+ cnt := 7;
+ m := 1<<7;
+ code = 0;
+ sr := v;
+ i = 1;
+ for(;;i++){
+ if(sr & m)
+ code |= 1;
+ if(code <= maxcode[i])
+ break;
+ code <<= 1;
+ m >>= 1;
+ if(m == 0){
+ t.shift[v] = 0;
+ t.value[v] = -1;
+ continue Bytes;
+ }
+ cnt--;
+ }
+ t.shift[v] = 8-cnt;
+ t.value[v] = t.val[t.valptr[i]+(code-t.mincode[i])];
+ }
+
+ return (nsize, nil);
+}
+
+quanttables(h: ref Header, b: array of byte): string
+{
+ if(h.qt == nil)
+ h.qt = array[4] of array of int;
+ err: string;
+ n: int;
+ for(l:=0; l<len b; l+=1+n){
+ (n, err) = quanttable(h, b[l:]);
+ if(err != nil)
+ return err;
+ }
+ return nil;
+}
+
+quanttable(h: ref Header, b: array of byte): (int, string)
+{
+ (pq, tq) := nibbles(b[0]);
+ if(pq > 1)
+ return (0, sys->sprint("ReadJPG: unknown quantization table class %d", pq));
+ if(tq > 3)
+ return (0, sys->sprint("ReadJPG: unknown quantization table index %d", tq));
+ q := array[64] of int;
+ h.qt[tq] = q;
+ for(i:=0; i<64; i++){
+ if(pq == 0)
+ q[i] = int b[1+i];
+ else
+ q[i] = int2(b, 1+2*i);
+ }
+ return (64*(1+pq), nil);
+}
+
+zig := array[64] of {
+ 0, 1, 8, 16, 9, 2, 3, 10, 17, # 0-7
+ 24, 32, 25, 18, 11, 4, 5, # 8-15
+ 12, 19, 26, 33, 40, 48, 41, 34, # 16-23
+ 27, 20, 13, 6, 7, 14, 21, 28, # 24-31
+ 35, 42, 49, 56, 57, 50, 43, 36, # 32-39
+ 29, 22, 15, 23, 30, 37, 44, 51, # 40-47
+ 58, 59, 52, 45, 38, 31, 39, 46, # 48-55
+ 53, 60, 61, 54, 47, 55, 62, 63 # 56-63
+};
+
+decodescan(h: ref Header): (ref Rawimage, string)
+{
+ ss := h.ss;
+ Ns := int ss[0];
+ if((Ns!=3 && Ns!=1) || Ns!=h.Nf)
+ return (nil, "ReadJPG: can't handle scan not 3 components");
+
+ image := ref Rawimage;
+ image.r = ((0, 0), (h.X, h.Y));
+ image.cmap = nil;
+ image.transp = 0;
+ image.trindex = byte 0;
+ image.fields = 0;
+ image.chans = array[h.Nf] of array of byte;
+ if(Ns == 3)
+ image.chandesc = CRGB;
+ else
+ image.chandesc = CY;
+ image.nchans = h.Nf;
+ for(k:=0; k<h.Nf; k++)
+ image.chans[k] = array[h.X*h.Y] of byte;
+
+ # build per-component arrays
+ Td := array[Ns] of int;
+ Ta := array[Ns] of int;
+ data := array[Ns] of array of array of real;
+ H := array[Ns] of int;
+ V := array[Ns] of int;
+ DC := array[Ns] of int;
+
+ # compute maximum H and V
+ Hmax := 0;
+ Vmax := 0;
+ for(comp:=0; comp<Ns; comp++){
+ if(h.comp[comp].H > Hmax)
+ Hmax = h.comp[comp].H;
+ if(h.comp[comp].V > Vmax)
+ Vmax = h.comp[comp].V;
+ }
+
+ # initialize data structures
+ allHV1 := 1;
+ for(comp=0; comp<Ns; comp++){
+ # JPEG requires scan components to be in same order as in frame,
+ # so if both have 3 we know scan is Y Cb Cr and there's no need to
+ # reorder
+ cs := int ss[1+2*comp];
+ (Td[comp], Ta[comp]) = nibbles(ss[2+2*comp]);
+ H[comp] = h.comp[comp].H;
+ V[comp] = h.comp[comp].V;
+ nblock := H[comp]*V[comp];
+ if(nblock != 1)
+ allHV1 = 0;
+ data[comp] = array[nblock] of array of real;
+ DC[comp] = 0;
+ for(m:=0; m<nblock; m++)
+ data[comp][m] = array[8*8] of real;
+ }
+
+ ri := h.ri;
+
+ h.buf[0] = byte 16rFF; # see nextbyte()
+ h.cnt = 0;
+ h.sr = 0;
+ nacross := ((h.X+(8*Hmax-1))/(8*Hmax));
+ nmcu := ((h.Y+(8*Vmax-1))/(8*Vmax))*nacross;
+ zz := array[64] of real;
+ err := "";
+ for(mcu:=0; mcu<nmcu; ){
+ for(comp=0; comp<Ns; comp++){
+ dcht := h.dcht[Td[comp]];
+ acht := h.acht[Ta[comp]];
+ qt := h.qt[h.comp[comp].Tq];
+
+ for(block:=0; block<H[comp]*V[comp]; block++){
+ # F-22
+ t := decode(h, dcht);
+ diff := receive(h, t);
+ DC[comp] += diff;
+
+ # F-23
+ zz[0:] = zeroreals;
+ zz[0] = real (qt[0]*DC[comp]);
+ k = 1;
+ for(;;){
+ rs := decode(h, acht);
+ (rrrr, ssss) := nibbles(byte rs);
+ if(ssss == 0){
+ if(rrrr != 15)
+ break;
+ k += 16;
+ }else{
+ k += rrrr;
+ z := receive(h, ssss);
+ zz[zig[k]] = real (z*qt[k]);
+ if(k == 63)
+ break;
+ k++;
+ }
+ }
+
+ idct(zz, data[comp][block]);
+ }
+ }
+
+ # rotate colors to RGB and assign to bytes
+ if(Ns == 1) # very easy
+ colormap1(h, image, data[0][0], mcu, nacross);
+ else if(allHV1) # fairly easy
+ colormapall1(h, image, data[0][0], data[1][0], data[2][0], mcu, nacross);
+ else # miserable general case
+ colormap(h, image, data[0], data[1], data[2], mcu, nacross, Hmax, Vmax, H, V);
+
+ # process restart marker, if present
+ mcu++;
+ if(ri>0 && mcu<nmcu-1 && mcu%ri==0){
+ restart := mcu/ri-1;
+ rst, nskip: int;
+ nskip = 0;
+ do{
+ do{
+ rst = nextbyte(h, 1);
+ nskip++;
+ }while(rst>=0 && rst!=16rFF);
+ if(rst == 16rFF){
+ rst = nextbyte(h, 1);
+ nskip++;
+ }
+ }while(rst>=0 && (rst&~7)!=int RST);
+ if(nskip != 2)
+ err = sys->sprint("skipped %d bytes at restart %d\n", nskip-2, restart);
+ if(rst < 0)
+ return (nil, readerror());
+ if((rst&7) != (restart&7))
+ return (nil, sys->sprint("ReadJPG: expected RST%d got %d", restart&7, int rst&7));
+ h.cnt = 0;
+ h.sr = 0;
+ for(comp=0; comp<Ns; comp++)
+ DC[comp] = 0;
+ }
+ }
+ return (image, err);
+}
+
+colormap1(h: ref Header, image: ref Rawimage, data: array of real, mcu, nacross: int)
+{
+ pic := image.chans[0];
+ minx := 8*(mcu%nacross);
+ dx := 8;
+ if(minx+dx > h.X)
+ dx = h.X-minx;
+ miny := 8*(mcu/nacross);
+ dy := 8;
+ if(miny+dy > h.Y)
+ dy = h.Y-miny;
+ pici := miny*h.X+minx;
+ k := 0;
+ for(y:=0; y<dy; y++){
+ for(x:=0; x<dx; x++){
+ r := clamp[int (data[k+x]+128.)+CLAMPOFF];
+ pic[pici+x] = r;
+ }
+ pici += h.X;
+ k += 8;
+ }
+}
+
+colormapall1(h: ref Header, image: ref Rawimage, data0, data1, data2: array of real, mcu, nacross: int)
+{
+ rpic := image.chans[0];
+ gpic := image.chans[1];
+ bpic := image.chans[2];
+ minx := 8*(mcu%nacross);
+ dx := 8;
+ if(minx+dx > h.X)
+ dx = h.X-minx;
+ miny := 8*(mcu/nacross);
+ dy := 8;
+ if(miny+dy > h.Y)
+ dy = h.Y-miny;
+ pici := miny*h.X+minx;
+ k := 0;
+ for(y:=0; y<dy; y++){
+ for(x:=0; x<dx; x++){
+ Y := data0[k+x]+128.;
+ Cb := data1[k+x];
+ Cr := data2[k+x];
+ r := int (Y+1.402*Cr);
+ g := int (Y-0.34414*Cb-0.71414*Cr);
+ b := int (Y+1.772*Cb);
+ rpic[pici+x] = clamp[r+CLAMPOFF];
+ gpic[pici+x] = clamp[g+CLAMPOFF];
+ bpic[pici+x] = clamp[b+CLAMPOFF];
+ }
+ pici += h.X;
+ k += 8;
+ }
+}
+
+colormap(h: ref Header, image: ref Rawimage, data0, data1, data2: array of array of real, mcu, nacross, Hmax, Vmax: int, H, V: array of int)
+{
+ rpic := image.chans[0];
+ gpic := image.chans[1];
+ bpic := image.chans[2];
+ minx := 8*Hmax*(mcu%nacross);
+ dx := 8*Hmax;
+ if(minx+dx > h.X)
+ dx = h.X-minx;
+ miny := 8*Vmax*(mcu/nacross);
+ dy := 8*Vmax;
+ if(miny+dy > h.Y)
+ dy = h.Y-miny;
+ pici := miny*h.X+minx;
+ H0 := H[0];
+ H1 := H[1];
+ H2 := H[2];
+ for(y:=0; y<dy; y++){
+ t := y*V[0];
+ b0 := H0*(t/(8*Vmax));
+ y0 := 8*((t/Vmax)&7);
+ t = y*V[1];
+ b1 := H1*(t/(8*Vmax));
+ y1 := 8*((t/Vmax)&7);
+ t = y*V[2];
+ b2 := H2*(t/(8*Vmax));
+ y2 := 8*((t/Vmax)&7);
+ x0 := 0;
+ x1 := 0;
+ x2 := 0;
+ for(x:=0; x<dx; x++){
+ Y := data0[b0][y0+x0++*H0/Hmax]+128.;
+ Cb := data1[b1][y1+x1++*H1/Hmax];
+ Cr := data2[b2][y2+x2++*H2/Hmax];
+ if(x0*H0/Hmax >= 8){
+ x0 = 0;
+ b0++;
+ }
+ if(x1*H1/Hmax >= 8){
+ x1 = 0;
+ b1++;
+ }
+ if(x2*H2/Hmax >= 8){
+ x2 = 0;
+ b2++;
+ }
+ r := int (Y+1.402*Cr);
+ g := int (Y-0.34414*Cb-0.71414*Cr);
+ b := int (Y+1.772*Cb);
+ rpic[pici+x] = clamp[r+CLAMPOFF];
+ gpic[pici+x] = clamp[g+CLAMPOFF];
+ bpic[pici+x] = clamp[b+CLAMPOFF];
+ }
+ pici += h.X;
+ }
+}
+
+# decode next 8-bit value from entropy-coded input. chart F-26
+decode(h: ref Header, t: ref Huffman): int
+{
+ maxcode := t.maxcode;
+ if(h.cnt < 8)
+ nextbyte(h, 0);
+ # fast lookup
+ code := (h.sr>>(h.cnt-8))&16rFF;
+ v := t.value[code];
+ if(v >= 0){
+ h.cnt -= t.shift[code];
+ return v;
+ }
+
+ h.cnt -= 8;
+ if(h.cnt == 0)
+ nextbyte(h, 0);
+ h.cnt--;
+ cnt := h.cnt;
+ m := 1<<cnt;
+ sr := h.sr;
+ code <<= 1;
+ i := 9;
+ for(;;i++){
+ if(sr & m)
+ code |= 1;
+ if(code <= maxcode[i])
+ break;
+ code <<= 1;
+ m >>= 1;
+ if(m == 0){
+ sr = nextbyte(h, 0);
+ m = 16r80;
+ cnt = 8;
+ }
+ cnt--;
+ }
+ h.cnt = cnt;
+ return t.val[t.valptr[i]+(code-t.mincode[i])];
+}
+
+#
+# load next byte of input
+# we should really just call h.fd.getb(), but it's faster just to use Bufio
+# to load big chunks and manage our own byte-at-a-time input.
+#
+nextbyte(h: ref Header, marker: int): int
+{
+ b := int h.buf[h.bufi++];
+ if(b == 16rFF){
+ # check for sentinel at end of buffer
+ if(h.bufi >= h.nbuf){
+ underflow := (h.bufi > h.nbuf);
+ h.nbuf = h.fd.read(h.buf, NBUF);
+ if(h.nbuf <= 0){
+ h.ch <-= (nil, readerror());
+ exit;
+ }
+ h.buf[h.nbuf] = byte 16rFF;
+ h.bufi = 0;
+ if(underflow) # if ran off end of buffer, just restart
+ return nextbyte(h, marker);
+ }
+ if(marker)
+ return b;
+ b2 := h.buf[h.bufi++];
+ if(b2 != byte 0){
+ if(b2 == DNL){
+ h.ch <-= (nil, "ReadJPG: DNL marker unimplemented");
+ exit;
+ }else if(b2<RST && RST7<b2){
+ h.ch <-= (nil, sys->sprint("ReadJPG: unrecognized marker %x", int b2));
+ exit;
+ }
+ # decode is reading into restart marker; satisfy it and restore state
+ if(h.bufi < 2){
+ # misery: must shift up buffer
+ h.buf[1:] = h.buf[0:h.nbuf+1];
+ h.nbuf++;
+ h.buf[0] = byte 16rFF;
+ h.bufi -= 1;
+ }else
+ h.bufi -= 2;
+ b = 16rFF;
+ }
+ }
+ h.cnt += 8;
+ h.sr = (h.sr<<8)|b;
+ return b;
+}
+
+# return next s bits of input, MSB first, and level shift it
+receive(h: ref Header, s: int): int
+{
+ while(h.cnt < s)
+ nextbyte(h, 0);
+ v := h.sr >> (h.cnt-s);
+ m := (1<<s);
+ v &= m-1;
+ h.cnt -= s;
+ # level shift
+ if(v < (m>>1))
+ v += ~(m-1)+1;
+ return v;
+}
+
+# IDCT based on Arai, Agui, and Nakajima, using flow chart Figure 4.8
+# of Pennebaker & Mitchell, JPEG: Still Image Data Compression Standard.
+# Remember IDCT is reverse of flow of DCT.
+
+a0: con 1.414;
+a1: con 0.707;
+a2: con 0.541;
+a3: con 0.707;
+a4: con 1.307;
+a5: con -0.383;
+
+# scaling factors from eqn 4-35 of P&M
+s1: con 1.0196;
+s2: con 1.0823;
+s3: con 1.2026;
+s4: con 1.4142;
+s5: con 1.8000;
+s6: con 2.6131;
+s7: con 5.1258;
+
+# overall normalization of 1/16, folded into premultiplication on vertical pass
+scale: con 0.0625;
+
+idct(zin: array of real, zout: array of real)
+{
+ x, y: int;
+
+ r := array[8*8] of real;
+
+ # transform horizontally
+ for(y=0; y<8; y++){
+ eighty := y<<3;
+ # if all non-DC components are zero, just propagate the DC term
+ if(zin[eighty+1]==0.)
+ if(zin[eighty+2]==0. && zin[eighty+3]==0.)
+ if(zin[eighty+4]==0. && zin[eighty+5]==0.)
+ if(zin[eighty+6]==0. && zin[eighty+7]==0.){
+ v := zin[eighty]*a0;
+ r[eighty+0] = v;
+ r[eighty+1] = v;
+ r[eighty+2] = v;
+ r[eighty+3] = v;
+ r[eighty+4] = v;
+ r[eighty+5] = v;
+ r[eighty+6] = v;
+ r[eighty+7] = v;
+ continue;
+ }
+
+ # step 5
+ in1 := s1*zin[eighty+1];
+ in3 := s3*zin[eighty+3];
+ in5 := s5*zin[eighty+5];
+ in7 := s7*zin[eighty+7];
+ f2 := s2*zin[eighty+2];
+ f3 := s6*zin[eighty+6];
+ f5 := (in1+in7);
+ f7 := (in5+in3);
+
+ # step 4
+ g2 := f2-f3;
+ g4 := (in5-in3);
+ g6 := (in1-in7);
+ g7 := f5+f7;
+
+ # step 3.5
+ t := (g4+g6)*a5;
+
+ # step 3
+ f0 := a0*zin[eighty+0];
+ f1 := s4*zin[eighty+4];
+ f3 += f2;
+ f2 = a1*g2;
+
+ # step 2
+ g0 := f0+f1;
+ g1 := f0-f1;
+ g3 := f2+f3;
+ g4 = t-a2*g4;
+ g5 := a3*(f5-f7);
+ g6 = a4*g6+t;
+
+ # step 1
+ f0 = g0+g3;
+ f1 = g1+f2;
+ f2 = g1-f2;
+ f3 = g0-g3;
+ f5 = g5-g4;
+ f6 := g5+g6;
+ f7 = g6+g7;
+
+ # step 6
+ r[eighty+0] = (f0+f7);
+ r[eighty+1] = (f1+f6);
+ r[eighty+2] = (f2+f5);
+ r[eighty+3] = (f3-g4);
+ r[eighty+4] = (f3+g4);
+ r[eighty+5] = (f2-f5);
+ r[eighty+6] = (f1-f6);
+ r[eighty+7] = (f0-f7);
+ }
+
+ # transform vertically
+ for(x=0; x<8; x++){
+ # step 5
+ in1 := scale*s1*r[x+8];
+ in3 := scale*s3*r[x+24];
+ in5 := scale*s5*r[x+40];
+ in7 := scale*s7*r[x+56];
+ f2 := scale*s2*r[x+16];
+ f3 := scale*s6*r[x+48];
+ f5 := (in1+in7);
+ f7 := (in5+in3);
+
+ # step 4
+ g2 := f2-f3;
+ g4 := (in5-in3);
+ g6 := (in1-in7);
+ g7 := f5+f7;
+
+ # step 3.5
+ t := (g4+g6)*a5;
+
+ # step 3
+ f0 := scale*a0*r[x];
+ f1 := scale*s4*r[x+32];
+ f3 += f2;
+ f2 = a1*g2;
+
+ # step 2
+ g0 := f0+f1;
+ g1 := f0-f1;
+ g3 := f2+f3;
+ g4 = t-a2*g4;
+ g5 := a3*(f5-f7);
+ g6 = a4*g6+t;
+
+ # step 1
+ f0 = g0+g3;
+ f1 = g1+f2;
+ f2 = g1-f2;
+ f3 = g0-g3;
+ f5 = g5-g4;
+ f6 := g5+g6;
+ f7 = g6+g7;
+
+ # step 6
+ zout[x] = (f0+f7);
+ zout[x+8] = (f1+f6);
+ zout[x+16] = (f2+f5);
+ zout[x+24] = (f3-g4);
+ zout[x+32] = (f3+g4);
+ zout[x+40] = (f2-f5);
+ zout[x+48] = (f1-f6);
+ zout[x+56] = (f0-f7);
+ }
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-12-10 17:25:45 +0000