Introduction

Convex hull is a fundamental problem in computational geometry. It involves finding the smallest convex polygon that completely encloses a set of points in the plane. Convex hulls have various applications in fields like image processing, navigation, and robotics. Andrew’s algorithm is a well-known method for constructing the convex hull of a set of points in Euclidean space.

Implementation

(*
 * Convex hulls by Andrew's monotone chain algorithm.
 *
 * For a description of the algorithm, see
 * https://en.wikibooks.org/w/index.php?title=Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain&stableid=40169
 *)

(*------------------------------------------------------------------*)
(* Just enough plane geometry for our purpose.  *)

module Plane_point =
  struct
    type t = float * float

    let make (xy : t) = xy
    let to_tuple ((x, y) : t) = (x, y)

    let x ((x, _) : t) = x
    let y ((_, y) : t) = y

    let equal (u : t) (v : t) = (x u = x v && y u = y v)

    (* Impose a total order on points, making it one that will work
       for Andrew's monotone chain algorithm. *)
    let order (u : t) (v : t) =
      (x u < x v) || (x u = x v && y u < y v)

    (* The Array module's sort routines expect a "cmp" function. *)
    let cmp u v =
      if order u v then
        (-1)
      else if order v u then
        (1)
      else
        (0)

    (* Subtraction is really a vector or multivector operation. *)
    let sub (u : t) (v : t) = make (x u -. x v, y u -. y v)

    (* Cross product is really a multivector operation. *)
    let cross (u : t) (v : t) = (x u *. y v) -. (y u *. x v)

    let to_string ((x, y) : t) =
      "(" ^ string_of_float x ^ " " ^ string_of_float y ^ ")"
  end
;;

(*------------------------------------------------------------------*)
(* We want something akin to array_delete_neighbor_dups of Scheme's
   SRFI-132. *)

let array_delete_neighbor_dups equal arr =
  (* Returns a Seq.t rather than an array. *)
  let rec loop i lst =
    (* Cons a list of non-duplicates, going backwards through the
       array so the list will be in forwards order. *)
    if i = 0 then
      arr.(i) :: lst
    else if equal arr.(i - 1) arr.(i) then
      loop (i - 1) lst
    else
      loop (i - 1) (arr.(i) :: lst)
  in
  let n = Array.length arr in
  List.to_seq (if n = 0 then [] else loop (n - 1) [])
;;

(*------------------------------------------------------------------*)
(* The convex hull algorithm. *)

let cross_test pt_i hull j =
  let hull_j = hull.(j)
  and hull_j1 = hull.(j - 1) in
  0.0 < Plane_point.(cross (sub hull_j hull_j1)
                       (sub pt_i hull_j1))

let construct_lower_hull n pt =
  let hull = Array.make n (Plane_point.make (0.0, 0.0)) in
  let () = hull.(0) <- pt.(0)
  and () = hull.(1) <- pt.(1) in
  let rec outer_loop i j =
    if i = n then
      j + 1
    else
      let pt_i = pt.(i) in
      let rec inner_loop j =
        if j = 0 || cross_test pt_i hull j then
          begin
            hull.(j + 1) <- pt_i;
            j + 1
          end
        else
          inner_loop (j - 1)
      in
      outer_loop (i + 1) (inner_loop j)
  in
  let hull_size = outer_loop 2 1 in
  (hull_size, hull)

let construct_upper_hull n pt =
  let hull = Array.make n (Plane_point.make (0.0, 0.0)) in
  let () = hull.(0) <- pt.(n - 1)
  and () = hull.(1) <- pt.(n - 2) in
  let rec outer_loop i j =
    if i = (-1) then
      j + 1
    else
      let pt_i = pt.(i) in
      let rec inner_loop j =
        if j = 0 || cross_test pt_i hull j then
          begin
            hull.(j + 1) <- pt_i;
            j + 1
          end
        else
          inner_loop (j - 1)
      in
      outer_loop (i - 1) (inner_loop j)
  in
  let hull_size = outer_loop (n - 3) 1 in
  (hull_size, hull)

let construct_hull n pt =

  (* Side note: Construction of the lower and upper hulls can be done
     in parallel. *)
  let (lower_hull_size, lower_hull) = construct_lower_hull n pt
  and (upper_hull_size, upper_hull) = construct_upper_hull n pt in

  let hull_size = lower_hull_size + upper_hull_size - 2 in
  let hull = Array.make hull_size (Plane_point.make (0.0, 0.0)) in

  begin
    Array.blit lower_hull 0 hull 0 (lower_hull_size - 1);
    Array.blit upper_hull 0 hull (lower_hull_size - 1)
      (upper_hull_size - 1);
    hull
  end

let plane_convex_hull points =
  (* Takes an arbitrary sequence of points, which may be in any order
     and may contain duplicates. Returns an ordered array of points
     that make up the convex hull. If the sequence of points is empty,
     the returned array is empty. *)
  let pt = Array.of_seq points in
  let () = Array.fast_sort Plane_point.cmp pt in
  let pt = Array.of_seq
             (array_delete_neighbor_dups Plane_point.equal pt) in
  let n = Array.length pt in
  if n <= 2 then
    pt
  else
    construct_hull n pt
;;

(*------------------------------------------------------------------*)

let example_points =
  [Plane_point.make (16.0, 3.0);
   Plane_point.make (12.0, 17.0);
   Plane_point.make (0.0, 6.0);
   Plane_point.make (-4.0, -6.0);
   Plane_point.make (16.0, 6.0);
   Plane_point.make (16.0, -7.0);
   Plane_point.make (16.0, -3.0);
   Plane_point.make (17.0, -4.0);
   Plane_point.make (5.0, 19.0);
   Plane_point.make (19.0, -8.0);
   Plane_point.make (3.0, 16.0);
   Plane_point.make (12.0, 13.0);
   Plane_point.make (3.0, -4.0);
   Plane_point.make (17.0, 5.0);
   Plane_point.make (-3.0, 15.0);
   Plane_point.make (-3.0, -9.0);
   Plane_point.make (0.0, 11.0);
   Plane_point.make (-9.0, -3.0);
   Plane_point.make (-4.0, -2.0);
   Plane_point.make (12.0, 10.0)]
;;

let hull = plane_convex_hull (List.to_seq example_points)
;;

Array.iter
  (fun p -> (print_string (Plane_point.to_string p);
             print_string " "))
  hull;
print_newline ()
;;

(*------------------------------------------------------------------*)

Andrew’s monotone chain algorithm is a method for computing the convex hull of a set of points in the plane. The algorithm proceeds by sorting the points by their x-coordinates, then using two lower and upper “chains” of points to iteratively construct the convex polygon. The chains start with the two leftmost and rightmost points, respectively, and extend in an alternating fashion as new points are added. At each step, the algorithm determines whether a newly added point creates a convex turn, and if so, whether it creates a left or right turn.

Step-by-step Explanation

  1. Begin by taking a set of points in the plane.
  2. Sort the points by their x-coordinates. If two points have the same x-coordinate, sort them by their y-coordinates.
  3. Create two empty stacks, one for the lower chain of points and one for the upper chain of points.
  4. Push the two leftmost points onto the lower chain stack, and the two rightmost points onto the upper chain stack.
  5. Iterate over the remaining points, adding each to the lower or upper chain as appropriate.
  6. To add a point to a chain, first remove all points from the chain that create a right turn with the new point. These removed points are now not part of the convex hull and can be safely discarded.
  7. Once all right-turning points have been removed, add the new point to the chain.
  8. Repeat steps 6-7 for the other chain.
  9. The points remaining in the lower and upper chains, plus the first point of each chain, form the convex hull of the original set of points.

Complexity Analysis

Andrew’s algorithm has a time complexity of O(n log n) for sorting the points, and O(n) for constructing the hull. The space complexity is also O(n), since the convex hull may contain all the original points. In practice, the algorithm has been shown to be one of the fastest methods for computing the convex hull of a set of points in the plane, and is widely used in applications requiring this computation.