Master computational geometry fundamentals — cross product for orientation testing, point-in-polygon detection, convex hull algorithms (Graham scan and Jarvis march), closest pair of points, line segment intersection, and applications in graphics and competitive programming.
Introduction
Computational geometry is the branch of computer science that studies algorithms for solving geometric problems — problems involving points, lines, polygons, and other geometric objects. It is fundamental to computer graphics, robotics, GIS (geographic information systems), and competitive programming.
The most powerful tool in computational geometry is the cross product, which tells you the orientation of three points (clockwise, counterclockwise, or collinear) in O(1) time.
The 2D cross product of vectors (a, b) and (c, d) is: a×d - b×c
For three points P, Q, R:
| If result > 0 | R is to the LEFT of PQ (counterclockwise turn) |
| If result < 0 | R is to the RIGHT of PQ (clockwise turn) |
| If result = 0 | P, Q, R are COLLINEAR |
def cross_product(O: Point, A: Point, B: Point) -> float:
"""
2D cross product of vectors OA and OB.
Returns:
> 0 if OAB makes a counterclockwise (left) turn
< 0 if OAB makes a clockwise (right) turn
= 0 if O, A, B are collinear
"""
return (A[0] - O[0]) * (B[1] - O[1]) - (A[1] - O[1]) * (B[0] - O[0])
# Test
P, Q, R = (0, 0), (1, 0), (0, 1)
print(cross_product(P, Q, R)) # 1 > 0 → counterclockwise
P, Q, R = (0, 0), (1, 0), (0, -1)
print(cross_product(P, Q, R)) # -1 < 0 → clockwise
P, Q, R = (0, 0), (1, 0), (2, 0)
print(cross_product(P, Q, R)) # 0 → collinear
Convex Hull — Graham Scan: O(n log n)
The convex hull is the smallest convex polygon that contains all the given points. Think of stretching a rubber band around all the points and letting it snap — the shape it forms is the convex hull.
Graham Scan algorithm:
- Find the bottom-most point (lowest y, leftmost if tie) — call it P0
- Sort all other points by polar angle relative to P0
- Process each point: maintain a stack where each turn is counterclockwise
- Pop clockwise turns — the remaining stack is the convex hull
def convex_hull(points: List[Point]) -> List[Point]:
"""
Graham Scan convex hull algorithm.
Time: O(n log n) — dominated by sorting
Space: O(n)
Returns: List of hull vertices in counterclockwise order
"""
n = len(points)
if n <= 2:
return points
# Find bottom-most point (leftmost if tie)
pivot = min(points, key=lambda p: (p[1], p[0]))
def polar_angle_key(p):
"""Sort key: polar angle from pivot, then distance."""
if p == pivot:
return (float('-inf'), 0)
cp = cross_product(pivot, p, pivot) # Not useful here
# Compute angle using atan2
import math
angle = math.atan2(p[1] - pivot[1], p[0] - pivot[0])
dist = distance_sq(pivot, p)
return (angle, dist)
# Sort by polar angle
sorted_pts = sorted(points, key=polar_angle_key)
# Graham scan: build lower and upper hull
stack = []
for p in sorted_pts:
# While last three points make a clockwise or collinear turn, pop
while len(stack) >= 2 and cross_product(stack[-2], stack[-1], p) <= 0:
stack.pop()
stack.append(p)
return stack
# Alternative: using the monotone chain algorithm (more numerically stable)
def convex_hull_monotone(points: List[Point]) -> List[Point]:
"""
Monotone chain convex hull.
Builds lower and upper hulls separately.
Time: O(n log n), more numerically stable than Graham scan.
"""
points = sorted(set(points)) # Sort lexicographically, remove duplicates
if len(points) <= 1:
return points
# Build lower hull
lower = []
for p in points:
while len(lower) >= 2 and cross_product(lower[-2], lower[-1], p) <= 0:
lower.pop()
lower.append(p)
# Build upper hull
upper = []
for p in reversed(points):
while len(upper) >= 2 and cross_product(upper[-2], upper[-1], p) <= 0:
upper.pop()
upper.append(p)
# Concatenate (remove duplicate endpoints)
return lower[:-1] + upper[:-1]
points = [(0,0),(1,1),(2,2),(0,2),(2,0),(1,0),(0,1)]
hull = convex_hull_monotone(points)
print("Convex hull:", hull) # [(0,0),(2,0),(2,2),(0,2)]
Point in Polygon — Ray Casting: O(n)
Problem: Given a polygon and a point P, determine if P is inside the polygon.
Ray casting algorithm: Cast a ray from P in any direction (e.g., right). Count how many polygon edges it crosses. If odd → inside; if even → outside.
def point_in_polygon(point: Point, polygon: List[Point]) -> bool:
"""
Ray casting algorithm.
Casts a horizontal ray to the right from point.
Counts intersections with polygon edges.
Time: O(n), Space: O(1)
"""
x, y = point
n = len(polygon)
inside = False
j = n - 1
for i in range(n):
xi, yi = polygon[i]
xj, yj = polygon[j]
# Does the ray from (x,y) going right intersect edge (xi,yi)-(xj,yj)?
if ((yi > y) != (yj > y)) and (x < (xj - xi) * (y - yi) / (yj - yi) + xi):
inside = not inside
j = i
return inside
# Test
square = [(0,0), (4,0), (4,4), (0,4)]
print(point_in_polygon((2, 2), square)) # True — inside
print(point_in_polygon((5, 2), square)) # False — outside
print(point_in_polygon((2, -1), square)) # False — below
triangle = [(0,0), (4,0), (2,4)]
print(point_in_polygon((2, 1), triangle)) # True
print(point_in_polygon((0, 3), triangle)) # False
Line Segment Intersection: O(1)
Problem: Do two line segments intersect?
def segments_intersect(p1: Point, p2: Point, p3: Point, p4: Point) -> bool:
"""
Check if segment p1-p2 intersects segment p3-p4.
Uses cross product orientation tests.
Time: O(1)
"""
d1 = cross_product(p3, p4, p1)
d2 = cross_product(p3, p4, p2)
d3 = cross_product(p1, p2, p3)
d4 = cross_product(p1, p2, p4)
if ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and \
((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0)):
return True
# Collinear cases
def on_segment(p: Point, q: Point, r: Point) -> bool:
return (min(p[0], r[0]) <= q[0] <= max(p[0], r[0]) and
min(p[1], r[1]) <= q[1] <= max(p[1], r[1]))
if d1 == 0 and on_segment(p3, p1, p4): return True
if d2 == 0 and on_segment(p3, p2, p4): return True
if d3 == 0 and on_segment(p1, p3, p2): return True
if d4 == 0 and on_segment(p1, p4, p2): return True
return False
print(segments_intersect((0,0),(4,4),(0,4),(4,0))) # True — X-shaped
print(segments_intersect((0,0),(2,2),(3,3),(5,5))) # False — parallel/same line no overlap
print(segments_intersect((0,0),(1,1),(2,0),(3,1))) # False — parallel
Closest Pair of Points: O(n log n)
Problem: Given n points, find the pair with minimum distance.
Divide and conquer beats O(n²) brute force:
def closest_pair(points: List[Point]) -> Tuple[float, Point, Point]:
"""
Find closest pair of points using divide and conquer.
Time: O(n log n), Space: O(n)
"""
def brute_force(pts):
min_dist = float('inf')
pair = None
for i in range(len(pts)):
for j in range(i+1, len(pts)):
d = distance(pts[i], pts[j])
if d < min_dist:
min_dist, pair = d, (pts[i], pts[j])
return min_dist, pair
def closest_rec(pts_sorted_x, pts_sorted_y):
n = len(pts_sorted_x)
if n <= 3:
return brute_force(pts_sorted_x)
mid = n // 2
mid_point = pts_sorted_x[mid]
left_x = pts_sorted_x[:mid]
right_x = pts_sorted_x[mid:]
left_y = [p for p in pts_sorted_y if p[0] <= mid_point[0]]
right_y = [p for p in pts_sorted_y if p[0] > mid_point[0]]
d_left, pair_left = closest_rec(left_x, left_y)
d_right, pair_right = closest_rec(right_x, right_y)
if d_left < d_right:
d, best_pair = d_left, pair_left
else:
d, best_pair = d_right, pair_right
# Check strip around the dividing line
strip = [p for p in pts_sorted_y if abs(p[0] - mid_point[0]) < d]
for i in range(len(strip)):
j = i + 1
while j < len(strip) and (strip[j][1] - strip[i][1]) < d:
dist = distance(strip[i], strip[j])
if dist < d:
d, best_pair = dist, (strip[i], strip[j])
j += 1
return d, best_pair
pts_x = sorted(points, key=lambda p: p[0])
pts_y = sorted(points, key=lambda p: p[1])
return closest_rec(pts_x, pts_y)
points = [(2,3),(12,30),(40,50),(5,1),(12,10),(3,4)]
dist, (p1, p2) = closest_pair(points)
print(f"Closest pair: {p1} and {p2}, distance: {dist:.4f}")
def polygon_area(polygon: List[Point]) -> float:
"""
Calculate area of polygon using the Shoelace formula.
Works for any simple (non-self-intersecting) polygon.
Returns positive area regardless of vertex order.
Time: O(n)
"""
n = len(polygon)
area = 0.0
for i in range(n):
j = (i + 1) % n
area += polygon[i][0] * polygon[j][1]
area -= polygon[j][0] * polygon[i][1]
return abs(area) / 2.0
square = [(0,0), (4,0), (4,4), (0,4)]
print(polygon_area(square)) # 16.0
triangle = [(0,0), (4,0), (0,3)]
print(polygon_area(triangle)) # 6.0
Applications
| Application | Algorithm | Complexity |
|---|
| Convex hull | Graham Scan / Monotone Chain | O(n log n) |
| Point in polygon | Ray casting | O(n) |
| Segment intersection | Cross product | O(1) per pair |
| Closest pair | Divide and conquer | O(n log n) |
| Polygon area | Shoelace formula | O(n) |
| All intersections | Sweep line (Bentley-Ottmann) | O((n+k) log n) |
Real-world uses:
- Computer graphics: collision detection, clipping, visibility
- GPS/Maps: polygon containment, route planning
- Robotics: obstacle avoidance, motion planning
- GIS: geographic queries, map overlays
Summary
Computational geometry problems are solved using a small set of fundamental tools:
- Cross product: Determines orientation (left/right/straight turn) of three points — the most important tool
- Convex hull: Smallest convex polygon containing all points — Graham Scan O(n log n)
- Point in polygon: Ray casting algorithm O(n)
- Segment intersection: Cross product orientation tests O(1)
- Closest pair: Divide and conquer O(n log n) — better than O(n²) brute force
- Polygon area: Shoelace formula O(n)
Master the cross product and orientation test — they underlie virtually every computational geometry algorithm.
*Advanced Topics section complete. All 8 files written.*