Polygons & Proximity

The first three lessons built the planar toolkit: the orientation primitive, the convex hull, and the plane sweep. This lesson applies those tools to four problems every geometry library ships, each fast and exact by a counting or sign argument: the closest pair of a point set, the point-in-polygon test, the area of a polygon, and the diameter and width of a convex polygon. Three of them need nothing beyond the cross product; the fourth turns the hull into a linear scan by walking two pointers in lock-step.

Closest pair of points

The brute force compares all $\left(\genfrac{}{}{0ex}{}{n}{2}\right)$ pairs in $\Theta (n^2)$. The classic improvement is divide-and-conquer, and it is the geometry that makes it work: a packing bound caps the cross-boundary work at a constant per point. This is the full treatment — the Selection lesson introduces it as a divide-and-conquer classic alongside sorting, and the sweep-line lesson gives a one-pass alternative on a dynamic strip.

Divide. Sort by $x$ and split at the median $x$ into a left half $L$ and a right half $R$ by a vertical line $x=m$. Conquer. Recurse on each half; let ${\delta }_L$ and ${\delta }_R$ be the best distances found, and $\delta =\min ({\delta }_L,{\delta }_R)$. Combine. The closest pair overall is either one of these two or a pair that straddles the line — one point in $L$, one in $R$ — at distance $<\delta$. The whole difficulty is bounding that straddling case.

A straddling pair closer than $\delta$ must have both endpoints within horizontal distance $\delta$ of the line, so both lie in a vertical strip of width $2\delta$ centred on $x=m$. Discard everything outside the strip; sort the survivors by $y$ (in practice, filter the already-$y$-sorted list); then walk up the strip. The decisive fact is that each point need only be compared against a constant number of its successors in $y$-order.

The combine step therefore makes $O(1)$ distance checks per strip point, for $O(n)$ work. The recurrence is the familiar one,

which the Master Theorem solves in its balanced case ($a=b^d$ with $a=2,b=2,d=1$), giving $T(n)=\Theta (n\log n)$. The one implementation subtlety is the inner sort: a naive $y$-sort of the strip each level would add a $\log n$ factor. Sorting once by $y$ up front and filtering that order at each level keeps the combine linear, so the bound holds.¹

The distance comparisons can stay exact: comparing $\mathrm{dist}(p,q)$ against $\delta$ is comparing squared distances $(p_x-q_x{)}^2+(p_y-q_y{)}^2$ against ${\delta }^2$, which are integers for integer inputs — no square root, no rounding, the same discipline as the orientation primitive.

from typing import NamedTuple

class Point(NamedTuple):
  """
    A planar point with integer x and y components.\n
    Integer coordinates keep every derived primitive exact, so the sign\n
    tests built on them never misfire from floating-point rounding.\n
  """
  x: int
  y: int

def cross(origin: Point, first: Point, second: Point) -> int:
  """
    Orientation of the ordered triple origin -> first -> second.\n
    Returns (first - origin) x (second - origin): > 0 left turn,\n
    < 0 right turn, 0 collinear. Twice the signed triangle area.\n
  """
  return (
    (first.x - origin.x) * (second.y - origin.y)
    - (first.y - origin.y) * (second.x - origin.x)
  )

def squared_distance(first: Point, second: Point) -> int:
  """
    Squared Euclidean distance between two points.\n
    Squared so the comparison stays exact integer arithmetic — taking the\n
    root is never needed to decide which of two pairs is farther apart.\n
  """
  delta_x: int = first.x - second.x
  delta_y: int = first.y - second.y
  return delta_x * delta_x + delta_y * delta_y

def on_segment(start: Point, end: Point, query: Point) -> bool:
  """
    Whether `query` lies on the closed segment from `start` to `end`.\n
    The point must be collinear with the segment and within its bounding\n
    box; both checks are exact integer tests.\n
  """
  # must be collinear first, then inside the segment's bounding box.
  if cross(start, end, query) != 0:
    return False
  return (
    min(start.x, end.x) <= query.x <= max(start.x, end.x)
    and min(start.y, end.y) <= query.y <= max(start.y, end.y)
  )

