matrices.pas

ezw的pascal源码

{
MATRICES.PAS

Unit for matrix and vector manipulation.
A matrix is treated as a set of vectors.

Don't even try to understand the pointer operations.

(C) C. Valens, <[e-mail]>

Created     : 10/04/1998
Last update : 29/10/1998
}


unit matrices;

{******************************************************}
{*                     INTERFACE                      *}
{******************************************************}
interface

type
{ Define the data type to be stored in the matrix and vectors.}
  element_type = integer;

{ Vector type definitions. }
  element_ptr = ^element_type;
  vector = record
             v: element_ptr;
             l: integer;
           end;

{ Matrix type definitions. }
  vector_ptr = ^vector;
  matrix = record
             m: vector_ptr;
             r: integer;
             c: integer;
           end;

  statistics = record
                 N: real;
                 mean: real;
                 variance: real;
               end;

const
  min_element_type = -32768;
  max_element_type =  32767;

{ Vector operations. }
procedure create_vector(var v: vector; l: integer; var error: boolean);
procedure destroy_vector(var v: vector);
function get_vector_element(v: vector; i: integer): element_type;
procedure put_vector_element(var v: vector; i: integer; d: element_type);
function vector_max(v: vector): element_type;
function abs_vector_max(v: vector): element_type;
function vector_min(v: vector): element_type;
procedure copy_vector(src: vector; var dest: vector);
procedure write_vector(v: vector);
procedure clear_vector(var v: vector);
procedure absolute_clip_vector(var v: vector; clip_level: element_type);
procedure clip_vector(var v: vector; clip_level: element_type);
procedure absolute_threshold_vector(var v: vector; threshold: element_type);
procedure threshold_vector(var v: vector; threshold: element_type);
function vector_mean(v: vector): real;
procedure vector_statistics(v: vector; var stats: statistics);

{ Matrix operations. }
procedure create_matrix(var m: matrix; r, c: integer; var error: boolean);
procedure destroy_matrix(var m: matrix);
function get_matrix_element(m: matrix; r, c: integer): element_type;
procedure put_matrix_element(var m: matrix; r, c: integer; d: element_type);
procedure get_row(m: matrix; r: integer; var v: vector);
procedure put_row(var m: matrix; r: integer; v: vector);
procedure get_column(m: matrix; c: integer; var v: vector);
procedure put_column(var m: matrix; c: integer; v: vector);
function matrix_max(m: matrix): element_type;
function abs_matrix_max(m: matrix): element_type;
function matrix_min(m: matrix): element_type;
procedure copy_matrix(src: matrix; var dest: matrix);
procedure write_matrix(m: matrix);
procedure normalize_matrix(var m: matrix; norm: real);
procedure unnormalize_matrix(var m: matrix; min, max: element_type);
procedure reorder_matrix(src: matrix; var dest: matrix; max_level: integer);
procedure deorder_matrix(src: matrix; var dest: matrix; max_level: integer);
procedure test_matrix(m: matrix);
procedure clear_matrix(var m: matrix);
procedure absolute_clip_matrix(var m: matrix; clip_level: element_type);
procedure clip_matrix(var m: matrix; clip_level: element_type);
procedure absolute_threshold_matrix(var m: matrix; threshold: element_type);
procedure threshold_matrix(var m: matrix; threshold: element_type);
function matrix_psnr(m, ref: matrix): real;
function matrix_mean(m: matrix): real;
procedure matrix_statistics(m: matrix; var stats: statistics);
procedure matrix_shift(m: matrix; shift: element_type);


{******************************************************}
{*                  IMPLEMENTATION                    *}
{******************************************************}
implementation

{******************************************************}
{*                Vector manipulations                *}
{******************************************************}
procedure create_vector(var v: vector; l: integer; var error: boolean);
begin
  error := FALSE;
(*  GetMem(v.v,l*SizeOf(element_ptr));*)
  GetMem(v.v,l*SizeOf(element_type));
  if v.v=NIL then begin
    error := TRUE;
    v.l := 0;
  end
  else v.l := l;
