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 twodimensional plane. One of the easiest to implement is a boundingbox, which is a rectangle that spans the set from its minimum and maximum points in the X and Y planes.
Creating a boundingbox 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 easytocompute 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 minimumenergy 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
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 Xaxis and Yaxis. 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 Xaxis 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_{n1} is above or below the line formed by p_{n2} 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 testandremove 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 counterclockwise 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:

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:

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:

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

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:

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 nexttolast 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 2D 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, 132133 (1972).
[2] A. M. Andrew. Another efficient algorithm for convex hulls in two dimensions. Inform. Process. Lett., 9(5):216219, 1979. (A note about the this algorithm can be found here.)
24 users commented in " Building the Convex Hull "
Followup comment rss or Leave a Trackback[...] » Building the Convex Hull Mark Nelson: Programming, mostly. Finding the Convex Hull of a set of points is an interesting problem in computational geometry. (tags: algorithms c++ graphics math) [...]
Awsome work!!
Very useful. Best exposition of an algorithm for finding a convex hull that I've been able to find.
Hi! Nice article. You left out one valuable thing in both your description and code: testing whether a point is inside an established hull. Maybe just because it's reasonably trivial =)
Perhaps there's a more clever way to do it, but I'd do a O(n) search calling direction() for each line segment on the hull boundary, testing that the point is always on the righthand side.
@geoff:
I don't know if there is a sublinear solution to this problem, but people love their convex hulls, so it should be easy to find.
My solution: be the point, turn in a circle 360 degrees. If all you see is hull all the way around, you're inside.
Very Good explanation
[...] Points[Len] & Points[Len1] respectively. There is a great animation demonstrating this process here. After the upper and lower partitions are processed (giving us 2 sets of points that we could [...]
Very useful. Best exposition of an algorithm for finding a convex hull that I've been able to find [2]. Really saved the day ! Thank u !
Thanks a lot for this article.
The article "[1] R.L. Graham, An efficient algorithm for determining the convex hull of a finite planar set, Info. Proc. Lett. 1, 132133 (1972)." has moved to:
http://www.math.ucsd.edu/~ronspubs/72_10_convex_hull.pdf
Antoon
@van Breda,
Thanks, I updated the link!
 Mark
Nice article.
To elaborate on the pointinconvexhull test. This can be done in O(log n):
Given a point (x,y) look whether it is in the lower or upper halfspace, using the same test as applied in partition_points().
Then lookup the index 'a' of the first coordinate with xvalue larger than x in lower_partition_points if (x,y) is in the lower halfspace.
Lookup the index 'a' of the first coordinate with xvalue larger than x
in upper_partition_points if (x,y) is in the upper halfspace.
This lookup can be done using a binary search (http://en.wikipedia.org/wiki/Binary_search_algorithm) as the arrays are sorted by xcoordinate, So in O(log n) time.
Then get the vector from the coordinate at index 'a1' to 'a'. If (x,y) is in the lower halfspace and on or left of the vector then it is inside the convex hull.
Similarly if (x,y) is is in the upper halfspace and on or right of the vector then it is inside the convex hull.
You can use the direction() function from the c++ code for the vector side test.
If index 'a' or 'a1' does not exist then (x,y) is outside the convex hull or on the hull at rightmost hull point. To check against the latter case, test if x equals 'right' (from the c++ class) and if true,
set 'a' = last index of lower_partition_points (or upper_partition_points depending on whether (x,y) is in the lower or upper halfspace) and test if (x,y) is on the vector from 'a1' to 'a' using direction()
How about a convex hull in dimensions above 2?
Well, this problem is so important in Windows update for the screen shots, fast games, and 3D cg scenes!
Well...
that which you call "a convex hull", I apparently know how (in logic) to create in "an infinite number of dimensions" (i.e. SPACE (itself) becomes "relative" (unless there are _________________ (censored) in which situation space is computed (= maximum efficiency (in theory) (noise considerations can be factored in via a more sophisticated technique I reckon) of storing data (as "realtime "objects")(the whole thing becomes a _______________ (the _____________________)(curiously it is a "perpetual motion machine" in that everything (can) rotate around everything else ("perpetual motion") within internally (or sometimes externally in other versions e.g. teleportation dynamic scenario ..) generated "limits" ("machine" like..)  the whole thing becomes like a solid object, ... an algorithm partial 'interferometer' i.e. you don't need an algorithm to store the data, the data itself " becomes " an ____________(quantum dynamical; interspace (a between space betwen space: a hyperspace bypass!!!!!!!! )
Thanks. Very useful and explanatory.
Question : In class GrahamScan, rather than generating N random points, how can one read data from some file i.e. 2 columns and n rows file with maximum value 1 and minimum value 0, then use the data to compute the convex hull.
Just started programming in C++ so if you could direct me on which files of the graham.zip i should modify and how they should be modified, i'll be very grateful.
If you are just learning C++, I would say the modifications you are asking about will be a great learning exercise.
 Mark
one of the best I found on the internet...
Thanks for the article. I implemented the system in another programming environment and found that very hard to do. It worked for almost all cases for me, but sometimes produced a concave polygon. Probably because of a (my) coding error.
I found a much simpler way to do it, which follows similiar pattern.
1.) Sort for x like you did
2.) Build the upper and lower list like you did
3.) then build another list with all the points going clockwise, that is copy the upper list into it and reverse add the lower list without the endpoints. You can think of a tube now that contains all points in clockwise order. Lets call the list 'inList'.
5.) add the leftmost point to the new list 'outList'.
4.) go through inList drawing a line between point 0 and point 2 and test if point 1 is on the left side. If it is, add it to outList. If not, increase a counter. Pseudocode:
Call this function until it returns 0. Swap the lists for each call (make the inList the outList and vice versa).
After some 'rounds' your convex hull is in outList. This works perfectly for me. Most time (with my input data around 2050 points) the algorithm is done in three rounds, eliminating around 90% of the points in the first round. Even for highly complex scenarios it doesn't take more than 4 rounds. I think it is easier to understand and more straightforward to code.
uups, sorry, code is messed up by formatting. here the next unformatted try:
I hope you're a teacher or a professor somewhere. This is simply the most clear and coherent explanation of a geometry based problem I have ever seen on the internet. Right on down to the animations.
Bravo, sir. Bravo.
Hello this is a wonderful article. Please can you let me know the license on use of the code in this article?
I am looking to use a convex hull algorithm in a closed source game.
I can write my own script but would find it easier to convert your code to C# for use in the project.
Regards,
Daniel Vernall
@Daniel:
Should be no problem, see: http://marknelson.us/codeusepolicy/
Hello,
thank you for the wonderful description. Please correct me if I'm wrong, but I think there is a mistake in your description:
You wrote:
'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.'
However, I believe it is the other way round, if the determinant is positive, it is in the upper set with P0 being the leftmost point, P2 the rightmostpoint and P1 being the point to be sorted?
Leave A Reply
You can insert source code in your comment without fear of too much mangling by tagging it to use the iG:Syntax Hiliter plugin. A typical usage of this plugin will look like this:[c]
int main()
{
printf( "Hello, world!\n" );
return 1;
}
[/c]
Note that tags are enclosed in square brackets, not angle brackets. Tags currently supported by this plugin are: as (ActionScript), asp, c, cpp, csharp, css, delphi, html, java, js, mysql, perl, python, ruby, smarty, sql, vb, vbnet, xml, code (Generic).
If you post your comment and you aren't happy with the way it looks, I will do everything I can to edit it to your satisfaction.