Unit Sterren;
{
This software has been designed and it is
CopyLefted by Han de Bruijn:

                          (===)
                         @-O^O-@
                          #/_\#
                           ###

DELAUNAY TESSELLATION FIELD ESTIMATOR
=====================================
}
INTERFACE

Uses Algemeen;

type
  dubbel = array of double;

procedure histogram(rij : punten; mesh : drietal; var DTFE : dubbel); { DTFE }
procedure Makkelijk(rij : punten; mesh : drietal; var DTFE : dubbel); { DTFE }
procedure voorbeeld(M,N : integer; var rij : punten; var mesh : drietal); { Example }

IMPLEMENTATION

procedure histogram(rij : punten; mesh : drietal; var DTFE : dubbel);
{
  Delaunay Tessellation Field Estimator
  with sites (rij), connectivity (mesh)
  giving normed output (DTFE) = 1/area
  with the area of the Voronoi Regions
}
var
  i,j,k,L,M,N,t,D : integer;
  nr,grens : integers;
  Opp : array of double;
  niet,rand : array of boolean;
  a,b,c,MP,PM : punt;
  uit : drietal;
  h : drie;
  max : double;
begin
{ Calculate size of (uit) }
  N := Length(rij);
  SetLength(nr,N);
  for k := 0 to N-1 do
  begin
    nr[k] := 0;
  end;
  M := Length(mesh);
  for L := 0 to M-1 do
  begin
    i := mesh[L].A; j := mesh[L].B; k := mesh[L].C;
    nr[i] := nr[i]+1; nr[j] := nr[j]+1; nr[k] := nr[k]+1;
  end;
  SetLength(grens,N+1);
  grens[0] := 0;
  for k := 1 to N do
  begin
    grens[k] := 0;
    for t := 0 to k-1 do
    begin
      grens[k] := grens[k] + nr[t];
    end;
  end;
  L := grens[N];
  SetLength(uit,L);

{ All triangles around a site
  ordered to site numbering }
  for k := 0 to N-1 do
  begin
    nr[k] := 0;
  end;
  for L := 0 to M-1 do
  begin
    i := mesh[L].A; j := mesh[L].B; k := mesh[L].C;
    uit[grens[i]+nr[i]].A := i; uit[grens[i]+nr[i]].B := L;
    uit[grens[j]+nr[j]].A := j; uit[grens[j]+nr[j]].B := L;
    uit[grens[k]+nr[k]].A := k; uit[grens[k]+nr[k]].B := L;
    nr[i] := nr[i]+1; nr[j] := nr[j]+1; nr[k] := nr[k]+1;
  end;
  for L := 0 to grens[N]-1 do
  begin
    t := uit[L].A; D := uit[L].B;
    i := mesh[D].A; j := mesh[D].B; k := mesh[D].C;
    if t = i then
    begin
      uit[L].B := j; uit[L].C := k;
    end;
    if t = j then
    begin
      uit[L].B := k; uit[L].C := i;
    end;
    if t = k then
    begin
      uit[L].B := i; uit[L].C := j;
    end;
  end;

{ Sort triangles in counterclockwise order }
  for k := 0 to N-1 do
  begin
    for L := grens[k]+1 to grens[k+1]-1 do
    begin
      if uit[L].B = uit[L-1].C then Continue;
      for t := L+1 to grens[k+1]-1 do
      begin
        if uit[t].B = uit[L-1].C then
        begin
          h := uit[L];
          uit[L] := uit[t];
          uit[t] := h;
          Break;
        end;
      end;
    end;
  end;

{ Sites at boundary
  shall be not done }
  SetLength(niet,N);
  SetLength(rand,N);
  for k := 0 to N-1 do
  begin
    niet[k] := false;
    rand[k] := false;
  end;
  for k := 0 to N-1 do
  begin
    for L := grens[k] to grens[k+1]-1 do
    begin
      niet[uit[L].B] := not niet[uit[L].B];
      niet[uit[L].C] := not niet[uit[L].C];
    end;
    for L := grens[k] to grens[k+1]-1 do
    begin
      if niet[uit[L].B] or niet[uit[L].C]
      then rand[k] := true;
      niet[uit[L].B] := false;
      niet[uit[L].C] := false;
    end;
  end;

