Building the Convex Hull

Test Code

A demo program called graham.cpp (source code available for download here, along with Visual Studio 2003 and 2005 solutions, plus a g++ Makefile) implements this algorithm fairly faithfully. 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 I present 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:

1. Construct the GrahamScan object. This also creates the random set of points.
2. Partition the points into the upper and lower hull sets.
3. 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() creates a gnuplot command file that produces output like that 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 ] );
    }
}

By storing the points in an std::pair, with x in the first member and y in the second member, you 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 in Figure 4.