Approximated MinAreaCircle

program new;
{$I SRL/osr.simba}

var
  bmp: TMufasaBitmap;

//technically this can be solved in n log n time.. but meh.. so much work
procedure FurthestPoints(const TPA:TPointArray; out A,B:TPoint);
var
  i,j: Integer;
  dist,tmp: Single;
begin
  for i:=0 to High(TPA) do
    for j:=i+1 to High(TPA) do
    begin
      tmp := Hypot(TPA[i].x-TPA[j].x, TPA[i].y-TPA[j].y);
      if tmp > dist then
      begin
        dist := tmp;
        A := TPA[i];
        B := TPA[j];
      end;
    end;
end;

(*
  Approximate smallest bounding circle algorithm by Slacky
  This formula makes 3 assumptions of outlier points that defines the circle.

  p1 and p2 are the furthest possible points from one another
  p3 is the furthest point from the center of the line of p1 and p2 at a 90 degree angle.

  Already now we have a circle that is very close to encapsulating.

  We assume that p3 is nearly always correct, and if there is one wrong assumption it's either p1, or p2.
  We move p1 and/or p2 until the criteria is reached, odds are one of them is correct already alongside p3.
  If we actually didn't solve it so far, we move p3 to see if that does it.

  https://en.wikipedia.org/wiki/Smallest-circle_problem

  Time complexity is somewhere along the lines of
    O(max(n log n, h^2))

  Where `h` is the remainding points after convex hull is performed.
  `n` is the number of points, and `n log n` is assuming optimal convexhull algorithm.
*)
function MinAreaCircle(arr: TPointArray): TCircle;
  function Circle3(P1,P2,P3: TPoint): TCircle;
  var
    t,d,l: Single;
    p,pcs,q,qcs: record x,y: Single; end;
  begin
    t   := ArcTan2(p1.y-p2.y, p1.x-p2.x) - PI/2;
    p   := [(p1.x+p2.x) / 2, (p1.y+p2.y) / 2];
    pcs := [Cos(t), Sin(t)];

    t   := ArcTan2(p1.y-p3.y, p1.x-p3.x) - PI/2;
    q   := [(p1.x+p3.x) / 2, (p1.y+p3.y) / 2];
    qcs := [Cos(t), Sin(t)];

    d := pcs.x * qcs.y - pcs.y * qcs.x;
    if Abs(d) < 0.001 then Exit();
    l := (qcs.x * (p.y - q.y) - qcs.y * (p.x - q.x)) / d;
    Result.x := Round(p.x + l * pcs.x);
    Result.y := Round(p.y + l * pcs.y);
    Result.Radius := Ceil(Sqrt(Sqr((p.x + l * pcs.x)-p1.x) + Sqr((p.y + l * pcs.y)-p1.y)));
  end;

  function InCircle(arr: TPointArray; c: TCircle): Boolean;
  var i: Int32;
  begin
    Result := True;
    for i:=0 to High(arr) do
      if Hypot(arr[i].x-c.x, arr[i].y-c.Y) > c.Radius+1 then Exit(False);
  end;

  function pass(p1,p2,p3: TPoint; now:TCircle; var new: TPoint): TCircle;
  var i: Int32; c: TCircle;
  begin
    new := p2;
    Result := now;
    for i:=0 to High(arr) do
    begin
      if (arr[i] = p3) or (arr[i] = p1) then continue;
      c := Circle3(p1,arr[i],p3);
      if (c.Radius < Result.Radius) and InCircle(arr, c) then
      begin
        Result := c;
        new := arr[i];
      end;
    end;
  end;
var
  i: Int32;
  v,q,p1,p2,p3: TPoint;
begin
  arr := arr.ConvexHull();
  FurthestPoints(arr, p1,p2);

  v := [(p1.x + p2.x) div 2, (p1.y + p2.y) div 2];
  p3 := arr[0];
  for i:=0 to High(arr) do
    if Sqrt(Sqr(v.x-arr[i].x)+Sqr(v.y-arr[i].y)) > Sqrt(Sqr(v.x-p3.x)+Sqr(v.y-p3.y)) then
      p3 := arr[i];

  if (p3 = p2) or (p3 = p1) then
    Exit(TCircle([Round(v.x), Round(v.y), Ceil(Distance(p1,p2) / 2)]));

  Result := Circle3(p1,p2,p3);
  if Result.Radius = 0 then
    Exit(TCircle([Round(v.x), Round(v.y), Ceil(Distance(p1,p2) / 2)]));

  if InCircle(arr, Result) then Exit(Result);

  Result.Radius := $FFFF;
  Result := pass(p1,p2,p3, Result, v);
  if v <> p2 then q := v;
  Result := pass(p2,p1,p3, Result, v);
  if v <> p1 then Result := pass(v,q,p3, Result, v);

  if InCircle(arr, Result) and (Result.Radius <> $FFFF) then Exit(Result);

  Result := pass(p2,p3,p1, Result, v); // very rarely is the p3 assumption is wrong.

  if (Result.Radius = $FFFF) then Result.Radius := Ceil(Distance(p1,p2) / 2);
end;

function ToString(s: TPoint): string; override;
begin
  Result := '['+ToString(s.x)+','+ToString(s.y)+']';
end;

function InCircle(arr: TPointArray; c: TCircle): Boolean;
var i: Int32;
begin
  Result := True;
  for i:=0 to High(arr) do
    if Hypot(arr[i].x-c.x, arr[i].y-c.Y) > c.Radius+1 then Exit(False);
end;


var
  TPA: TPointArray;
  a: T2DPointArray;
  e1,e2,i: Int32;
  t1,t2,t1_sum,t2_sum: Double;
  c1,c2: TCircle;
begin
  bmp.Init();
  bmp.SetSize(900,900);
  bmp.Debug();

  while True do
  begin
    TPA := RandomTPA(Random(5,150), [250,250,550,550]);

    c1 := MinAreaCircle(TPA);
    c2 := TPA.MinAreaCircle();

    c1.Radius += 1;
    (*
    if (not InCircle(TPA.ConvexHull(),c1)) then
    begin
      WriteLn(c1);
      WriteLn(TPA.ConvexHull());
      bmp.DrawCircle(c1.x, c1.y, c1.Radius, $FF);
      bmp.DrawCircle(c2.x, c2.y, c2.Radius, $00FF00);
      bmp.DrawTPA(TPA.ConvexHull().Connect(), $FF00FF);
      bmp.Debug();
      bmp.Clear();

      WaitUntil(isKeyDown(VK_ENTER), 50, 60000);
      Wait(20)
    end;
    *)

    if (Abs(c1.Radius - c2.Radius) > 10) then
    begin
      WriteLn(c1);
      WriteLn(TPA.ConvexHull());
      bmp.DrawCircle(c1.x, c1.y, c1.Radius, $FF);
      bmp.DrawCircle(c2.x, c2.y, c2.Radius, $00FF00);
      bmp.DrawTPA(TPA.ConvexHull().Connect(), $FF00FF);
      bmp.Debug();
      bmp.Clear();

      WriteLn(TPA.ConvexHull());
      WaitUntil(isKeyDown(VK_ENTER), 50, 60000);
      Wait(20);

      if c1.Radius > c2.Radius then Inc(e1);
      if c1.Radius < c2.Radius then Inc(e2);
    end;

    Inc(i);
    if i mod 1000 = 0 then
    begin
      WriteLn('Tested: ', i, '. Error count (new/old): ', e1, '/', e2);
    end;
  end;

  //bmp.Debug();
end.