Error message

  • Notice: Trying to get property of non-object in filter_default_format() (line 532 of /home/ntroutman/webapps/nt_drupal/modules/filter/filter.module).
  • Notice: Undefined variable: options in filter_process_format() (line 911 of /home/ntroutman/webapps/nt_drupal/modules/filter/filter.module).

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[5] and points[13]. 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[0]] - points[edge1_idx[1]]
        for idx2, edge2_index in enumerate(edges[(idx+1):]):            
            e2 = points[edge2_index[0]] - points[edge2_index[1]]
            
            # 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[0] + v * mins[1] + w * mins[2]
    corner += shift
    
    # create the sides which are vectors with the magnitude the length
    # of that side
    v1 = u * (maxes[0] - mins[0])
    v2 = v * (maxes[1] - mins[1])
    v3 = w * (maxes[2] - mins[2])
    
    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

Comments

Hello Nathaniel,

 

I am writing because I too am interested in calculating the minimum bounding volume for a set of points in 3D space.  My MS work is developing spatial statistics in GIS for 3D bones at a fossil site and volume is extremely important in the development.  I haven't been able to find anything on calculating this volume til I came across your site.  Would it be ok if I implemented your python script into my own?  Also, as I am still a newb at python, where would I insert my own set of points into your script?  Any help will be greatly appreciated.

 

Thanks,

Winn

Yes, you can implement my algorithm. To cite it I would request citing my master's thesis: http://www.nathanieltroutman.net/content/masters-thesis

As for getting your own points in. You will have to extract the convex hull, which I did using qhull. Then calling findOBBEdge(edges, points) will return the AABB where edges is the list of edges of the convex hull as a list of tuples which index into the points array.

Perhaps I will writeup an article that goes through the whole process form a list of 3D points to the AABB.

 

 

 

 

Add new comment