Building the Convex Hull
Dr. Dobb's Portal
September, 2007 Article on DDJ site |
Finding the Convex Hull of a set of points is an interesting problem in computational geometry. It is useful as a building block for a diverse set of applications, including thing such as:
- Collision detection in video games, providing a useful replacement for bounding boxes.
- Visual pattern matching/object detection
- Mapping
- Path determination
Just as an example, if one of the Mars rovers has to chart a path around a boulder, the convex hull can be used to provide the shortest path past the obstacle, given a map that shows the points where the boulder abuts the ground.
This article will go over the definition of the 2D convex hull, describe Graham’s efficient algorithm for finding the convex hull of a set of points, and present a sample C++ program that can be used to experiment with the algorithm.
The Convex Hull
There are many ways to draw a boundary around a set of points in a two-dimensional plane. One of the easiest to implement is a bounding-box, which is a rectangle that spans the set from its minimum and maximum points in the X and Y planes.
Creating a bounding-box is easy, but it doesn’t form as tight a wrapper as we might like around a set of points. Consider the bounding box around the three points shown in Figure 1.
Figure 1 - A standard bounding box around three points
We can certainly wrap those points much more tightly using easy-to-compute straight lines, and Figure 2 shows an example that is significantly more compact:
Figure 2 - A convex hull around the three points from Figure 1
As it happens, Figure 2 is a convex hull.
So what is the definition of a convex hull? The common visualization analogy for a 2D convex hull is to imagine the set of points on the plane as nails pounded into a board. If you wrap the entire set in an appropriately sized rubber band, the band will snap into place, forming a convex hull, which is the minimum-energy wrapper that encloses all the points.
An informal definition that has a little more precision but is still easy to understand might say that the convex hull meets the following properties:
- The hull is a cycle graph whose vertices are composed of a subset of the the points in set S.
- No points in S lie outside the graph.
- All interior angles in the graph are less than 180 degrees.
Computing the Convex Hull
So given a set of points, how do we compute the convex hull without benefit of hammer, nails, and rubber bands?
For some problems, a brute force solution is adequate. In the case of a convex hull, a reasonable brute force algorithm might look like this:
for all points p in S
for all points q in S
if p != q
draw a line from p to q
if all points in S except p and q lie to the left of the line
add the directed vector pq to the solution set
if all points in S except p and q lie to the right of the line
add the directed vector pq to the solution set
After running this algorithm, you’ve got a list of point pairs that compose the solution set, and you simply have to put them together in the correct order.
This solution will work, but with just a quick look at the code you can see a big problem - a triply nested loop that runs over the magnitude of N, making this an O(n^{3}) algorithm. That’s not going to scale up as well as we might like.
Fortunately, a little searching will show you that there are algorithms that calculate a convex hull in a 2D space considerably faster - in O(n·lgn) time, as a matter of fact.
The Graham Scan
The algorithm I’ll demonstrate here is referred to as the Andrew’s variant of the Graham scan. Ronald Graham’s 1972 paper [1] proposed a convex hull construction algorithm that ran in O(n·lgn) time, and Andrews variation is a simplification that requires a bit less computation. First I’ll give the terse definition of the algorithm, then explain each step in more detail.
Algorithm ConvexHull( S ) Sort all points in S based on their position on the X axis Designate point left as the leftmost point Designate point right as the rightmost point Remove left and right from S While there are still points in S remove Point from S if Point is above the line from left to right add Point to the end of array upper else add Point to the end of array lower // // Construct the lower hull // Add left to lower_hull While lower is not empty add lower[0] to the end of lower_hull remove lower[0] from lower while size(lower_hull >= 3 and the last 3 points lower_hull are not convex remove the next to last element from lower_hull // // Construct the upper hull // Add left to upper_hull While upper is not empty add upper[0] to the end of upper_hull remove upper[0] from upper while size(upper_hull >= 3 and the last 3 points upper_hull are not convex remove the next to last element from upper_hull // Merge upper_hull and lower_hull to form hull return hull
The details
The algorithm starts with an unordered set of points defined by cartesian coordinates - each point has a position on the X-axis and Y-axis. To illustrate the algorithm we’ll start with the points shown in Figure 3.
Figure 3 - The raw points used as input to the algorithm
Before construction of the upper and lower hull can take place, we have to first sort the input data based on its X-axis value, then partition the resulting set into a leftmost point, a rightmost point, and the sets of points above and below the line between the leftmost and rightmost point.
Obviously once the data is sorted it’s trivial to find the leftmost and rightmost points, and remove them from the set of points - these are just the first and last members of the sorted array.
Sorting the remaining points into the upper and lower sets requires that we have some function that determines whether a point is above or below a line. I accomplish this using the following strategy. Given a set of points on a line, p0, p1, and p2, I first perform a coordinate translation so that p1 is at 0,0. I then take the determinant of p0 and p2. The resulting value will be negative if p2 angled off in the left direction, positive if it has moved to the right, and 0 if it is colinear with the first two points.
The matrix that is used to get this determinant looks like this:
A= | (p0_{x}-p1_{x}) | (p2_{x}-p1_{x}) |
(p0_{y}-p1_{y}) | (p2_{y}-p1_{y}) |
For this 2x2 matrix, the formula for the determinant will be:
det = ((p0_{x}-p1_{x})x(p2_{y}-p1_{y})) - ((p2_{x}-p1_{x})x(p0_{y}-p1_{y}))
Partitioning set S into upper and lower is simply a matter of iterating over each point, calculating the determinant, and moving it into lower for a value ≥ 0, or upper for a value < 0.
Once this is complete, the resulting partitions will look like the ones shown in Figure 4.
Figure 4 - The set of points after partitioning
Green points are in the upper partition, red in the lower
Now that things are partitioned, the actual construction of the hull can begin. From the algorithm description, you can see that the upper and lower hull construction steps are symmetrical. Both proceed from right to left, starting at the leftmost point and moving to the rightmost point. The only difference is the source of their input points, and the direction they check to insure convexity.
The upper or lower hull is started by simply adding left to the output hull. Points are then added from the correct input source. As each point is added, if the number of points in the working hull is equal to 3 or more, a test is made to see if the last three points have created a convex angle.
Testing for the convex angle is done using the same determinant formula as shown above. If the hull has n points, we simply test to see if p_{n-1} is above or below the line formed by p_{n-2} and p_{n}. When constructing the lower hull, if we see that the point is above the line, we have violated convexity, and the middle point is removed from the hull. The opposite test is made when constructing the upper hull. This test-and-remove process is repeated until the last three points are convex, or there are fewer than 3 points in the working hull.
An animation of this process is shown in Figure 5.
Figure 5 - An animation of the hull being built
The final step, merging the upper and lower hulls, is a trivial matter of appending one hull to the other, and removing the extra copy of right. Once that is done, the actual convex hull definition can be given as a list of points, starting with left and moving counter-clockwise around the hull. Figure 6 is the last frame in the animation, which shows the complete hull.
Figure 6 - The finished convex hull
Test Code
A small demo program called graham.cpp implements this algorithm fairly faithfully. In order to make the experimentation a little more interesting, graham provides two forms of output showing the results of the program:
- Standard text output listing the various data sets used in the program
- A command file (gnuplot.cmd) that can be used with gnuplot to visualize the process
The images shown in this article were collected using gnuplot 4, which is a fully featured 2D and 3D plotting program, with excellent multiplatform support. Seeing the algorithm operate visually in real time is very helpful in gaining a good understanding of how it works.
The C++ program utilizes a class called GrahamScan
that takes care of all these
details. By creating the object and then calling its methods, creation and display of the convex
hull is easy to follow. In my test program, main()
executes the entire operation by
creating the object and then calling its methods:
If you ignore method calls that start with log_ or plot_, the real work takes place in only three steps:
- Construct the
GrahamScan
object. This also creates the random set of points. - Partition the points into the upper and lower hull sets.
- Build the convex hull.
The constructor
When calling the GrahamScan
constructor, you will pass in five numbers: the number
of points, and the min and max values for the X and Y axis. The min and max values not only bound
the range of the randomly generated points, they also determine the scope of the axes that will
be displayed when the values are shown in gnuplot.
The bulk of the work in the constructor is in these few lines of code:
The data set is stored as std::pair<int,int>
objects in a vector aptly called
raw_points
. After calling the constructor, the log_raw_points()
method
can be called, which will product output like this:
Creating raw points: (97,90) (27,10) (59,8) (58,19) (85,90) (62,91) (94,42) (84,68) (16,21) (49,14) (31,84) (40,25) (59,95) (55,89) (81,95) (22,46) (27,80) (18,90) (59,37) (38,45) |
Calling plot_raw_points()
will create a gnuplot command file that produces output
like that shown in Figure 3.
The Partitioning Code
The algorithm definition tells us that the partition step needs to identify the leftmost point, the rightmost point, and the two sets of points above and below the line between leftmost and rightmost.
This is all accomplished in a straightforward manner:
Note that by storing the points in an std::pair<int,int>
, with x in the
first
member and y in the second
member, we can use the standard
library sort()
routine to order the array.
After sorting and then moving each point into one of the two partitions, we have the following
members of GrahamScan
defined and ready to use in building the hull:
left
, the leftmost point in the set of pointsright
, the rightmost point in the set of pointsupper_partition_points
, the points above the line betweenleft
andright
.lower_partition_points
, the points below or on the line betweenleft
andright
.
If you call the log_partitioned_points()
method, you’ll get output that looks
something like this after partitioning:
Partitioned set: Left : (16,21) Right : (97,90) Lower partition: (27,10)(40,25)(49,14)(58,19)(59,8)(59,37)(84,68) (94,42) Upper partition: (18,90)(22,46)(27,80)(31,84)(38,45)(55,89)(59,95) (62,91)(81,95) (85,90) |
Calling the plot_partitioned_points() creates a gnuplot command sequence that will display a plot like that shown in Figure 4.
Creating the Hull
The actual creation of the hull is done by method build_hull
, which calls method
build_half_hull
twice, once with the points on the lower hull, and once with the
points on the upper hull. The output of build_half_hull
is sent on the first call to
array lower_hull
and next to array upper_hull
:
So the bulk of the work is actually done in build_half_hull
, which looks like this:
The main loop in this routine simply pulls points out of the input
array, adds them
to the output
array, and then performs the check to make sure that convexity has not
been violated. If it has, points are removed until it is again correct. Figure 5 gives a nice
animated view of the process.
Once this is done, calling the log_hull routine produces output that looks like this:
Lower hull: (16,21)(27,10)(59,8)(94,42)(97,90) Upper hull: (16,21)(18,90)(59,95)(81,95)(97,90) Convex hull: (16,21) (27,10) (59,8) (94,42) (97,90) (81,95) (59,95) (18,90) (16,21) |
The plot_hull()
method can then be called to create a gnuplot command file that will
display the convex hull, as shown in Figure 6.
There are quite a few variations on the 2-D convex hull building process, and this program ought to be amenable to trying out many of them. If you use the existing data structures and just change the algorithms, you can use the existing gnuplot routines to animate your work and get a good feel for how it is working. Enjoy!
References and Links
Source code:
graham.zip, which contains VS 2003 and 2005 solutions, plus a g++ Makefile.
[1] R.L. Graham, An efficient algorithm for determining the convex hull of a finite planar set, Info. Proc. Lett. 1, 132-133 (1972).
[2] A. M. Andrew. Another efficient algorithm for convex hulls in two dimensions. Inform. Process. Lett., 9(5):216-219, 1979. (A note about the this algorithm can be found here.)