ref: d3d1b3b47ff80f451c6c9f8b7f1262fef6545f2b
dir: /appl/lib/readjpg.b/
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);
}
}