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

                          (===)
                         @-O^O-@
                          #/_\#
                           ###
  
          COMPUTATIONAL GEOMETRY IMPROVED
          ===============================
}
INTERFACE

Uses Algemeen;

procedure Delaunay(rij : punten; var mesh : drietal; var MP : punten); { Delaunay Triangulation brute force }
procedure convex_hull(rij : punten; var omheen : drietal); { Convex Hull brute force }
procedure Voronoi(rij : punten; BIG : double; var p,q : punten); { Make Voronoi Regions out of Delaunay Triangulation }
procedure Duaal(mesh : drietal; var buren : drietal); { Make dual (buren) out of triangle (mesh) }

IMPLEMENTATION

procedure Delaunay(rij : punten; var mesh : drietal; var MP : punten);
{
  Delaunay Triangulation brute force
  ----------------------------------
  out of sites (rij) as input,
  giving connectivity (mesh) and
  midpoints of circumscribes circles (MP)
}
var
  N,i,j,k,L,tel : integer;
  x,y,RR,det : double;
  OK : boolean;
  p,q,m : punt;

function minus(a,b : punt) : punt;
var
  h : punt;
begin
  h.x := a.x-b.x;
  h.y := a.y-b.y;
  minus := h;
end;

begin
  N := Length(rij);
  SetLength(mesh,2*N);
  SetLength(MP,2*N);
  tel := 0;
  for i := 0 to N-1 do
  begin
  { Three random points }
    for j := i+1 to N-1 do
    begin
      for k := j+1 to N-1 do
      begin
      { Circumscribed circles }
        omcirkel(rij[i],rij[j],rij[k],m);
        x := rij[i].x; y := rij[i].y;
        RR := sqr(x-m.x)+sqr(y-m.y);
        OK := true;
      { No points inside circles }
        for L := 0 to N-1 do
        begin
          if (L=i) or (L=j) or (L=k) then Continue;
          x := rij[L].x ; y := rij[L].y;
          if sqr(x-m.x)+sqr(y-m.y) < RR then OK := false;
        end;
        if not OK then Continue;
      { Accompanying triangles }
        mesh[tel].A := i; mesh[tel].B := j; mesh[tel].C := k;
        MP[tel] := m;
        tel := tel + 1;
      end;
    end;
  end;
  SetLength(mesh,tel);
  SetLength(MP,tel);

{ Same orientation for all triangles }
  for k := 0 to tel-1 do
  begin
    p := minus(rij[mesh[k].B],rij[mesh[k].A]);
    q := minus(rij[mesh[k].C],rij[mesh[k].A]);
    det := p.x*q.y-p.y*q.x;
    if det < 0 then
    begin
      L := mesh[k].A;
      mesh[k].A := mesh[k].B;
      mesh[k].B := L;
    end;
  end;
end;

procedure convex_hull(rij : punten; var omheen : drietal);
{
  Convex Hull brute force
  -----------------------
  out of sites (rij) as input,
  giving connectivity (omheen)
}
var
  N,i,j,k,tel : integer;
  p,q,r : punt;
  L : double;
  min,plus : boolean;
begin
  N := Length(rij);
  tel := 0;
  SetLength(omheen,N);
  for i := 0 to N-1 do
  begin
    p := rij[i];
    for j := i+1 to N-1 do
    begin
      q := rij[j];
      plus := true; min := true;
      for k := 0 to N-1 do
      begin
        if (k = i) or (k = j) then Continue;
        r := rij[k];
        L := (q.y-p.y)*(r.x-p.x)-(q.x-p.x)*(r.y-p.y);
        plus := plus and ((L > 0) or (L = 0));
        min := min and ((L < 0) or (L = 0));
      end;
      if not (plus or min) then Continue;
    { Ensure same orientation for all line segments }
      if plus then
      begin
        omheen[tel].A := i; omheen[tel].B := j;
      end;
      if min then
      begin
        omheen[tel].B := i; omheen[tel].A := j;
      end;
      tel := tel + 1;
    end;
  end;
  SetLength(omheen,tel);
end;

function normaal(BIG : double; p,q : punt) : punt;
{
  Normal at Line Segment (p,q)
}
var
  x1,x2,y1,y2,xm,ym,L : double;
  n : punt;
begin
  x1 := p.x; y1 := p.y;
  x2 := q.x; y2 := q.y;
  xm := (x1+x2)/2; ym := (y1+y2)/2;
{ Form1.Image1.Canvas.MoveTo(x2i(xm),y2j(ym)); }
  L := sqrt(sqr(x2-x1)+sqr(y2-y1));
  n.x := xm - (y2-y1)/L*BIG;
  n.y := ym + (x2-x1)/L*BIG;
{ Form1.Image1.Canvas.LineTo(x2i(n.x),y2j(n.y)); }
  normaal := n;
end;

