Algorithm 11
[Home] [Overview] [History] [Algorithms] [Books] [Web Sites] [Gift Shop]

A Fast Approximate
2D Convex Hull Algorithm

by Dan Sunday

Last month we looked at some of the fastest O(n log n) algorithms for computing The Convex Hull of a Planar Point Set or Polygon.  This month, we present an algorithm that gives an approximation for the 2D convex hull.  The advantage of this algorithm is that it is much faster with just an O(n) runtime.  In many applications, an approximate hull suffices, and the gain in speed can be significant for very large point sets or polygons.

The BFP Approximate Hull Algorithm

Sometimes a good approximation for the convex hull is sufficient for an application, especially where the n data points have a sampling error.  The main advantage of an approximation algorithm would be that it is significantly more efficient.  In 2D, a very efficient approximate convex hull algorithm is the one of [Bentley-Faust-Preparata, 1982] which runs in O(n) time.

It is known that the speed of an algorithm for the convex hull of a 2D point set S is dominated by the need to initially sort the n points of the set, which takes O(n log n) time.  After that, the best algorithms, such as the Graham Scan or the Monotone Chain, only require O(n) time.  So, to improve efficiency, one must replace the sort operation with something else.  The Bentley-Faust-Preparata (BFP) algorithm does this by dividing the plane into narrow vertical strips, and then (in a single linear pass) finding the maximum and minimum points of S in each strip.  Then, they use the minimum points to get the lower hull chain, and the maximum points to get the upper hull chain.  These chains are computed using the monotone chain algorithm [Andrew, 1979] that we previously presented.  The "sorting" of the maximum and minimum point sets comes for free in the natural ordering of the vertical strip bins used by the algorithm.  So, we have replace sorting with a linear O(n) procedure.  But, one pays for this by possibly throwing away points which are non-extreme in their vertical strip bin, yet are still on the convex hull boundary.  However, they can't be too far from that convex hull, and thus we get a good approximation to the actual convex hull.  Further, the approximation gets better as we make the strips narrower.  The extreme points in these vertical strips are circled in the following diagram that also shows the approximate hull as a bold dashed contour.  In this example, several points are outside the approximate hull.

More accurately, the BFP algorithm finds the extrema of the vertical strip bins as follows:

1. Scan S once to find, in O(n) time:

   1. the minimum and maximum x-coordinates, x_min and x_max.

   2. P_-- and P_-+ be the points with P.x = x_min first and then min y or max y respectively.

   3. P_+- and P₊₊ be the points with P.x = x_max first and then min y or max y respectively.

   4. L_min is the lower line joining the two points P_-- and P_+- .

   5. L_max is the upper line joining the two points P_-+ and P₊₊ .

2. Divide the range [x_min, x_max] into k equal subranges, B₁ to B_k, each of width w_k = (x_max-x_min) / k .  That is:  B_i = [x_i-1 , x_i ] where x_i = x_min + i w_k .

3. Define a (k+2)-element array of bins B[] to record the points in each B_i with a y-min and y-max values.  Let:

   1. B[0] be the bin for P.x = x_min.  Thus, B[0].min = P_--, and B[0].max = P_-+ .

   2. B[i] (i=1,k) be the bin for the subrange B_i .

   3. B[k+1] be the bin for P.x = x_max.  Thus, B[k+1].min = P_+-, and B[k+1].max = P₊₊ .

4. Scan S again to find, in O(n) time:

   1. B[i].min = the y-minimum point P below L_min with P.x in B_i .

   2. B[i].max = the y-maximum point P above L_max with P.x in B_i .

Note that for some ranges B_i these points, B[i].min and B[i].max, will not exist.  But that is ok because all we are trying to do is select points that are valid candidates for the lower or upper convex hulls.  That is why we only consider minimum points below L_min, since the lower hull must be below this line.  Similarly, only maximum points above L_max are considered since the upper hull is above this line.  These selections are shown in the following diagram.  Notice that there are fewer extreme points due to filtering with respect to L_min and L_max.  This will make the algorithm run slightly faster.

The algorithm now proceeds by computing the lower and upper hulls using a (simplified) variation of Andrew's monotone chain hull algorithm.  Putting it all together in pseudo-code, we have:

Pseudo-Code: BFP Approximate Hull Algorithm

Our implementation nearHull_2D() of this algorithm in C++ is given below.

Error Analysis

As for how close the approximate hull is to the exact convex hull, it can be shown that:

Thus, the approximation can be made arbitrarily close by picking a large enough k.  The proof of this fact is easily illustrated by the following diagram.

which demonstrates that a point P below the lower hull and above the bin minimum, has an error d < d(P,Q) < (x_max-x_min) / k .

Implementation

Here is a "C++" implementation of the approximate hull algorithm.

// Copyright 2001, softSurfer (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.

// Assume that classes are already given for the objects:
//    Point with coordinates {float x, y;}
//===================================================================

// isLeft(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2 on the line
//            <0 for P2 right of the line
//    See: the January 2001 Algorithm on Area of Triangles
inline float
isLeft( Point P0, Point P1, Point P2 )
{
    return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
}
//===================================================================
 

#define NONE  (-1)

typedef struct range_bin Bin;
struct range_bin {
    int    min;    // index of min point P[] in bin (>=0 or NONE)
    int    max;    // index of max point P[] in bin (>=0 or NONE)
};
 

