Figure 10: Computing the lower tangent.
The asymptotic running time of the algorithm can be expressed by a
recurrence. Given an input of size n, consider the time needed to
perform all the parts of the procedure, ignoring the recursive
calls. This includes the time to partition the point set, compute the
two tangents, and return the final result. Clearly the first and third
of these steps can be performed in O(n) time, assuming a linked list
representation of the hull vertices. Below we will show that the
tangents can be computed in O(n) time. Thus, ignoring constant
factors, we can describe the running time by the following recurrence.
T(n) =
1 if n <= 3
n + 2T (n/2) otherwise.
This is the same recurrence that arises in Mergesort. It is easy to
show that it solves to T (n) 2 O(n log n) (see CLR).
All that remains is showing how to compute the two tangents. One
thing that simplifies the process of computing the tangents is that
the two point sets A and B are separated from each other by a vertical
line (assuming no duplicate x coordinates). Let's concentrate on the
lower tangent, since the upper tangent is symmetric. The algorithm
operates by a simple ``walking'' procedure. We initialize a to be the
rightmost point of HA and b is the leftmost point of HB . (These can
be found in linear time.) Lower tangency is a condition that can be
tested locally by an orientation test of the two vertices and
neighboring vertices on the hull. (This is a simple exercise.) We
iterate the following two loops, which march a and b down, until they
reach the points lower tangency.
Finding the Lower Tangent
LowerTangent(HA ; HB ) :
(1) Let a be the rightmost point of HA .
(2) Let b be the leftmost point of HB .
(3) While ab is not a lower tangent for HA and HB do
(a) While ab is not a lower tangent to HA do a = a - 1 (move a clockwise).
(b) While ab is not a lower tangent to HB do b = b + 1 (move b counterclockwise).
(4) Return ab.
Proving the correctness of this procedure is a little tricky, but not
too bad. Check O'Rourke's book out for a careful proof. The important
thing is that each vertex on each hull can be visited at most once by
the search, and hence its running time is O(m), where
m = |HA| + |HB| <= |A| + |B|.
This is exactly what we needed to get the overall O(n log n) running time.
QuickHull: If the divide and conquer algorithm can be viewed as a sort
of generalization of Merge Sort, one might ask whether there is
corresponding generalization of other sorting algorithm for computing
convex hulls. In particular, the next algorithm that we will consider
can be thought of as a generalization of the QuickSort sorting
procedure. The resulting algorithm is called QuickHull. Like
QuickSort, this algorithm runs in O(n log n) time for favorable inputs
but can take as long as O(n^2) time for unfavorable inputs. However,
unlike QuickSort, there is no obvious way to convert it into a
randomized algorithm with O(n log n) expected running
time. Nonetheless, QuickHull tends to perform very well in practice.
The intuition is that in many applications most of the points lie in
the interior of the hull. For example, if the points are uniformly
distributed in a unit square, then it can be shown that the expected
number of points on the convex hull is O(log n). The idea behind
QuickHull is to discard points that are not on the hull as quickly as
pos sible. QuickHull begins by computing the points with the maximum
and minimum, x and y coordinates. Clearly these points must be on the
hull. Horizontal and vertical lines passing through these points are
support lines for the hull, and so define a bounding rectangle, within
which the hull is contained. Furthermore, the convex quadrilateral
defined by these four points lies within the convex hull, so the
points lying within this quadrilateral can be eliminated from further
consideration. All of this can be done in O(n) time.
Figure 11: QuickHull's initial quadrilateral. Points inside are discarded.
To continue the algorithm, we classify the remaining points into the
four corner triangles that remain. In general, as this algorithm
executes, we will have an inner convex polygon, and associated with
each edge we have a set of points that lie ``outside'' of that
edge. (More formally, these points are witnesses to the fact that this
edge is not on the convex hull, because they lie outside the half
plane defined by this edge.) When this set of points is empty, the
edge is a final edge of the hull. Consider some edge ab. Assume that
the points that lie ``outside'' of this hull edge have been placed in
a bucket that is associated with ab. Our job is to find a point c
among these points that lies on the hull, discard the points in the
triangle abc, and split the remaining points into two subsets, those
that lie outside ac and those than lie outside of cb. We can classify
each point by making two orientation tests.
Figure 12: QuickHull elimination procedure.
How should c be selected? There are a number of possible selection
criteria that one might think of. The method that is most often
proposed is to let c be the point that maximizes the perpendicular
distance from the line ab. (For example, another possible choice might
be the point that maximizes the angle cba or cab. It turns out that
these can be are very poor choices because they tend to produce
imbalanced partitions of the remaining points.) We replace the edge ab
with the two edges ac and cb, and classify the points as lying in one
of three groups: those that lie in the triangle abc, which are
discarded, those that lie outside of ac, and those that lie outside of
cb. We put these points in buckets for these edges, and recurse. (We
claim that it is not hard to classify each point p, by computing the
orientations of the triples acp and cbp.)
The running time of Quickhull, as with QuickSort, depends on how
evenly the points are split at each stage. Let T (n) denote the
running time on the algorithm assuming that n points remain outside of
some edge. In O(n) time we can select a candidate splitting point c
and classify the points in the bucket in O(n) time. Let n1 and n2
denote the number of remaining points, where n1 + n2 <= n. Then the
running time is given by the recurrence:
T(n) =
1 if n = 1
T(n1) + T (n2) where n1 + n2 <= n. otherwise
In order to solve this recurrence, it would be necessary to determine
the ``reasonable'' values for n1 and n2 . If we assume that the
points are ``evenly'' distributed, in the sense that
max(n1 ; n2 ) <= a * n for some constant a < 1,
then by applying the same analysis as that used in QuickSort (see
Cormen, Leiserson, Rivest) the running time will solve to O(n log n),
where the constant factor depends on a. On the other hand, if the
splits are not balanced, then the running time can easily increase to
O(n^2).
Does QuickHull outperform Graham's scan? This depends to a great
extent on the distribution of the point set. There are variations of
QuickHull that are designed for specific point distributions
(e.g. points uniformly distributed in a square) and their authors
claim that they manage to eliminate almost all of the points in a
matter of only a few iterations.
Gift Wrapping and Jarvis's March: The next algorithm that we will
consider is a variant on an O(n^2 ) sorting algorithm called
SelectionSort. For sorting, this algorithm repeatedly finds the next
element to add to the sorted order from the remaining items. The
corresponding convex hull algorithm is called Jarvis's march. which
builds the hull in O(nh) time by a process called ``gift
wrapping''. The algorithm operates by considering any one point that
is on the hull, say, the lowest point. We then find the ``next'' edge
on the hull in counterclockwise order. Assuming that p(k) and p(k-1)
were the last two points added to the hull, compute the point q that
maximizes the angle [ p(k-1) p(k) q ]. Thus, we can find the point q in
O(n) time. After repeating this h times, we will return back to the
starting point and we are done. Thus, the overall running time is
O(nh). Note that if h is o(log n) (asymptotically smaller than log n)
then this is a better method than Graham's algorithm.
One technical detail is how we to find an edge from which to
start. One easy way to do this is to let p(1) be the point with the
lowest y coordinate, and let p(0) be the point (Inf, 0), which is
infinitely far to the right. The point p(0) is only used for computing
the initial angles, after which it is discarded.
An Output Sensitive Algorithm for Convex Hulls
Reading: Today's material is not covered in our text. Chan's
algorithm can be found in T. Chan, ``Optimal output sensitive convex
hull algorithms in two and three dimensions'', Discrete and
Computational Geometry, 16, 1996, 361--368.
Figure 13: Jarvis's march.
The Omega(n log h) lower bound on convex hulls is given in the paper,
D. G. Kirkpatrick and R. Seidel, ``The ultimate planar convex hull
algorithm,'' SIAM J. Comput., 15, 1986, 287--299. The O(logn) bound
on the complexity of the convex hull is folklore.
Output Sensitive Convex Hull Algorithms: It turns out that in the
worst case, convex hulls cannot be computed faster than in O( n log n)
time. One way to see this intuitively is to observe that the convex
hull itself is sorted along its boundary, and hence if every point
lies on the hull, then computing the hull requires sorting of some
form. Yao proved the much harder result that determining which points
are on the hull (without sorting them along the boundary) still
requiresO( n log n) time. However both of these results rely on the
fact that all (or at least a constant fraction) of the points lie on
the convex hull. This is often not true in practice.
The QuickHull and Jarvis's March algorithms that we saw last time
suggest the question of how fast can convex hulls be computed if we
allow the running time to be described in terms of both the input size
n and the output size h. Many geometric algorithms have the property
that the output size can be a widely varying function of the input
size, and worst case output size may not be a good indicator of what
happens typically. An algorithm which attempts to be more efficient
for small output sizes is called an output sensitive algorithm, and
running time is described as a asymptotic function of both input size
and output size.
Chan's Algorithm: Given than any convex hull algorithm must take at
least O(n) time, and given that ``log n'' factor arises from the fact
that you need to sort the at most n points on the hull, if you were
told that there are only h points on the hull, then a reasonable
target running time is O(n log h). (Below we will see that this is
optimal.) Kirkpatrick and Seidel discovered a relatively complicated
O(n log h) time algorithm, based on a clever pruning method in 1986.
The problem was considered closed until around 10 years later when
Timothy Chan came up with a much simpler algorithm with the same
running time. One of the interesting aspects of Chan's algorithm is
that it involves combining two slower algorithms (Graham's scan and
Jarvis's March) together to form an algorithm that is faster than
either one.
The problem with Graham's scan is that it sorts all the points, and
hence is doomed to having an Omega( n log n) running time, irrespective of
the size of the hull. On the other hand, Jarvis's march can perform
better if you have few vertices on the hull, but it takes O( n) time
for each hull vertex.
Chan's idea was to partition the points into groups of equal
size. There are m points in each group, and so the number of groups is
r = dn=me. For each group we compute its hull using Graham's scan,
which takes O(m log m) time per group, for a total time of O(rm log m)
= O(n log m). Next, we run Jarvis's march on the groups. Here we take
advantage of the fact that you can compute the tangent between a point
and a convex m gon in O(logm) time. (We will leave this as an
exercise.) So, as before there are h steps of Jarvis's march, but
because we are applying it to r convex hulls, each step takes only
O(r log m) time, for a total of O(hr log m) = ((hn/m) log m)
time. Combining these two parts, we get a total of
O((n + (hn/m)) log m) time.
Observe that if we set m = h then the total running time will be O(n
log h), as desired. There is only one small problem here. We do not
know what h is in advance, and therefore we do not know what m should
be when running the algorithm. We will see how to remedy this
later. For now, let's imagine that someone tells us the value of
m. The following algorithm works correctly as long as m – h. If we
are wrong, it returns a special error status.
Chan's Partial Convex Hull Algorithm (Given m)
PartialHull(P; m) :
(1) Let r = ceil(n/m). Partition P into disjoint subsets P(1),P(2),... P(r),
each of size at most m.
(2) For i = 1 to r do:
(a) Compute Hull(P(i)) using Graham's scan and store the vertices in
an ordered array.
(3) Let p0 = (-Inf; 0) and let p1 be the bottommost point of P .
(4) For k = 1 to m do:
(a) For i = 1 to r do:
Compute point q in P(i) that maximizes the angle
p(k-1) p(k) q
(b) Let p(k+1) be the point q in q(1),q(2),...q(r) than maximizes the angle
p(k-1) p(k) q
(c) If p(k+1) = p(1) then return {p(1), p(2), ... p(k)}.
(5) Return ``m was too small, try again.''
Figure 14: Chan's Convex Hull Algorithm.
We assume that we store the convex hulls from step (2a) in an ordered
array so that the step inside the for loop of step (4a) can be solved
in O(logm) time using binary search. Otherwise, the analysis follows
directly from the comments made earlier.
The only question remaining is how do we know what value to give to m?
The last trick, is a common trick used in algorithms to guess the
value of a parameter that affects the running time. We could try m =
1, 2, 3, ... , until we luck out and have m >= h, but this would take
too long. Binary search would be a more efficient option, but if we
guess to large a value for m (e.g. m = n=2) then we are immediately
stuck with O(n log n) time, and this is too slow. Instead, the trick
is to start with a small value of m and increase it rapidly. Since the
dependence on m is only in the log term, as long as our value of m is
within a polynomial of h, that is, m = h^c for some constant c, then
the running time will still be O(n log h). So, our approach will be to
guess successively larger values of m, each time squaring the previous
value, until the algorithm returns a successful result. This trick is
often called doubling search (because the unknown parameter is
successively doubled), but in this case we will be squaring rather
than doubling.
Chan's Complete Convex Hull Algorithm
Hull(P ) :
(1) For t = 1; 2; : : : do:
(a) Let m = min(2^(2^t),n)
(b) Invoke PartialHull(P, m), returning the result in L.
(c) If L != ``try again'' then return L.
Note that 2^2^t has the effect of squaring the previous value of
m. How long does this take? The t-th iteration takes
O(n log 2^2^t ) = O(n* 2^t ) time.
We know that it will stop as soon as 2^2^t >= h, that is if
t = ceil(lg lg n). (We will use lg to denote logarithm base 2.)
So, the total running time (ignoring the constant factors) is:
Sum(t = 1 ... lg lg h) n*2^t =
n * Sum(t = 1 ... lg lg h) 2^t <=
n * 2^(1+lg lg h) =
2n lg h =
O(n log h).
which is just what we want.
Lower Bound: Next we will show that Chan's result is asymptotically
optimal in the sense that any algorithm for computing the convex hull
of n points with h points on the hull requires O( n log h) time. The
proof is a generalization of the proof that sorting a set of n numbers
requires O( n log n) comparisons. If you recall the proof that
sorting takes at least O( n log n) comparisons, it is based on the
idea that any sorting algorithm can be described in terms of a
decision tree. Each comparison has at most 3 outcomes (!, =, or
?). Each such comparison corresponds to an internal node in the
tree. The execution of an algorithm can be viewed as a traversal along
a path in the resulting 3-ary tree. The height of the tree is a lower
bound on the worst case running time of the algorithm. There are at
least n! different possible inputs, each of which must be reordered
differently, and so you have a 3-ary tree with at least n leaves. Any
such tree must have
O(log(base 3) n!) height. Using Stirling's approximation which says
(n! is approximately sqrt(2*PI*n)*(n/e)^n),
this solves to O( n log n) height.
We will give an O( n log h) lower bound for the convex hull
problem. In fact, we will give an O( n log h) lower bound on the
following simpler decision problem, whose output is either yes or no.
Convex Hull Size Verification Problem (CHSV):
Given a point set P and integer h, does the convex hull of P have h
distinct vertices?
Clearly if this takes O( n log h) time, then computing the hull must
take at least as long. As with sorting, we will assume that the
computation is described in the form of a decision tree. Assuming
that the algorithm uses only comparisons is too restrictive for
computing convex hulls, so we will generalize the model of computation
to allow for more complex functions. We assume that we are allowed to
compute any algebraic function of the point coordinates, and test the
sign of the resulting function. The result is called a algebraic
decision tree. The input to the CHSV problem is a sequence of 2n = N
real numbers. We can think of these numbers as forming a vector
(z(1), z(2), ... z(N) = vector(z) in R^N, (a vector z in N-dimensional space),
which we will call a configuration.
Each node of the decision tree is associated with a multivariate
algebraic formula of degree at most d. for example:
f(z) = z(1)*z(4) - 2z(3)z(6) + 5(z(6)^2).
would be an algebraic function of degree 2. The node branches in one
of three ways, depending on whether the result is negative, zero, or
positive. (Computing orientations and dot products both fall into this
category.) Each leaf of the resulting tree corresponds to a possible
answer that the algorithm might give.
For each input vector z to the CHSV problem, the answer is either
``yes'' or ``no''. The set of all ``yes'' points is just a subset of
points Y in R^N , that is a region in this space. Given an arbitrary
input z the purpose of the decision tree is to tell us whether this
point is in Y or not. This is done by walking down the tree,
evaluating the functions on z and following the appropriate branches
until arriving at a leaf, which is either labeled ``yes'' (meaning z
in Y) or ``no''. An abstract example (not for the convex hull problem)
of a region of configuration space and a possible algebraic decision
tree (of degree 1) is shown in the following figure. (We have
simplified it by making it a binary tree.) In this case the input is
just a pair of real numbers.
Figure 15: The geometric interpretation of an algebraic decision tree.
We say that two points u, v in Y are in the same connected component
of Y if there is a path in R^N from u to v such that all the points
along the path are in the set Y . (There are two connected components
in the figure.) We will make use of the following fundamental result
on algebraic decision trees, due to Ben Or. Intuitively, it states
that if your set has M connected components, then there must be at
least M leaves in any decision tree for the set, and the tree must
have height at least the logarithm of the number of leaves.
Theorem: Let Y in R^N be any set and let T be any d th order algebraic
decision tree that determines membership in W . If W has M disjoint
connected components, then T must have height at least Omega(log M - N).
We will begin our proof with a simpler problem.
Multiset Size Verification Problem (MSV):
Given a multiset of n real numbers and an integer k, confirm that the
multiset has exactly k distinct elements.
Lemma: The MSV problem requires O( n log k) steps in the worst case in
the d th order algebraic decision tree
Proof: In terms of points in R^n , the set of points for which the
answer is ``yes'' is
Y = {(z(1), z(2),... z(n) ) in R^n such that |{z(1), z(2), ... ,
z(n)}| = k}
It suffices to show that there are at least k!*k^(n-k) different
connected components in this set, because by Ben Or's result it would
follow that the time to test membership in Y would be
Omega(log(k!*k^( n-k )) - n) =
Omega(k log k + (n - k) log k - n) =
Omega( n log k)
Consider the all the tuples (z(1) ... z(n) ) with z(1), ..., z(k) set
to the distinct integers from 1 to k, and z(k+1) ... z(n) each set to
an arbitrary integer in the same range. Clearly there are k! ways to
select the first k elements and k^(n-k) ways to select the
remaining elements. Each such tuple has exactly k distinct items, but
it is not hard to see that if we attempt to continuously modify one of
these tuples to equal another one, we must change the number of
distinct elements, implying that each of these tuples is in a
different connected component of Y .
To finish the lower bound proof, we argue that any instance of MSV can
be reduced to the convex hull size verification problem (CHSV). Thus
any lower bound for MSV problem applies to CHSV as well.
Theorem: The CHSV problem requires O( n log h) time to solve. Proof:
Let Z = (z(1), ..., z(n) ) and k be an instance of the MSV problem. We
create a point set p(1), ...,p(n) in the plane where p(i) = (z(i),
z(i)^), and set h = k.
(Observe that the points lie on a parabola, so that all the points are
on the convex hull.) Now, if the multiset Z has exactly k distinct
elements, then there are exactly h = k points in the point set (since
the others are all duplicates of these) and so there are exactly h
points on the hull. Conversely, if there are h points on the convex
hull, then there were exactly h = k distinct numbers in the multiset
to begin with in Z.
Thus, we cannot solve CHSV any faster than O( n log h) time, for
otherwise we could solve MSV in the same time.
The proof is rather unsatisfying, because it relies on the fact that
there are many duplicate points. You might wonder, does the lower
bound still hold if there are no duplicates? Kirkpatric and Seidel
actually prove a stronger (but harder) result that the O( n log h)
lower bound holds even you assume that the points are distinct.
Expected Number of Hull Points
For different distributions a different number of points is ex pected
to lie on the hull. One of the simplest cases to analyze is that of
points uniformly distributed in a unit square.
Theorem: If n points are sampled from a uniform distribution in a unit
square, then the expected number of points on the convex hull is O(log
n). Proof: First, we will break the hull into four parts, the upper
left, upper right, lower left, and lower right hulls. This is done by
breaking the hull at its leftmost, rightmost, topmost and bottommost
points. We will show that each has O(log n) points. By symmetry, we
may consider one, say the upper right hull.
Rather than bounding the number of points on the hull, we will bound a
larger quantity. A point p(i) = (x(i), y(i) ) is said to be dominated
by another point p(j) if x(j) >= x(i) and y(j) >= y(i) . The maxima
of a point set P are the points that are not dominated by any other
point in P . Observe that each vertex of the convex hull is a maxima
(but not vice versa), and therefore there are at least as many maxima
points as convex hull vertices. So it suffices to prove that the
expected number of maxima is O(log n).
Suppose that the points are sorted in decreasing order of x
coordinate. We will assume that all points have distinct coordinates
to simplify things. Clearly p(i) is in the maxima if and only if p(i)
has the largest y coordinate among the subset {p(1), p(2),
... p(i)}. What is the probability that this is true? The y
coordinates of these points are independent and identically
distributed. Therefore each of these points is equally likely to be
the maximum, and so the probability that p(i) has the largest y
coordinate is just 1/i. Thus, each point p(i) is a maxima with
probability 1/i.
Figure 16: Upper right hull and maxima.
The expected number of maxima is the sum of these expectations:
sum(i = 1..n) 1/i
The sum is the harmonic series and it is well known that it is very
close to ln n (see Cormen, Leiserson, and Rivest). From the comments
made earlier it follows that the expected number of points on all the
hulls is at most four times this quantity, which is O(log n).
This course is modeled on the Computational Geometry Course taught by Dave Mount, at the University of Maryland. These notes are modifications of his Lecture Notes, which are copyrighted as follows: Copyright, David M. Mount, 2000, Dept. of Computer Science, University of Maryland, College Park, MD, 20742. These lecture notes were prepared by David Mount for the course CMSC 754, Computational Geometry, at the University of Maryland, College Park. Permission to use, copy, modify, and distribute these notes for educational purposes and without fee is hereby granted, provided that this copyright notice appear in all copies.