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

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

   MORE EFFICIENT DELAUNAY TRIANGULATION
   =====================================
   Proceeding inwards from (convex) hull
}
INTERFACE

Uses Algemeen, Franco;

procedure Delaunay(row : punten; var mesh : drietal);

IMPLEMENTATION

{ EDGE BASED INSERT AND LOOKUP
  ============================ }
type
  bijeen = record
    nr : integer;
    x,y : double;
  end;
  edge = record
    A,B : integer; { Endpoints }
    g : Longint; { for sorting }
  end;
  edges = array of edge;
{
  o---->----o   Orientation of the edges
  A         B   is essential for further
  o----<----o   processing in Zijdelings
  B         A                 ----------
}
var
  zijde : edges; { Main structure }
  boven : integer;
  nr,grens : integers;
  gedaan : array of boolean;
  vrij,nieuw : array of bijeen;
  klein : integer; { Optimization }
  rond : edges; { Convex Hull }
  rij : punten; { Random Sites }
  bekijk : array of integer; { Inspect }

function BinZoek(T : edge; var n : integer) : boolean;
{
  Binary Search (Wikipedia)
  -------------------------
  Edge present (BinZoek) in
  sorted array (zijde) at (n)
}
var
  E,L,R,m : integer;
  OK : boolean;
begin
  L := 0;
  R := boven-1;
  OK := false;
  E := Length(zijde);
  T.g := E*T.A+T.B;
  while L <= R do
  begin
    m := ((L + R) div 2);
    if zijde[m].g < T.g then
      L := m + 1
    else if zijde[m].g > T.g then
      R := m - 1
    else begin
      OK := true; n := m; Break;
    end;
  end;
  BinZoek := OK;
end;

function invoegen(i,j : integer) : boolean;
{
  Sorted edges (i,j) Insert
}
var
  E,k,m,L,nop: integer;
  g,h : Longint;
  z : edge;
begin
  E := Length(zijde);
  z.A := i; z.B := j; z.g := i*E+j;
{ First entry }
  if boven = 0 then
  begin
    zijde[0] := z; boven := 1;
    invoegen := true; Exit;
  end;
  g := E*i + j;
  h := E*zijde[boven-1].A + zijde[boven-1].B;
{ On top is larger }
  if g > h then
  begin
    boven := boven + 1;
    zijde[boven-1] := z;
    invoegen := true; Exit;
  end;
{ Equal entry not done }
  if BinZoek(z,nop) then
  begin
    invoegen := false; Exit;
  end;
{ On top is smaller }
  m := boven; L := 0;
  for k := 1 to m do
  begin
    h := E*zijde[m-k].A + zijde[m-k].B;
    if g < h then
      zijde[m-k+1] := zijde[m-k]
    else Break;
    L := k;
  end;
{ Finishing touches }
  zijde[m-L+1] := zijde[m-L];
  zijde[m-L] := z;
  boven := boven + 1;
  invoegen := true;
end;

procedure test;
{
  Just a Test
}
const
  E : integer = 27;
var
  k,i,j,m : integer;
  ff : edge;
  OK : boolean;
begin
  SetLength(zijde,E);
  for k := 0 to E-1 do
  begin
    zijde[k].A := 0; zijde[k].B := 0;
    zijde[k].g := 0;
  end;
  boven := 0;
  for k := 0 to E-1 do
  begin
    i := Round(10*Random); j := Round(10*Random);
    if i = j then Continue;
    Writeln(i:3,j:3);
    invoegen(i,j);
    for m := 0 to boven-1 do
    begin
      Writeln(zijde[m].A:3,zijde[m].B:3,zijde[m].g:6);
    end;
    Writeln;
  { Readln; }
  end;
  Writeln(boven,' < ',E);
  Writeln;
  ff.A := 7;
  for k := 0 to 9 do
  begin
    ff.B := k;
    Write(ff.A:3,ff.B:3);
    m := -1;
    OK := BinZoek(ff,m);
    Writeln(OK:6,m:3);
  end;
end;

function gevonden(i,j,k : integer) : boolean;
{
  Circumscribed circle of (i,j,k)
  with no other points (L) inside
}
var
  L,n : integer;
  a,b,c,m : punt;
  RR,x,y : double;
  OK : boolean;
