Quickhull

/*
 * Convex Hull (Quickhull)
 *   The convex hull of a set of points is the smallest convex region
 *   containing the points
 */

/* quickhull(Points, ConvexHullVertices) is true if Points is a list of      */
/*   points in the form p(X,Y), and ConvexHullVertices are the vertices in   */
/*   the form p(X,Y) of the convex hull of the Points, in clockwise order,   */
/*   starting and ending at the smallest point (as determined by X-values,   */
/*   and by Y-values to resolve ties).                                       */
/* e.g. quickhull([p(0,6),p(3,7),p(4,6),p(4,5),p(3,4),p(2,4),p(5,0)],        */
/*                [p(0,6),p(3,7),p(4,6),p(5,0),p(0,6)]).                     */
quickhull(Points, [P1|Vertices]):-
  min_max_list(Points, P1, P2),
  points_to_left(Points, P1, P2, LeftMost, Ls),
  points_to_left(Points, P2, P1, RightMost, Rs),
  \+ both_empty(Ls, Rs), !,
  quickhull_1(Rs, P2, RightMost, P1, Vertices0, []),
  quickhull_1(Ls, P1, LeftMost, P2, Vertices, Vertices0).
quickhull(Points, [P1,P1]):-
  min_max_list(Points, P1, P2),
  P1=P2, /* All points identical */
  !.
quickhull(Points, [P1,P2,P1]):-
  /* All points collinear */
  min_max_list(Points, P1, P2).

/* quickhull_1(Ps, P1, P2, P3, Ys, []) is true if Ys is that part of the     */
/*   convex hull to the left of the line from P1 to P3 of those points in Ps */
/*   outside the triangle defined by the points P1, P2 and P3, where P2 is   */
/*   that element of Ps furthest to the left of the line from P1 to P3.      */
quickhull_1([], _, _, P3, [P3|Ys], Ys):-!.
quickhull_1(Ps, P1, P2, P3, Ys, Zs):-
  points_to_left(Ps, P1, P2, LeftMost12, Ps1),
  points_to_left(Ps, P2, P3, LeftMost23, Ps2),
  quickhull_1(Ps2, P2, LeftMost23, P3, Ws, Zs),
  quickhull_1(Ps1, P1, LeftMost12, P2, Ys, Ws).

both_empty([], []).

/* min_max_list(List, Min, Max) is true if Min is the minimum value in the   */
/*   List, and Max is the maximum value.                                     */
min_max_list([X|Xs], Min, Max):-min_max_list_1(Xs, X, X, Min, Max).

min_max_list_1([], Min, Max, Min, Max).
min_max_list_1([X|Xs], OldMin, OldMax, Min, Max):-
  min_max(X, OldMin, OldMax, NewMin, NewMax),
  min_max_list_1(Xs, NewMin, NewMax, Min, Max).

min_max(p(X,Y), p(X0,_), Max, p(X,Y), Max):-X < X0, !.
min_max(p(X,Y), p(X,Y0), Max, p(X,Y), Max):-Y < Y0, !. /* Resolve ties */
min_max(p(X,Y), Min, p(X0,_), Min, p(X,Y)):-X > X0, !.
min_max(p(X,Y), Min, p(X,Y0), Min, p(X,Y)):-Y > Y0, !. /* Resolve ties */
min_max(_, Min, Max, Min, Max).

/* points_to_left(Ps, P1, P2, LeftMost, Ls) is true if Ls are the points of  */
/*   Ps strictly to the left of the line through P1 and P2, and LeftMost is  */
/*   the point furthest to the left of the line through P1 and P2.           */
points_to_left(Ps, P1, P2, LeftMost, Ls):-
  points_to_left_1(Ps, P1, P2, a(0.0,void,LeftMost), Ls).

points_to_left_1([], _, _, a(_,LeftMost,LeftMost), []).
points_to_left_1([P|Ps], P1, P2, State, Ls):-
  area(P1, P2, P, Area),
  points_to_left_2(Area, Ps, P, P1, P2, State, Ls).

points_to_left_2(Area, Ps, P, P1, P2, a(Area0,_,LeftMost), [P|Ls]):-
  Area > Area0, !, /* P is to left, and further than P0 */
  points_to_left_1(Ps, P1, P2, a(Area,P,LeftMost), Ls).
points_to_left_2(Area, Ps, P, P1, P2, a(Area,P0,LeftMost), [P|Ls]):-
  /* P is to left, and same distance as P0 */
  area(P1, P0, P, Area1),
  Area1 > 0, !,     /* P1 is nearer to P than P0 */
  points_to_left_1(Ps, P1, P2, a(Area,P,LeftMost), Ls).
points_to_left_2(Area, Ps, P, P1, P2, State, [P|Ls]):-
  Area > 0.0, !, /* P is to left and nearer than P0 */
  points_to_left_1(Ps, P1, P2, State, Ls).
points_to_left_2(_, Ps, _, P1, P2, State, Ls):-
  /* P is to right or collinear */
  points_to_left_1(Ps, P1, P2, State, Ls).

/* area(Pa, Pb, Pc, Area) is true if Area is twice the signed area of the    */
/*   triangle determined by Pa, Pb and Pc. The area is positive if Pa, Pb    */
/*   and Pc are oriented anti-clockwise, negative if clockwise, and zero if  */
/*   the points are collinear.                                               */
area(p(Xa,Ya), p(Xb,Yb), p(Xc,Yc), Area):-
  Area is (Xb-Xa) * (Yc-Ya) - (Xc-Xa) * (Yb-Ya).

LPA Index     Home Page

Valid HTML 4.01 Strict