end;

procedure destroy_vector(var v: vector);
begin
(*  if v.v<>NIL then FreeMem(v.v,v.l*SizeOf(element_ptr));*)
  if v.v<>NIL then FreeMem(v.v,v.l*SizeOf(element_type));
  v.l := 0;
end;

function get_vector_element(v: vector; i: integer): element_type;
begin
  get_vector_element := element_ptr(longint(v.v)+i*SizeOf(element_type))^;
end;

procedure put_vector_element(var v: vector; i: integer; d: element_type);
begin
  element_ptr(longint(v.v)+i*SizeOf(element_type))^ := d;
end;

function vector_max(v: vector): element_type;
var
  i: integer;
  max, temp: element_type;
begin
  max := get_vector_element(v,0);
  for i:=0 to v.l-1 do begin
    temp := get_vector_element(v,i);
    if max<temp then max := temp;
  end;
  vector_max := max;
end;

function abs_vector_max(v: vector): element_type;
var
  i: integer;
  max, temp: element_type;
begin
  max := 0;
  for i:=0 to v.l-1 do begin
    temp := Abs(get_vector_element(v,i));
    if max<temp then max := temp;
  end;
  abs_vector_max := max;
end;

function vector_min(v: vector): element_type;
var
  i: integer;
  min, temp: element_type;
begin
  min := get_vector_element(v,0);
  for i:=0 to v.l-1 do begin
    temp := get_vector_element(v,i);
    if min>temp then min := temp;
  end;
  vector_min := min;
end;

procedure copy_vector(src: vector; var dest: vector);
var
  i: integer;
begin
  for i:=0 to src.l-1 do begin
    put_vector_element(dest,i,get_vector_element(src,i));
  end;
end;

procedure write_vector(v: vector);
var
  i: integer;
begin
  for i:=0 to v.l-1 do begin
    write(get_vector_element(v,i):3);
  end;
end;

procedure clear_vector(var v: vector);
var
  i: integer;
begin
  for i:=0 to v.l-1 do begin
    put_vector_element(v,i,0);
  end;
end;

procedure absolute_clip_vector(var v: vector; clip_level: element_type);
var
  i: integer;
begin
  for i:=0 to v.l-1 do begin
    if Abs(get_vector_element(v,i))>clip_level then begin
      put_vector_element(v,i,clip_level);
    end
  end;
end;

procedure clip_vector(var v: vector; clip_level: element_type);
var
  i: integer;
begin
  for i:=0 to v.l-1 do begin
    if get_vector_element(v,i)>clip_level then begin
      put_vector_element(v,i,clip_level);
    end
  end;
end;

procedure absolute_threshold_vector(var v: vector; threshold: element_type);
var
  i: integer;
begin
  for i:=0 to v.l-1 do begin
    if Abs(get_vector_element(v,i))<=threshold then begin
      put_vector_element(v,i,0);
    end
  end;
end;

procedure threshold_vector(var v: vector; threshold: element_type);
var
  i: integer;
begin
  for i:=0 to v.l-1 do begin
    if get_vector_element(v,i)<=threshold then begin
      put_vector_element(v,i,0);
    end
  end;
end;

function vector_mean(v: vector): real;
var
  i: integer;
  mean, temp: real;
begin
  mean := 0;
  for i:=0 to v.l-1 do begin
    mean := mean + get_vector_element(v,i);
  end;
  temp := i;
  if temp<>0 then vector_mean := mean/temp
  else vector_mean := -1;
end;

procedure vector_statistics(v: vector; var stats: statistics);
var
  i: integer;
  data, mean, mean2: real;
begin
  mean := 0;
  mean2 := 0;
  stats.N := v.l;
  if stats.N<>0 then begin
    for i:=0 to v.l-1 do begin
      data := get_vector_element(v,i);
      mean := mean + data;
      mean2 := mean2 + Sqr(data);
    end;
    mean := mean/stats.N;
    mean2 := mean2/stats.N;
    stats.mean := mean;
    stats.variance := mean2 - Sqr(mean);
  end
  else begin
    stats.mean := 0;
    stats.variance := -1;
  end;
