834: Convex Hull — Graham Scan

Difficulty: 4 Level: Advanced Find the smallest convex polygon enclosing a set of 2D points in O(n log n) — the gateway algorithm to computational geometry.

The Problem This Solves

The convex hull is the "rubber band stretched around all points" — the minimal convex polygon containing the entire point set. It's needed in collision detection (the hull of a mesh), geographic analysis (the boundary region of a GPS track), robotics (workspace reachability), and computer vision (object shape descriptors). It's also the preprocessing step for many harder geometry problems: rotating calipers for diameter, Minkowski sum, half-plane intersection. Graham scan is O(n log n) — dominated by the sorting step. Once sorted, the sweep is O(n): each point is pushed and popped at most once. This is optimal: any algorithm must sort the points (or implicitly sort them), which requires Ω(n log n) comparisons. Jarvis march is simpler to code but O(n × h) where h is hull size — worse when the hull is large. The cross product is the fundamental geometric primitive underlying every 2D algorithm: three points, one determinant, three cases: left turn (keep), right turn (pop), collinear (depends on goal). Learning to think in cross products unlocks the rest of 2D computational geometry.

The Intuition

Sort all points by polar angle from the bottom-left pivot. Then sweep: maintain a stack where every consecutive triple makes a left turn (counterclockwise). When a new point makes a right turn with the top two stack elements, pop — that middle point is interior to the hull. The CCW invariant is maintained throughout. Cross product of (B-A) × (C-A): positive = C is left of line AB (left turn), negative = right turn, zero = collinear. This single determinant is all you need for any turn test in 2D.

How It Works in Rust

#[derive(Clone, Copy, Debug, PartialEq)]
struct Point { x: f64, y: f64 }

// Cross product of vectors (b-a) and (c-a)
// > 0: left turn (CCW), < 0: right turn (CW), = 0: collinear
fn cross(a: &Point, b: &Point, c: &Point) -> f64 {
 (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)
}

fn convex_hull(mut points: Vec<Point>) -> Vec<Point> {
 let n = points.len();
 if n <= 1 { return points; }

 // Step 1: Find bottom-left pivot (lowest y, then leftmost x)
 let pivot_idx = points.iter().enumerate()
     .min_by(|(_, a), (_, b)| a.y.partial_cmp(&b.y).unwrap()
         .then(a.x.partial_cmp(&b.x).unwrap()))
     .map(|(i, _)| i).unwrap();
 points.swap(0, pivot_idx);
 let pivot = points[0];

 // Step 2: Sort by polar angle from pivot — O(n log n)
 // Collinear points: keep closer one first
 points[1..].sort_by(|a, b| {
     let c = cross(&pivot, a, b);
     if c > 0.0 { std::cmp::Ordering::Less }     // a comes before b: left turn
     else if c < 0.0 { std::cmp::Ordering::Greater }
     else { pivot.dist2(a).partial_cmp(&pivot.dist2(b)).unwrap() }
 });

 // Step 3: Stack sweep — O(n)
 let mut stack: Vec<Point> = Vec::with_capacity(n);
 for &p in &points {
     // Pop points that make a right turn (non-left turn)
     while stack.len() >= 2
         && cross(&stack[stack.len()-2], &stack[stack.len()-1], &p) <= 0.0
     {
         stack.pop();
     }
     stack.push(p);
 }
 stack
}

`partial_cmp` on floats returns `Option<Ordering>` — use `.unwrap()` when you know inputs are non-NaN, or `.unwrap_or(Ordering::Equal)` for defensive code. For integer points, use integer cross products to avoid floating-point precision issues entirely.

What This Unlocks

Collision detection: Hull of a mesh provides the cheapest broad-phase test; two objects can't collide if their hulls don't overlap (test with GJK or SAT).
Rotating calipers: After computing the hull, rotating calipers finds the diameter (farthest pair), minimum bounding rectangle, and Minkowski sum in O(n) — building directly on Graham scan output.
Competitive geometry: Most computational geometry contest problems start with "given n points…" and the first step is usually the convex hull.

Key Differences

Concept	OCaml	Rust
Point type	`type point = { x: float; y: float }`	`struct Point { x: f64, y: f64 }`
Sort by angle	`List.sort` with comparison function	`sort_by` with closure — same semantics
Stack	Recursive list as stack	`Vec<Point>` with `push`/`pop`
Cross product	Pure function `cross a b c`	Function or method on `Point`
Float comparison	Raises exception on NaN	`partial_cmp` returns `Option<Ordering>`

