Subscribe
to eAlerts
for Geom Site
Updates

CLICK HERE

by Dan Sunday

In Algorithm 10, we looked at some of the fastest algorithms for computing The Convex Hull of a Planar Point Set. We now present an algorithm that gives a fast approximation for the 2D convex hull. The advantage of this algorithm is that it is much faster with just an 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] (BFP) which runs in 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 time. After that, the best algorithms, such as the Graham Scan or the Monotone Chain, only require 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 (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 a variation of the monotone chain algorithm [Andrew, 1979] that we presented in Algorithm 10. 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 replace sorting with a linear 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 time:
1. the minimum and maximum x-coordinates, x_min and x_max.
2. and be the points with first and then min y or max y respectively.
3. and be the points with first and then min y or max y respectively.
4. is the lower line joining the two points and .
5. is the upper line joining the two points and .
   
2. Divide the range [x_min, x_max] into k equal subranges, B₁ to B_k, each of width . That is: where .
   
3. Define a (k+2)-element array of bins B[] to record the points B[].min and B[].max in each B_i with the minimum and maximum y values in that bin. Let:
1. B[0] be the bin for . Thus, B[0].min = , and B[0].max = .
2. B[i] (i=1,k) be the bin for the subrange B_i .
3. B[k+1] be the bin for . Thus, B[k+1].min = , and B[k+1].max = .
   
4. Scan S again to find, in time:
1. B[i].min = the y-minimum point P below with P.x in B_i.
2. B[i].max = the y-maximum point P above 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 , since the lower hull must be below this line. Similarly, only maximum points above 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 and . 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. This is easily illustrated by the following diagram.

which demonstrates that a point P below the lower hull and above the bin y-minimum, has an error: .

Implementation

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

// Copyright 2001 softSurfer, 2012 Dan Sunday
// 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: Algorithm 1 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 2012 Dan Sunday, 2001 softSurfer