Convex Hull Algorithms

Convex Hull Algorithms

The convex hull of a set of points is the smallest convex set that contains all the points. In two dimensions, this shape is easily visualized by imagining a rubber band stretched around the entire point set and allowed to snap tight. The vertices of the resulting polygon are a subset of the original points, and the region enclosed by these vertices forms the convex hull. Finding this boundary computationally is a fundamental problem in computational geometry, acting as a prerequisite for more complex spatial analyses, collision detection systems, and pattern recognition tasks.

The problem of computing the convex hull is mathematically well defined across any number of dimensions. Given a finite set of points in a Euclidean space, the convex hull is the intersection of all convex sets containing those points, or equivalently, the set of all convex combinations of the points. The internal structure of this space is governed by Carathéodory's theorem, which states that for any set of points in a $d$-dimensional space, any point within the convex hull can be expressed as a convex combination of at most $d+1$ points from the original set Carathéodory's theorem - Wikipedia^[1]. While the continuous geometric definition applies universally, algorithms designed to compute the hull focus on identifying the discrete subset of points that form the boundary vertices and the faces connecting them. The complexity of this task scales with both the number of input points and the dimensionality of the space.

Research into convex hull algorithms began in earnest during the early 1970s, driven by the need to optimize routing and spatial queries in early computer systems. Since then, the field has produced a variety of algorithms optimized for different constraints, ranging from simple two-dimensional planar methods to highly complex structures capable of handling multidimensional data. These algorithms are classified based on their time and space complexity, with modern optimal methods achieving bounds that depend not only on the size of the input but also on the size of the output.

Mathematical Definition and Problem Statement

In Euclidean geometry, a set is convex if, for every pair of points within the set, the straight line segment connecting them is also entirely within the set. The convex hull of a discrete set of points is formally defined as the set of all convex combinations of those points. A convex combination is a linear combination of the points where all coefficients are non-negative and sum to exactly one. Alternatively, the convex hull can be defined as the intersection of all half-spaces that contain the point set. This alternative definition highlights a fundamental concept in projective geometry known as point-line duality.^[31] Computing the convex hull of a set of points is mathematically dual to computing the intersection of a set of half-spaces, meaning an algorithm designed to solve one problem can be adapted to solve the other.

From an algorithmic perspective, the input is typically a finite array of $n$ points in a $d$-dimensional space. The required output is a representation of the boundary of the convex hull. In two dimensions, this output is an ordered list of vertices that form a convex polygon, usually sorted in counterclockwise order. In three dimensions, the output is a convex polyhedron, represented by its vertices, edges, and planar faces. In higher dimensions, the boundary consists of facets, which are $(d-1)$-dimensional simplices.

The computational complexity of finding the convex hull is bounded from below. In 1975, Michael Shamos demonstrated that computing the convex hull of a planar set of points requires $\Omega(n \log n)$ time in the worst case Geometric complexity (Shamos 1975)^[2]. This lower bound is established by reducing the problem of sorting $n$ real numbers to the problem of computing a convex hull. Given $n$ numbers, one can map them to $n$ points on a parabola defined by the equation $y = x^2$. Because all points on a parabola are strictly convex, every point will become a vertex of the convex hull. If an algorithm computes the convex hull of these points, it will return them in sequential order along the boundary, effectively sorting the original numbers. If an algorithm could compute the convex hull in less than $O(n \log n)$ time, it could also sort numbers in less than $O(n \log n)$ time, which violates the theoretical lower bound for comparison-based sorting. Therefore, any convex hull algorithm that operates in $O(n \log n)$ time is considered optimal in the worst case for a general point set.^[31]

Early History and Foundational Algorithms

Before the formalization of computational geometry as a distinct discipline, problems involving spatial boundaries were often solved using application-specific heuristics. The systematic study of convex hull algorithms began in 1972 when Ronald Graham, a mathematician at Bell Labs, published the first efficient algorithm for the planar convex hull problem ^[3].^[32] Graham was motivated by a practical problem involving the routing of a robotic arm, which required computing the shortest path around a set of obstacles. His method, now known as the Graham scan, achieved a time complexity of $O(n \log n)$ and established a baseline for future research.