from typing import Sequence

from geometry import Point, squared_distance

# A candidate result: the squared distance and the pair achieving it.
Pair = tuple[int, Point, Point]

def _brute_force(points: Sequence[Point]) -> Pair:
  """
    Closest pair among at most three points by checking every pair.\n
    The recursion's base case, where divide-and-conquer is not worth it.\n
  """
  # seed with the first pair, then check every remaining pair.
  best: Pair = (squared_distance(points[0], points[1]), points[0], points[1])
  for left_index in range(len(points)):
    for right_index in range(left_index + 1, len(points)):
      first: Point = points[left_index]
      second: Point = points[right_index]

      distance: int = squared_distance(first, second)
      if distance < best[0]:
        best = (distance, first, second)

  return best

def _strip_closest(strip: Sequence[Point], best: Pair) -> Pair:
  """
    Improve `best` using pairs that straddle the dividing line.\n
    `strip` holds the survivors within delta of the line, sorted by y. By\n
    the strip-packing lemma each point need only be compared against its\n
    next seven successors in y-order, so this scan is linear.\n
  """
  result: Pair = best
  for lower_index in range(len(strip)):
    upper_index: int = lower_index + 1

    # compare each point with at most its next seven y-successors.
    while upper_index < len(strip) and upper_index <= lower_index + 7:
      lower: Point = strip[lower_index]
      upper: Point = strip[upper_index]

      # the y-gap alone already exceeds the bound — no later point can help.
      gap_y: int = upper.y - lower.y
      if gap_y * gap_y >= result[0]:
        break

      # keep the closer straddling pair.
      distance: int = squared_distance(lower, upper)
      if distance < result[0]:
        result = (distance, lower, upper)
      upper_index += 1

  return result

def _closest(by_x: Sequence[Point], by_y: Sequence[Point]) -> Pair:
  """
    Recursive driver: `by_x` is sorted by x, `by_y` the same points by y.\n
    Returns the squared distance of the closest pair and the pair itself.\n
  """
  # small enough to settle directly.
  count: int = len(by_x)
  if count <= 3:
    return _brute_force(by_x)

  # split by x at the median; ties on the divider stay in the left half.
  middle: int = count // 2
  divider_x: int = by_x[middle].x
  left_by_x: Sequence[Point] = by_x[:middle]
  right_by_x: Sequence[Point] = by_x[middle:]

  # partition the y-order to each side while preserving it, so neither
  # recursive call has to re-sort. Identity tags carry the x-split's tie-break.
  left_labels: set[int] = {id(point) for point in left_by_x}
  left_by_y: list[Point] = [p for p in by_y if id(p) in left_labels]
  right_by_y: list[Point] = [p for p in by_y if id(p) not in left_labels]

  # recurse on each half and keep the better of the two.
  left_best: Pair = _closest(left_by_x, left_by_y)
  right_best: Pair = _closest(right_by_x, right_by_y)
  best: Pair = left_best if left_best[0] <= right_best[0] else right_best

  # the strip: y-sorted points within sqrt(best) of the dividing line.
  delta_squared: int = best[0]
  strip: list[Point] = [
    point
    for point in by_y
    if (point.x - divider_x) * (point.x - divider_x) < delta_squared
  ]
  return _strip_closest(strip, best)

def closest_pair(points: Sequence[Point]) -> tuple[Point, Point]:
  """
    The two points of `points` at smallest Euclidean distance.\n
    Requires at least two points. Duplicate coordinates are allowed and\n
    yield a zero-distance pair. Runs in Theta(n log n).\n
  """
  if len(points) < 2:
    raise ValueError("closest_pair needs at least two points")

  # pre-sort once by x and by y; the recursion reuses both orders.
  by_x: list[Point] = sorted(points)
  by_y: list[Point] = sorted(points, key=lambda point: (point.y, point.x))

  _, first, second = _closest(by_x, by_y)
  return (first, second)

Point in polygon

Given a simple polygon $P$ (vertices $V_0,\dots ,V_{n-1}$ in order) and a query point $q$, is $q$ inside, outside, or on the boundary? Two tests dominate; both were named in the primitives lesson, and we make them precise here, boundary caveats included.