begin
  a := rij[i]; b := rij[j]; c := rij[k];
  if determinant(a,b,c) > 0 then OK := true else OK := false;
{ Thus trying to avoid this inner loop }
  if OK then
  begin
    omcirkel(a,b,c,m);
    RR := sqr(a.x-m.x)+sqr(a.y-m.y);
    for n := 0 to klein-1 do
    begin
      L := vrij[n].nr;
      if (L = i) or (L = j) or (L = k) then Continue;
      x := vrij[n].x; y := vrij[n].y;
      if sqr(x-m.x) + sqr(y-m.y) < RR then OK := false;
      if not OK then Break;
    end;
  end;
  gevonden := OK;
end;

procedure minder;
{
  Define search boundaries (grens)
  and which sites are done (gedaan)
}
var
  k,N,i,j,t,nop : integer;
  OK : boolean;
  W : edge;
begin
  N := Length(rij);
  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;
  for k := 0 to N-1 do
  begin
    if gedaan[k] then Continue;
    OK := (nr[k] > 0);
    if OK then
    for t := grens[k] to grens[k+1]-1 do
    begin
      i := zijde[t].A; j := zijde[t].B;
      W.A := j; W.B := i;
      OK := OK and BinZoek(W,nop);
    end;
    gedaan[k] := OK;
  end;
end;

procedure krimpen;
{
  Reduce the number (klein) 
  of relevant sites (vrij)
}
var
  k,L,tel : integer;
begin
  tel := 0;
  for k := 0 to klein-1 do
  begin
    L := vrij[k].nr;
    if gedaan[L] then Continue;
    nieuw[tel] := vrij[k];
    tel := tel + 1;
  end;
  klein := tel;
  for k := 0 to klein-1 do
  begin
    vrij[k] := nieuw[k];
  end;
end;

function toestaan(i,j,k : integer) : boolean;
{
  Are there triangles
  inside any triangle?
}
var
  L,m : integer;
  OK : boolean;
  a,b,c : punt;
begin 
  OK := true;
  a := rij[i]; b := rij[j]; c := rij[k];
  for L := grens[i] to grens[i+1]-1 do
  begin
    m := zijde[L].B;
    bekijk[m] := bekijk[m] + 1;
  end;
  for L := grens[j] to grens[j+1]-1 do
  begin
    m := zijde[L].B;
    bekijk[m] := bekijk[m] + 1;
  end;
  for L := grens[k] to grens[k+1]-1 do
  begin
    m := zijde[L].B;
    bekijk[m] := bekijk[m] + 1;
    if bekijk[m] = 3 then
    begin
      OK := OK and (not binnen(a,b,c,rij[m]));
    end;
  end;
  for L := grens[i] to grens[i+1]-1 do
  begin
    m := zijde[L].B;
    bekijk[m] := 0;
  end;
  for L := grens[j] to grens[j+1]-1 do
  begin
    m := zijde[L].B;
    bekijk[m] := 0;
  end;
  for L := grens[k] to grens[k+1]-1 do
  begin
    m := zijde[L].B;
    bekijk[m] := 0;
  end;
  toestaan := OK;
end;

procedure vanouds(var mesh : drietal);
{
  Produce triangle connectivity (mesh)
  out of sorted edges structure (zijde)
}
var
  E,N,R,F,p,q : integer;
  i,j,k,L,m,tel : integer;
  V : edge;
  a,b,c : punt;
  OK : boolean;
{ Nothing to do anymore: }
  gereed : array of boolean;
begin
  N := Length(rij);
  SetLength(bekijk,N);
  for k := 0 to N-1 do
  begin
    bekijk[k] := 0;
  end;
  E := Length(zijde);
  SetLength(gereed,E);
  for L := 0 to E-1 do
  begin
    gereed[L] := false;
  end;
{ R := 3*N - (E div 2) - 3; }
  R := Length(rond);
  F := 2*N-R-2;
  SetLength(mesh,F);