The following year, in 1973, R. A. Jarvis at the Australian National University published a different approach ^[4]. Jarvis conceptualized the problem similarly to wrapping a piece of string around the points, a technique that became known as the gift wrapping algorithm or Jarvis march. While its worst-case time complexity was $O(n^2)$, its performance depended heavily on the number of points that actually formed the hull, making it highly efficient for specific types of datasets.

As the decade progressed, researchers sought methods that were faster in practice or easier to implement. In 1977, William Eddy published a recursive algorithm A new convex hull algorithm for planar sets (Eddy 1977)^[5], which was independently discovered and published by A. Bykat in 1978 Convex hull of a finite set of points in two dimensions (Bykat 1978)^[6]. This method, dubbed Quickhull due to its conceptual similarity to the Quicksort algorithm, partitioned the point set recursively. Also in 1977, Franco P. Preparata and Se June Hong introduced a divide and conquer approach capable of solving the problem in both two and three dimensions ^[7]. In 1978, Selim Akl and Godfried Toussaint published another fast heuristic approach that discarded interior points early in the computation to speed up the average case A fast convex hull algorithm (Akl and Toussaint 1978)^[8].

While general point sets require $O(n \log n)$ time, researchers discovered that specific input configurations allow for faster computation. In 1987, Avraham Melkman published an algorithm that computes the convex hull of a simple polyline or a simple polygon in strictly linear $O(n)$ time On-line construction of the convex hull of a simple polyline (Melkman 1987)^[9].^[33] Melkman's algorithm processes the vertices sequentially using a double-ended queue, avoiding the sorting step entirely because the geometric non-self-intersecting nature of the simple polygon provides an implicit ordering.

Mechanism of Action: Primary Two-Dimensional Algorithms

The foundational algorithms developed in the 1970s rely on different geometric principles to identify the boundary vertices. The three most prominent methods are the Graham scan, the Jarvis march, and Quickhull.

Graham Scan and Monotone Chain

The Graham scan operates by first identifying the point with the lowest y-coordinate, which is guaranteed to be on the convex hull. It then sorts the remaining points based on the polar angle they make with this anchor point. The algorithm iterates through the sorted points, maintaining a stack of potential hull vertices. For each new point, it checks whether the turn formed by the last two points in the stack and the new point is a left turn or a right turn using the two-dimensional cross product.

The cross product of three points $A(x_1, y_1)$, $B(x_2, y_2)$, and $C(x_3, y_3)$ is calculated as $(x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1)$. If the result is positive, the points form a counterclockwise turn. If negative, they form a clockwise turn. If zero, the points are collinear. If the turn is clockwise, the previous point cannot be on the convex hull and is popped from the stack. This popping process continues iteratively until the top two points on the stack and the new point form a counterclockwise left turn. The new point is then pushed onto the stack. The sorting step takes $O(n \log n)$ time, and the stack evaluation takes $O(n)$ time, resulting in an overall time complexity of $O(n \log n)$ Graham scan - Wikipedia^[10].

In 1979, A. M. Andrew published a variant of the Graham scan known as the Monotone Chain algorithm Another efficient algorithm for convex hulls in two dimensions (Andrew 1979)^[11]. Instead of sorting by polar angle, Andrew's algorithm sorts the points lexicographically by their x and y coordinates. It then constructs the upper and lower halves of the hull separately using the same cross product logic. This avoids the use of trigonometric functions required for polar angle calculations, making it faster and less prone to floating-point inaccuracies in practice.

function MonotoneChain(points):
    Sort points lexicographically (first by x, then by y)
    lower_hull = empty list
    for point in points:
        while length(lower_hull) >= 2 and cross_product(lower_hull[-2], lower_hull[-1], point) <= 0:
            lower_hull.pop()
        lower_hull.push(point)

    upper_hull = empty list
    for point in reversed(points):
        while length(upper_hull) >= 2 and cross_product(upper_hull[-2], upper_hull[-1], point) <= 0:
            upper_hull.pop()
        upper_hull.push(point)

    Remove the last point of each half because it is repeated
    return lower_hull + upper_hull

Jarvis March (Gift Wrapping)