/// Convex Hull: Graham Scan O(n log n).
///
/// Sort by polar angle from bottom-left pivot, then stack-sweep
/// keeping only left turns (CCW orientation).

#[derive(Clone, Copy, Debug, PartialEq)]
struct Point {
    x: f64,
    y: f64,
}

impl Point {
    fn new(x: f64, y: f64) -> Self { Point { x, y } }

    fn dist2(&self, other: &Point) -> f64 {
        let dx = self.x - other.x;
        let dy = self.y - other.y;
        dx * dx + dy * dy
    }
}

/// Cross product of vectors (b-a) and (c-a).
/// > 0: left turn (CCW), < 0: right turn (CW), = 0: collinear.
fn cross(a: &Point, b: &Point, c: &Point) -> f64 {
    (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)
}

/// Graham scan: returns CCW convex hull.
fn convex_hull(mut points: Vec<Point>) -> Vec<Point> {
    let n = points.len();
    if n <= 1 { return points; }

    // Find bottom-most (then left-most) pivot
    let pivot_idx = points
        .iter()
        .enumerate()
        .min_by(|(_, a), (_, b)| {
            a.y.partial_cmp(&b.y).unwrap()
             .then(a.x.partial_cmp(&b.x).unwrap())
        })
        .map(|(i, _)| i)
        .unwrap();
    points.swap(0, pivot_idx);
    let pivot = points[0];

    // Sort by polar angle w.r.t. pivot
    points[1..].sort_by(|a, b| {
        let c = cross(&pivot, a, b);
        if c > 0.0 { std::cmp::Ordering::Less }
        else if c < 0.0 { std::cmp::Ordering::Greater }
        else {
            // Collinear: closer point first
            pivot.dist2(a).partial_cmp(&pivot.dist2(b)).unwrap()
        }
    });

    // Stack-based sweep
    let mut stack: Vec<Point> = Vec::with_capacity(n);
    for p in &points {
        // Pop while last three make a right turn or are collinear
        while stack.len() >= 2
            && cross(&stack[stack.len() - 2], &stack[stack.len() - 1], p) <= 0.0
        {
            stack.pop();
        }
        stack.push(*p);
    }
    stack
}

/// Hull area via shoelace formula.
fn hull_area(hull: &[Point]) -> f64 {
    let n = hull.len();
    if n < 3 { return 0.0; }
    let sum: f64 = (0..n)
        .map(|i| {
            let j = (i + 1) % n;
            hull[i].x * hull[j].y - hull[j].x * hull[i].y
        })
        .sum();
    sum.abs() / 2.0
}