{ Finish edges at Hull }
  for k := 0 to R-1 do
  begin
    BinZoek(rond[k],p);
    gereed[p] := true;
  end;

  tel := 0;
  for L := 0 to E-1 do
  begin
    if gereed[L] then Continue;
    i := zijde[L].A; j := zijde[L].B;
    gereed[L] := true;
    for m := grens[j] to grens[j+1]-1 do
    begin
      if gereed[m] then Continue;
      k := zijde[m].B;
      V.A := k; V.B := i;
      OK := BinZoek(V,p);
    { Triangle (i,j,k) not found }
      if not OK then Continue;
      a := rij[i]; b := rij[j]; c := rij[k];
    { All triangles positive orientation }
      OK := (determinant(a,b,c) > 0);
      if not OK then Continue;
    { No triangles allowed inside triangle }
      if not toestaan(i,j,k) then Continue;
    { Search succesfully completed }
      gereed[p] := true;
      V.A := j; V.B := k;
      BinZoek(V,q);
      gereed[q] := true;
      mesh[tel].A := i; mesh[tel].B := j; mesh[tel].C := k;
      tel := tel + 1;
      Break;
    end;
  end;
end;

procedure Delaunay(row : punten; var mesh : drietal);
{
  Delaunay Triangulation
  starting with Convex Hull
}
const
  d : integer = 2;
var
  N,R,E,nop : integer;
  k,m,i,j,t,tel : integer;
  OK : boolean;
  rand : drietal; { Hull }
  ruimte : edges;
  V,W : edge;
  nrs : integers; { Hull }
begin
{ Backup needed }
  N := Length(row);
  SetLength(rij,N);
  for k := 0 to N-1 do
  begin
    rij[k] := row[k];
  end;

{ Convex Hull processing }
  Franco.convex_hull(rij,nrs);
  R := Length(nrs);
  SetLength(rand,R);
  SetLength(rond,R);
  rand[0].B := nrs[R-1];
  rand[0].A := nrs[0];
  for k := 1 to R-1 do
  begin
    rand[k].B := nrs[k-1];
    rand[k].A := nrs[k];
  end;
  for k := 0 to R-1 do
  begin
    rond[k].A := rand[k].A;
    rond[k].B := rand[k].B;
  end;

{ Edges initialization }
  N := Length(rij);
  SetLength(nr,N);
  SetLength(grens,N+1);
  SetLength(gedaan,N);
  for k := 0 to N-1 do
  begin
    nr[k] := 0;
    gedaan[k] := false;
  end;
  R := Length(rand);
  E := 2*(3*N-R-3);
  SetLength(zijde,E);
  boven := 0;
  for k := 0 to E-1 do
  begin
    zijde[k].A := 0; { start point }
    zijde[k].B := 0; { stop point }
    zijde[k].g := 0; { for sorting }
  end;
  for k := 0 to R-1 do
  begin
    i := rand[k].A; j := rand[k].B;
    OK := invoegen(i,j); { insert }
    if OK then nr[i] := nr[i] + 1;
  end;
{ Needed for efficiency }
  klein := N;
  SetLength(vrij,klein);
  for k := 0 to N-1 do
  begin
    vrij[k].x := rij[k].x;
    vrij[k].y := rij[k].y;
    vrij[k].nr := k;
  end;
  SetLength(nieuw,klein);

  SetLength(ruimte,E);
  SetLength(rand,E);
  tel := R;
{ Main iterations loop }
  while (boven < E) do
  begin
    t := 0;
    for k := 0 to tel-1 do
    begin
      j := rand[k].A; i := rand[k].B;
      OK := invoegen(i,j); { insert }
      if OK then nr[i] := nr[i] + 1;
      for m := 0 to N-1 do
      begin
        if (m = i) or (m = j) then Continue;
        if gedaan[m] then Continue;
      { Delaunay triangle found }
        OK := gevonden(i,j,m);
        if not OK then Continue;
        ruimte[t].A := j; ruimte[t].B := m;
        t := t + 1;
        ruimte[t].A := m; ruimte[t].B := i;
        t := t + 1;
      end;
    end;
    for k := 0 to t-1 do
    begin
      i := ruimte[k].A; j := ruimte[k].B;
      OK := invoegen(i,j); { insert }
      if OK then nr[i] := nr[i] + 1;
    end;
    tel := 0;
    for k := 0 to boven-1 do
    begin
      V := zijde[k];
      W.A := V.B; W.B := V.A;
    { Single edge: not at inside boundary }
      if BinZoek(W,nop) then Continue;
      i := V.A; j := V.B;
      rand[tel].A := i; rand[tel].B := j;
      tel := tel + 1;
    end;
    minder; krimpen;
  end;
{ Finishing touch }
  vanouds(mesh);
end;

END.