The Jarvis march begins by finding the leftmost point in the dataset. From this starting vertex, the algorithm searches for the next vertex by evaluating all other points to find the one that makes the largest right hand turn, or equivalently, the smallest polar angle relative to the current vertex. To avoid expensive trigonometric operations, the algorithm uses the cross product to compare the relative angles of any two candidate points. This search requires iterating through all $n$ points. Once the next vertex is found, it becomes the new current vertex, and the process repeats until the algorithm wraps back to the starting point. If the final convex hull has $h$ vertices, the algorithm performs $n$ comparisons $h$ times, resulting in a time complexity of $O(nh)$ Gift wrapping algorithm - Wikipedia^[12]. This makes the Jarvis march highly efficient when $h$ is very small, but it degrades to $O(n^2)$ if all points lie on the boundary.

Quickhull

Quickhull uses a divide and conquer strategy. It begins by finding the points with the minimum and maximum x-coordinates, which must lie on the hull. The line connecting these two points divides the dataset into two subsets. For each subset, the algorithm finds the point furthest from the dividing line. This point forms a triangle with the line segment, and any points inside this triangle are discarded, as they cannot be on the hull. The algorithm then recursively applies this process to the two new outward-facing line segments formed by the triangle. The expected average time complexity is $O(n \log n)$, but in the worst case, such as when points are arranged on the perimeter of a circle, it degrades to $O(n^2)$ Quickhull - Wikipedia^[13].

Output-Sensitive Algorithms and Optimal Bounds

While the Graham scan guarantees $O(n \log n)$ performance, and the Jarvis march achieves $O(nh)$, researchers sought an algorithm that could guarantee optimal performance relative to both the input size $n$ and the output size $h$. Such an algorithm would need a worst-case time complexity of $O(n \log h)$.

In 1986, David Kirkpatrick and Raimund Seidel published an algorithm that achieved this theoretical lower bound The ultimate planar convex hull algorithm (Kirkpatrick and Seidel 1986)^[14]. Their approach introduced a "marriage before conquest" strategy. Traditional divide and conquer algorithms compute the hulls of the left and right halves of the data and then merge them. Kirkpatrick and Seidel reversed this logic. Their algorithm uses a linear time median finding algorithm to divide the points by a vertical line, then finds the "bridge" edge of the convex hull that crosses this line before computing the sub-hulls. Points beneath the bridge are discarded, and the algorithm recurses on the remaining points. While theoretically optimal, the Kirkpatrick Seidel algorithm is highly complex to implement and carries a large constant factor, making it less practical for standard applications Introduction to Algorithms (CLRS)^[15].

Ten years later, in 1996, Timothy Chan introduced a much simpler algorithm that also achieved the optimal $O(n \log h)$ time complexity Optimal output-sensitive convex hull algorithms in two and three dimensions (Chan 1996)^[16]. Chan's algorithm elegantly combines the Graham scan and the Jarvis march. It divides the $n$ points into blocks of size $m$, computes the convex hull of each block using the Graham scan in $O(n \log m)$ time, and then uses a modified Jarvis march to wrap around the sub-hulls. Instead of checking every single point, the modified Jarvis march treats each sub-hull as a single entity and finds the extreme tangent line to each sub-hull in $O(\log m)$ time using binary search. The wrapping phase takes $O(h (n/m) \log m)$ time.

Because the final number of vertices $h$ is unknown at the start, Chan's algorithm guesses the value of $m$. If the guess $m$ is too small, the algorithm aborts and tries again with a larger $m$. To maintain the $O(n \log h)$ time complexity, $m$ is not simply incremented. Instead, it is squared at each iteration, setting $m = 2^{2^t}$ for iteration $t$. The time spent on each iteration is proportional to $n \log m$. Because the sequence of guesses grows doubly exponentially, the sum of the time spent on all failed iterations is bounded by a geometric series that is strictly less than the time spent on the final, successful iteration. This ensures the total time remains $O(n \log h)$ Chan's algorithm - Wikipedia^[17].

Higher Dimensions and the Upper Bound Theorem

Computing the convex hull in three or more dimensions introduces significant mathematical and computational complexity. In two dimensions, the boundary is a sequence of edges. In three dimensions, the boundary is a polyhedron composed of vertices, edges, and planar faces, typically triangulated. In $d$ dimensions, the boundary consists of facets, which are $(d-1)$-dimensional simplices.

