Naive 2D Delaunay Triangulation

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 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).

LPA Index     Home Page