MinAreaCircle - December 2024 version

program new;

var
  img: TImage;
  TPA: TPointArray;
  c: TCircle;


function TPointArray.FurthestPoint(from: TPoint): TPoint;
var i: Int32;
begin
  if Length(Self) = 0 then Exit;
  Result := Self[0];
  for i:=1 to High(Self) do
    if from.DistanceTo(Self[i]) > from.DistanceTo(Result)  then
      Result := Self[i];
end;

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;


(*
  Smallest bounding circle algorithm by Slacky
  https://en.wikipedia.org/wiki/Smallest-circle_problem

  Time complexity: O(max(n log n, h^2))

  Where `h` is the length of the convex hull, and `n` is the data size.
  Can be guaranteed O(n log n) with rotating calipers algorithm for FurthestPoints,
  dont think it's worth the added complexity though.
*)
function TPointArray.MinAreaCircle2(): TCircle;
var
  poly: TPointArray;
  p1,p2,p3,p,q,t: TPoint;

  function Circle3(const 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;

begin
  poly := Self.ConvexHull();  // O(n log n)
  poly.FurthestPoints(p1,p2); // O(m^2) call - ugh

  p := [(p1.x+p2.x) div 2, (p1.y+p2.y) div 2];
  t := poly.FurthestPoint(p);

  // Two Point Circle
  if (p1 = t) or (p2 = t) then
    Exit([p.x,p.y, Ceil( p1.DistanceTo(p2) / 2 )]);

  // Three Point Circle
  p := Circle3(p1, p2, t).Center;
  q := poly.FurthestPoint(p);

  if Sign(CrossProduct(q,  p1, p2)) = Sign(CrossProduct(t, p1, p2)) then
    p3 := q
  else
    p3 := t;

  p := Circle3(p1,p2,p3).Center;
  q := poly.FurthestPoint(p);
  if (p1 <> q) and (p2 <> q) and (p3 <> q) then
  begin
    p := [(p3.x+q.x) div 2, (p3.y+q.y) div 2];
    if p1.DistanceTo(p) < p2.DistanceTo(p) then
    begin
      p1 := q;
      Swap(p1, p2);
    end else
      p2 := q;

    p := Circle3(p1, p2, p3).Center;
    q := poly.FurthestPoint(p);
    if (p1 <> q) and (p2 <> q) and (p3 <> q) then
      p1 := q;
  end;

  Result := Circle3(p1, p2, p3);
end;

function TPointArray.ToString(): string;
var p: TPoint;
begin
  Result += '[';
  for p in self do
    Result += '['+ToString(p.x)+','+ToString(p.y)+'],';
  SetLength(Result, Length(Result)-1);
  Result += ']';
end;

var
  i: Int32;
  poly: TPointArray;
  c1,c2: TCircle;
  p: TPoint;
begin
  img := TImage.Create(900,900);
  img.Show(False);
  RandSeed := 12345;
  Writeln(RandSeed);

  i := 0;
  while True do
  begin
    TPA := RandomTPA(50, [250,250,550,550]);

    img.DrawColor := $FF;
    img.DrawPolygon(TPA.ConvexHull());

    c2 := TPA.MinAreaCircle2();
    img.DrawColor := $FFFF00;
    img.DrawCircle(c2.Center, c2.Radius);

    img.DrawColor := $FF;
    img.DrawPolygon(TPA.ConvexHull());


    img.Show(False);
    img.Clear();

    Inc(i);
  end;

  img.Show(False);
end.