The complexity of the convex hull in higher dimensions is dictated by the Upper Bound Theorem, proven by Peter McMullen in 1970 The 1970 paper by Peter McMullen is titled "The maximum number of faces of a convex polytope"^[18]. The theorem states that the maximum number of faces of a $d$-dimensional convex polytope with $n$ vertices is proportional to $n^{\lfloor d/2 \rfloor}$. The worst-case scenario is represented by the cyclic polytope, formed by choosing points on the moment curve, defined parametrically as $(t, t^2, t^3, ..., t^d)$. For a cyclic polytope, every subset of $\lfloor d/2 \rfloor$ vertices forms a face of the convex hull, maximizing the number of faces. Consequently, any algorithm computing the full convex hull in $d$ dimensions must have a worst-case time complexity of at least $\Omega(n^{\lfloor d/2 \rfloor})$.

The divide and conquer algorithm proposed by Preparata and Hong in 1977 successfully computes the three-dimensional convex hull in $O(n \log n)$ time, which is optimal for three dimensions. However, extending algorithms to generalized $d$-dimensional space required further innovation. In 1993, Bernard Chazelle published an algorithm that computes the convex hull in any fixed dimension $d$ in $O(n \log n + n^{\lfloor d/2 \rfloor})$ time, achieving the theoretical optimum for all dimensions An optimal convex hull algorithm in any fixed dimension (Chazelle 1993)^[19].

For practical implementations in higher dimensions, the Quickhull algorithm was generalized by C. Bradford Barber, David P. Dobkin, and Hannu Huhdanpaa in 1996 The Quickhull algorithm for convex hulls (Barber, Dobkin, Huhdanpaa 1996)^[20]. Their implementation, known as Qhull, handles floating-point inaccuracies by merging coplanar facets, ensuring that the resulting hull remains strictly convex even with precision errors. Qhull operates efficiently in practice for dimensions up to eight or nine, though it succumbs to the exponential growth of facets in higher dimensions.

Real-World Applications

Convex hull algorithms are embedded in numerous software systems, acting as the mathematical foundation for spatial analysis and physical simulations.

Collision Detection and Physics Engines

In video games, robotics, and physics engines, determining whether two complex three-dimensional objects intersect is computationally expensive. To optimize this, objects are often enclosed in bounding volumes, such as convex hulls. The Gilbert-Johnson-Keerthi distance algorithm, published in 1988, is the standard method for computing the minimum distance between two convex shapes ^[21]. This algorithm relies heavily on the concept of the Minkowski difference. The Minkowski difference of two shapes is the set of all vectors formed by subtracting a point in the second shape from a point in the first shape Minkowski addition - Wikipedia^[22]. If the two original shapes intersect, their Minkowski difference will contain the origin. The algorithm uses support mappings to iteratively search for a simplex inside the Minkowski difference that encloses the origin. If the convex hulls of two objects do not intersect, the complex objects inside them cannot intersect, allowing the system to bypass detailed polygon-level collision checks.

Geographic Information Systems

Geographic Information Systems use convex hulls to estimate the territorial range of a dataset, such as the habitat of an animal species based on sighting coordinates or the coverage area of a network of cell towers. Spatial database extensions provide native functions to compute these boundaries. For example, the PostGIS extension for PostgreSQL includes the ST_ConvexHull function, which computes the smallest convex geometry that encloses a given spatial object or collection of objects PostGIS ST_ConvexHull^[23].

Computer Vision and Image Processing

In computer vision, convex hulls are used for shape analysis and feature extraction. When analyzing the silhouette of an object, such as a human hand, the convex hull forms a boundary around the fingertips. The spaces between the actual hand contour and the convex hull are known as convexity defects. By calculating the depth and position of these defects, algorithms can recognize specific hand gestures, count the number of extended fingers, and track hand movements in real time.

Machine Learning and Operations Research

