                    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, xmin and xmax.
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 [xmin, xmax] into k equal subranges, B1 to Bk, 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 Bi with the minimum and maximum y values in that bin. Let:
1. B be the bin for . Thus, B.min = , and B.max = .
2. B[i] (i=1,k) be the bin for the subrange Bi .
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 Bi.
2. B[i].max = the y-maximum point P above 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 , 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

 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 = the join of the  lower and upper hulls.    Output: = 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 of it.

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.x, xmax = P.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];
}

// Next, get the max and min points in the k range bins
Bin*   B = new Bin[k+2];   // first allocate the bins
B.min = minmin;          B.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
}

## 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)   