7 min read
On this page

Geometric Primitives

Points, Lines, Segments, and Polygons

Point: (x,y)R2(x, y) \in \mathbb{R}^2. Represented as a pair of coordinates (integer or floating-point).

Line: defined by equation ax+by+c=0ax + by + c = 0, or parametrically as p+tdp + t \cdot d for direction dd and tRt \in \mathbb{R}. Two points determine a unique line.

Line segment: the portion of a line between two endpoints p1,p2p_1, p_2: {p1+t(p2p1):t[0,1]}\{p_1 + t(p_2 - p_1) : t \in [0,1]\}.

Ray: {p+td:t0}\{p + t \cdot d : t \geq 0\}.

Polygon: ordered sequence of vertices (v0,v1,,vn1)(v_0, v_1, \ldots, v_{n-1}) defining a closed region via edges vivi+1modnv_i v_{i+1 \mod n}.

  • Simple polygon: edges intersect only at shared vertices (Jordan curve theorem applies: divides plane into interior and exterior)
  • Convex polygon: all interior angles <180°< 180°; equivalently, any line segment between two interior points lies entirely inside
  • Monotone polygon: any vertical line intersects the boundary at most twice

Orientation Test (Cross Product)

The cross product of vectors u=(ux,uy)\vec{u} = (u_x, u_y) and v=(vx,vy)\vec{v} = (v_x, v_y):

u×v=uxvyuyvx\vec{u} \times \vec{v} = u_x v_y - u_y v_x

This equals the signed area of the parallelogram spanned by u\vec{u} and v\vec{v}.

Orientation of three points p,q,rp, q, r: compute the cross product (qp)×(rp)(q - p) \times (r - p):

orient(p,q,r)=(qxpx)(rypy)(qypy)(rxpx)\text{orient}(p, q, r) = (q_x - p_x)(r_y - p_y) - (q_y - p_y)(r_x - p_x)

  • >0> 0: counterclockwise (left turn)
  • =0= 0: collinear
  • <0< 0: clockwise (right turn)

This is the fundamental predicate of computational geometry. Nearly every algorithm reduces to orientation tests. Equivalent to computing the sign of a 3×33 \times 3 determinant:

orient(p,q,r)=pxpy1qxqy1rxry1\text{orient}(p, q, r) = \begin{vmatrix} p_x & p_y & 1 \\ q_x & q_y & 1 \\ r_x & r_y & 1 \end{vmatrix}

Segment Intersection

Two segments p1p2\overline{p_1 p_2} and p3p4\overline{p_3 p_4} intersect iff:

  1. General position: orient(p1,p2,p3)\text{orient}(p_1, p_2, p_3) and orient(p1,p2,p4)\text{orient}(p_1, p_2, p_4) have different signs, AND orient(p3,p4,p1)\text{orient}(p_3, p_4, p_1) and orient(p3,p4,p2)\text{orient}(p_3, p_4, p_2) have different signs
  2. Collinear case: if any orientation is zero, check if the point lies on the segment (bounding box test)

For the intersection point (if it exists), solve the parametric system:

p1+t(p2p1)=p3+s(p4p3),t,s[0,1]p_1 + t(p_2 - p_1) = p_3 + s(p_4 - p_3), \quad t, s \in [0, 1]

This yields parameters via Cramer's rule in O(1)O(1) time.

Proper intersection: segments cross transversally (not at endpoints). Improper: endpoint touches the other segment or segments overlap.

Area Computation (Shoelace Formula)

For a simple polygon with vertices (x0,y0),(x1,y1),,(xn1,yn1)(x_0, y_0), (x_1, y_1), \ldots, (x_{n-1}, y_{n-1}) in order:

A=12i=0n1(xiyi+1xi+1yi)A = \frac{1}{2}\left|\sum_{i=0}^{n-1}(x_i y_{i+1} - x_{i+1} y_i)\right|

where indices are taken modulo nn. The signed version (without absolute value) is positive for counterclockwise vertex ordering.

Derivation: sum of signed areas of triangles formed with the origin. Each edge contributes a trapezoid. The formula is exact for integer/rational coordinates.

Centroid of a simple polygon:

xˉ=16Ai=0n1(xi+xi+1)(xiyi+1xi+1yi)\bar{x} = \frac{1}{6A}\sum_{i=0}^{n-1}(x_i + x_{i+1})(x_i y_{i+1} - x_{i+1} y_i)

and analogously for yˉ\bar{y}.

Point-in-Polygon

Ray Casting

Cast a horizontal ray from the query point to ++\infty. Count the number of edge crossings:

  • Odd crossings: inside
  • Even crossings: outside

Edge cases: ray passes through a vertex or is collinear with an edge. Handle by treating the ray as passing slightly above: count an edge crossing only if one endpoint is strictly above and the other is at or below the ray's yy-coordinate.

Complexity: O(n)O(n) per query for an nn-vertex polygon. Can be preprocessed for O(logn)O(\log n) queries via trapezoidal decomposition.

Winding Number

The winding number counts how many times the polygon boundary winds around the query point:

w=12πi=0n1angle(q,vi,vi+1)w = \frac{1}{2\pi}\sum_{i=0}^{n-1} \text{angle}(q, v_i, v_{i+1})

  • w=0w = 0: outside
  • w0w \neq 0: inside (for simple polygons, w=±1w = \pm 1)

