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

Bounding Containers for Polygons,
Polyhedra, and Point Sets (2D & 3D)

by Dan Sunday

It is often useful to have a bounding container (BC), such as a bounding box or ball, for a finite geometric object.  Using a BC can significantly speed up software for ray tracing, collision avoidance, hidden object detection, etc.  Before invoking a computationally expensive intersection or containment algorithms for a complicated object, often a simple test with an uncomplicated bounding container can exclude the possibility of intersection or containment, and no further computation is needed.  If a point is outside a bounding container, then it is also outside the object inside the container.  If two bounding containers are disjoint, then so are the objects they contain, and thus they cannot intersect.

We will restrict attention to finite linear objects, such as points, segments, triangles, polygons, and polyhedra.  These objects, and sets of them, are specified by linear combinations of their vertices, and their complexity can be measured by the total number n of vertices they have.  All containers we discuss can be computed directly from the set of vertices, and so we just have to consider algorithms for sets of points. 

The usefulness of such containers, especially "the bounding box", is sometimes only given passing mention as an obvious preprocessing test to make before attempting other algorithms.  However, there are a number of different choices for bounding containers that a programmer should be familiar with.  Important characteristics of a good bounding container are:

We consider two basic types of BC's: linear containers (such as rectangular boxes and convex polygons), and quadratic containers (such as spheres or ellipses).  In general, linear containers excel for criteria (C) and (D), whereas quadratic containers excel at (B) but can take more time to compute, and do not approximate the object as well.  An example of different containers for the same object, a polygon, is shown in the diagram:

Linear Containers

A linear container is one whose interior is specified by a finite number of linear inequalities.  This implies that a bounded linear container is either a convex polygon (2D) or a convex polyhedron (3D).  For example, in 2D, a container W could be specified by k inequalities: f_i(x,y) = a_ix + b_iy + c_i £ 0 (i=1,k), all of which would have to be true for a point (x,y) to be in the region.  If any one of the inequalities failed, then the test point would be outside W.  So, when a point is outside W, this can be discovered on the average by testing half of the inequalities.  Also, each inequality defines a (convex) half-space H_i bounded by the line L_i: f_i(x,y) = 0.  The region W is the intersection of all these half-spaces, and this implies that it is also convex.  When bounded, W is a convex polygon with segments of the lines L_i as edges.

Similarly in 3D, a linear container Y can be specified by k linear inequalities: f_i(x,y) = a_ix + b_iy + c_iz + d_i £ 0 (i=1,k), each of which defines a half-space H_i bounded by a plane p_i.  The region inside Y is now a 3D convex polyhedron with convex polygons W_i (contained in the planes p_i) as faces.

The Bounding Box

A "box" is a rectangular region whose edges are parallel to the coordinate axes, and is thus defined by its maximum and minimum extents for all axes.  So, a 2D box is given by all (x,y) coordinates satisfying x_min £ x £ x_max and y_min £ y £ y_max, and is specified by the extreme points (x_min, y_min) and (x_max, y_max) which are its bottom-left and top-right corners.  A 3D box is specified by the extreme corners (x_min, y_min, z_min) and (x_max, y_max, z_max).  Inclusion of a point P=(x,y) in a box is tested by verifying that all inequalities are true; and if any one of them fails, then the point is not inside, but rather outside the box.  In 2D there are 4 inequalities, and in 3D there are 6.  But, on the average, a point will be rejected after 2 tests in 2D and 3 tests in 3D.  Similarly, one can test the disjointness of two boxes B₁ and B₂ by comparing their minimum and maximum extents.  With respect to the x-axis, if either x_max1 < x_min2 or x_max2 < x_min1, then B₁ and B₂ are disjoint.  There are equivalent disjointness tests with respect to the y-axis and z-axis.  If any of the tests is true, then B₁ and B₂ are disjoint.  If all of them are false, then the two boxes are intersect.

