Calculating Minimum Volume Bounding Box
In my work I ended up needing to know the minimum volume bounding box (MVBB) of a given of points in 3d. This is also called the oriented bounding box (OBB). Finding the axis-aligned bounding box (AABB) is a simple task, however it can have a volume much larger than needed. Come to find out, finding the minimum volume bounding box is a very difficult problem with several techniques for finding it. However, I didn't find any useful code for finding the MVBB.
One technique is brute force the search by rotating the standard basis for R3 through the full range of angles until you find the minimum. This is ugly, but simple. I didn't try this method as I figured it would be way to slow.
The next technique if ran across used some combination of the eigen vectors of the covariance matrices of the faces, or something like that. I couldn't find any clear information or example code, so I kept searching.
Somewhere I ran across someone saying that the MVBB would share a face with convex hull of the points. So I ran with this as I had all of the convex hull information already. Well, my first attempt worked quite well on my test set of points. However, when I began to try additional random point clouds problems began to show up. Turned out that I made a bad assumption in my code by mistake. So I searched around to figure out how to correct things. In my search I can across a paper that said that it has been proven that the MVBB is not required to share a face with the convex hull. Shame, all the code wasted.
So that paper I found was "Finding Minimal Enclosing Boxes" by Joseph O'Rourke from 1985, in it was a new bit of information, two edges of the convex hull will lie on two adjacent faces (one on each) of the MVBB. Using this I wrote some code that finds the MVBB of a set of points. The algorithm takes in the points of the convex hull and a list of the edges of the hull. The edges are stored as indexes into the points, hence the edge (5, 13) is between the two points and points. Below should be working code (it was pulled from my project so I've cleaned it removing references to the project, but it hasn't been run in its current form). You will need numpy to use the code as its math heavy and doing it all in plain python would have been really slow.
import sys import numpy as np def findOBBEdge(edges, points): # shift the points such that the minimum x, y, z values # in the entire set of points is 0. shift = points.min(axis=0) points = points - shift min_volume = sys.maxint # try every pair of edges (ordering is not important) for idx, edge1_idx in enumerate(edges): e1 = points[edge1_idx] - points[edge1_idx] for idx2, edge2_index in enumerate(edges[(idx+1):]): e2 = points[edge2_index] - points[edge2_index] # transform the two edges into a orthogonal basis w = vec_cross(e1, e2) # returns normalized vector u = vec_cross(w, edge1) v = vec_cross(u, w) # project all the points on to the basis u1, u2 u3 p = calcProjections(points, u1, u2, u3) volume, mins, maxes = calcVolume(p) # we are looking for the minimum volume box if volume <= min_volume: min_volume = volume specs = u, v, w, mins, maxes, volume u, v, w, mins, maxes, volume = specs # get the corner by using our projections, then shift it to move # it back into the same origin as the original set of points corner = u * mins + v * mins + w * mins corner += shift # create the sides which are vectors with the magnitude the length # of that side v1 = u * (maxes - mins) v2 = v * (maxes - mins) v3 = w * (maxes - mins) return corner, v1, v2, v3 def calcVolume(p): """Calculates the volume of the box that would encompass the given points using the given projection. projection is sized (NxM) where N is the number of points and M is the number of vectors they were projected onto. Also return the minimum and maximum bounds along each of those vectors.""" # the minimum and maximum projection of each basis vector mins = p.min(axis=0) maxes = p.max(axis=0) # the volume product of each difference between the maximum and # minimum values from the projection onto each basis vector volume = np.prod(maxes - mins) return volume, mins, maxes def calcProjections(points, *vectors): """Calculates the projection of points (NxD) onto the vectors (MxD) and return the projections p which is a matrix sized (N, M) where N is the number of points and M is the number of vectors. p[i][j], is the projection of points[i] onto vectors[j] (which is between 0 and 1).""" u = np.array(vectors) # project the points onto the vectors into on fell swoop d = np.dot(points, u.T) # this is the dot product of each vector with itself v2 = np.diag(np.inner(u, u)) p = d / v2 return p def vec_cross(u, v): """Return the normalized cross product of u and v.""" w = np.cross(u, v) w /= np.sqrt(np.sum(v**2)) return w