Efficient computation: avoid trigonometry by tracking sign changes in the yy-coordinate relative to the query point and using orientation tests. Same O(n)O(n) complexity as ray casting.

Advantage over ray casting: correctly handles self-intersecting polygons (where "inside" is ambiguous for ray casting but well-defined by winding number convention).

Convexity Testing

A polygon is convex iff all cross products of consecutive edge vectors have the same sign:

sign((vi+1vi)×(vi+2vi+1)) is constant for all i\text{sign}((v_{i+1} - v_i) \times (v_{i+2} - v_{i+1})) \text{ is constant for all } i

O(n)O(n) time.

Point-in-convex-polygon: O(logn)O(\log n) by binary search on the angle from a fixed interior point (e.g., the first vertex), then a single orientation test.

Distance Computations

Point to line (line through a,ba, b):

d(p,)=(ba)×(pa)bad(p, \ell) = \frac{|(b - a) \times (p - a)|}{|b - a|}

Point to segment ab\overline{ab}: project pp onto the line; if the projection parameter t=(pa)(ba)ba2t = \frac{(p-a) \cdot (b-a)}{|b-a|^2} is in [0,1][0,1], use the point-to-line distance. Otherwise, use the distance to the nearer endpoint.

Segment to segment: check if projections overlap; otherwise, minimum of four point-to-segment distances. Important for collision detection.

Polygon Triangulation

Every simple polygon can be triangulated into n2n - 2 triangles (for nn vertices).

Ear-clipping: an "ear" is a triangle formed by three consecutive vertices whose diagonal lies inside the polygon. Every simple polygon with 4\geq 4 vertices has at least two ears (two-ears theorem). Repeatedly clip ears: naive O(n2)O(n^2); optimized O(n)O(n) via careful bookkeeping (though the O(n)O(n) algorithm by Chazelle is complex).

Monotone polygon triangulation: O(n)O(n) after O(nlogn)O(n \log n) monotone decomposition (sweep-line). This is the standard practical approach.

Fan triangulation for convex polygons: O(n)O(n), all triangles share a common vertex.

Convex Polygon Operations

Intersection of two convex polygons: O(n+m)O(n + m) via rotating calipers or simultaneous traversal (O'Rourke's algorithm).

Minkowski sum of two convex polygons: O(n+m)O(n + m) by merging sorted edge vectors.

Tangent lines, supporting lines: O(logn)O(\log n) via binary search on the polygon.

Exact Arithmetic and Robustness

The Robustness Problem

Floating-point orientation tests can give incorrect signs due to rounding, leading to:

  • Topological inconsistencies (a point reported inside and outside)
  • Infinite loops in sweep-line algorithms
  • Invalid geometric constructions (self-intersecting "convex" hulls)

A classic failure: three nearly collinear points with coordinates like (0,0),(1,ϵ),(2,2ϵ+δ)(0, 0), (1, \epsilon), (2, 2\epsilon + \delta) where δ\delta is below machine precision.

Shewchuk's Exact Predicates

Jonathan Shewchuk's adaptive precision arithmetic (1997) provides exact geometric predicates using only floating-point hardware:

Key idea: represent intermediate results as non-overlapping expansions (sums of floating-point numbers with decreasing magnitude). The orientation predicate becomes:

  1. Compute with standard double precision (fast filter)
  2. If the result is clearly positive/negative, return immediately
  3. If near zero, progressively refine using error-free transformations (Dekker/Knuth two-product, two-sum)
  4. Final exact result uses O(1)O(1) machine words for 2D orientation

Performance: the fast filter succeeds for >99% of non-degenerate cases, costing only ~2x standard floating-point. Exact path is 10-50x slower but rarely needed.

Predicates provided: orient2d, orient3d, incircle, insphere. These suffice for Delaunay triangulation, convex hull, and most fundamental algorithms.

Alternative Approaches

  • Exact integer arithmetic: use integer coordinates and multiprecision integers for determinants. GMP library. Exact but slow
  • Interval arithmetic: compute with intervals; if the interval contains zero, refine. CGAL uses this approach
  • Simulation of simplicity (SoS): symbolic perturbation to eliminate degeneracies. Edelsbrunner and Mucke (1990). Perturb input by infinitesimal amounts ϵi\epsilon^i and evaluate the sign using lexicographic rules
  • Controlled perturbation: actually perturb coordinates by small random amounts, then verify geometric consistency

CGAL's Exact Computation Paradigm

The Exact Geometric Computation (EGC) paradigm: all geometric decisions (predicates) are computed exactly; only numerical constructions (coordinates of intersection points) may be approximate. Since algorithms branch only on predicates, EGC guarantees topological correctness.

CGAL implements this via:

  • Lazy exact arithmetic: compute with doubles, fall back to exact (MPFR/GMP) only when needed
  • Algebraic number types for constructions involving square roots (circular arcs)
  • Separation bounds to determine when floating-point precision suffices

Degeneracies

Common degeneracies: three collinear points, four cocircular points, coincident points, vertical segments in sweep-line algorithms.

Handling strategies:

  1. General position assumption: state and prove algorithms assuming non-degeneracy, then handle special cases
  2. Symbolic perturbation (SoS): automatically resolves all degeneracies consistently
  3. Explicit case analysis: handle each degenerate case in the code (error-prone but sometimes necessary)

Correct handling of degeneracies is often where the majority of implementation complexity lies. It is the primary reason computational geometry libraries (CGAL, S2 Geometry) are strongly preferred over hand-rolled implementations for production use.