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