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:

PLAIN TEXT

CODE:

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

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³) 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:

PLAIN TEXT

C++:

int main(int argc, char* argv[])
{
std::ofstream gnuplot_file( "gnuplot.cmd" );
const int N = 20;
GrahamScan g( N, 0, 100, 0, 100 );
g.log_raw_points( std::cout );
g.plot_raw_points( gnuplot_file );
g.partition_points();
g.log_partitioned_points( std::cout );
g.plot_partitioned_points( gnuplot_file );
//
// Build the hull
//
g.build_hull( gnuplot_file );
g.log_hull( std::cout );
g.plot_hull( gnuplot_file, "complete" );
return 0;
}

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:

PLAIN TEXT

C++:

srand( static_cast<unsigned int>( time(NULL) ) );
for ( size_t i = 0 ; i <N ; i++ ) {
int x = ( rand() % ( x_range.second - x_range.first + 1 ) ) + x_range.first;
int y = ( rand() % ( y_range.second - y_range.first + 1 ) ) + y_range.first;
raw_points.push_back( std::make_pair( x, y ) );
}

The data set is stored as std::pair 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:

PLAIN TEXT

C++:

void partition_points()
{
//
// Step one in partitioning the points is to sort the raw data
//
std::sort( raw_points.begin(), raw_points.end() );
//
// The the far left and far right points, remove them from the
// sorted sequence and store them in special members
//
left = raw_points.front();
raw_points.erase( raw_points.begin() );
right = raw_points.back();
raw_points.pop_back();
//
// Now put the remaining points in one of the two output sequences
//
for ( size_t i = 0 ; i <raw_points.size() ; i++ )
{
int dir = direction( left, right, raw_points[ i ] );
if ( dir <0 )
upper_partition_points.push_back( raw_points[ i ] );
else
lower_partition_points.push_back( raw_points[ i ] );
}
}

Note that by storing the points in an std::pair, 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 points
right, the rightmost point in the set of points
upper_partition_points, the points above the line between left and right.
lower_partition_points, the points below or on the line between left and right.

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:

PLAIN TEXT

C++:

void build_hull( std::ofstream &f )
{
build_half_hull( f, lower_partition_points, lower_hull, 1 );
build_half_hull( f, upper_partition_points, upper_hull, -1 );
}

So the bulk of the work is actually done in build_half_hull, which looks like this:

PLAIN TEXT

C++:

void build_half_hull( std::ostream &f,
std::vector<std::pair<int,int>> input,
std::vector<std::pair<int,int>> &output,
int factor )
{
// The hull will always start with the left point, and end with the right
// point. Initialize input and output accordingly
//
input.push_back( right );
output.push_back( left );
while ( input.size() != 0 ) {
//
// Repeatedly add the leftmost point to the null, then test to see
// if a convexity violation has occurred. If it has, fix things up
// by removing the next-to-last point in the output sequence until
// convexity is restored.
//
output.push_back( input.front() );
input.erase( input.begin() );
while ( output.size()>= 3 ) {
size_t end = output.size() - 1;
if ( factor * direction( output[ end - 2 ],
output[ end ],
output[ end - 1 ] ) <= 0 ) {
output.erase( output.begin() + end - 1 );
} else
break;
}
}
}

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.)