// nearHull_2D(): the BFP fast approximate 2D convex hull algorithm
//     Input:  P[] = an (unsorted) array of 2D points
//              n = the number of points in P[]
//              k = the approximation accuracy (large k = more accurate)
//     Output: H[] = an array of the convex hull vertices (max is n)
//     Return: the number of points in H[]
int
nearHull_2D( Point* P, int n, int k, Point* H )
{
    int    minmin=0,  minmax=0;
    int    maxmin=0,  maxmax=0;
    float  xmin = P[0].x,  xmax = P[0].x;
    Point* cP;                 // the current point being considered
    int    bot=0, top=(-1);  // indices for bottom and top of the stack

    // Get the points with (1) min-max x-coord, and (2) min-max y-coord
    for (int i=1; i<n; i++) {
        cP = &P[i];
        if (cP->x <= xmin) {
            if (cP->x < xmin) {        // new xmin
                xmin = cP->x;
                minmin = minmax = i;
            }
            else {                      // another xmin
                if (cP->y < P[minmin].y)
                    minmin = i;
                else if (cP->y > P[minmax].y)
                    minmax = i;
            }
        }
        if (cP->x >= xmax) {
            if (cP->x > xmax) {        // new xmax
                xmax = cP->x;
                maxmin = maxmax = i;
            }
            else {                      // another xmax
                if (cP->y < P[maxmin].y)
                    maxmin = i;
                else if (cP->y > P[maxmax].y)
                    maxmax = i;
            }
        }
    }
    if (xmin == xmax) {      // degenerate case: all x-coords == xmin
        H[++top] = P[minmin];           // a point, or
        if (minmax != minmin)           // a nontrivial segment
            H[++top] = P[minmax];
        return top+1;                   // one or two points
    }

    // Next, get the max and min points in the k range bins
    Bin*   B = new Bin[k+2];   // first allocate the bins
    B[0].min = minmin;         B[0].max = minmax;        // set bin 0
    B[k+1].min = maxmin;       B[k+1].max = maxmax;      // set bin k+1
    for (int b=1; b<=k; b++) { // initially nothing is in the other bins
        B[b].min = B[b].max = NONE;
    }
    for (int b, i=0; i<n; i++) {
        cP = &P[i];
        if (cP->x == xmin || cP->x == xmax) // already have bins 0 and k+1
            continue;
        // check if a lower or upper point
        if (isLeft( P[minmin], P[maxmin], *cP) < 0) {  // below lower line
            b = (int)( k * (cP->x - xmin) / (xmax - xmin) ) + 1;  // bin #
            if (B[b].min == NONE)       // no min point in this range
                B[b].min = i;           // first min
            else if (cP->y < P[B[b].min]->y)
                B[b].min = i;           // new min
            continue;
        }
        if (isLeft( P[minmax], P[maxmax], *cP) > 0) {  // above upper line
            b = (int)( k * (cP->x - xmin) / (xmax - xmin) ) + 1;  // bin #
            if (B[b].max == NONE)       // no max point in this range
                B[b].max = i;           // first max
            else if (cP->y > P[B[b].max]->y)
                B[b].max = i;           // new max
            continue;
        }
    }

    // Now, use the chain algorithm to get the lower and upper hulls
    // the output array H[] will be used as the stack
    // First, compute the lower hull on the stack H
    for (int i=0; i <= k+1; ++i)
    {
        if (B[i].min == NONE)  // no min point in this range
            continue;
        cP = &P[ B[i].min ];   // select the current min point

        while (top > 0)        // there are at least 2 points on the stack
        {
            // test if current point is left of the line at the stack top
            if (isLeft( H[top-1], H[top], *cP) > 0)
                break;         // cP is a new hull vertex
            else
                top--;         // pop top point off stack
        }
        H[++top] = *cP;        // push current point onto stack
    }

    // Next, compute the upper hull on the stack H above the bottom hull
    if (maxmax != maxmin)      // if distinct xmax points
        H[++top] = P[maxmax];  // push maxmax point onto stack
    bot = top;                 // the bottom point of the upper hull stack
    for (int i=k; i >= 0; --i)
    {
        if (B[i].max == NONE)  // no max point in this range
            continue;
        cP = &P[ B[i].max ];   // select the current max point

        while (top > bot)      // at least 2 points on the upper stack
        {
            // test if current point is left of the line at the stack top
            if (isLeft( H[top-1], H[top], *cP) > 0)
                break;         // current point is a new hull vertex
            else
                top--;         // pop top point off stack
        }
        H[++top] = *cP;        // push current point onto stack
    }
    if (minmax != minmin)
        H[++top] = P[minmin];  // push joining endpoint onto stack

    delete B;                  // free bins before returning
    return top+1;              // # of points on the stack
}

References

A.M. Andrew, "Another Efficient Algorithm for Convex Hulls in Two Dimensions", Info. Proc. Letters 9, 216-219 (1979)

Jon Bentley, G.M. Faust & Franco Preparata, "Approximation Algorithms for Convex Hulls", Comm. ACM 25, 64-68 (1982)

Franco Preparata & Michael Shamos, Computational Geometry: An Introduction, Section 4.1.2 "Approximation Algorithms for Convex Hull" (1985)

Copyright © 2001-2006 softSurfer.  All rights reserved.
Email comments and suggestions to feedback@softsurfer.com