Ray casting (the parity test)

Shoot a ray from $q$ in a fixed direction — conventionally along $+x$ — and count how many polygon edges it crosses. The Jordan curve theorem guarantees that a ray from an interior point crosses the boundary an odd number of times, and a ray from an exterior point an even number. So an odd crossing count means inside.

The arithmetic per edge $\overline{V_i{V}_{i+1}}$ is a two-line test: the edge can cross the horizontal ray at $y=q_y$ only if one endpoint is strictly above $q_y$ and the other not above it, and then the crossing $x$ must lie to the right of $q_x$. Both conditions reduce to sign comparisons; the second avoids dividing by cross-multiplying.

The pitfalls are all degenerate rays — those grazing a vertex or running along an edge — and they are what make naive implementations report nonsense.

from enum import Enum
from typing import Sequence
from geometry import Point, cross, on_segment


class Location(Enum):
  """
    Where a query point sits relative to a polygon.\n
  """
  INSIDE = "inside"
  OUTSIDE = "outside"
  BOUNDARY = "boundary"


def _ray_crossings(polygon: Sequence[Point], query: Point) -> int:
  """
    Parity of +x-ray crossings: 1 if odd (inside), 0 if even (outside).\n
    Each edge is half-open in y — counted only when exactly one endpoint is\n
    strictly above the query — so a grazed vertex counts once and horizontal\n
    edges are ignored. The crossing-to-the-right check cross-multiplies\n
    instead of dividing, keeping it exact.\n
  """
  vertex_count: int = len(polygon)
  is_inside: bool = False
  for index in range(vertex_count):
    start: Point = polygon[index]
    end: Point = polygon[(index + 1) % vertex_count]

    # half-open rule: exactly one endpoint strictly above the ray's height.
    if (start.y > query.y) != (end.y > query.y):

      # edge's x at y = query.y vs query.x, cross-multiplied to avoid dividing.
      numerator: int = (query.y - start.y) * (end.x - start.x)
      denominator: int = end.y - start.y
      offset: int = (query.x - start.x) * denominator

      # the multiply flips the inequality when the denominator is negative.
      if denominator > 0:
        crosses_right: bool = numerator > offset
      else:
        crosses_right = numerator < offset

      if crosses_right:
        is_inside = not is_inside

  return 1 if is_inside else 0


def point_in_polygon(
  polygon: Sequence[Point],
  query: Point,
  boundary_is_inside: bool = True,
) -> Location:
  """
    Classify `query` against a simple `polygon` via ray casting.\n
    Vertices are listed in order (either orientation). A boundary point is\n
    reported as BOUNDARY; set `boundary_is_inside` False only changes how a\n
    caller may treat that label, not the classification. Inside is decided\n
    by an odd +x-ray crossing count.\n
  """
  if on_boundary(polygon, query):
    return Location.BOUNDARY
  if _ray_crossings(polygon, query) == 1:
    return Location.INSIDE
  return Location.OUTSIDE

Winding number

The parity test is blind to how the boundary wraps around $q$, which is why it breaks on self-intersecting polygons. The winding number $w(q)$ counts the net number of full turns the boundary makes around $q$ as it is traversed once. For a simple polygon $w(q)=\pm 1$ inside (sign by orientation) and $0$ outside; for a self-intersecting polygon the magnitude can exceed $1$, and the nonzero rule ($q$ inside iff $w(q)\ne 0$) is the one graphics systems use.

Computed naively this would sum signed angles — irrational, floating point. The exact version sums signed edge crossings of the $+x$ ray instead, weighting by direction: an edge crossing the ray upward contributes $+1$, downward$-1$. The sum is $w(q)$, computed entirely from orientation signs.

Both tests are $O(n)$ per query. When the polygon is convex and many queries are expected, the primitives lesson's $O(\log n)$ fan binary search wins after an $O(n)$ preprocess.

from enum import Enum
from typing import Sequence
from geometry import Point, cross, on_segment


