Logo_iSurf_org

Subscribe
to eAlerts

for Geom Site
Updates

CLICK HERE

 

A Fast  Convex Hull Approximation

by Dan Sunday

 

 

In Algorithm 10, we looked at some of the fastest O(nlog-n) 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 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] (BFP) 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(nlog-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 (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 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, xmin and xmax.
    2. Pminmin and Pminmax be the points with P.x=x-min first and then min y or max y respectively.
    3. Pmaxmin and Pmaxmax 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 Pminmin and Pmaxmin.
    5. L-max is the upper line joining the two points Pminmax and Pmaxmax.
  2. Divide the range [xmin, xmax] into k equal subranges, B1 to Bk, each of width wk=(xmax-xmin).div.k. That is: Bi=[xi-1,xi] where xi=xmin+iwk.
  3. Define a (k+2)-element array of bins B[] to record the points B[].min and B[].max in each Bi with the minimum and maximum y values in that bin. Let:
    1. B[0] be the bin for P.x=x-min. Thus, B[0].min = Pminmin, and B[0].max = Pminmax.
    2. B[i] (i=1,k) be the bin for the subrange Bi .
    3. B[k+1] be the bin for P.x=x-max. Thus, B[k+1].min = Pmaxmin, and B[k+1].max = Pmaxmax.
  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 Bi.
    2. B[i].max = the y-maximum point P above L-max with P.x in Bi.


Note that for some ranges Bi 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.

approx_hull-2

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

    Input: a set S = {P = (P.x,P.y)} of n points
           k = the approximation accuracy (large k = more accurate)

    Get the points with 1st x min or max and 2nd y min or max
        minmin = index of P with min x first and min y second
        minmax = index of P with min x first and max y second
        maxmin = index of P with max x first and min y second
        maxmax = index of P with max x first and max y second.
    Let L_min be the lower line joining P[minmin] with P[maxmin].
    Let L_max be the upper line joining P[minmax] with P[maxmax].

    Let B[k+2] be an array of subrange bins as described above.
    Compute the y-min and y-max in each subrange i=1,k to get:
        B[i].min = the P in the i-th subrange with min y below L_min
        B[i].max = the P in the i-th subrange with max y above L_max

    Compute the lower hull stack as follows:
    (1) Push P[minmin] onto the stack.
    (2) for (each bin B[i], i = 1 to k)
        {
            if (B[i].min does not exist)
                 continue;
            Let Pi = B[i].min.
            while (there are at least 2 points on the stack)
            {
                 Let PT1 = the top point on the stack.
                 Let PT2 = the second point on the stack.
                 if (Pi is strictly left of the line from PT2 to PT1)
                     break out of this while loop.
                 Pop the top point PT1 off the stack.
            }
            Push Pi onto the stack.
        }
    (3) Push P[maxmin] onto the stack.

    Similarly, compute the upper hull stack.

    Let OMEGA = the join of the  lower and upper hulls.

    Output: OMEGA = the convex hull of S.


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:

[Preparata & Shamos, 1985, p-156, Theorem  4.6
Any point P in S that is not inside the approximate convex hull is within distance (xmax-xmin).div.k of it.


Thus, the approximation can be made arbitrarily close by picking a large enough k. This is easily illustrated by the following diagram.

approx_error

which demonstrates that a point P below the lower hull and above the bin y-minimum, has an error: delta.lt.d(P,Q).lt.(xmax-xmin).div.k.


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