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

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

Uses Algemeen;

procedure convex_hull(rand : punten; var nrs : integers);

IMPLEMENTATION

procedure convex_hull(rand : punten; var nrs : integers);
{ 
AN EFFICIENT ALGORITHM FOR DETERMINING
THE CONVEX HULL OF A FINITE PLANAR SET

           L.R.Graham (1972)
}
type
{ Cartesian & polar coordinates
  and necessary bookkeeping }
  tmp = record
    x,y,c,r : double;
    yes : boolean;
    nr : integer;
  end;
  temp = array of tmp;
var
  rij : temp;
  N,k,kk,tel,L : integer;
  min,det : double;
  O,p,q,h : tmp;
  max : double;

function Hoek(v : tmp) : double;
{
  Angle equivalent of
  vector with x-axis
    -1 < Hoek < +1
}
var
  L,c,s : double;
begin
{ Length : }
  L := sqrt(sqr(v.x)+sqr(v.y));
{ Sine and cosine: }
  c := v.x/L ; s := v.y/L;
{ The triangle spanned by (c,s) and (1,0) is equilateral
  Hence the sum-vector (1+c,s) divides the top angle in
  two equal pieces. The accompanying tangens is s/(1+c).
  Must be doubled to obtain the original angle again. Do
  it in a safe manner & make the denominator non-zero:
  h := 2*arctan(s/(1+abs(c)));
  With positive cosine, result is between -Pi/2 & +Pi/2
  However, the equivalent shortcut below is sufficient. }
  Hoek := s/(1+abs(c));
end;

procedure DonaldKnuth(var rij : temp);
var
  L : integer;

  procedure KwikSort(Left, Right : integer);
  { 
    Quick sort.
    Described on page 114 of "The Art of Computer Programming"
    volume 3 / "Sorting and Searching" by Donald E. Knuth.
    Programmed as in: Turbo Pascal 7, The Complete Reference,
    page 705
  }
  var 
    i,j : integer;
    k : double;
  
    procedure Switch(var A, B : tmp);
    var
      C : tmp;
    begin
      if (A.c <> B.c) then
      begin
        C := A;
        A := B;
        B := C;
      end;
    end;
  
  begin
    k := (rij[left].c + rij[right].c)/2;
    i := left;
    j := right;
    repeat
      while rij[i].c < k do
        Inc(i,1);
      while k < rij[j].c do
        Dec(j,1);
      if i <= j then
      begin
        Switch(rij[i],rij[j]);
        Inc(i,1);
        Dec(j,1);
      end;
    until i > j;
    if left < j then
     KwikSort(left, j);
    if i < right then
     KwikSort(i, right);
  end;

begin
  L := Length(rij)-1;
  KwikSort(0,L);
end;

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

begin
{ Make copy of points }
  N := Length(rand);
  SetLength(rij,N);
  for k := 0 to N-1 do
  begin
    rij[k].x := rand[k].x; rij[k].y := rand[k].y;
    rij[k].yes := true; rij[k].nr := k;
  end;
  min := rij[0].x;
{ Determine Origin }
  kk := 0;
  for k := 1 to N-1 do
  begin
    if rij[k].x < min then
    begin
      min := rij[k].x;
      kk := k;
    end;
  end;
  O := rij[kk];

{ Compute polar coordinates }
  tel := 0;
  for k := 0 to N-1 do
  begin
    if k = kk then Continue;
    rij[tel].nr := rij[k].nr;
    rij[tel].x := rij[k].x-O.x;
    rij[tel].y := rij[k].y-O.y;
    rij[tel].c := Hoek(rij[tel]);
    rij[tel].r := sqrt(sqr(rij[tel].x)+sqr(rij[tel].y));
    tel := tel + 1;
  end;
  SetLength(rij,tel);
  DonaldKnuth(rij);

{ Eliminate Special Cases }
  N := Length(rij);
  tel := 0;
  for k := 0 to N-1 do
  begin
    if rij[k].r = 0 then Continue;
    rij[tel] := rij[k];
    tel := tel + 1;
  end;
  SetLength(rij,tel);
{ Especially equal angles }
  N := Length(rij);
  for k := 0 to N-1 do
  begin
    max := rij[k].r;
    kk := k;
    for L := k+1 to N-1 do
    begin
      if rij[L].c > rij[k].c then Break;
      if rij[L].r > max then
      begin
        kk := L;
        max := rij[L].r;
      end;
    end;
    rij[k] := rij[kk];
  end;
  tel := 1;
  h := rij[0];
  for k := 1 to N-1 do
  begin
    if rij[k].c = h.c then Continue;
    rij[tel] := rij[k];
    h := rij[k];
    tel := tel + 1;
  end;
  SetLength(rij,tel);

{ Determine Convex Hull }
  while true do
  begin
    N := Length(rij);
    for k := 1 to N-2 do
    begin
      p := minus(rij[k],rij[k-1]);
      q := minus(rij[k+1],rij[k]);
      det := p.x*q.y-p.y*q.x;
      if det < 0 then rij[k].yes := false;
    end;
    tel := 0;
    for k := 0 to N-1 do
    begin
      if not rij[k].yes then Continue;
      rij[tel] := rij[k];
      tel := tel + 1;
    end;
    if tel = N then Break;
    SetLength(rij,tel);
  end;

{ Finishing Touches }
  tel := tel + 1;
  SetLength(rij,tel);
  rij[tel-1] := O;
  SetLength(nrs,tel);
  for k := 0 to tel-1 do
  begin
    nrs[k] := rij[k].nr;
  end;
end;

END.