Convex Hull: Graham Scan Algorithm Explained with Code

Have you ever needed to find the boundary of a set of points? Whether you're building a game with collision detection, analyzing geographic data, or simplifying shapes in computer graphics, the convex hull is a fundamental concept. The convex hull of a set of points is the smallest convex polygon that contains all the points. Think of it as the shape formed by a rubber band stretched around the outermost points.

Among the many algorithms to compute the convex hull, the Graham Scan is one of the most elegant and efficient. It runs in O(n log n) time, thanks to an initial sorting step, and uses a simple stack-based approach to build the hull. In this tutorial, you'll learn how the Graham Scan works, step by step, with clear Python code. We'll also cover common pitfalls, debugging strategies, and real-world applications. By the end, you'll be able to implement the algorithm confidently and understand when to use it.

What is the Convex Hull?

The convex hull of a set of points is the smallest convex polygon that contains all the points. In other words, it's the shape you get if you wrap a rubber band around the outermost points. The convex hull is a fundamental concept in computational geometry with applications in pattern recognition, image processing, statistics (e.g., the convex hull of a data set), and game development (e.g., collision detection).

Formally, a set S is convex if for any two points in S, the line segment connecting them lies entirely in S. The convex hull of a set of points P is the intersection of all convex sets containing P. For a finite set of points, the convex hull is a convex polygon whose vertices are a subset of the original points.

# Example: Convex hull of a set of points
points = [(0,0), (1,1), (2,2), (2,0), (1,0), (0,2)]
# The convex hull would be [(0,0), (2,0), (2,2), (0,2)]

Graham Scan Algorithm Overview

The Graham Scan algorithm computes the convex hull of a set of points in O(n log n) time. It works in three main steps:

1. Find the point with the lowest y-coordinate (and leftmost if ties). This point is guaranteed to be on the hull.
2. Sort all other points by polar angle relative to this base point. If two points have the same angle, keep only the farthest one (since the closer one is inside the hull edge).
3. Build the hull using a stack: Iterate through the sorted points. For each point, while the last two points on the stack and the new point make a non-left turn (i.e., clockwise or collinear), pop the stack. Then push the new point.

The algorithm uses the orientation test (cross product) to determine the turn direction. A left turn means the point is part of the hull; a right turn means the previous point is not on the hull.

def orientation(p, q, r):
    # Returns positive if counter-clockwise, negative if clockwise, 0 if collinear
    return (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0])

def graham_scan(points):
    # Step 1: Find lowest point (and leftmost)
    start = min(points, key=lambda p: (p[1], p[0]))
    # Step 2: Sort by polar angle
    sorted_pts = sorted(points, key=lambda p: (math.atan2(p[1]-start[1], p[0]-start[0]), p))
    # Step 3: Build hull
    stack = []
    for p in sorted_pts:
        while len(stack) >= 2 and orientation(stack[-2], stack[-1], p) <= 0:
            stack.pop()
        stack.append(p)
    return stack

Implementing the Orientation Function

The orientation function is the heart of the Graham Scan. It determines whether three points make a left turn, right turn, or are collinear. Given points p, q, r, the cross product (q - p) × (r - p) gives: - Positive: counter-clockwise (left turn) - Negative: clockwise (right turn) - Zero: collinear

In Python, we can implement this as:

``python def orientation(p, q, r): return (q[0] - p[0]) (r[1] - p[1]) - (q[1] - p[1]) (r[0] - p[0]) ``

This function is used in the while loop to pop points that are not part of the convex hull. When orientation is <= 0, we pop because we want to keep only left turns (counter-clockwise). If we want to include collinear points on the hull, we can change the condition to < 0 (strict left turn) and handle collinear points separately.

def orientation(p, q, r):
    val = (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0])
    if val == 0:
        return 0  # collinear
    return 1 if val > 0 else -1  # 1: ccw, -1: cw

# Example
p = (0,0)
q = (1,1)
r = (2,2)
print(orientation(p, q, r))  # 0 (collinear)

Handling Collinear Points

Collinear points on the hull boundary are tricky. If multiple points lie on the same line segment that is part of the hull, we only want the endpoints. The standard Graham Scan pops points when orientation <= 0, which removes collinear points. However, if we want to include all collinear points on the hull (e.g., for a more detailed boundary), we need to adjust.

One approach: When sorting, if two points have the same angle relative to the start, keep only the farthest one. Then, during the stack building, use orientation < 0 (strict left turn) to pop, and handle collinear points by keeping them if they are farther from the start. Alternatively, you can first compute the hull with strict left turns, then add collinear points that lie on the hull edges.

In practice, for most applications, you want the minimal convex hull, so you can safely remove collinear points. But if you need all boundary points, you must modify the algorithm.

import math

def graham_scan_collinear(points):
    start = min(points, key=lambda p: (p[1], p[0]))
    # Sort by angle, and for same angle, by distance (farther first)
    def angle_key(p):
        return (math.atan2(p[1]-start[1], p[0]-start[0]), 
                (p[0]-start[0])**2 + (p[1]-start[1])**2)
    sorted_pts = sorted(points, key=angle_key)
    # Remove duplicates with same angle (keep farthest)
    unique = []
    for p in sorted_pts:
        if unique and orientation(start, unique[-1], p) == 0:
            # Same angle: replace if farther
            if (p[0]-start[0])**2 + (p[1]-start[1])**2 > (unique[-1][0]-start[0])**2 + (unique[-1][1]-start[1])**2:
                unique[-1] = p
        else:
            unique.append(p)
    # Build hull with strict left turn (orientation > 0)
    stack = []
    for p in unique:
        while len(stack) >= 2 and orientation(stack[-2], stack[-1], p) <= 0:
            stack.pop()
        stack.append(p)
    return stack

