## Thursday, March 31, 2011

### Fitting Oriented Bounding Boxes

It can be useful to fit oriented bounding boxes (OBBs) to sets of geometric primitives for a number of applications. This post describes one way to do this, based on the PhD thesis of Stefan Gottschalk. I have put together a sample implementation that you can use and adapt as you wish.

The method of fitting is based around the covariance matrix of the items to be fit. The covariance matrix is a generalization of variance to higher dimensions. The use of the covariance matrix in fitting an oriented bounding box is that if the covariance matrix is constructed carefully, its eigenvectors determine the rotation required to obtain a tightly fitting box.

The overall process for fitting a bounding box will then be:
• Compute the mean position and covariance matrix of the inputs
• Extract the eigenvectors of the covariance matrix
• Build a transformation matrix using the mean position and eigenvectors
• Size the box to fit the input

Given a decent linear algebra library (such as gmm++), only the first task is particularly complicated, so I'll discuss it first.

In the case of fitting a bounding box to a collections of objects, there are three dimensions X, Y and Z, so the covariance matrix is 3x3. The covariance matrix for this case is shown below:

$\mathbf{C} = \left[ \begin{array}{ccc} E[(x-E[x])(x-E[x])] & E[(x-E[x])(y-E[y])] & E[(x-E[x])(z-E[z])] \\ E[(y-E[y])(x-E[x])] & E[(y-E[y])(y-E[y])] & E[(y-E[y])(z-E[z])] \\ E[(z-E[z])(x-E[x])] & E[(z-E[z])(y-E[y])] & E[(z-E[z])(z-E[z])] \end{array} \right]$

In the above equation x, y and z are the three global coordinates, and E[arg] is the expected value of the argument arg, which in this case is the mean of arg. If the mean values of x, y and z are denoted by $\hat{x}$, $\hat{y}$ and $\hat{z}$, $\mathbf{C}$ can be written as

$\mathbf{C} = \left[ \begin{array}{ccc} E[x x] - \hat{x} \hat{x} & E[x y] - \hat{x} \hat{y} & E[x z] - \hat{x} \hat{z} \\ E[y x] - \hat{y} \hat{x} & E[y y] - \hat{y} \hat{y} & E[y z] - \hat{y} \hat{z} \\ E[z x] - \hat{z} \hat{x} & E[z y] - \hat{z} \hat{y} & E[z z] - \hat{z} \hat{z} \\ \end{array} \right]$

where the simplification is due to the properties of the expected value. Armed with this definition, we can begin computing covariance matrices for sets of points. This is the simplest case, and can produce decent bounding boxes if the points are well distributed.

We start by computing the mean positions:

$\hat{x}=\frac{1}{N}\sum_{i=1}^{N} x_i \hspace{1cm} \hat{y}=\frac{1}{N}\sum_{i=1}^{N} y_i \hspace{1cm} \hat{z}=\frac{1}{N}\sum_{i=1}^{N} z_i$

where N is the number of points in the point set. This leaves us needing to compute the $E[xx]$, $E[xy]$, ..., $E[zz]$ terms. Since we are dealing with points, these turn out to be very simple:

$E[xx] = \frac{1}{N} \sum_{i=1}^{N} x_i x_i \hspace{1cm} E[xy] = \frac{1}{N} \sum_{i=1}^{N} x_i y_i, \hspace{1cm}..., \hspace{1cm} E[xz] = \frac{1}{N} \sum_{i=1}^{N} x_i z_i$

and the analogous expressions for the other terms. With these expressions, we can now build the covariance matrix $\matbf{C}$ and extract its eigenvectors $\mathbf{v_0}$, $\mathbf{v_1}$ and $\mathbf{v_2}$ using our math library. Each of the three eigenvectors will be a three-component vector, and will be orthogonal to the other two eigenvectors since $\mathbf{C}$ is symmetric.

To build the OBB transformation matrix we compute

$\mathbf{r}=\frac{\mathbf{v_0}}{||\mathbf{v_0}||} \hspace{1cm} \mathbf{u}=\frac{\mathbf{v_1}}{||\mathbf{v_1}||} \hspace{1cm} \mathbf{f}=\frac{\mathbf{v_2}}{||\mathbf{v_2}||}$

and then use them as the columns of the rotation portion of the OBB transformation matrix $\mathbf{T}$, i.e.:

$\mathbf{T} = \left[ \begin{array}{cccc} r_x & u_x & f_x & 0\\ r_y & u_y & f_y & 0\\ r_z & u_z & f_z & 0\\ 0 & 0 & 0 & 1 \end{array} \right]$

All that remains is to figure out how big to make the box and where to put it so that it contains all the points. To do this, we transform each point $\mathbf{p}_i=[x_i, y_i, z_i]^T$ into the local frame of the box using $\mathbf{T}^{-1}$ and keep track of the minimum and maximum coordinates in each direction. Since the rotation part of the matrix is orthogonal and the translation part is zero, a shortcut for transforming the points without having to invert $\mathbf{T}$ is as follows

