We are not aware to any previously-published algorithm that solves this problem.

Three-dimensional boxes which enclose sets of points are used for maintaining hierarchical partitioning of sets of points. These data structures have important applications in computer graphics (e.g., for fast rendering of a scene or for collision detection), statistics (for storing and performing range-search queries on a large database of samples), etc. From a top-down viewpoint, the problem in such applications is divided into two (admittedly, related) problems of splitting a given set of points into two (or more) subsets, and of computing a nearly-optimal box (or another generic shape) that encloses each subset. In this paper we concentrate on the second problem.

Numerous heuristics have been proposed for computing a box which encloses a given set of points. The simplest heuristic was naturally to compute the axis-aligned bounding box (AABB) of the point set. Twodimensional variants of this heuristic include the well-known R-tree, the packed R-tree [ RL85 ], the R + -tree [ SRF87 ], the R * -tree [ BKSS90 ], etc. [ HKM95 ] use a minimum-volume AABB trimmed in a fixed number of directions for speeding up collision detection. [ GLM96 ] implement in their RAPID system the OBB-tree (a tree of arbitrarilyoriented bounding boxes), where each box which encloses a set of polygons is aligned with the principal components of the distribution of polygon vertices. A similar idea is used by [ BCG+96 ] for the BOXTREE. The latter work suggests a few variants, in which the computed box is aligned with only one principal component of the point distribution (e.g., the one that corresponds to the smallest or largest eigenvalue of the covariance matrix of the point coordinates), and the other two directions are determined by another method (e.g., by computing

We first present the notation used in this paper.

We denote the origin of coordinates by o and the unit cube by C. Throughout the paper we denote by Q a point set in two dimensions, and use S to denote a similar set in three dimensions. Unless specified otherwise, the set is assumed to contain n (or a comparable number of) points. We denote by CH(Q) (resp., CH(S)) the 2-dimensional (resp., 3-dimensional) convex hull of Q (resp., S). The symbol P is used for denoting a convex polyhedron in three-space.

The symbols R and B are used for rectangles and boxes, respectively We denote by R opt (Q) the minimum-area bounding rectangle of Q and by B opt (S) the minimum-volume bounding box of S.1 Let V be a set of orthogonal vectors in R 3 . We denote by B opt (S, V ) for the minimumvolume bounding box of S whose set of directions contains V . We denote some constant approximation of the minimum-volume bounding box of S (for some predefined positive constant) by B * (S). (Such a box is computed in Lemma 3.4 and Lemma 3.6 -see Section 3 .)

Finally, given a box

In this section we present our main approximation algorithm to the minimum-volume bounding box of a set of points in three dimensions.

First we need a good-enough approximation of the diameter of the point set. In particular, we can approximate in linear time the diameter of a point set in three dimensions by a factor of 1/ √ 3. (In fact, one can find in linear time a (1/ √ 3)-approximation of the diameter in any dimension; see [ EK89 ].)

Actually, we can find any arbitrarily-good approximation of the diameter. Since we were not able to find any reference with a proof of the following folklore lemma (see [ Har99 ] for a similar result), we include an easy proof of it here. Lemma 3.3. Given a set S of n points in R d (for a fixed d) and ε > 0, one can compute in O(n + 1/ε2(d-1) ) time a pair of points s, t ∈ S such that |st| ≥ (1 -ε)D(S).

Proof: Let B be the minimum axis-parallel box containing S, and let G = Grid((ε/(2 √ d))B). For a point x ∈ R d , denote by x G the closest point of G to x. Define S G = {x G | x ∈ S}. Finally, let l be the length of the longest diagonal of a cell of G. Clearly, l ≤ (ε/2)D(S). For every pair of points x, y ∈ S we have

Thus, D(S G ) ≥ D(S) -l. Let s, t be the two points in S for which s G , t G realize the diameter of S G . Now,

The set S G can be computed in O(n) time (where the hidden constant of proportionality contains a factor d). The cardinality of S G is O(1/ε d ). Thus, we can compute the diameter of S G (in a brute force manner) in O(1/ε 2d ) time. The two points of S that correspond to the diameter of S G are, by the above analysis, the sought (1 -ε)-approximation of the diameter of S.

Any point of S G that lies between two other points of S G (along one of the axes) can not correspond to a point of S that defines D(S). By removing all such points in O(n) time we can consider only O(1/ε d-1 ) points. Hence the running time of the algorithm is improved to

The two points of S G that realize its diameter must be vertices of the convex-hull of S G . Thus, the running time of the algorithm of Lemma 3.3 can be further improved by first computing the set of vertices of CH(S G ), denoted by

For higher dimensions we use an output-sensitive algorithm [ Cha96 ]. Let m denote the number of points in S G . Clearly,

Computing the diameter of

time. Overall, the running time of the algorithm of Lemma 3.3 can be improved to O (1/ε)

.

We can do even better in three dimensions if we are willing to sacrifice simplicity. In this case we compute the exact diameter of

The bottleneck of the cited algorithm is the computation of CH(Q); when the latter is given in advance, R opt (Q) can be computed in optimal Θ(n) time. Proof: By using the algorithm of Lemma 3.2 we compute in O(n) time two points s, t ∈ S which form a (1/ √ 3)-approximation of the diameter of S. Let H be a plane perpendicular to st, and let Q be the orthogonal projection of S into H. Now, by using the algorithm of Lemma 3.2 again, we compute in O(n) time two points s ′ , t ′ ∈ Q (see Figure 3.1 ) for which

Let µ be a direction perpendicular to st and s ′ t ′ . We claim that the box B * = B opt (S, {st, s ′ t ′ , µ}) is a (6 √ 6)-approximation of B opt (S). Indeed, let R be the minimum-area bounding rectangle of Q in the directions s ′ t ′ and µ, let ω be the length of the edge of R in the direction µ, and let u, v be the two points of Q lying on the two edges of R parallel to

On the other hand, the quadrilateral

, and its area is

. Let s H be the orthogonal projection of s into H, and let T 1 , . . . , T m be a triangulation of CH(Q) (within H) in which all the triangles share the vertex s H .

We next show that Vol(S) ≥ Area(CH(Q))|st|/3. The outer edge p ′ i q ′ i of T i is a projection of a pair of points p i , q i ∈ S (for i = 1, . . . , m) into H. Thus, each triangle T i corresponds to a tetrahedron T i = CH({p i , q i , s, t}) (for i = 1, . . . , m). These tetrahedra are pairwise disjoint in their interiors, and

where we use the fact that the volume of a tetrahedron does not change when we translate one of its vertices u in a direction parallel to an edge (of the tetrahedron) that does not coincide with u (or, in more generality,

The convex hull of the projection of S contains a large quadrangle parallel to the opposite face of u). We thus have

On the other hand,

Since S ⊆ B * , we have Vol(S) ≤ Vol(B * ) ≤ 6 √ 6 Vol(S), as required.

Note that the approximation factor for the diameter obtained in Lemma 3.4 can be improved by using better approximation schemes, e.g., the algorithm of Lemma 3.3 . In any case a better constant approximation factor will be reflected by a higher constant hidden in the big-Oh notation of the running time.

The algorithm of Lemma 3.4 can be extended to a d-dimensional space by choosing the direction of the exact diameter v of the point set S as one direction of the bounding box, projecting the points to a hyperplane H ∈ R d-1 perpendicular to v, solving recursively the (d -1)-dimensional problem in H, and outputing the Cartesian product of the (d -1)-dimensional solution and v. The volume of the computed box is at most d! times the volume of the optimal (minimum-volume) bounding box of S.3 Moreover, the side lengths of the bounding box are found in decreasing order. [ FL95 ] use a similar method for visualizing a set of points in a high-dimensional space. They compute in a similar manner the first two (or three) directions, project the points into the subspace spanned by these directions, and display the projected points in this subspace.

We next show that a large-enough convex polyhedron contained in the unit cube C must contain a largeenough axis-parallel cube.

Lemma 3.5. Let P be a convex polyhedron with the properties that P ⊆ C and Vol(P ) ≥ 1/15. Then there exists a translation v ∈ R 3 for which 1 107 C + v ⊆ P .

Proof: It is easy to verify that the area of the intersection of C with any plane is at most 3π/4. (Indeed, the intersection area is maximized when the plane passes through the center of C, and such a cross-section is contained in a disk of radius √ 3/2.) Hence,

,

where Width(P ) is the minimum distance between two parallel planes supporting P . Let r(P ) be the radius of the largest ball K inscribed in P . It is known [ GK92 ] that r(P ) ≥ Width(P )/(2 √ 3). This implies r(P ) ≥ 2/(45 √ 3π). Finally, K inscribes an axis-parallel cube C ′ whose side is of length (2/ √ 3)r(P ) ≥ 1/107, and C ′ ⊆ K ⊆ P , as asserted.

Let A and B be two sets in R 3 . The Minkowski sum of A and B is the set

The idea of our main algorithm is approximating CH(S) by a low-complexity convex polyhedron P ⊇ CH(S), followed by computing (exactly) B opt (P ). The polyhedron P is chosen such that B opt (P ) is an (1 + ε)approximation of B opt (S).

First we need a combined version of Lemma 3.4 and Lemma 3. Note that this lemma is true even if S is a planar set, in which case the minimum-area bounding rectangle provided by Lemma 3.4 degenerates to a segment, and the volume of the computed bounding box is 0. This guarantees that the approximation algorithm (described below) produces a degenerate box (of volume 0) for a degenerate (planar) input point set.

We are now ready to present the approximation algorithm for B opt (S). Let B = B * (S) be the bounding box of S computed by Lemma 3.6 , and let B ε be a translated copy of ε 428 B centered at o. In addition, define S ⊕ = CH(S) ⊕ B ε and G = Grid( 1 2 B ε ). We approximate S on G. For each point p ∈ S let G(p) be the set of eight vertices of the cell of G that contains p, and let S G = ∪ p∈S G(p). Define P = CH(S G ). Clearly, CH(S) ⊆ P ⊆ S ⊕ . Moreover, one can compute P in O(n + (1/ε 2 ) log (1/ε)) time. On the other hand, P ⊆ B ⊕ B ε . The latter term is a box which contains at most k = 860/ε + 4 grid points along each of the directions set by B, so k is also an upper bound for the number of grid points contained by P in each direction. [ And63 ] showed that the complexity of P is O(k 3/2 ) = O(1/ε 3/2 ). We exploit this result in the analysis of the running time of the algorithm. Finally, we compute B opt (P ) exactly.

It remains to show that B opt (P ) is a (1 + ε)-approximation of B opt (S). Let B ε opt be a translation of ε 4 B opt (S) that contains B ε . (The existence of B ε opt is guaranteed by Lemma 3.6 .) Thus, tells us that the constant c may be chosen small enough so as to make the interior of every cell of M contain at least one grid point of G.

Is it fairly easy to prove that:

Lemma 4.1. Let R and H be a rectangle and a plane, respectively, in R 3 . There exists a rectangle R ′ ⊂ H that contains the orthogonal projection of R into H, such that Area(R ′ ) ≤ Area(R).

Our next goal is to show that there exist two G-grid points, where the minimum-volume bounding box of S perpendicular to the direction defined by these points is a good approximation of B opt (S). (Note that in the grid-search algorithm the first point is actually o.) Lemma 4.2. Let p and q be points of G ∩ B M (0,0,5) and G ∩ B M (0,0,20m+15) , respectively. Then B opt (S, {pq}) is a (1 + ε)-approximation of B opt (S).

Proof: Refer to Figure 4.2 . Denote by H p and H q the two planes that pass through p and q, respectively, and are perpendicular to pq. Also denote the angle between pq and the z + axis by α. By construction, the direction pq is almost vertical: A simple calculation shows that tan(α) ≤ √ 2/(10m). Let H xy p be the xy-parallel plane that passes through p, let r ∈ H p be an arbitrary point for which |pr| ≤ 2|pq|, and let r ′ be the orthogonal projection of r into H xy p . Then,

where we use the fact that 3/(10m) ≤ 4/(10m + 20) if and only if m ≥ 6, which is guaranteed. This implies that H p cannot intersect B opt (S) since the vertical distance between H p and H xy p (for such a point r) is less than 4δ|b 1 | (four cells of M), whereas the vertical distance between H xy p and the bottom of B opt (S) is at least 4δ|b 1 |. Similarly, the plane H q cannot intersect the box B opt (S). The same argument shows that H p and H q intersect only the vertical edges of B but do not intersect its top and bottom faces.

Let A be the orthogonal projection of B opt (S) into H p (along pq), and let V be the product of A and pq (that is, a prism whose base is A). It is easy to verify that B opt (S) ⊆ V ⊆ B. Obviously, the rectangle R = B ∩ H xy p contains the set V ∩H xy p . Moreover, according to Lemma 4.1 the orthogonal projection of R into

Let B be the box with R ′ as its base and the opposite face lying on H q . Clearly, S ⊆ B opt (S) ⊆ V ⊆ B. Finally, we have

Since Vol( B) ≥ Vol(B opt (S, {pq})), the latter box B opt (S, {pq}) is also a (1 + ε)-approximation of B opt (S).

We complete this discussion by showing that ⃗ pq is indeed “short,” namely, that |pq| ∞ (measured by G-grid units) is small. Intuitively, this follows from the fact that the grid sizes of M and G are comparable up to a multiplicative constant. By Lemma 3.5 we have that a copy of the unit box of M (that is, δB), scaled down by a constant factor and translated, is contained by the unit box of G (that is, cεB * (S)). In particular, this implies that every grid box of M is covered by a constant number of grid boxes of G. In addition, the segment pq is contained in ∪ 20m+15 i=5 B M (0,0,i) , whose height (size along z) is O(1/ε). Thus, the segment pq can be covered by O(1/ε) grid boxes of G. Hence, all the coordinates of q -p (in G units) are O(1/ε), where the constant of proportionality hidden in the big-Oh notation is c (the constant used by the algorithm). We thus establish the following theorem: Theorem 4.3. Let S be a set of n points in R 3 , and let 0 < ε ≤ 1 be a parameter. One can compute in O(n log n + n/ε 3 ) time a bounding box B(S) with Vol(B(S)) ≤ (1 + ε) Vol(B opt (S)).

The algorithm described above may be too slow to use in practice because the constant of proportionality hidden in the big-Oh notation (affected by the value of c) may be too large. However, it suggests the heuristic of computing the bounding boxes B opt (S, {v}) induced by directions defined by grid points of G “close” to o. Theorem 4.3 implies that the higher the bound on the “length” of v is, the better the approximation is.

We have implemented software that computes the exact 2-dimensional minimum-area bounding rectangle of a planar point set (by computing its convex hull and then applying a rotating-calipers algorithm). Based on this tool we have implemented software that computes the exact 3-dimensional minimum-volume bounding box of a spatial point set, one of whose directions is given. We have used the latter tool for implementing several approximation heuristics for the minimum-volume bounding box, and report here on several of them. The entire software was implemented in plain C and it runs under any Unix-like operating system. The running times reported here were measured in a Linux environment on a 200-MHz Pentium-Pro machine. The software consists of about 1,500 lines of code.

It was easy to observe that for any bounding box B(S), one can always “locally” improve (decrease) the volume of the box by projecting it into a plane perpendicular to one of the directions of B(S), followed by computing the minimum-area bounding rectangle of the projected set in that plane, and by using this rectangle as the base of an improving bounding box of S. Our experimental results revealed many examples in which this procedure converges to a local (but not a global) minimum-volume bounding box. Nevertheless, this procedure improves (by a small amount) the solutions produced by the approximation algorithms described in this paper. We performed each experiment twice, without and with this solution-improving step.

Here are three examples (out of the many examples we experimented with) of the performance of the heuristics. Figure 5.1 (a.1) shows a rotated version of the set S = {(-1, -0.1, 0), (-1, 0.1, 0), (1, 0, -0.1), (1, 0, 0.1)}.

(The points are displayed as small triangles.) Figure 5.1 (a.2) and Figure 5.1 (a.3) show B * (S) (which is also optimal among all boxes aligned with some diagonal of CH(S)), and the improved B * (S)-G(20) (see below), respectively. box of S among those aligned with some hull diagonal. Figure 5.1 (c) shows the same types of boxes bounding a set of 100 points which were randomly and uniformly selected on the unit sphere. Table 5.1 shows the volumes of several bounding boxes of the three spatial sets, and the corresponding running times of our software. The box B * (S) is the minimum-volume bounding box aligned with the diameter of the set S. The “all-pairs” box is obtained by minimizing the volume of all the boxes aligned with directions which connect some two points of S. The suffix “-G(k)” stands for checking all the boxes aligned with directions obtained by connecting the origin o with a grid point whose L ∞ norm is at most k. The table reports results for Grid(B * (S)) and for the regular Cartesian grid. The column entitled “MVBB(v)” details the number of calls to the function that computes B opt (S, {v}).

Note that although the more laborious heuristics require much time to run, some of the faster heuristics perform reasonably well in practice. For example, the B * (S)-G(5) variant combined with the improvement step produces good approximating boxes and runs fast. Moreover, this heuristic performs considerably better than the uniform-grid heuristic for “long and skinny” sets of points.

In this paper we present an efficient algorithm for approximating the minimum-volume bounding box of a point set in R 3 . We also present a simpler algorithm which we implemented and experimented with on numerous three-dimensional point sets.

Jeff Erickson pointed out in a personal communication that it was possible to reduce the O(1/ε 4.5 ) term in the running time of the first algorithm to O(1/ε 3 ) (in the cost of adding an O(n log n) term) by using Dudley’s method [ Dud74 ,AHSV97 ]. The main idea is to scale down the space so as to transform B * (S) to a unit cube. There one computes a (cε)-approximation (for a suitable constant c) of CH(S) and scales up this approximating polyhedron back to the original space. Similarly to our argumentation in Section 3 , one can show that the minimum-volume bounding box of the scaled-up polyhedron is a (1 + ε)-approximation of the minimum-volume bounding box of S. This version of the algorithm requires O(n log n + 1/ε 3 ) time.

The implementation of the algorithm described in this paper is available online [ Har00 ]. We conclude by mentioning two open problems:

• Can one maintain dynamically a (1 + ε)-approximation of the minimum-volume bounding box of a moving point set in R 3 ?

• Can one compute efficiently and by a simple algorithm a (1 + ε)-approximation of the minimum-volume bounding ellipsoid of a point set in R 3 ?

Since we are interested in the area (or volume) of the object, we may pick an arbitrary minimum rectangle (or box) if several exist.

Note that we approximate the diameter from below, while we approximate the minimum-volume bounding box from above.

Here we can trade time for approximation quality. By investing less time we can compute a (1/ √ 2)-(resp., (1/ √ 3)-) approximation of the diameter in two (resp., k, for 3 ≤ k ≤ d) dimensions, and obtain a ( √ 2 • 3 (d-2)/2 d!)-approximation of B opt (S). By investing even less time we can compute a (1/ √ k)-approximation of the diameter in k dimensions (for 2 ≤ k ≤ d), and obtain a (d!) 3/2 -approximation of B opt (S).