/appl/lib/readjpg.b
Brainfuck | 973 lines | 885 code | 88 blank | 0 comment | 67 complexity | 1a454a5e87f1da8d21170a805059b2ed MD5 | raw file
Possible License(s): GPL-2.0, LGPL-2.1, MPL-2.0-no-copyleft-exception
- 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);
- }
- }