In a Delaunay triangulation of a set of two-dimensional points, the circumcircle of each triangle contains no points other than the vertices of that triangle.
/* delaunay(Points, Triangles) is true if the elements t(P1,P2,P3) of */ /* Triangles comprise the Delaunay triangulation of the Points, the */ /* vertices of the triangles going anticlockwise. */ /* e.g. delaunay([p(1,5),p(9,3),p(0,9),p(5,5),p(4,2)], X) */ /* gives X=[t(p(1,5),p(5,5),p(0,9)), t(p(1,5),p(4,2),p(5,5)), */ /* t(p(9,3),p(0,9),p(5,5)), t(p(9,3),p(5,5),p(4,2))] */ /* This is a naive O(n^4) algorithm. It may fail or give incorrect results */ /* if there are repeated points, or three collinear points, or four */ /* co-circular points. */ delaunay(Points, Triangles):- findall(Triangle, delaunay_1(Points, Triangle), Triangles). delaunay_1(Points, Triangle):- rest(P1, Points, Points1), rest(P2, Points1, Points2), member(P3, Points2), circle_3(P1, P2, P3, Centre, RadiusSquared), all_outside(Points, [P1,P2,P3], Centre, RadiusSquared), orient(t(P1,P2,P3), Triangle). /* 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 * (A*(Y3 - Y2) - B*(X3 - X2)), G =\= 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. /* all_outside(Points, Ignore, Centre, RadiusSquared) is true if all the */ /* Points not in Ignore are outside the boundary of the circle with the */ /* given Centre and RadiusSquared. */ all_outside([], _, _, _). all_outside([P|Points], Ignore, Centre, RadiusSquared):- member(P, Ignore), !, all_outside(Points, Ignore, Centre, RadiusSquared). all_outside([P|Points], Ignore, Centre, RadiusSquared):- distance_squared(P, Centre, DistanceSquared), DistanceSquared > RadiusSquared, all_outside(Points, Ignore, Centre, RadiusSquared). orient(t(P1,P2,P3), t(P1,P3,P2)):- strictly_to_right(P3, P1, P2), !. orient(Triangle, Triangle). /* strictly_to_right(Pa, Pb, Pc) is true if Pa is strictly to the right of */ /* the directed line from Pb to Pc. */ strictly_to_right(p(Xa,Ya), p(Xb,Yb), p(Xc,Yc)):- (Xb-Xa) * (Yc-Ya) - (Xc-Xa) * (Yb-Ya) < 0. /* 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. /* rest(X, Ys, Zs) is true if X is a member of the list Ys, and the list Zs */ /* is the rest of the list following X. On backtracking, all solutions */ /* will be found. */ rest(X, [X|Ys], Ys). rest(X, [_|Ys], Zs):-rest(X, Ys, Zs). /* member(X, Xs) is true if the element X is contained in the list Xs. On */ /* backtracking, all solutions will be found. */ %member(X, [X|_]). %member(X, [_|Xs]):-member(X, Xs).