Optimal Output-Sensitive Convex Hull Algorithms in Two and Three Dimensions
Timothy M. Chan* Department of Computer Science
University of British Columbia Vancouver, B.C. V6T 1Z4, Canada
Abstract
We present simple output-sensitive algorithms that construct the convex hull of a set of n points in two or three dimensions in worst-case optimal O(nlogh) time and O(n) space, where h denotes the number of vertices of the convex hull.
1 Introduction
Given a set P of n points in the Euclidean plane E2 or Euclidean space E3, we consider the problem of computing the convex hull of P, conv(P), which is defined as the smallest convex set containing P. The convex hull problem has received considerable attention in computational geometry [ff, 2f, 23, 25]. In E2, an algorithm known as Graham's scan [15] achieves O(ralogra) running time, and in E3, an algorithm by Preparata and Hong [24] has the same complexity. These algorithms are optimal in the worst case, but if h, the number of hull vertices, is small, then it is possible to obtain better time bounds. For example, in E2, a simple algorithm called Jarvis's march [19] can construct the convex hull in O(nh) time. This bound was later improved to O(nlogh) by an algorithm due to Kirkpatrick and Seidel [20], who also provided a matching lower bound; a simplification of their algorithm has been recently reported by Chan, Snoeyink, and Yap [2]. In E3, one can obtain an 0(n/i)-time algorithm using the gift-wrapping method, an extension of Jarvis's march originated by Chand and Kapur [3]. A faster but more involved algorithm in E3 was discovered by Edelsbrunner and Shi [13], having running time 0(ralog2 h). Finally, by derandomizing an algorithm of Clarkson and Shor [8], Chazelle and Matousek [7] succeeded in attaining optimal 0(n\ogh) time in E3. These algorithms, with complexity measured as a function of both n and the "output size" h, are said to be output-sensitive.
In this note, we point out a simple output-sensitive convex hull algorithm in E2 and its extension in E3, both running in optimal O(nlogh) time. Previous optimal (deterministic) methods, including the algorithm by Kirkpatrick and Seidel and its improvement by Chan, Snoeyink, and Yap, all rely on the existence of a linear-time procedure for finding medians. In Chazelle and Matousek's three-dimensional algorithm, even more complex tools for derandomization, such as e-approximations, are used.   Our algorithms avoid median-finding and derandomization altogether; Dobkin-Kirkpatrick
* Supported by a Killam Predoctoral Fellowship and an NSERC Postgraduate Scholarship.
1
Figure 1: Wrapping a set of \n/m] convex polygons of size to.
hierarchies [9, 10] are the only data structures used in the three-dimensional case. Our idea is to speed up Jarvis's march and the gift-wrapping method by using a very simple grouping trick.
2   An Output-Sensitive Algorithm in Two Dimensions
Let P C E2 be a set of n > 3 points. For simplicity, we assume that the points of P are in general position, i.e., no three points are collinear; see Section 4 for how to deal with degenerate point sets.
Recall that Jarvis's march [19, 23, 25] computes the h vertices of the convex hull one at a time, in counterclockwise (ccw) order, by a sequence of h wrapping steps: if Pk-i and pk are the previous two vertices computed, then the next vertex Pk+i is set to be the point p £ P that maximizes the angle Lpk-iPkP with p / p^. One wrapping step can obviously be done in O(n) time by scanning all n points; with an appropriate initialization the method constructs the entire convex hull in O(nh) time.
We observe that a wrapping step can be done faster if we preprocess the points. Choose a parameter to between 1 and n and partition P into \n/m] groups each of size at most to. Compute the convex hull of each group in 0(mlogto) time by, say, Graham's scan [15]. This gives us [ra/m] possibly overlapping convex polygons each with at most to vertices, after a preprocessing time of 0(^(m log to)) = O(ralogTO). Now, a wrapping step can be done by scanning all \n/m] polygons and computing tangents or supporting lines of the polygons through the current vertex p^, as shown in Figure 1. Since tangent finding takes logarithmic time for a convex polygon by binary or Fibonacci search [5, 25] (the dual problem is to intersect a convex polygon with a ray), the time required for a wrapping step is then O(^MogTO). As h wrapping steps are needed to compute the hull, the total time of the algorithm becomes 0(n log to +       log to)) = 0(n(l + h/m) log to).
The following is a pseudocode of the algorithm just described. The procedure always runs within 0(n(l + H/m) log to) time and successfully returns the list of vertices of conv(P) in ccw order when
2
H>h.
Algorithm Hull2D(P, to, H), where P C E2, 3 < to < n, and H > 1
1. partition P into subsets Pi,..., P\njm-\ each of size at most to
2. for i = 1,..., [n/m] do
3. compute conv(P8) by Graham's scan and store its vertices in an array in ccw order
4. p0 <- (0, -oo)
5. pi <— the rightmost point of P
6. for k = 1,..., H do
7. for i = 1,..., \n/m] do
8. compute the point g8- £ P4- that maximizes Lp^-xpyn      / p^) by performing a binary search on the vertices of conv(P8)
9. Pk+i <— the point q from {qi,.. .,^|-„/m]} that maximizes Apk-ip^q
10. if Pk+i = Pi then return the list (pi,.. .,pk)
11. return incomplete
By choosing m = H, the complexity of the algorithm is then 0(ra(l + H/m) log to) = O(nlogH). Since the value of h is not known in advance, we use a sequence of if's to "guess" its value as shown below (the same strategy is used in Chazelle and Matousek's algorithm):
Algorithm Hull2D(P), where P C E2
1. for t= 1,2,... do
2. L*r- Hull2D(P, to, H), where to = H = min{22', n}
3. if L / incomplete then return L
The procedure stops with the list of hull vertices as soon as the value of H in the for-loop reaches or exceeds h. The number of iterations in the loop is [loglog h] (using base-2 logarithms), and the t-th iteration takes O(nlogH) = 0(n2t) time. Therefore, the total running time of the algorithm is 0(Et=°iglOg/l1 n2t) = O(ra2rios1°s/ll+1) = O(nlogh). The storage requirement is clearly linear.
3   An Output-Sensitive Algorithm in Three Dimensions
Let P C E3 be a set of n > 4 points. Again we assume general position, i.e., no four points are coplanar (see Section 4). It suffices to construct the 2h — 4 facets (triangular faces) of the convex hull; with the aid of a dictionary, we can easily produce the set of h vertices and 3h — 6 edges together with their adjacency and order information in additional O(hlogh) time.
The higher-dimensional analogue of Jarvis's march is Chand and Kapur's gift-wrapping method [3, 25, 26], which computes the hull facets one at a time as follows: from a given facet /, we generate its three adjacent facets fj by performing a wrapping step about each of the three edges ej of / (j = 1, 2, 3). Here, a wrapping step about ej is to compute a point pj £ P that maximizes the angle between / and conv(ej U{pj}) with pj (j£ ej. Since such a step can be done in O(n) time, we can find the facets adjacent to / in O(n) time. Assuming an initial facet fo is given (which can be found in two wrapping steps), a breadth-first or depth-first search can then generate all facets of the convex
3
hull. Using a dictionary to detect duplication, we can ensure that each facet is processed once. This implies that the algorithm performs 3(2h — 4) wrapping steps and thus runs in O(nh) time.
We can use the same grouping idea from the previous section to improve the time complexity to optimal O(nlogh) while maintaining linear space. The calls to Graham's scan (line 3 of Hull2D(P, to, H)) are now replaced by calls to Preparata and Hong's three-dimensional convex hull algorithm [24], which has the same complexity. To make line 8 work in E3, we need to calculate tangents or supporting planes of convex polyhedra through a given line (or, in dual space, intersect convex polyhedra with a ray). If we use the hierarchical representation of Dobkin and Kirkpatrick [9, 10] to store these polyhedra (which requires only linear-time preprocessing), then the tangents can be computed in logarithmic time each, as before. The analysis is thus identical to that of the two-dimensional algorithm. The pseudocode is as follows:
Algorithm Hull3D(P, to, H), where P C E3, 4 < to < n, and H > 1
1. partition P into subsets Pi,..., P\n/m-\ each of size at most to
2. for i = 1,..., \n/m] do
3. compute conv(P8) by Preparata and Hong's algorithm and store it in a Dobkin-Kirkpatrick hierarchy
4. F, Q <— {/o}j where fo is some initial facet of conv(P)
5. for k = 1,.. .,2H - 4 do
6. if Q = 0 then return F
7. pick some / £ Q and set Q <— Q — {/}
8. let e,j be the edges of / (j = 1, 2, 3)
9. for j = 1,2,3 do
10. for i = 1,..., \n/m] do
11. compute the point g8- £ Pi that maximizes the angle between / and conv(ej U {qi}) by searching the hierarchy of conv(P8)
12. pj <— the point q from {qi,..., q^n/m-j} that maximizes the angle between
/ and conv(ej U {q}) (q (j£ ej)
13. /_,• <- conv(ej U {p3})
14. if f3 £ F then
15. F^-FU {/,•}, Q^QU{fj}
16. return incomplete
We can use a queue or a stack to implement Q and a dictionary to implement F. As there are only 0(h) dictionary operations, they can be carried out in O(hlogh) time. In fact, more clever implementations of the gift-wrapping method via a shelling order replace the need for dictionaries with just a priority queue.
As before, we choose the group size to = H and guess the value of h with a sequence of if's:
Algorithm Hull3D(P), where P C E3
1. for t= 1,2,... do
2. L*r- Hull3D(P, to, H), where to = H = min{22', n}
3. if L / incomplete then return L
4
4 Refinements
In this section, we suggest ideas on possible improvements that may speed up our algorithms in practice; we also discuss how degenerate cases can be handled.
Idea 1. First, points found to be in the interior of conv(P8) in line 3 of Hull2D(P, to, H) or Hull3D(P, to, H) can be eliminated from further consideration. This may potentially save work during future iterations of the algorithm, although it does not affect the worst-case complexity.
Idea 2. In Hull2D(P) and Hull3D(P), we choose the group size to = H so as to balance the O(ralogTO) preprocessing cost and the 0(H(^ log to)) cost for the 0(H) wrapping steps. Alternatively, we can choose to = m'm{H log if, n} (or set H = to/log to). This choice of to does not affect the former cost except in the lower-order terms, but it reduces the latter cost from O(nlogH) to O(n) and thus results in a smaller constant factor overall.
Idea 3. With Idea 2, the dominant cost of algorithm Hull2D(P, to, H) lies in the preprocessing, i.e., the computation of the convex hulls of the groups in line 3. To reduce this cost, we may consider reusing hulls computed from the previous iteration and merging them as the group size is increased. Suppose to' is the previous group size. Since the convex hull of two convex polygons can be computed in linear time (the dual problem is to intersect two convex polygons), we can compute the convex hull of [to/to'] convex m'-gons in 0(m log (to/to') ) time by the standard "mergehull" divide-and-conquer algorithm [25]. Thus, the \n/m] hulls in line 3 can be constructed in 0(n log (to/to') ) rather than O(ralogTO) time. The same can be said for the three-dimensional case, but merging two convex polyhedra, though possible in linear time [4], is more complicated.
Idea 4. In Hull2D(P), we use the sequence of group sizes to = 22', t = 1,2,..., to guess h. The improvements from Ideas 2 and 3 in fact permit us to choose slower growing sequences and still retain optimal O(nlogh) complexity. For example, one possible sequence is simply to = 2f, £ = 2, 3,..., which corresponds to doubling the group size after each iteration. Note that a coarser sequence approximates h less well while a denser sequence requires more iterations. We may try to optimize the worst-case constant factor and lower-order terms using sequences with different growth rates. We suggest the sequence to = 2f ,£ = 2,3,...
Idea 5. E. Welzl has observed that the binary search in line 8 of algorithm Hull2D(P, to, H) can be replaced by a simpler linear search without changing the time complexity of the algorithm. The following monotonicity property provides the justification: during the course of the algorithm, the variable g8- in line 8 can only advance in the ccw direction along conv(P8) for each fixed i. As a result, the h-vertex convex hull of p convex polygons with a total of n vertices can be computed in 0(n-\-hp) time by gift-wrapping; the two-polygon (p = 2) version of the algorithm is in fact the dual of an intersection algorithm by O'Rourke et al. [22] (see also [23, 25]). The total cost of Hull2D(P, to, H) can then be reduced to 0(n\ogm-\-H(n/m)) time, which is a logm factor saving in the second term. Although the overall constant factor is unaffected by the saving if Idea 2 is employed (as the first term is the dominant one), the linear search is easier to implement. There does not seem to be an analogous simplification in three dimensions.
5
Degeneracies. In both algorithms Hull2D(P, to, H) and Hull3D(P, to, H), we have assumed that the points of P are in general position. One way to cope with degenerate point sets is to apply general perturbation methods such as [12, 14]; however, these methods may cause the output size h to increase, as a point that is not a hull vertex but lies on the hull boundary may become a vertex after perturbation. Thus, it is better to handle the degenerate cases directly. For algorithm Hull2D(P, to, H), this is not difficult to do: when there is more than one point q that maximizes the angle Lpk-iPkq in line 9, pick the point q that is farthest from p^; use the same rule to break ties in line 8.
For algorithm Hull3D(P, to, H), we can do the following: In line 8, let ej = cijbj with cij and bj oriented in ccw order around /. When there is more than one point q that maximizes the angle between / and conv(ej U {q}) in line 12, pick the point q that maximizes the angle /.bjCijq; and if there is still more than one q that achieves the maximum, pick the one farthest from cij. Use the same rule to break ties in line 11. For degenerate point set, it is easier to keep track of edges rather than facets, since facets can be convex polygons rather than triangles. So, make F and Q sets of
--► __y
edges instead, and in line 15, add the oriented edges bjCij and ajij to F and Q. Although we may not have a complete description of the facet incident to these two edges, we know the equation of the plane containing the facet; this equation is sufficient to perform wrapping about these edges.
5 Extensions
We have presented new optimal output-sensitive convex hull algorithms in E2 and E3. The algorithms are simpler than previous O(nlogh) algorithms, particularly in the three-dimensional case, and the constant factors behind the big-Oh are likely to be smaller than those of the previous algorithms (in the worst case).
Besides its simplicity, our approach has the advantage that it is applicable to a variety of other problems. As an illustration, consider the problem of computing the lower envelope C(S) of a set S of n line segments in the plane, which we define as the boundary of Uses * where s denotes the unbounded trapezoid conv(s U {(0, +oo)}) for a given segment s. (Convex hulls correspond to lower envelopes of lines in the dual.) Let h be the output size, i.e., the number of edges in the envelope; it is known that h is at most 0(na(n)) [16]. Hershberger [17] has given a worst-case optimal algorithm that computes lower envelopes in O(ralogra) time. We now describe how his algorithm can be made output-sensitive with our technique.
First, observe that we can trace the h edges in C(S) from left to right by performing h ray shooting operations, where a ray shooting operation is: given a ray p emanating from a point on or beneath C(S), find the first trapezoid s (s £ S) that p crosses. As such an operation can be done in O(n) time, this gives us a naive O(nh) method, like Jarvis's march. To improve the running time, partition S into \n/m] groups each of at most to segments and compute the lower envelope of each group by Hershberger's algorithm; this takes O(ralogTO) time in total. Using known data structures such as [6, 18], we can perform ray shooting under each of these \n/m] envelopes in O(logra) time after 0(ma(m)) preprocessing (the ray shooting methods can be simplified in our case since envelopes are monotone). This implies that the h ray shooting operations on C(S) can be done in 0(h(^logm)) time. Choosing an appropriate group size to and guessing the output size h give us an optimal output-sensitive O(nlogh) algorithm for computing the lower envelope.
6
Other applications of our technique can be found in [1], including the output-sensitive construction of higher-dimensional convex hulls and ^-levels. In many cases, our grouping idea, combined with appropriate data structures, can be used to obtain optimal O(nlogh) algorithms if the output size h is sufficiently small, i.e., if h = o(na) for a suitable constant a.
Acknowledgements
I wish to thank Jack Snoeyink for his great support and encouragement as well as for his many valuable suggestions.
References
[1] T. M. Chan. Output-sensitive results on convex hulls, extreme points, and related problems. Discrete Comput. Geom., this issue.
[2] T. M. Chan, J. Snoeyink, and C.-K. Yap. Output-sensitive construction of polytopes in four dimensions and clipped Voronoi diagrams in three. In Proc. 6th ACM-SIAM Sympos. on Discrete Algorithms, 282291, 1995.
[3] D. R. Chand and S. S. Kapur. An algorithm for convex polytopes. J. ACM, 17:78-86, 1970.
[4] B. Chazelle. An optimal algorithm for intersecting three-dimensional convex polyhedra. SIAM J. Corn-put., 21:671-696, 1992.
[5] B. Chazelle and D. P. Dobkin. Intersection of convex objects in two and three dimensions. J. ACM, 34:1-27, 1987.
[6] B. Chazelle, H. Edelsbrunner, M. Grigni, L. Guibas, J. Hershberger, M. Sharir, and J. Snoeyink. Ray shooting in polygons using geodesic triangulations. In Proc. 18th Int. Colloq. Automata, Languages, and Programming, Lect. Notes in Comput. Sei., vol. 510, 661-673, 1991.
[7] B. Chazelle and J. Matoušek. Derandomizing an output-sensitive convex hull algorithm in three dimensions. Comput. Geom. Theory AppL, 5:27-32, 1995.
[8] K. L. Clarkson and P. W. Shor. Applications of random sampling in computational geometry, II. Discrete Comput. Geom., 4:387-421, 1989.
[9] D. P. Dobkin and D. G. Kirkpatrick. Fast detection of polyhedral intersection. Theoretical Comput. Sei., 27:241-253, 1983.
[10] D. P. Dobkin and D. G. Kirkpatrick. Determining the separation of preprocessed polyhedra: a unified approach. In Proc. 17th Int. Colloq. Automata, Languages, and Programming, Lect. Notes in Comput. Sei., vol. 443, 440-413, 1990.
[11] H. Edelsbrunner. Algorithms in Combinatorial Geometry. Springer-Verlag, Berlin, 1987.
[12] H. Edelsbrunner and E. P. Mücke. Simulation of simplicity: A technique to cope with degenerate cases in geometric algorithms. ACM Trans. Graphics, 9:66-104, 1990.
[13] H. Edelsbrunner and W. Shi. An 0(nlog2 h) time algorithm for the three-dimensional convex hull problem. SIAM J. Comput, 20:259-277, 1991.
[14] I. Emiris and J. Canny. An efficient approach to removing geometric degeneracies. In Proc. 8th ACM Sympos. Comput. Geom., 74-82, 1992.
7
[15] R. L. Graham. An efficient algorithm for determining the convex hull of a finite planar set. Inform. Process. Lett., 1:132-133, 1972.
[16] S. Hart and M. Sharir. Nonlinearity of Davenport-Schinzel sequences and of generalized path compression schemes. Combinatorial, 6:151-177, 1986.
[17] J. Hershberger. Finding the upper envelope of n line segments in O(nlogn) time. Inform. Process. Lett., 33:169-174, 1989.
[18] J. Hershberger and S. Suri. A pedestrian approach to ray shooting: shoot a ray, take a walk. In Proc. 4-th ACM-SIAM Sympos. on Discrete Algorithms, 54-63, 1993.
[19] R. A. Jarvis. On the identification of the convex hull of a finite set of points in the plane. Inform. Process. Lett., 2:18-21, 1973.
[20] D. G. Kirkpatrick and R. Seidel.   The ultimate planar convex hull algorithm?    SIAM J. Comput., 15:287-299, 1986.
[21] K. Mulmuley.   Computational Geometry: An Introduction Through Randomized Algorithms. Prentice Hall, New York, 1993.
[22] J. O'Rourke, C.-B. Chien, T. Olson, and D. Naddor. A new linear algorithm for intersecting convex polygons. Comput. Graph. Image Process., 19:384-391, 1982.
[23] J. O'Rourke. Computational Geometry in C. Cambridge University Press, 1994.
[24] F. P. Preparata and S. J. Hong.   Convex hulls of finite sets of points in two and three dimensions. Commun. ACM, 20:87-93, 1977.
[25] F. P. Preparata and M. I. Shamos. Computational Geometry: An Introduction. Springer-Verlag, New York, 1985.
[26] G. F. Swart. Finding the convex hull facet by facet. J. Algorithms, 6:17-48, 1985.
8