fn main() {
    let points = vec![
        Point::new(0.0, 0.0), Point::new(1.0, 1.0), Point::new(2.0, 2.0),
        Point::new(0.0, 2.0), Point::new(2.0, 0.0), Point::new(1.0, 0.0),
    ];
    let hull = convex_hull(points);
    println!("Hull ({} points):", hull.len());
    for p in &hull { println!("  ({:.1}, {:.1})", p.x, p.y); }

    // Unit square + interior point
    let square = vec![
        Point::new(0.0, 0.0), Point::new(1.0, 0.0),
        Point::new(1.0, 1.0), Point::new(0.0, 1.0),
        Point::new(0.5, 0.5), // interior
    ];
    let hull2 = convex_hull(square);
    println!("\nSquare hull ({} points):", hull2.len());
    println!("Area: {:.1} (expected 1.0)", hull_area(&hull2));
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_square_hull() {
        let pts = vec![
            Point::new(0.0, 0.0), Point::new(1.0, 0.0),
            Point::new(1.0, 1.0), Point::new(0.0, 1.0),
            Point::new(0.5, 0.5),
        ];
        let hull = convex_hull(pts);
        assert_eq!(hull.len(), 4); // Interior point excluded
    }

    #[test]
    fn test_triangle_hull() {
        let pts = vec![
            Point::new(0.0, 0.0), Point::new(2.0, 0.0),
            Point::new(1.0, 2.0), Point::new(1.0, 0.5),
        ];
        let hull = convex_hull(pts);
        assert_eq!(hull.len(), 3);
    }

    #[test]
    fn test_collinear_points() {
        // All collinear: hull should have 2 or more points
        let pts = vec![
            Point::new(0.0, 0.0), Point::new(1.0, 1.0), Point::new(2.0, 2.0),
        ];
        let hull = convex_hull(pts);
        assert!(hull.len() >= 2);
    }

    #[test]
    fn test_single_point() {
        let pts = vec![Point::new(1.0, 2.0)];
        let hull = convex_hull(pts);
        assert_eq!(hull.len(), 1);
    }

    #[test]
    fn test_all_hull_points_are_extreme() {
        let pts = vec![
            Point::new(0.0, 0.0), Point::new(4.0, 0.0),
            Point::new(4.0, 4.0), Point::new(0.0, 4.0),
            Point::new(2.0, 2.0), Point::new(1.0, 1.0),
        ];
        let hull = convex_hull(pts.clone());
        // Every hull point should have been in the original set
        for h in &hull {
            assert!(pts.contains(h), "{h:?} not in original set");
        }
    }

    #[test]
    fn test_area_unit_square() {
        let pts = vec![
            Point::new(0.0, 0.0), Point::new(1.0, 0.0),
            Point::new(1.0, 1.0), Point::new(0.0, 1.0),
        ];
        let hull = convex_hull(pts);
        let area = hull_area(&hull);
        assert!((area - 1.0).abs() < 1e-9, "area = {area}");
    }

    #[test]
    fn test_cross_product() {
        let a = Point::new(0.0, 0.0);
        let b = Point::new(1.0, 0.0);
        let c = Point::new(0.0, 1.0);
        assert!(cross(&a, &b, &c) > 0.0); // CCW / left turn
        assert!(cross(&a, &c, &b) < 0.0); // CW / right turn
    }
}

(* Convex Hull: Graham Scan in OCaml *)

type point = { x: float; y: float }

let dist2 a b =
  let dx = a.x -. b.x and dy = a.y -. b.y in
  dx *. dx +. dy *. dy

(* Cross product of vectors (b-a) and (c-a) *)
(* > 0: CCW (left turn), < 0: CW (right turn), = 0: collinear *)
let cross a b c =
  (b.x -. a.x) *. (c.y -. a.y) -. (b.y -. a.y) *. (c.x -. a.x)

(* Graham scan: returns convex hull in CCW order *)
let convex_hull (points : point list) : point list =
  let n = List.length points in
  if n <= 1 then points
  else begin
    let pts = Array.of_list points in
    (* Find bottom-most (then left-most) point as pivot *)
    let pivot = Array.fold_left (fun best p ->
      if p.y < best.y || (p.y = best.y && p.x < best.x) then p else best
    ) pts.(0) pts in
    (* Sort by polar angle with respect to pivot *)
    let sorted = Array.copy pts in
    Array.sort (fun a b ->
      if a = pivot then -1
      else if b = pivot then 1
      else
        let c = cross pivot a b in
        if c > 0.0 then -1
        else if c < 0.0 then 1
        else compare (dist2 pivot a) (dist2 pivot b)
    ) sorted;
    (* Stack-based Graham scan *)
    let stack = Array.make n {x=0.0; y=0.0} in
    let top = ref (-1) in
    Array.iter (fun p ->
      while !top >= 1 && cross stack.(!top - 1) stack.(!top) p <= 0.0 do
        decr top
      done;
      incr top;
      stack.(!top) <- p
    ) sorted;
    Array.to_list (Array.sub stack 0 (!top + 1))
  end

let () =
  let points = [
    {x=0.0;y=0.0}; {x=1.0;y=1.0}; {x=2.0;y=2.0};
    {x=0.0;y=2.0}; {x=2.0;y=0.0}; {x=1.0;y=0.0};
  ] in
  let hull = convex_hull points in
  Printf.printf "Hull (%d points):\n" (List.length hull);
  List.iter (fun p -> Printf.printf "  (%.1f, %.1f)\n" p.x p.y) hull;

  let square = [
    {x=0.0;y=0.0}; {x=1.0;y=0.0}; {x=1.0;y=1.0}; {x=0.0;y=1.0};
    {x=0.5;y=0.5}  (* interior point *)
  ] in
  let h2 = convex_hull square in
  Printf.printf "\nSquare hull (%d points, interior excluded):\n" (List.length h2);
  List.iter (fun p -> Printf.printf "  (%.1f, %.1f)\n" p.x p.y) h2