{
EZW.PAS

Unit for Embedded Zero Tree coding.

Based on "Embedded Image Coding Using Zerotrees of Wavelet Coefficients"
by Jerome M. Shapiro, IEEE Transactions on Signal Processing, Vol.41, No.12,
December 1993, pp 3445-3462.

A fifo is used in the dominant pass which results in a so-called Morton order
scan instead of Shapiro's raster scan (see figure 2 in "Analysis Based Coding
of Image Transform and Subband Coefficients" by V. Ralph Algazi and Robert
R. Estes, Jr.).

Morton order scan:
==================

   1 | 2 |  5   6 | 17  18  21  22
  ---+---|        |
   3 | 4 |  7   8 | 19  20  23  24
  -------+--------|
   9  10 | 13  14 | 25  26  29  30
         |        |
  11  12 | 15  16 | 27  28  31  32
  ----------------+---------------
  33  34   37  38 | 49  50  53  54
                  |
  35  36   39  40 | 51  52  55  56
                  |
  41  42   45  46 | 57  58  61  62
                  |
  43  44   47  48 | 59  60  63  64


Raster scan:
============

   1 | 2 |  5   6 | 17  18  19  20
  ---+---|        |
   3 | 4 |  7   8 | 21  22  23  24
  -------+--------|
   9  10 | 13  14 | 25  26  27  28
         |        |
  11  12 | 15  16 | 29  30  31  32
  ----------------+---------------
  33  34   35  36 | 49  50  51  52
                  |
  37  38   39  40 | 53  54  55  56
                  |
  41  42   43  44 | 57  58  59  60
                  |
  45  46   47  48 | 61  62  63  64


Subband distribution:
=====================

  LL | HL | HL   HL | HL   HL   HL   HL
  ---+--- |         |
  LH | HH | HL   HL | HL   HL   HL   HL
  --------+---------|
  LH   LH | HH   HH | HL   HL   HL   HL
          |         |
  LH   LH | HH   HH | HL   HL   HL   HL
  ------------------+------------------
  LH   LH   LH   LH | HH   HH   HH   HH
                    |
  LH   LH   LH   LH | HH   HH   HH   HH
                    |
  LH   LH   LH   LH | HH   HH   HH   HH
                    |
  LH   LH   LH   LH | HH   HH   HH   HH


(C) C. Valens, <c.valens@mindless.com>

Created    : 01/05/1998
Last update: 25/05/1998
}

program ezwdos;

{$define debug}
{$define disk}

uses
{$ifdef disk}
  matrices, diskfifo, disklist;
{$else}
  matrices, fifo, list;
{$endif}

type
  ezw_file_header = record
    height, width: integer;
    threshold: element_type;
  end;

