14.4. Convex Hulls¶
- File:
ConvexHulls.ml
In the last section of our brief introduction to problems, techniques, and algorithms of computational geometry, we take a look at a very simple yet practically relevant problem: constructing a convex hull of a set of points. A convex hull for a set of points \(S\) is a smalles polygon (in terms of are) such that all points from \(S\) are contained within it, or are lying on its boundary.
The definition implies that the vertices of the polygons are some points from \(S\).
The image below shows an example of a convex hull for a set of 50 points:
14.4.1. Plane-sweeping algorithm¶
Let us study the following elegant algorithm, known as Graham scan, for computing a convex hull fora set of two-dimensional points.
The algorithm relies on sorting and implements a “plane-sweeping” intuition by considering all points in a certain sequence, making sure to include only those to the hull-in-construction, that do not disrupt the convexity property.
Our construction of a convex hull will rely on imperative stacks, which we will extend with the following definitions:
include Polygons
open Stacks
module StackX (S: AbstractStack) = struct
include S
let top s = match pop s with
| None -> None
| Some x ->
push s x;
Some x
let next_to_top s = match pop s with
| None -> None
| Some x ->
let y = top s in
push s x;
y
let list_of_stack s =
let res = ref [] in
while not (is_empty s) do
let e = Util.get_exn @@ pop s in
res := e :: !res
done;
!res
end
The crux of the algorithm is to sort the set of points twice, using different comparisons:
- The first sort is done by \(y\)-coordinate;
- The second sort is done by a radial angle from the point \(p_0\) with a smallest \(y\) coordinate (and leftmost of points with the smallest \(y\)-coordinate).
For these sortings, we will need the following two procedures:
(* Sort by axis Y *)
let axis_y_sorter (Point (x1, y1)) (Point (x2, y2)) =
if y1 < y2 then -1 else if y1 > y2 then 1
else if x1 < x2 then -1 else if x1 > x1 then 1
else 0
(* Sort by polar angle wrt p0 *)
let polar_angle_sorter p0 p1 p2 =
let Polar (r1, a1) = p1 -- p0 |> polar_of_cartesian in
let Polar (r2, a2) = p2 -- p0 |> polar_of_cartesian in
if a1 < a2 then -1 else if a1 > a2 then 1
else if r1 < r2 then -1 else if r1 > r2 then 1
else 0
The main Graham algorithm is as follows:
(* Graham's Scan *)
let convex_hull points =
(* At least three points *)
assert (List.length points >= 3);
let y_sorted = List.sort axis_y_sorter points in
let p0 = y_sorted |> List.hd in
match List.tl y_sorted |>
List.sort (polar_angle_sorter p0) with
| p1 :: p2 :: rest ->
let open CHStack in
let s = mk_stack 0 in
push s p0;
push s p1;
push s p2;
let non_left_turn p =
let q1 = next_to_top s |> get_exn in
let q2 = top s |> get_exn in
direction q1 q2 p >= 0
in
(* Main algorithm *)
List.iter (fun p ->
while non_left_turn p do
ignore (pop s)
done;
push s p) rest;
list_of_stack s
| _ -> error "Cannot happen"
The main loop checks all the vertices in the order of the increasing
angle from p0
, making sure that the currently build convex hull
for a subset is convex. It removes the points that violat this
invariants from the
Question: What is the complexity of the procedure in terms of the size n
of the set of points?
14.4.2. Graham scan invariant¶
A moment from the middle of Graham scan is shown on a picture below:
It is easy to see that the invariant of the main loop is that all
points in the stack s
always correspond to a convex polygon
containing all the points observed so far. Hence, by induction, the
resulting polygon is also convex.
Since we have ways to check both these properties (convexity and containment of a point in a polygon, we can engineer the following tests for Graham scan).
First, let us generate a set of random points of a fixed size n
:
let gen_random_points ?dim:(dim = 550.) n =
let res = ref [] in
for _ = 0 to n - 1 do
let p = gen_random_point dim in
res := p :: !res
done;
!res
Second, let us use it in a randomised test:
open ConvexHulls
let test_random_ch n =
let ps = gen_random_points n in
let ch = convex_hull ps in
assert (is_convex ch);
assert (List.for_all (point_within_polygon ch) ps)
let%test _ =
for _ = 0 to 100 do
test_random_ch 50
done;
true