#### 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),
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 */
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,

/* 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([], _, _, _).
member(P, Ignore), !,
distance_squared(P, Centre, DistanceSquared),

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