The "bounding box" of a finite geometric object is the box with minimal area (in 2D), or minimal volume (in 3D or higher dimensions), that contains the object.   For any collection of linear objects (points, segments, polygons, and polyhedra), their bounding box is given by the minimum and maximum coordinate values for the point set S of the object's vertices: (x_min, x_max, y_min, y_max) in 2D, plus (z_min, z_max) in 3D.  These values are easily computed in O(n) time with a single scan of all the vertex points, sometimes while the object's vertices are being read or computed. 

The bounding box is the computationally simplest of all linear bounding containers, and the one most frequently used in many applications.  At runtime, the inequalities do not involve any arithmetic, and only compare raw coordinates with the precomputed min and max constants.

The Bounding Diamond

When one introduces some arithmetic, the simplest nontrivial expressions are those that just add and subtract raw coordinates.  In 2D, we have the expressions p=(x+y) and q=(x-y) which correspond to lines with slopes of (-1) and 1.  For a 2D set S of points, we can compute the minimum and maximum over S of each of these two expressions to get (p_min, p_max, q_min, q_max).  Then, the "bounding diamond" for the set S is the region given by coordinates (x,y) satisfying the inequalities:  p_min £ (x+y) £ p_max and q_min £ (x-y) £ q_max.  Geometrically, it is a rectangle rotated by 45° and resembles a diamond.

The bounding diamond involves a small bit more computation than the bounding box, and a priori they are equivalent approximations of the sets they contain.  However, after computing them both, one may be found to be better than the other.  All the minimums and maximums involved (8 of them) can be computed in O(n) time with a single scan of the set S, and then the areas of the bounding box B and the bounding diamond D can be compared, and the smaller container can be used if one wants.  Since everything is rectangular, we have:

Area(B) = (x_max- x_min)*(y_max- y_min)
Area(D) = (p_max- p_min)*(q_max- q_min)

Further, one could use both containers, the box and the diamond, to get an even smaller combined bounding octagon defined by all 8 inequalities.  In fact, this is probably the preferred usage.  Test for inclusion or intersection with the bounding box first, and then the bounding octagon second.  If one wants to test a point P=(x,y) for inclusion in a polygon W, the only overhead is computing the expressions (x+y) and (x-y) just before testing the diamond inequalities after P was found to be inside the bounding box.  To test disjointness of two polygons, no further expressions need to be evaluated since one only has to test the opposing parallel extremes of the bounding octagon,

In 3D, one can also use maximums and minimums of the four expressions (x ± y ± z) to get 8 inequalities that define a 3D bounding diamond or bounding octahedron.  Specifically, the inequalities are:

(x + y + z)_min £ x + y + z £ (x + y + z)_max
(x + y - z)_min £ x + y - z £ (x + y - z)_max
(x - y + z)_min £ x - y + z £ (x - y + z)_max
(x - y - z)_min £ x - y - z £ (x - y - z)_max

If the coordinates of a point fail to satisfy any of these, then it is outside the bounding octahedron.  Also, two bounding octahedrons are disjoint if the max of an expression for one is less than the min of the same expression for the other.  This gives 8 inequalities to test for disjointness.

One could also combine the 8 octahedron inequalities with the 6 inequalities for the 3D bounding box, and get a set of 14 inequalities which give a 3D bounding cuboctahedron.  Whether it is efficient to test this many inequalities depends on the complexity of the geometric object inside the bounding container.

It should be noted though, that in k-dimensions, with coordinates (x₁,x₂,...,x_k), the k-D bounding box is given by 2k inequalities.  But the k-D bounding diamond is given by 2^k min and max inequalities with the 2^k-1 expressions (x₁±x₂...±x_k).  Thus, this method does not generalize efficiently to high dimensions.

The Convex Hull