{ Area of Voronoi Regions }
  SetLength(Opp,N);
  for k := 0 to N-1 do
  begin
    Opp[k] := 0;
    if rand[k] then Continue;
    L := grens[k];
    a := rij[uit[L].A];
    b := rij[uit[L].B];
    c := rij[uit[L].C];
    omcirkel(a,b,c,MP);
    for L := grens[k]+1 to grens[k+1]-1 do
    begin
      PM := MP;
      a := rij[uit[L].A];
      b := rij[uit[L].B];
      c := rij[uit[L].C];
      omcirkel(a,b,c,MP);
      Opp[k] := Opp[k] + determinant(a,PM,MP);
    end;
    PM := MP;
    L := grens[k];
    a := rij[uit[L].A];
    b := rij[uit[L].B];
    c := rij[uit[L].C];
    omcirkel(a,b,c,MP);
    Opp[k] := Opp[k] + determinant(a,PM,MP);  
  end;

{ Delaunay Tessellation
  Field Estimator normed }
  SetLength(DTFE,N);
  for k := 0 to N-1 do
  begin
    DTFE[k] := 0;
    if rand[k] then Continue;
    DTFE[k] := 1/Opp[k];
  end;
  max := 0;
  for k := 0 to N-1 do
  begin
    if DTFE[k] > max then max := DTFE[k];
  end;
  for k := 0 to N-1 do
  begin
    DTFE[k] := DTFE[k]/max;
  end;
end;

procedure Makkelijk(rij : punten; mesh : drietal; var DTFE : dubbel);
{
  Area of Delaunay Triangles
  instead of Voronoi Regions
}
var
  i,j,k,L,F,N : integer;
  a,b,c : punt;
  Opp : dubbel;
  det,max : double;
begin
  N := Length(rij);
  SetLength(Opp,N);
  for k := 0 to N-1 do
  begin
    Opp[k] := 0;
  end;
  F := Length(mesh);
  for L := 0 to F-1 do
  begin
    i := mesh[L].A; j := mesh[L].B; k := mesh[L].C;
    a := rij[i]; b := rij[j]; c := rij[k];
    det := Determinant(a,b,c);
    Opp[i] := Opp[i] + det;
    Opp[j] := Opp[j] + det;
    Opp[k] := Opp[k] + det;
  end;
{ for k := 0 to N-1 do
  begin
    Writeln(k:5,Opp[k]);
  end; }

{ Delaunay Tessellation
  Field Estimator normed }
  SetLength(DTFE,N);
  for k := 0 to N-1 do
  begin
    DTFE[k] := 1/Opp[k];
  end;
  max := 0;
  for k := 0 to N-1 do
  begin
    if DTFE[k] > max then max := DTFE[k];
  end;
  for k := 0 to N-1 do
  begin
    DTFE[k] := DTFE[k]/max;
  end;
end;

procedure voorbeeld(M,N : integer; var rij : punten; var mesh : drietal);
{
  Wheel shape
}
var
  dp,dr,x,y : double;
  i,j,k,F,t : integer;
begin
{ Coordinates }
  SetLength(rij,M*N+1);
  dp := 2*Pi/M;
  dr := 1/N;
  rij[0].x := 0; rij[0].y := 0;
  for i := 0 to N-1 do
  begin
    for j := 0 to M-1 do
    begin
      k := j*N + i+1;
      x := (i+1)*dr*cos(j*dp); y := (i+1)*dr*sin(j*dp);
      rij[k].x := x; rij[k].y := y;
    end;
  end;

{ Connectivity }
  F := 2*M*N-M;
  SetLength(mesh,F);
  for j := 1 to M-1 do
  begin
    mesh[j-1].A := 0;
    mesh[j-1].B := (j-1)*N + 1;
    mesh[j-1].C := j*N + 1;
  end;
  mesh[M-1].A := 0;
  mesh[M-1].B := (M-1)*N + 1;
  mesh[M-1].C := 1;
  t := M;
  for i := 1 to N-1 do
  begin
    for j := 1 to M-1 do
    begin
      mesh[t].A := (j-1)*N + i;
      mesh[t].B := (j-1)*N + i+1;
      mesh[t].C := j*N + i+1;
      t := t + 1;
      mesh[t].A := j*N + i+1;
      mesh[t].B := j*N + i;
      mesh[t].C := (j-1)*N + i;
      t := t + 1;
    end;
    mesh[t].A := (M-1)*N + i;
    mesh[t].B := (M-1)*N + i+1;
    mesh[t].C := i+1;
    t := t + 1;
    mesh[t].A := i+1;
    mesh[t].B := i;
    mesh[t].C := (M-1)*N + i;
    t := t + 1;
  end;
end;

END.