Minimal Enclosing Circle

/*
 * Minimal Enclosing Circle
 *   The minimal enclosing circle of a set of points is the circle with the
 *   smallest radius containing all the points
 */

/* mec(Points, p(Xc,Yc), RadiusSquared) is true if the circle whose centre   */
/*   is at (Xc,Yc) and whose radius squared is RadiusSquared is the minimum  */
/*   enclosing circle of the set of Points.                                  */
/*   This is Emo Welzl's O(n) time algorithm.                                */
/* e.g. mec([p(3,1),p(4,8),p(3,10),p(8,8)], C, R) gives C=p(4.1,5.5),R=21.46 */
mec([Point], Point, 0.0):-!.
mec(Points, Centre, RadiusSquared):-
  mec_1(Points, [], Centre, RadiusSquared).

/* mec_1(OtherPoints, BoundaryPoints, Centre, RadiusSquared) is true if the  */
/*   BoundaryPoints are on the circle with the given Centre and              */
/*   RadiusSquared, and the OtherPoints are on or inside this same circle.   */
mec_1([], [R1,R2], Centre, RadiusSquared):-!,
  circle_2(R1, R2, Centre, RadiusSquared).
mec_1(_, [R1,R2,R3], Centre, RadiusSquared):-
  circle_3(R1, R2, R3, Centre, RadiusSquared), !.
mec_1(_, [_,_,_], p(0,0), 0.0):-!,
  write(`Two identical points`), nl,
  abort.
mec_1([P|Ps], Rs, Centre, RadiusSquared):-
  mec_1(Ps, Rs, Centre, RadiusSquared),
  distance_squared(P, Centre, DistanceSquared),
  DistanceSquared =< RadiusSquared, !.
mec_1([P|Ps], Rs, Centre, RadiusSquared):-
  mec_1(Ps, [P|Rs], Centre, RadiusSquared).

/* circle_2(p(X1,Y1), p(X2,Y2), p(Xc,Yc), RadiusSquared) is true if the      */
/*   centre of the smallest circle passing through the two points (X1,Y1)    */
/*   and (X2,Y2) is (Xc,Yc), and RadiusSquared is the square of its radius.  */
/* e.g. circle_2(p(-6,8), p(-2,-4), p(X,Y), R) gives X=-4, Y=2, R=40         */
circle_2(p(X1,Y1), p(X2,Y2), p(Xc,Yc), RadiusSquared):-
  Xc is (X1+X2) / 2.0,
  Yc is (Y1+Y2) / 2.0,
  X1_X2 is X1 - X2,
  Y1_Y2 is Y1 - Y2,
  RadiusSquared is (X1_X2*X1_X2 + Y1_Y2*Y1_Y2) / 4.0.

/* circle_3(p(X1,Y1), p(X2,Y2), p(X3,Y3), p(Xc,Yc), RadiusSquared) is true   */
/*   if the centre of circle passing through the three points (X1,Y1),       */
/*   (X2,Y2) and (X3,Y3) is (Xc,Yc), and RadiusSquared is the square of its  */
/*   radius. Fails if the three points are collinear.                        */
/*   This algorithm is borrowed from the comp.graphics.algorithms FAQ.       */
/* e.g. circle_3(p(-6,5), p(-3,-4), p(2,1), p(X,Y), R) gives X=-3, Y=1, R=25 */
circle_3(p(X1,Y1), p(X2,Y2), p(X3,Y3), p(Xc,Yc), RadiusSquared):-
  A is X2 - X1,
  B is Y2 - Y1,
  G is 2.0 * (A*(Y3 - Y2) - B*(X3 - X2)),
  G =\= 0.0,
  C is X3 - X1,
  D is Y3 - Y1,
  E is A*(X1 + X2) + B*(Y1 + Y2),
  F is C*(X1 + X3) + D*(Y1 + Y3),
  Xc is (D*E - B*F) / G,
  Yc is (A*F - C*E) / G,
  H is X1 - Xc,
  I is Y1 - Yc,
  RadiusSquared is H*H + I*I.

/* distance_squared(p(X1,Y1), p(X2,Y2), DistanceSquared) is true if          */
/*   DistanceSquared is the square of the distance between the points        */
/*   (X1,Y1) and (X2,Y2).                                                    */
distance_squared(p(X1,Y1), p(X2,Y2), DistanceSquared):-
  X1_X2 is X1 - X2,
  Y1_Y2 is Y1 - Y2,
  DistanceSquared is X1_X2*X1_X2 + Y1_Y2*Y1_Y2.

LPA Index     Home Page