The "convex hull" H of a point set S is the smallest convex region containing the points.  Any other convex container for S must contain the convex hull H as a subset.  This means that the convex hull is the bounding container that is the closest and least area approximation for the object it contains.  It is easy to show that H is a polygon in 2D or a polyhedron in 3D with vertices from the set S.  A good way to visualize the 2D convex hull is to imagine that the points of S are pegs stuck in a plane (or nails in a piece of wood), and that H is formed by an elastic band stretched around the outside of all the pegs.  The band contracts to enclose the peg set as tightly as possible.  This is the boundary of the convex hull.  Each edge of the boundary can be expressed as a linear equation, and the half space containing the hull is given by an inequality: (ax+by+c) £ 0 in 2D or (ax+by+cz+d) £ 0 where a, b, c, d are floating point numbers.  The region inside the hull is defined by the collection of all these inequalities.  For example:

The downside of using the convex hull as a container is that it may have a lot of edges, and to check for the inclusion of a point in the hull by testing many independent inequalities takes a lot of computation.  However, there are faster methods to test for inclusion of a point in a convex polygon (see Algorithm 3 about Inclusion of a Point  in a Polygon).  But further, it is not straightforward to test that two different noncongruent convex hulls are disjoint.  This is because the hull does not have opposed parallel faces which can be tested to be outside each other.   So, although the convex hull gives a lower bound for how small a bounding container can be, it is not really practical to use it for this purpose.