function zelfde_zijde(mesh : drietal; i,j : integer) : boolean;
{
  Triangles (i) and (j) share the same edge
}
var
  geval : array of boolean;
  OK : boolean;
  k : integer;
begin
  SetLength(geval,9);

  geval[0] := (mesh[i].A = mesh[j].B) and (mesh[i].B = mesh[j].A);
  geval[1] := (mesh[i].A = mesh[j].A) and (mesh[i].B = mesh[j].C);
  geval[2] := (mesh[i].A = mesh[j].C) and (mesh[i].B = mesh[j].B);

  geval[3] := (mesh[i].B = mesh[j].C) and (mesh[i].C = mesh[j].B);
  geval[4] := (mesh[i].B = mesh[j].B) and (mesh[i].C = mesh[j].A);
  geval[5] := (mesh[i].B = mesh[j].A) and (mesh[i].C = mesh[j].C);

  geval[6] := (mesh[i].C = mesh[j].A) and (mesh[i].A = mesh[j].C);
  geval[7] := (mesh[i].C = mesh[j].C) and (mesh[i].A = mesh[j].B);
  geval[8] := (mesh[i].C = mesh[j].B) and (mesh[i].A = mesh[j].A);

  OK := false;
  for k := 0 to 8 do
  begin
    if not geval[k] then Continue;
    OK := true; Break;
  end;
  zelfde_zijde := OK;
end;

function aan_rand(mesh,nrs : drietal; i,j : integer) : boolean;
{
  Triangle is at Boundary
}
var
  geval : array of boolean;
  OK : boolean;
  k : integer;
begin
  SetLength(geval,3);
  geval[0] := (nrs[i].A = mesh[j].B) and (nrs[i].B = mesh[j].A);
  geval[1] := (nrs[i].A = mesh[j].C) and (nrs[i].B = mesh[j].B);
  geval[2] := (nrs[i].A = mesh[j].A) and (nrs[i].B = mesh[j].C);
  OK := false;
  for k := 0 to 2 do
  begin
    if not geval[k] then Continue;
    OK := true; Break;
  end;
  aan_rand := OK;
end;

procedure Voronoi(rij : punten; BIG : double; var p,q : punten);
{
  Make Voronoi Regions out
  of Delaunay Triangulation
}
var
  mesh,nrs : drietal;
  MP : punten;
  N,F,R,i,j,k,tel : integer;
  OK : boolean;
begin
  N := Length(rij);
  convex_hull(rij,nrs);
  R := Length(nrs);
  Delaunay(rij,mesh,MP);
  F := Length(mesh);
{ Triangles at Hull }
  for i := 0 to R-1 do
  begin
    nrs[i].C := -1;
    for j := 0 to F-1 do
    begin
      if aan_rand(mesh,nrs,i,j) then nrs[i].C := j;
    end;
  end;
{ Voronoi Regions partially }
  Setlength(p,F+N-1);
  Setlength(q,F+N-1);
  tel := 0;
  for i := 0 to F-1 do
  begin
    for j := i+1 to F-1 do
    begin
      OK := zelfde_zijde(mesh,i,j);
      if not OK then Continue;
      p[tel] := MP[i]; q[tel] := MP[j];
      tel := tel + 1;
    end;
  end;
{ Outward Normals at Convex Hull }
  for k := 0 to R-1 do
  begin
    p[tel] := MP[nrs[k].C];
    q[tel] := normaal(BIG,rij[nrs[k].A],rij[nrs[k].B]);
    tel := tel + 1;
  end;
end;

procedure Duaal(mesh : drietal; var buren : drietal);
{
  Make dual (buren) out of triangle (mesh)
}
var
  F,i,j,k : integer;
  OK : boolean;
  geval : array[0..3] of boolean;
begin
  F := Length(mesh);
  SetLength(buren,F);
  for k := 0 to F-1 do
  begin
    buren[k].A := -1;
    buren[k].B := -1;
    buren[k].C := -1;
  end;
  for i := 0 to F-1 do
  begin
    for j := 0 to F-1 do
    begin
      if i = j then Continue;
      OK := zelfde_zijde(mesh,i,j);
      if not OK then Continue;

      geval[0] := not (mesh[i].A = mesh[j].A);
      geval[1] := not (mesh[i].A = mesh[j].B);
      geval[2] := not (mesh[i].A = mesh[j].C);
      if geval[0] and geval[1] and geval[2]
      then buren[i].A := j;

      geval[0] := not (mesh[i].B = mesh[j].A);
      geval[1] := not (mesh[i].B = mesh[j].B);
      geval[2] := not (mesh[i].B = mesh[j].C);
      if geval[0] and geval[1] and geval[2]
      then buren[i].B := j;

      geval[0] := not (mesh[i].C = mesh[j].A);
      geval[1] := not (mesh[i].C = mesh[j].B);
      geval[2] := not (mesh[i].C = mesh[j].C);
      if geval[0] and geval[1] and geval[2]
      then buren[i].C := j;
    end;
  end;
end;

END.