class Location(Enum):
  """
    Where a query point sits relative to a polygon.\n
  """
  INSIDE = "inside"
  OUTSIDE = "outside"
  BOUNDARY = "boundary"


def point_in_polygon_winding(
  polygon: Sequence[Point],
  query: Point,
) -> Location:
  """
    Classify `query` by the nonzero-winding rule.\n
    Reports BOUNDARY for points on an edge, INSIDE when the winding number\n
    is nonzero, OUTSIDE otherwise. Stays correct for self-intersecting\n
    polygons, where the parity test can misclassify.\n
  """
  if on_boundary(polygon, query):
    return Location.BOUNDARY
  if winding_number(polygon, query) != 0:
    return Location.INSIDE
  return Location.OUTSIDE

Polygon area: the shoelace formula

The primitives lesson introduced the shoelace formula for the area of a simple polygon with vertices $V_0,\dots ,V_{n-1}$ in order (indices mod $n$):

Here we prove it. The cleanest argument is the trapezoid sum: each edge, with the $x$-axis, bounds a trapezoid whose signed area is added as the boundary is walked, and the trapezoids below exterior edges cancel against those below interior edges, leaving exactly the enclosed region.

The sum is an integer for integer vertices, so the area is an exact half-integer.² Dropping the absolute value keeps the sign, which gives the winding direction for free: positive means the vertices are counterclockwise, negative clockwise — the cheapest orientation test for a whole polygon.

from fractions import Fraction
from typing import Sequence

from geometry import Point

def twice_signed_area(polygon: Sequence[Point]) -> int:
  """
    The shoelace sum, equal to twice the signed area of `polygon`.\n
    Positive for counterclockwise vertices, negative for clockwise, zero\n
    for a degenerate (collinear) polygon. Exact integer arithmetic.\n
  """
  # sum each edge's cross product around the closed boundary (wraps at the end).
  vertex_count: int = len(polygon)
  return sum(
    polygon[index].x * polygon[(index + 1) % vertex_count].y
    - polygon[(index + 1) % vertex_count].x * polygon[index].y
    for index in range(vertex_count)
  )

def signed_area(polygon: Sequence[Point]) -> Fraction:
  """
    The signed area of `polygon` as an exact half-integer.\n
    Sign carries the winding direction: positive counterclockwise, negative\n
    clockwise.\n
  """
  return Fraction(twice_signed_area(polygon), 2)

def polygon_area(polygon: Sequence[Point]) -> Fraction:
  """
    The (unsigned) area enclosed by `polygon`, an exact half-integer.\n
  """
  return abs(signed_area(polygon))

def is_counterclockwise(polygon: Sequence[Point]) -> bool:
  """
    Whether `polygon`'s vertices are listed counterclockwise.\n
    The cheapest orientation test for a whole polygon: the sign of the\n
    shoelace sum, with no per-triple turn checks.\n
  """
  return twice_signed_area(polygon) > 0

Rotating calipers: diameter and width

The convex hull lesson closed by promising that the hull turns many extremal queries into linear scans. The technique is rotating calipers, and the cleanest instance is the diameter — the greatest distance between any two points of a set.

Checking all $\left(\genfrac{}{}{0ex}{}{h}{2}\right)$ hull-vertex pairs is $\Theta (h^2)$. Rotating calipers do it in $O(h)$ by exploiting monotonicity: imagine two parallel lines (calipers) squeezing the hull, and rotate them together through $180°$. As the calipers turn, each touches a hull vertex; the contact points advance around the hull in step, never backtracking, because the support direction is monotone. Every antipodal pair is realized as some caliper orientation, so walking the two contact pointers once around the hull visits all candidate pairs.

Concretely, for hull vertices in counterclockwise order, the area of the triangle $V_i{V}_{i+1}{V}_j$ (a cross product) is maximized, over $j$, at the vertex farthest from edge $\overline{V_i{V}_{i+1}}$. As $i$ advances by one edge, the farthest vertex $j$ only moves forward, so a single linear walk of $j$ keeps pace with $i$. Each antipodal pair encountered is a diameter candidate.