In machine learning, the convex hull is a core component of Support Vector Machines for classification tasks. When separating two classes of data points, the algorithm essentially finds the shortest line segment connecting the convex hull of the first class to the convex hull of the second class. The hyperplane that bisects this line segment orthogonally is the optimal decision boundary that maximizes the margin between the two classes. If the convex hulls of the two classes intersect, the data is not linearly separable, and the algorithm must apply a kernel trick to map the points into a higher-dimensional space where the hulls no longer overlap. In operations research, the convex hull defines the feasible region of linear programming problems. The simplex algorithm navigates the vertices of the convex hull of the feasible region to find the optimal solution for resource allocation and facility location problems.

Software Implementations and Libraries

Due to the complexity of handling degenerate cases, such as collinear points or floating-point rounding errors, developers rarely write custom convex hull algorithms for production systems. Instead, they rely on established, heavily tested libraries.

The Qhull library, written in C, is the most widely used implementation for computing convex hulls, Delaunay triangulations, and Voronoi diagrams in multiple dimensions Qhull^[24]. It is the underlying engine for convex hull functions in major scientific computing environments. For instance, the SciPy library in Python exposes Qhull through the scipy.spatial.ConvexHull class scipy.spatial.ConvexHull^[25]. MATLAB and R also use Qhull for their native spatial functions.

In the domain of computational geometry, the Computational Geometry Algorithms Library provides rigorous C++ implementations of convex hull algorithms CGAL Convex Hull 2D^[26]. This library is designed to handle exact geometric computation. It separates the combinatorial algorithms from the geometric predicates, such as the orientation test. By using arbitrary-precision arithmetic for the predicates, the library guarantees that the combinatorial structure of the convex hull is always perfectly correct, eliminating the catastrophic topological failures that plague standard floating-point implementations Computational Geometry: Algorithms and Applications^[27].

For Java applications, the JTS Topology Suite provides a complete implementation of two-dimensional spatial algorithms, including convex hulls, which forms the basis for many open-source mapping applications JTS Topology Suite^[28].

Significance and Tradeoffs

The selection of a convex hull algorithm depends entirely on the constraints of the application, including the dimensionality of the data, the expected number of boundary points, and the tolerance for floating-point errors.

For general purpose two-dimensional applications, Andrew's Monotone Chain is the standard choice. Its $O(n \log n)$ time complexity is optimal for unsorted data, and its reliance on simple cross products makes it highly stable and easy to implement. When the output size $h$ is known to be extremely small compared to $n$, Chan's algorithm or the Jarvis march offers superior performance by scaling with the output size.

In three dimensions, the Preparata Hong algorithm or Quickhull are the standard approaches. Quickhull is particularly favored in practice because it extends cleanly to higher dimensions and handles degenerate data gracefully through facet merging.

However, convex hull algorithms face severe limitations in high-dimensional spaces due to the curse of dimensionality. Because the number of facets grows exponentially with the number of dimensions, computing the exact convex hull becomes computationally intractable for dimensions higher than ten. In such cases, exact algorithms are abandoned in favor of approximate convex hull methods or dimensionality reduction techniques.

Another tradeoff involves static versus dynamic data. The algorithms discussed above are designed for static point sets. If points are continuously added or removed, recomputing the entire hull from scratch is inefficient. Dynamic convex hull algorithms, such as those pioneered by Mark Overmars and Jan van Leeuwen in 1981, maintain the hull in a balanced search tree, allowing insertions and deletions in $O(\log^2 n)$ time Maintenance of configurations in the plane (Overmars and van Leeuwen 1981)^[29].

For objects in continuous motion, researchers developed kinetic data structures. Introduced by Julien Basch, Leonidas Guibas, and John Hershberger in 1997, kinetic data structures maintain the convex hull of moving points by scheduling events when the orientation of critical triangles changes Data structures for mobile data (Basch, Guibas, Hershberger 1997). These dynamic and kinetic structures are essential for real-time physics simulations where objects frequently deform or move, but they require significantly more memory and complex pointer management than their static counterparts Convex hull - Wikipedia^[30].

The study of convex hull algorithms encapsulates the broader evolution of computer science. It began with practical, heuristic-driven solutions in the 1970s, matured through mathematical proofs of optimality in the 1980s and 1990s, and stabilized into reliable, open-source libraries that now power modern spatial computing.