end;


{******************************************************}
{*                Matrix manipulations                *}
{******************************************************}
procedure create_matrix(var m: matrix; r, c: integer; var error: boolean);
var
  i: integer;
  temp: element_ptr;
begin
  error := FALSE;
  GetMem(m.m,r*SizeOf(vector));
  if m.m=NIL then error := TRUE
  else begin
    for i:=0 to r-1 do begin
      create_vector(vector_ptr(longint(m.m)+i*sizeof(vector))^,c,error);
    end;
  end;
  if error=FALSE then begin
    m.r := r;
    m.c := c;
  end
  else begin
    m.r := 0;
    m.c := 0;
  end;
end;

procedure destroy_matrix(var m: matrix);
var
  i: integer;
  temp: element_ptr;
begin
  for i:=0 to m.r-1 do begin
    destroy_vector(vector_ptr(longint(m.m)+i*sizeof(vector))^);
  end;
  if m.m<>NIL then FreeMem(m.m,m.r*SizeOf(vector));
  m.r := 0;
  m.c := 0;
end;

function get_matrix_element(m: matrix; r, c: integer): element_type;
begin
(*  get_matrix_element := get_vector_element(vector_ptr(longint(m.m)+r*sizeof(vector))^,c);*)
  get_matrix_element := element_ptr(longint(vector_ptr(longint(m.m)+r*sizeof(vector))^.v)+c*SizeOf(element_type))^;
end;

procedure put_matrix_element(var m: matrix; r, c: integer; d: element_type);
begin
(*  put_vector_element(vector_ptr(longint(m.m)+r*sizeof(vector))^,c,d);*)
  element_ptr(longint(vector_ptr(longint(m.m)+r*sizeof(vector))^.v)+c*SizeOf(element_type))^ := d;
end;

procedure get_row(m: matrix; r: integer; var v: vector);
var
  row: vector;
  i, N: integer;
begin
  row := vector_ptr(longint(m.m)+r*sizeof(vector))^;
  N := v.l;
  if v.l > m.c then N := m.c;
  for i:=0 to N-1 do begin
(*    put_vector_element(v,i,get_vector_element(row,i));*)
    put_vector_element(v,i,element_ptr(longint(row.v)+i*SizeOf(element_type))^);
  end;
end;

procedure put_row(var m: matrix; r: integer; v: vector);
var
  row: vector;
  i, N: integer;
begin
  row := vector_ptr(longint(m.m)+r*sizeof(vector))^;
  N := v.l;
  if v.l > m.c then N := m.c;
  for i:=0 to N-1 do begin
(*    put_vector_element(row,i,get_vector_element(v,i));*)
    put_vector_element(row,i,element_ptr(longint(v.v)+i*SizeOf(element_type))^);
  end;
end;

procedure get_column(m: matrix; c: integer; var v: vector);
var
  i, N: integer;
begin
  N := v.l;
  if v.l > m.r then N := m.r;
  for i:=0 to N-1 do begin
    put_vector_element(v,i,get_matrix_element(m,i,c));
  end;
end;

procedure put_column(var m: matrix; c: integer; v: vector);
var
  i, N: integer;
begin
  N := v.l;
  if v.l > m.r then N := m.r;
  for i:=0 to N-1 do begin
    put_matrix_element(m,i,c,get_vector_element(v,i));
  end;
end;

function matrix_max(m: matrix): element_type;
var
  row: vector;
  i: integer;
  max, temp: element_type;
  error: boolean;
begin
  create_vector(row,m.c,error);
  max := get_matrix_element(m,0,0);
  for i:=0 to m.r-1 do begin
    get_row(m,i,row);
    temp := vector_max(row);
    if max<temp then max := temp;