Because $j$ only ever advances and wraps once, the inner while does $O(h)$ total work across all $i$ — the pointer makes a single trip around the hull. The diameter scan is therefore $O(h)$, and $O(n\log n)$ end-to-end once the hull is built, against $\Theta (n^2)$ for all-pairs.³ The area calls are cross products, so distances enter only at the final max (and even there, comparing squared distances keeps it exact).

The same caliper sweep answers the dual question, the width — the smallest distance between two parallel supporting lines, i.e. the narrowest slab containing the polygon. The minimum-width slab always has one line flush with a hull edge, so rotate one caliper along each edge in turn and track the antipodal vertex's perpendicular distance (a cross product divided by the edge length); the minimum over edges is the width. It drives smallest-enclosing-rectangle and collision-clearance queries, again in $O(h)$ after the hull.

import math
from typing import Sequence

from geometry import Point, cross, squared_distance

def convex_hull(points: Sequence[Point]) -> list[Point]:
  """
    Convex hull of `points` in counterclockwise order (Andrew's monotone\n
    chain). Only true corners are kept; collinear interior points drop out.\n
    Fewer than three distinct points return the deduplicated, sorted points.\n
  """
  # too few distinct points to form a hull — return them sorted.
  ordered: list[Point] = sorted(set(points))
  if len(ordered) < 3:
    return ordered

  # sweep a half-hull, popping any vertex that makes a non-left turn.
  def build_chain(sweep: Sequence[Point]) -> list[Point]:
    chain: list[Point] = []
    for candidate in sweep:
      while len(chain) >= 2 and cross(chain[-2], chain[-1], candidate) <= 0:
        chain.pop()
      chain.append(candidate)
    return chain

  # stitch the lower and upper chains, dropping each one's shared endpoint.
  lower: list[Point] = build_chain(ordered)
  upper: list[Point] = build_chain(list(reversed(ordered)))
  return lower[:-1] + upper[:-1]

def _twice_triangle_area(first: Point, second: Point, apex: Point) -> int:
  """
    Twice the unsigned area of triangle first-second-apex.\n
    The distance of `apex` from line first-second, up to the fixed edge\n
    length — exactly what the caliper pointer maximizes.\n
  """
  return abs(cross(first, second, apex))

def diameter(points: Sequence[Point]) -> tuple[Point, Point]:
  """
    The farthest pair of `points` (largest Euclidean distance).\n
    Builds the hull, then advances a single antipodal pointer around it,\n
    checking each candidate antipodal pair. Requires at least two distinct\n
    points; runs in O(n log n).\n
  """
  # degenerate hulls: error on a single point, take the segment's endpoints.
  hull: list[Point] = convex_hull(points)
  hull_size: int = len(hull)
  if hull_size < 2:
    raise ValueError("diameter needs at least two distinct points")
  if hull_size == 2:
    return (hull[0], hull[1])

  best_squared: int = -1
  best_pair: tuple[Point, Point] = (hull[0], hull[1])
  antipodal: int = 1
  for edge_index in range(hull_size):
    edge_start: Point = hull[edge_index]
    edge_end: Point = hull[(edge_index + 1) % hull_size]

    # advance the antipodal pointer while the next vertex is farther from
    # the current edge; it never backtracks across the whole sweep.
    while (
      _twice_triangle_area(edge_start, edge_end, hull[(antipodal + 1) % hull_size])
      > _twice_triangle_area(edge_start, edge_end, hull[antipodal])
    ):
      antipodal = (antipodal + 1) % hull_size

    # both endpoints of the edge against the antipodal vertex are candidates.
    for endpoint in (edge_start, edge_end):
      distance: int = squared_distance(endpoint, hull[antipodal])
      if distance > best_squared:
        best_squared = distance
        best_pair = (endpoint, hull[antipodal])

  return best_pair