Complete Python Implementation

Here is a complete, production-ready implementation of the Graham Scan algorithm in Python. It includes error handling for edge cases and uses an epsilon for floating-point comparisons.

```python import math

def orientation(p, q, r, eps=1e-9): val = (q[0] - p[0]) (r[1] - p[1]) - (q[1] - p[1]) (r[0] - p[0]) if abs(val) < eps: return 0 return 1 if val > 0 else -1

def graham_scan(points): if len(points) < 3: return points[:] # Not enough points for a polygon # Find the point with the lowest y (and leftmost if tie) start = min(points, key=lambda p: (p[1], p[0])) # Sort by polar angle and distance def key_func(p): angle = math.atan2(p[1] - start[1], p[0] - start[0]) dist = (p[0] - start[0])2 + (p[1] - start[1])2 return (angle, dist) sorted_pts = sorted(points, key=key_func) # Remove duplicates of same angle (keep farthest) unique = [] for p in sorted_pts: if unique and orientation(start, unique[-1], p) == 0: # Same angle: keep the one farther from start if (p[0]-start[0])2 + (p[1]-start[1])2 > (unique[-1][0]-start[0])2 + (unique[-1][1]-start[1])2: unique[-1] = p else: unique.append(p) # Build hull using stack stack = [] for p in unique: while len(stack) >= 2 and orientation(stack[-2], stack[-1], p) <= 0: stack.pop() stack.append(p) return stack ```

This implementation returns the convex hull points in counter-clockwise order, starting from the lowest point.

import math

def orientation(p, q, r, eps=1e-9):
    val = (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0])
    if abs(val) < eps:
        return 0
    return 1 if val > 0 else -1

def graham_scan(points):
    if len(points) < 3:
        return points[:]
    start = min(points, key=lambda p: (p[1], p[0]))
    def key_func(p):
        angle = math.atan2(p[1] - start[1], p[0] - start[0])
        dist = (p[0] - start[0])**2 + (p[1] - start[1])**2
        return (angle, dist)
    sorted_pts = sorted(points, key=key_func)
    unique = []
    for p in sorted_pts:
        if unique and orientation(start, unique[-1], p) == 0:
            if (p[0]-start[0])**2 + (p[1]-start[1])**2 > (unique[-1][0]-start[0])**2 + (unique[-1][1]-start[1])**2:
                unique[-1] = p
        else:
            unique.append(p)
    stack = []
    for p in unique:
        while len(stack) >= 2 and orientation(stack[-2], stack[-1], p) <= 0:
            stack.pop()
        stack.append(p)
    return stack

# Example usage
points = [(0,0), (1,1), (2,2), (2,0), (1,0), (0,2)]
hull = graham_scan(points)
print(hull)

Testing and Validation

To ensure your Graham Scan implementation is correct, test with various cases:

1. Simple square: Points at (0,0), (0,1), (1,1), (1,0) should return all four.
2. Triangle with interior point: Points (0,0), (1,0), (0,1), (0.2,0.2) should return the triangle.
3. Collinear points: Points (0,0), (1,1), (2,2), (3,3) should return the endpoints (0,0) and (3,3).
4. Duplicate points: Should be handled gracefully (e.g., by removing duplicates first).
5. All points collinear: Should return the two farthest points.
6. Single point: Should return that point.
7. Two points: Should return both.

You can also validate by checking that all original points are inside or on the hull. Use a point-in-polygon test for verification.

def test_graham_scan():
    # Test 1: Square
    points = [(0,0), (0,1), (1,1), (1,0)]
    hull = graham_scan(points)
    assert len(hull) == 4
    # Test 2: Triangle with interior
    points = [(0,0), (1,0), (0,1), (0.2,0.2)]
    hull = graham_scan(points)
    assert len(hull) == 3
    # Test 3: Collinear
    points = [(0,0), (1,1), (2,2), (3,3)]
    hull = graham_scan(points)
    assert len(hull) == 2
    # Test 4: Single point
    points = [(1,1)]
    hull = graham_scan(points)
    assert hull == [(1,1)]
    print("All tests passed!")

test_graham_scan()

Real-World Applications

Convex hull algorithms are used in many fields:

• Collision Detection: In games and robotics, convex hulls approximate object shapes for fast collision checks.
• Image Processing: Hand gesture recognition uses convex hull of hand contour to detect fingers.
• Geographic Information Systems (GIS): Computing the boundary of a set of points (e.g., city limits).
• Statistics: The convex hull of a dataset can be used to identify outliers.
• Pattern Recognition: In machine learning, convex hulls help in support vector machines (SVMs) for finding support vectors.

For example, in a game, you might compute the convex hull of a player's character to detect collisions with obstacles. The hull simplifies the character's shape into a convex polygon, making collision checks faster.

# Example: Collision detection using convex hull
# Assume we have two convex hulls (list of points)
def check_collision(hull1, hull2):
    # Simple check: if any point of one hull is inside the other
    for point in hull1:
        if point_in_convex_polygon(point, hull2):
            return True
    for point in hull2:
        if point_in_convex_polygon(point, hull1):
            return True
    return False

# point_in_convex_polygon uses orientation to check if point is inside
# (implementation omitted for brevity)