const
(* Shapiro's example data *)
  example: array[0..7, 0..7] of integer = (
    ( 63,-34, 49, 10,  7, 13,-12,  7 ),
    (-31, 23, 14,-13,  3,  4,  6, -1 ),
    ( 15, 14,  3,-12,  5, -7,  3,  9 ),
    ( -9, -7,-14,  8,  4, -2,  3,  2 ),
    ( -5,  9, -1, 47,  4,  6, -2,  2 ),
    (  3,  0, -3,  2,  3, -2,  0,  4 ),
    (  2, -3,  6, -4,  3,  6,  3,  6 ),
    (  5, 11,  5,  6,  0,  3, -4,  4 )
  );

(* Code alphabet. *)
  zero = 0; (* binary 0 *)
  one  = 1; (* binary 1 *)
  ztr  = 2; (* binary 00 *)
  pos  = 3; (* binary 01 *)
  neg  = 4; (* binary 11 *)
  iz   = 5; (* binary 10 *)

var
  M: matrix;
  error: boolean;
  zeroes, ones: longint;
  ezw_file: file;
  output_byte, mask: char;
  header: ezw_file_header;

procedure load_data(var m: matrix);
var
  i, j: integer;
begin
  for i:=0 to 7 do begin
    for j:=0 to 7 do begin
      put_matrix_element(m,i,j,element_type(example[i,j]));
    end;
  end;
end;

(*
* Puts a bit in the output stream.
*)
procedure put_bit(bit: char);
begin
  if bit='1' then begin
    output_byte := Chr(Ord(output_byte) or Ord(mask));
    inc(ones);
  end
  else inc(zeroes);

  mask := Chr(Ord(mask) shr 1);
  if Ord(mask)=0 then begin
    BlockWrite(ezw_file,output_byte,1);
    output_byte := Chr(0);
    mask := Chr($80);
  end;
end;

(*
* Puts dominant-pass and subordinate-pass codes in the output stream.
*)
procedure output_code(code: integer);
begin
  case code of
    zero: begin
            put_bit('0');
{$ifdef debug}
            Write('0');
{$endif debug}
          end;
    one : begin
            put_bit('1');
{$ifdef debug}
            Write('1');
{$endif debug}
          end;
    pos : begin
            put_bit('0');
            put_bit('1');
{$ifdef debug}
            Write('p');
{$endif debug}
          end;
    neg : begin
            put_bit('1');
            put_bit('1');
{$ifdef debug}
            Write('n');
{$endif debug}
          end;
    ztr : begin
            put_bit('0');
            put_bit('0');
{$ifdef debug}
            Write('t');
{$endif debug}
          end;
    iz  : begin
            put_bit('1');
            put_bit('0');
{$ifdef debug}
            Write('i');
{$endif debug}
          end;
  end;
end;

(*
* Returns the largest value in a descendance tree.
*)
function max_descendant(m: matrix; x, y: integer): element_type;
var
  i, j, min_x, max_x, min_y, max_y: integer;
  temp, max: element_type;
begin
  if (x=0) and (y=0) then begin
    temp := get_matrix_element(m,0,0);
    put_matrix_element(m,0,0,min_element_type);
    max := abs_matrix_max(m);
    put_matrix_element(m,0,0,temp);
  end
  else begin
    min_x := x shl 1;
    min_y := y shl 1;
    max_x := (x+1) shl 1;
    max_y := (y+1) shl 1;
    if (min_x=m.c) or (min_y=m.r) then max_descendant := 0;

    max := 0;
    while (max_y<=m.r) and (max_x<=m.c) do begin
      for i:=min_y to max_y-1 do begin
        for j:=min_x to max_x-1 do begin
          temp := Abs(get_matrix_element(m,i,j));
          if temp>max then max := temp;
        end;
      end;
      min_x := min_x shl 1;
      max_x := max_x shl 1;
      min_y := min_y shl 1;
      max_y := max_y shl 1;
    end;

  end;

  max_descendant := max;

end;

(*
* Returns TRUE if descendance tree is a zerotree.
*)
function zerotree(m: matrix; x, y, threshold: integer): boolean;
var
  i, j, min_x, max_x, min_y, max_y: integer;
  temp, max: element_type;
  stop: boolean;
begin
  stop := FALSE;
  if (x=0) and (y=0) then begin
    temp := get_matrix_element(m,0,0);
    put_matrix_element(m,0,0,min_element_type);
    max := abs_matrix_max(m);
    put_matrix_element(m,0,0,temp);
    if max>=threshold then stop := TRUE;
  end
  else begin
    min_x := x shl 1;
    min_y := y shl 1;
    max_x := (x+1) shl 1;
    max_y := (y+1) shl 1;
    if (min_x=m.c) or (min_y=m.r) then zerotree := TRUE;

    max := 0;
    while (max_y<=m.r) and (max_x<=m.c) do begin
      for i:=min_y to max_y-1 do begin
        for j:=min_x to max_x-1 do begin
          temp := Abs(get_matrix_element(m,i,j));
          if temp>=threshold then begin
            stop := TRUE;
            break;
          end;
        end;
        if stop=TRUE then break;
      end;
      if stop=TRUE then break;
      min_x := min_x shl 1;
      max_x := max_x shl 1;
      min_y := min_y shl 1;
      max_y := max_y shl 1;
    end;
  end;
  if stop=TRUE then zerotree := FALSE
  else zerotree := TRUE;
end;

(*
* Returns a dominant-pass-code from the alphabet [pos,neg,ztr,iz].
*)
function code(m: matrix; x, y: integer; threshold: element_type): integer;
var
  temp: element_type;
begin
  temp := get_matrix_element(m,y,x);
  if Abs(temp)>=threshold then begin
    if temp>=0 then code := pos
    else code := neg;
  end
  else begin
(*    if (max_descendant(m,x,y)<threshold) then code := ztr*)
    if zerotree(m,x,y,threshold)=TRUE then code := ztr
    else code := iz;
  end;
end;

(*
* Appends a value to the subordinate list.
*)
procedure to_sub_list(value: element_type);
var
  d: list_type;
begin
(* Put only coefficient magnitude in list, sign is allready coded. *)
  d.x := Abs(value);
  d.y := 0;
  append_to_list(d);
end;

(*
* Builds a dominant pass EZW-element from a matrix element and a threshold.
*)
procedure process_element(m: matrix; threshold: element_type;
                          var s: ezw_element);
begin
  s.code := code(m,s.x,s.y,threshold);
  if (s.code=pos) or (s.code=neg) then begin
    to_sub_list(get_matrix_element(m,s.y,s.x));
    put_matrix_element(m,s.y,s.x,0);
  end;
end;

(*
* Performs one complete dominant pass. Dominant-pass-codes are sent to the
* output stream and the subordinate list is updated.
*)
procedure dominant_pass(m: matrix; threshold: element_type);
var
  s: ezw_element;
  min_x, max_x, min_y, max_y: integer;
  level: integer;
begin
  s.x := 0;
  s.y := 0;
  process_element(m,threshold,s);
  output_code(s.code);

  s.x := 1;
  s.y := 0;
  process_element(m,threshold,s);
  put_in_fifo(s);
  s.x := 0;
  s.y := 1;
  process_element(m,threshold,s);
  put_in_fifo(s);
  s.x := 1;
  s.y := 1;
  process_element(m,threshold,s);
  put_in_fifo(s);

  get_from_fifo(s);
  if fifo_empty=FALSE then output_code(s.code);

  while fifo_empty=FALSE do begin
    if s.code<>ztr then begin
      min_x := s.x shl 1;
      max_x := min_x+1;
      min_y := s.y shl 1;
      max_y := min_y+1;
      if (max_x<=m.c) and (max_y<=m.r) then begin
        for s.y:=min_y to max_y do begin
          for s.x:=min_x to max_x do begin
            process_element(m,threshold,s);
            put_in_fifo(s);
          end;
        end;
      end;
    end;
    get_from_fifo(s);
    if fifo_empty=FALSE then output_code(s.code);
  end;

end;

(*
* Performs one subordinate pass.
*)
procedure subordinate_pass(threshold: element_type);
var
  d: list_type;
  i: Longint;
  found: boolean;
begin
  if threshold>0 then begin
    for i:=0 to list_length-1 do begin
      get_list_element(d,i,found);
      if found=TRUE then begin
        if (d.x and threshold<>0) then output_code(one)
        else output_code(zero);
      end;
    end;
  end;
end;

(*
* EZW-codes matrix m, returns initial threshold.
*)
procedure EZW_code(m: matrix; threshold: element_type);
begin
  while threshold<>0 do begin
    dominant_pass(m,threshold);
    subordinate_pass(threshold shr 1);
    threshold := threshold shr 1;
  end;
end;

(*
* Main.
*)
begin
  Writeln;

  header.height := 8;
  header.width := 8;
  create_matrix(M,header.height,header.width,error);
  if error=TRUE then Exit;
  load_data(M);
  header.threshold := 1 shl Trunc((Ln(abs_matrix_max(M))/Ln(2)));

  Assign(ezw_file,'out.ezw');
  Rewrite(ezw_file,1);
  BlockWrite(ezw_file,header,SizeOf(header));

  zeroes := 0;
  ones := 0;
  output_byte := Chr(0);
  mask := Chr($80);
{$ifdef debug}
  write_matrix(M);
{$endif debug}
  EZW_code(M,header.threshold);
  if Ord(mask)<>0 then BlockWrite(ezw_file,output_byte,1);
{$ifdef debug}
  Writeln;
  Writeln(zeroes+ones,' bits: ',zeroes,' zeroes, ',ones,' ones.');
{$endif debug}

  Close(ezw_file);

  destroy_matrix(M);
  destroy_fifo;
  destroy_list;
end.