We will not pursue this topic further here, and will return to convex hulls some other time.  We just note that there are many algorithms for computing the convex hull in 2D, 3D and higher dimensions.  The fastest algorithms run in O(n log n) time for the 2D or 3D hull of a set of n points.  Many of these algorithms are described by [Preparata & Shamos, 1985, Chaps 3 & 4 ] and [O'Rourke, 1998, Chaps 3 & 4] both of whom have two whole chapters on computing convex hulls.  O'Rourke also gives explicit C code for computing both the 2D and 3D convex hull for a set of points which can be downloaded from his web site  .

The Minimal Rectangle

Nevertheless, it is useful to find minimal area bounding containers which only involve the evaluation of a few expressions to test for inclusion.  The minimal bounding rectangle is such a container.  It is the minimum area (or volume in 3D) rectangle with an arbitrary orientation that contains the vertex set S of a geometric object.  In 2D, it has two pairs of parallel edges, given by the maximum and minimum over S of two linear expressions f₁ = (a₁x + b₁y) and f₂ = (a₂x + b₂y).  In 3D, there are three expressions f_i = (a_ix + b_iy + c_iz) with i=1,3, defining 3 pairs of parallel planar faces.  Without the restriction of rectangularity, one could define the minimal bounding parallelogram (or parallelopied in 3D) which would have an even smaller area than the minimal rectangular box.  For such parallelotope regions, it is relatively easy to determine the inclusion of a point P=(x,y).  In 2D, the inequality tests are:

(a₁x + b₁y)_min £ a₁x + b₁y £ (a₁x + b₁y)_max
(a₂x + b₂y)_min £ a₂x + b₂y £ (a₂x + b₂y)_max

Again, if any one inequality fails, then the point is outside the parallelogram.  Although the more general parallelogram would be a useful container, most work in computational geometry has focused on algorithms for finding the minimal rectangle.  This is partly because these algorithms are slightly more efficient.

In 2D, there is a well-known algorithm for finding the minimum rectangle using a technique referred to as the "rotating calipers" ([Toussaint, 1983], [Pirzadeh, 1999]).  This technique applies to a convex polygon, and can find the minimal rectangle in O(n) time.  For nonconvex objects, one must first compute their convex hull which can be done in O(n log n) time.  The caliper method depends on the observation that one side of the minimal rectangle must coincide with one edge of the convex polygon it contains, as shown in the diagram:

Thus, one only has to consider the orientations given by edges of the convex polygon.  Further, there is an efficient method to determine the minimal rectangle containing a specific edge as one moves successively around the polygon.  This method can be imagined as having two pairs of parallel caliper lines pivoting about opposed vertices of the polygon until one caliper line meets an edge.  Then, the area of the rectangle for that orientation is noted, and the caliper line starts pivoting about the next vertex of that edge.  As a result, the algorithm keeps track of the other sides of the rectangle associated with each edge of the polygon.  So, a single O(n) scan (actually, a 90° rotation is enough) of the polygon edges suffices to find the minimal rectangle.

In 3D, it is considerably more complicated to find the minimum volume bounding box of a convex polyhedron.  The current fastest algorithm is by [O'Rourke, 1985], and it runs in O(n³) time.  But it may be possible to improve on this.  Regardless, since O(n³) algorithms are slow in practice for large point sets, there is considerable interest in fast approximations for the minimum volume box in 3D and higher dimensions.  Recently, [Barequet & Har-Peled, 1999] have described two randomized algorithms that run in O(n log n + 1/ε^4.5) and O(n log n + n/ε³) expected time, where Δ is the error of the approximation.  The deterministic variants of these algorithms run in O(n log² n + 1/ε^4.5) and O(n log² n + n/ε³) time.  Only the second of these algorithms is simple enough to implement in practice.

Quadratic Containers

There are two bounding containers whose boundary is given by a quadratic expression: the bounding ball (a disk or sphere), and the bounding ellipse or ellipsoid.  These are somewhat harder to determine than linear containers, but they are usually more efficient at runtime, especially in higher dimensions.

The Bounding Ball

The "bounding ball" (or "disk" or "sphere") of a geometric object is the smallest circle (in 2D space) or sphere (in 3D and k-D space) containing the object.  It is also called the "minimal spanning sphere" of the object.  In any dimension, the bounding ball of a linear geometric object (defined by its vertex set) is unique, and is specified by a center point C and a radius R.  Given a set S of points, there is a unique bounding ball containing S.

It is easy to test that a point is inside a bounding ball by checking that it is within distance R of the center C.  Also, two balls, say B₁ and B₂, are disjoint if the distance between their centers is greater than the sum of their radii, that is:  d(C₁,C₂) > R₁ + R₂.  So, these basic tests are very simple and efficient no matter what the dimension of the objects involved.  This is not true of linear containers where the number of inequality tests increases with the dimension of the space.  On the other hand, it is more difficult to precompute the bounding ball of a point set since quadratic expressions are involved.

There are several algorithms for exactly computing the minimal ball for a set S of n points.  Some authors have noted that the minimal ball can be derived directly from the "Furthest-point Voronoi Diagram" which can be computed in O(n log³ n) time.  This is a nice theoretical result, but is not easy to implement in practice.  On the other hand, it has been shown ([Welzl, 1991], [de Berg et al, 2000]) that the bounding ball can be computed using a randomized linear programming algorithm that runs in O(n) expected time.  The algorithm is incremental: it starts with a permutation of the point set S={p₁,p₂,...,p_n} and an initial ball B₂ containing the points p₁ and p₂.  It then incrementally adds the other points, constructing increasingly larger balls to contain them as needed.  So, B_k is the bounding ball for {p₁,p₂,...,p_k} and then we consider p_k+1.  If it is in B_k, then B_k+1=B_k.  Otherwise, it is outside B_k and we must expand B_k to a B_k+1 that includes p_k+1.  The trick is to make sure that the B_k+1 we end up with is minimal.  To do this, the incremental algorithm is applied recursively to the set {p_k+1,p₁,...,p_k} with the restriction that p_k+1 must always be on the boundary of all balls constructed.  In the process of doing this, there is yet another level of recursion where balls are constrained to have 2 specific points on their boundary.  After that, newly constructed balls must go through 3 non-collinear points, but in 2D three points uniquely determine a ball, and this determines B_k+1.  So, for a 2D ball, the algorithm has 2 levels of recursion.  In 3D, there are 3 levels to get 4 non-planar points that uniquely determine a 3D sphere. 

Despite sounding complicated, this algorithm is straightforward to implement (see [de Berg et al, 2000]).  Nevertheless, clever heuristics can speed it up significantly in practice.  An excellent fast implementation by [Bernd Gartner, 1999] is available on the web with downloadable source code.  His routine works for 2D, 3D, and higher dimensions up to about 30-D with only about 300 lines of C++ code.  It is numerically stable, and is very fast in low dimensions.

A Fast Approximate Bounding Ball

Another extremely fast O(n) algorithm by [Ritter, 1990] computes a good approximation for the bounding ball.  Although he presents a 3D algorithm, his method works efficiently in any k-dimensional space.  His deterministic incremental algorithm avoids doing any recursion, and just scans the vertex list twice.  For a point set S, it is as follows :

1. First, an initial good guess is made for a bounding ball B.  This is done by finding two points of S that are far from each other, and using the line between them as an initial diameter.  Then, the center of this diameter is the initial ball center, and half the length of the diameter is the initial ball radius.  One finds points of S that are far from each other by selecting ones on opposite extremes of the bounding box for S.
    

2. Next, each point p of S is tested for inclusion in the current ball (by checking that its distance from the center is less than or equal to the radius).  If the next point p_k+1 is in the current B_k, then B_k+1=B_k and one just proceeds to the next point.  But if  p_k+1 is outside B_k, then B_k is expanded just enough to include both itself as well as the point  p_k+1.  This is done by drawing a line from p_k+1 to the current center C_k of B_k and extending it further to intersect the far side of B_k.  This segment is then used as the new diameter for an expanded ball B_k+1.  As shown in the diagram, it clearly contains the prior ball B_k and (thus) all points of  S already considered, and no additional recursion is needed.

Each of the two stages of this algorithm does O(n) efficient computations.  Ritter estimates that this approximation is within 5% of the actual minimal bounding ball.  Below, we give a simplified 2D implementation fastBall().

The Bounding Ellipsoid

As you might expect by now, the "bounding ellipsoid" (or "minimal spanning ellipsoid") is the smallest volume ellipsoid (smallest area ellipse in 2D) that contains a vertex point set S.  Like the minimal rectangle, it can have any orientation.  Thus, it also is a very tight approximation for the object it contains, and is an excellent container.  Also, it is not difficult to test inclusion of a point in an ellipse, especially in 2D where the sum of the distances of a point from the two focal points is a constant on the boundary of the ellipse.

For a set of non-collinear points (non-planar in 3D) the bounding ellipsoid exists and is unique.  [Welzl, 1991] gives a fast randomized algorithm for computing the bounding ellipsoid in O(dd!n) expected time where d=k-1 for k-dimensional space.  A fast implementation has been developed by [Gartner & Schonherr, 1997].  The Welzl algorithm is basically the same as his randomized incremental bounding ball algorithm.  Having constructed the bounding ellipse E_k for the first k points, one then adds the next point p_k+1.  If p_k+1 is inside E_k, then E_k+1=E_k.  Otherwise the algorithm must inductively construct E_k+1.  Again, p_k+1 must be on the boundary of E_k+1, and so the algorithm can recursively construct the minimum ellipse for {p_k+1,p₁,...,p_k} with the restrict that  is on the boundary.  The recursion stops when there are enough constraints to uniquely define an ellipse.  Of course, this is somewhat more difficult to determine than for a spherical ball.  See [Gartner & Schonherr, 1997] for details.

Implementations

// 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 and Vector with
//        coordinates {float x, y;}
//        operators for:
//            Point  = Point ± Vector
//            Vector = Point - Point
//            Vector = Vector ± Vector
//            Vector = Scalar * Vector    (scalar product)
//            Vector = Vector / Scalar    (scalar division)
//    Ball with a center and radius {Point center; float radius;}
//===================================================================

// dot product which allows vector operations in arguments
#define dot(u,v)   ((u).x * (v).x + (u).y * (v).y)
#define norm2(v)   dot(v,v)        // norm2 = squared length of vector
#define norm(v)    sqrt(norm2(v))  // norm = length of vector
#define d(u,v)     norm(u-v)       // distance = norm of difference
 

// fastBall(): a fast approximation of the bounding ball for a point set
//               based on the algorithm given by [Jack Ritter, 1990]
//    Input:  an array V[] of n points
//    Output: a bounding ball = {Point center; float radius;}
void
fastBall( Point V[], int n, Ball* B)
{
    Point C;                           // Center of ball
    float rad, rad2;                   // radius and radius squared
    float xmin, xmax, ymin, ymax;      // bounding box extremes
    int   Pxmin, Pxmax, Pymin, Pymax;  // index of V[] at box extreme

    // find a large diameter to start with
    // first get the bounding box and V[] extreme points for it
    xmin = xmax = V[0].x;
    ymin = ymax = V[0].y;
    Pxmin = Pxmax = Pymin = Pymax = 0;
    for (int i=1; i<n; i++) {
        if (V[i].x < xmin) {
            xmin = V[i].x;
            Pxmin = i;
        }
        else if (V[i].x > xmax) {
            xmax = V[i].x;
            Pxmax = i;
        }
        if (V[i].y < ymin) {
            ymin = V[i].y;
            Pymin = i;
        }
        else if (V[i].y > ymax) {
            ymax = V[i].y;
            Pymax = i;
        }
    }
    // select the largest extent as an initial diameter for the ball
    Vector dVx = V[Pxmax] - V[Pxmin]; // diff of Vx max and min
    Vector dVy = V[Pymax] - V[Pymin]; // diff of Vy max and min
    float dx2 = norm2(dVx); // Vx diff squared
    float dy2 = norm2(dVy); // Vy diff squared
    if (dx2 >= dy2) {                     // x direction is largest extent
        C = V[Pxmin] + (dVx / 2.0);         // Center = midpoint of extremes
        rad2 = norm2(V[Pxmax] - C);         // radius squared
    }
    else {                                // y direction is largest extent
        C = V[Pymin] + (dVy / 2.0);         // Center = midpoint of extremes
        rad2 = norm2(V[Pymax] - C);         // radius squared
    }
    rad = sqrt(rad2);

    // now check that all points V[i] are in the ball
    // and if not, expand the ball just enough to include them
    Vector dV;
    float dist, dist2;
    for (int i=0; i<n; i++) {
        dV = V[i] - C;
        dist2 = norm2(dV);
        if (dist2 <= rad2)    // V[i] is inside the ball already
            continue;
        // V[i] not in ball, so expand ball to include it
        dist = sqrt(dist2);
        rad = (rad + dist) / 2.0;         // enlarge radius just enough
        rad2 = rad * rad;
        C = C + ((dist-rad)/dist) * dV;   // shift Center toward V[i]
    }
    B->Center = C;
    B->radius = rad;
    return;
}

References

Gill Barequet & Sariel Har-Peled, "Efficiently Approximating the Minimum-Volume Bounding Box of a Point Set in 3D", Proc. 10th ACM-SIAM Sympos. Discrete Algorithms (1999), 82-91

Mark de Berg et al, Computational Geometry : Algorithms and Applications,  Section 4.7 "Smallest Enclosing Disks" (2000)

Bernd Gartner, "Smallest Enclosing Balls - Fast and Robust in C++" Web Site (1999)

Bernd Gartner & Sven Schonherr, "Smallest enclosing ellipses - fast and exact", Tech. Report B 97-03, Free Univ. Berlin, Germany (1997)

Joseph O'Rourke, "Finding Minimal Enclosing Boxes", Int'l J. Comp. Info. Sci. 14 (1985), 183-199

Joseph O'Rourke, Computational Geometry in C (2nd Edition)  (1998)

Hormoz Pirzadeh, Rotating Calipers Web Site (1999)

Franco Preparata & Michael Shamos, Computational Geometry: An Introduction (1985)

Jack Ritter, "An Efficient Bounding Sphere" in Graphics Gems (1990)

Godfried Toussaint, "Solving geometric problems with the rotating calipers," Proc. IEEE MELECON'83 (1983)

Emo Welzl, "Smallest enclosing disks (balls and ellipsoids)" in New Results and New Trends in Computer Science, Lecture Notes in Computer Science, Vol. 555 (1991), 359-370

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