def width(points: Sequence[Point]) -> float:
  """
    The width of `points`: the narrowest slab (pair of parallel supporting\n
    lines) containing them all. One line is always flush with a hull edge,\n
    so rotate one caliper along each edge and track the antipodal vertex's\n
    perpendicular distance; the minimum over edges is the width. Returns 0.0\n
    when the hull is degenerate (a point or a segment). Runs in O(n log n).\n
  """
  # a point or segment has no width.
  hull: list[Point] = convex_hull(points)
  hull_size: int = len(hull)
  if hull_size < 3:
    return 0.0

  best_width: float = math.inf
  antipodal: int = 1
  for edge_index in range(hull_size):
    edge_start: Point = hull[edge_index]
    edge_end: Point = hull[(edge_index + 1) % hull_size]

    # carry the antipodal pointer forward to the vertex farthest from the edge.
    while (
      _twice_triangle_area(edge_start, edge_end, hull[(antipodal + 1) % hull_size])
      > _twice_triangle_area(edge_start, edge_end, hull[antipodal])
    ):
      antipodal = (antipodal + 1) % hull_size

    # perpendicular distance = (twice triangle area) / (edge length).
    twice_area: int = _twice_triangle_area(edge_start, edge_end, hull[antipodal])
    edge_length: float = math.sqrt(squared_distance(edge_start, edge_end))
    best_width = min(best_width, twice_area / edge_length)

  return best_width

Voronoi, Delaunay, and spatial trees

Closest-pair answers one proximity question for one query. Real systems ask it repeatedly — nearest hospital, nearest cell tower, nearest neighbor for a classifier — and want a structure, built once, that answers each query fast. The two canonical structures are duals of each other.

The Voronoi diagram of $n$ sites partitions the plane into one cell per site, the cell of site $s$ being every point closer to $s$ than to any other site. Its edges are the perpendicular bisectors between neighboring sites, and a nearest-site query becomes point location in this diagram. The Delaunay triangulation is its dual: connect two sites whenever their Voronoi cells share an edge. Delaunay triangulations have a defining property — the circumcircle of every triangle is empty of other sites — that makes them maximize the minimum angle, so they avoid the sliver triangles that wreck interpolation and mesh quality.⁴ Both are built in $O(n\log n)$, by Fortune's sweep from the previous lesson or by randomized incremental insertion, and either one yields the closest pair as a by-product: the two nearest sites are always joined by a Delaunay edge.

When the data is high-dimensional or the queries are approximate, Voronoi diagrams become impractical (their size grows as $O(n^{\lceil d/2\rceil })$), and the tool of choice is a spatial tree. A k-d tree recursively splits the points by alternating coordinate axes, giving $O(\log n)$ expected nearest-neighbor queries in low dimensions; R-trees group nearby objects into bounding rectangles for disk-resident spatial databases; and for the high-dimensional nearest-neighbor search behind modern embeddings, exact methods lose to locality-sensitive hashing and graph-based approximate indexes such as HNSW.⁵ All of these structures rest on the same primitive this module started with: every one decides which side and closer to which with orientation and distance comparisons, and keeps them exact where correctness of the combinatorics depends on it.

Takeaways

• Closest pair is divide-and-conquer: split at the median $x$, recurse, then scan a width-$2\delta$ strip. The strip-packing lemma caps the combine at $7$ comparisons per point, giving $T(n)=2T(n/2)+\Theta (n)=\Theta (n\log n)$ by the Master Theorem. Compare squared distances to stay exact.
• Point-in-polygon by ray casting: an odd count of $+x$-ray crossings means inside (Jordan curve), with half-open-edge tie-breaking for vertex/horizontal grazes and a separate on-segment test for the boundary.
• The winding number sums signed crossings ($+1$ up, $-1$ down); $w(q)\ne 0$ means inside and stays correct for self-intersecting polygons, where parity fails.
• The shoelace formula $\frac{1}{2}\left|{\sum }_i{V}_i\times V_{i+1}\right|$ is proved by the trapezoid sum: top edges add, bottom edges subtract, telescoping to the enclosed area; its sign gives the winding direction.
• Rotating calipers read the diameter (farthest pair, an antipodal hull pair) and the width (narrowest slab, flush with a hull edge) in $O(h)$ by advancing one monotone pointer around the hull — $O(n\log n)$ after the hull versus $\Theta (n^2)$ brute force.