$\mathbf{T}^{-1} \mathbf{p} = [ \mathbf{r} \cdot \mathbf{p_i}, \mathbf{u} \cdot \mathbf{p_i}, \mathbf{f} \cdot \mathbf{p_i} ]^T$

where $\cdot$ is the vector dot-product and everything else is defined as before. We then find the minimum and maximum of each coordinate of the transformed points for each coordinate, storing them in the 3D points $\mathbf{p_{min}}$ and $\mathbf{p_{max}}$. The half-extents of the bounding box will then be $\mathbf{\delta} = \frac{\mathbf{p_{max}}-\mathbf{p_{min}}}{2}$ and the center position $\mathbf{p_{cen}} =\frac{\mathbf{p_{max}}+\mathbf{p_{min}}}{2}$, both in the local frame of the box. The corresponding world-space coordinate of the center $\mathbf{t}$ can be found by multiplying $\mathbf{p_{cen}}$ by $\mathbf{T}$. Once $\mathbf{t}$ has been computed, it can be added into the transformation matrix to get the final transformation matrix for the OBB:

$\mathbf{T} = \left[ \begin{array}{cccc} r_x & u_x & f_x & t_x\\ r_y & u_y & f_y & t_y\\ r_z & u_z & f_z & t_z\\ 0 & 0 & 0 & 1 \end{array} \right]$

The result is an OBB that contains the input. I have used the Stanford bunny below as a sample model:

Unfortunately, if the points are not well distributed the box may not be a tight fit. The reason for this is that variations in point sampling density can skew the eigenvectors of the covariance matrix and result in a poor orientation being found.

One way to make this process more robust so that tighter bounding boxes can be found is to build the covariance matrix using the triangles of the model rather than the vertices, and integrate over the surface area of the model rather than simply summing the contributions to $\mathbf{C}$ from each vertex.

In this case, the formulas for the terms needed to build $\mathbf{C}$ are slightly more complicated than for points:

$\hat{x} = \frac{1}{A_m} \sum_{i=1}^{N} A_i \hat{x_i} \hspace{1cm} \hat{y} = \frac{1}{A_m} \sum_{i=1}^{N} A_i\hat{y_i} \hspace{1cm} \hat{z} = \frac{1}{A_m} \sum_{i=1}^{N} A_i \hat{z_i}$

where $A_m$ is the surface area of the entire model, $\hat{x_i}$, $\hat{y_i}$ and $\hat{z_i}$ are the coordinates of the centroid of triangle i and $A_i$ is the area of triangle i. The remaining terms $E[xx]$, $E[xy]$ (and analogous), are also a bit more complicated:

$E[xx] = \frac{1}{A_m} \sum_{i=1}^{N} \frac{A_i}{12} \left( 9 \hat{x_i} \hat{x_i} + p_x p_x + q_x q_x + r_x r_x \right)$

$E[xy] = \frac{1}{A_m} \sum_{i=1}^{N} \frac{A_i}{12} \left( 9 \hat{x_i} \hat{y_i} + p_x p_y + q_x q_y + r_x r_y \right)$

where p, q and r are the vertices of triangle i. The remaining expressions for $E[xz]$, $E[yy]$, $E[yz]$ and $E[zz]$ are analogous, and can be found by replacing x and y with the appropriate variables in the above expressions. With these expressions for the covariance matrix of the model, we can then build an OBB in the same manner as before. In most cases, this change will be enough to obtain a tight box. However, if the model has large interior features, these can again mess up the eigenvectors of the covariance matrix and result in a box that doesn't fit the object tightly.

Using triangles instead of points gives a slightly tighter result for the Stanford bunny:

There is one final improvement that can be made which will result in a method that is robust for all inputs. Instead of building the covariance matrix of the input itself we instead build the covariance matrix of the convex hull of the input. The convex hull approximates the shape of the input, and in the process removes interior surfaces that cause errors in fitting tight boxes. Computing the convex hull of a 3D object is somewhat involved, however there are libraries that do it which can be used (for example CGAL). Note however that the convex hull may well have a very non-uniform point distribution, so the triangle-based covariance matrix should be used and not the point-based version.

The convex hull version gives an even tighter box than the point version, but is substantially slower in my implementation (and should be in general) due to the complexity of building the convex hull of the input. Even though I have isotropically remeshed the bunny, there are still significant differences in the OBB orientations for the three methods, even though the point sampling and triangle sizes are relatively uniform.

 OBB fit using points OBB fit using triangles OBB fit using convex hull

I have put together a sample implementation of the three methods of fitting OBBs discussed in this posting. It depends on CGAL for computing convex hulls and gmm++ for linear algebra. Feel free to use it for whatever you'd like.