Likely, everyone has encountered bounding boxes as an image annotation method. Videos of how Tesla’s vehicles perceive their environment for example have already become quite familiar. I encountered rotated bounding boxes as a feature extraction method. Imagine a microscopy image showing cells. You have already determined for each pixel, whether it belongs to a cell or not. The bounding box can help you extract different features that might be useful for e.g. separating single cells from aggregates or other impurities using the ratio of the bounding box edges or the area.

The question I want to explore in this post is: How do you find this box? Let’s approach the problem from different angles to understand how it is solved algorithmically.

Given a set of points on a 2-dimensional plane, we want to get the smallest rectangle enclosing these points. Rotations are specifically allowed. Right off the bat, we can see that the only points relevant for this task are those on the convex hull of the object: Points “inside” of the structure cannot influence any vertex of the bounding box. Let’s reduce the set of points to the convex hull without going into detail of how this is done. A great explanation can be found here. In a nutshell: There are different methods to achieve this, but you can think of the convex hull as a tuple of edges which encloses all points and which you can traverse by only taking left turns (counter-clockwise traversal). If you are making (in this case) a right turn, you know that you can not be on the convex hull.

A closed-form solution?

To solve a problem, it is often a good idea to start thinking about an associated, simpler problem. Given that x and y values are centered at the origin, the area of the axis-aligned bounding box can be easily found by $(x_{max}-x_{min})*(y_{max}-y_{min})$.

The only step we have to take now is rotate the box and find the angle $\phi$ for which the area is at a minimum. After rotating the box around $90^\circ$, the box is in the starting position and thus must have passed through this minimum. Rotating further, this pattern will repeat, so with the angle $\phi$ given in radians, we know that $\phi \in [0,\frac{\pi}{2})$

The rotation matrix R for rotating a 2D-vector around the origin is

$$R= \begin{pmatrix} cos(\phi)&-sin(\phi)\\ sin(\phi)&cos(\phi)\\ \end{pmatrix}$$

Again, we only need a rotation and no translation because I centered the points at the origin, but we could do a similar thing using an affine transformation matrix, which would include the translation. Let’s put all of this in a function: the area of the bounding box is a function of the angle given a set of points $P$. The matrix of the rotated points $Z$ is

$$Z^T=RP^T= \begin{pmatrix} cos(\phi)&-sin(\phi)\\ sin(\phi)&cos(\phi)\\ \end{pmatrix} \begin{pmatrix} p_{x,1} & p_{x,2} & \cdots & p_{x,n} \\ p_{y,1} & p_{y,2} & \cdots & p_{y,n} \\ \end{pmatrix}$$

and the area is then

$$f(\phi)=(max(Z_{x,i})-min(Z_{x,i}))*(max(Z_{y,i})-min(Z_{y,i}))$$

Let’s have a look at the bounding box during rotation (left) and its associated area (right):

Now visually it is pretty simple to find the rotation angle for which the bounding box is at a minimum, but what about analytically? The minimum $\underset{\phi}{\operatorname{argmin}} f(\phi)$ can naively be found by setting the derivative to zero. But what is the derivative of a minimum?

Since we define the $min$ function to be

$$f(x, y) = \min(x,y) = \begin{cases} x & \text{if } x \le y \\ y & \text{if } x \gt y \end{cases}$$

the gradient would be

$$\nabla_f = \begin{cases} \begin{bmatrix} 1\\ 0\\ \end{bmatrix} & \text{if } x \lt y \\ \begin{bmatrix} 0\\ 1\\ \end{bmatrix} & \text{if } x \gt y \end{cases}$$

but that really does not help us as the function is not differentiable at x=y. The min/max operators can effectively be thought of as if/then statements which are bound to be present in the derivative as well. Thus, it is not defined at those points, but you can see an approximation to the derivative below. The real minimum must be one of the points where a (discontinuous) jump occurs from the negative to positive area in the slope.

Despite the fact that the function is effectively not very complex, we established that there is no useful closed-form solution. But there are much better options than to turn to general numerical optimizers. You may have noticed that rotation always results in a monotonous change of the bounding box area except when we hit an edge. Here is a visualization to show this more clearly:

On the left side you can see the area of the bounding box as a function of its rotation in degrees and on the right side an approximation to the derivative. Every time one of the edges of the bounding boxes is in parallel to the convex hull, I marked this with a red dot. Each of these represents a candidate for the rotation resulting in the minimum bounding box. (Not so) incidentally, those points occur when the bounding box aligns with an edge of the convex hull. From an intuitive perspective that also makes sense. So all that is left to do is find the red point with the lowest area.

The exhaustive approach was employed in the package shotGroups. On github you can find a nicely annotated version of this implementation. This requires the set of points to form a convex hull, because otherwise you would check edges which could not possibly be part of the bounding box. Let’s try it: We only need to check all rotations for which our bounding box aligns with an edge of the convex hull:

In practice, this is how it is done:

for each edge
   for each point on the hull
      project onto edge and its perpendicular
   compute the axis-aligned bounding box in that coordinate system
   save rotation and area of box to candidate list
select smallest box from candidate list

With this approach, we pass through all $n$ edges and then within this loop we have to do another $2$ passes through all other vertices. So regarding the computational complexity, the algorithm should be $\in O(n^2)$. In practice, it appears that the algorithm is sufficient for applications in 2 dimensions, but e.g. in three dimensions and a high number of edges, the computational demand might become noticeable.

Rotating calipers

First of all, the intuition we had before about the box sharing an edge with the convex hull is correct and even more so, there appears to be no way around checking all candidates. However is a more efficient way to check.

Toussaints rotating calipers algorithm also passes through the edges of the convex hull in the order of which they appear to form the polygon, but uses information about the angle between edges for this. You can step through the iterations in the figure below:

initialize the tight box along one support edge (red)
while the calipers have rotated less than 90 degrees
   for each of the four support vertices
      compute the angle from its rectangle side to the next hull edge
   rotate counter-clockwise by the smallest of these angles, reaching the next support edge (green)
   save rotation and area of the updated box to candidate list
select smallest box from candidate list

There are a couple of special cases here, that we won’t go further into. For example the same supporting vertex could appear after rotation, the bounding box might share two edges with the convex hull or multiple rotation angles might result in the minimum box area.

What is important: We only have to do a single pass this way and the computational complexity of this algorithm is $\in O(n)$.

Further points

This was a high-level overview which by no means begins to capture the intricacies implied in the mentioned methods. Finding the convex hull (which we glossed over) was a prerequisite in our case and this in itself can be numerically challenging due to rounding errors. If we are given a convex polygon to start with, the problem is passed on to the bounding box algorithms and can result in the two presented methods yielding slightly different results.

Employing rational arithmetic would solve this problem, but for the edge-aligned basis vectors to be unit-length, normalization has to occur, which is not possible with fractions. To find out how this problem is tackled, check out the great technical paper by David Eberly. Many thanks to Daniel Wollschläger for his help, especially explaining his approach in the package shotGroups.