An Algorithm for Computing the Approximate IoU Between Oriented 2D boxes

Notice that when the 2 boxes are too far apart such that there is no intersection, the union area is just the sum of the 2 areas.

When the 2 boxes do have intersections, the union area is the target/ground-truth box area, plus some extra bits from the predicted box that is outside of the target box, i.e. the area of \(A’ \setminus A^*\). This is represented as blue shading in the schematic in Figure 1.

Denote this \(A’ \setminus A^*\) area as \(A_{extra}\). Then we have an expression for the union area \(U\):

It is relatively easy to determine when \(I\) is guaranteed to be 0. We can compare the center distance between the 2 boxes with the sum of their diagonal lengths divided by 2. Even when the test has been too conservative (the 2 boxes have no intersection, when their center distance is smaller than half of the diagonal sum), the approximate algorithm still works.

Now the problem reduces to how to compute \(A_{extra}\).

Notice that the area of \(A_{extra}\) is formed by points along the edges of the blue prediction box, and the edge of the orange target box (see Figure 1).

How far away the individual points along the blue edges can be quantified using SDF to the target box.

However, the conventional definition of SDF to oriented boxes renders contour lines that have rounded corners. See Figure 2a (see also Figure 2a in Signed distance function to oriented 2D boxes). Such rounded corners are not helpful to our task of quantifying \(A_{extra}\).

This is where the modified version – SDF-L1 – comes into the equation. By replacing L2-norm distances with L1-norm, we remove the rounded corners from the distance contours. See Figure 2b.

This L1-norm formulation makes it easy to determine how far away an outside point \(P\) is from the infinite line along which the target box’s right edge is located, if \(P\) is located to the right of the target box.

Similarly, if \(P\) is located above the target, within the span of the target’s width, the distance between \(P\) and the target box is the distance between \(P\) and the infinite line that the top edge is located.

Having got the perpendicular distances between points along the prediction box’s edges, \(A_{extra}\) can then be approximated by sub-dividing the prediction box’s edges and forming a series of trapezoids. The heights of the trapezoids are the SDF-L1 values, and the bases are the small “lateral” steps \(dx\), or \(dy\), depending on whether the SDF-L1 values are horizontal or vertical.

Figure 3 below shows the process of sub-dividing the edge formed by vertices 1 and 2 of the prediction box, into 13 evenly spaced segments.

The SDF-L1 field is drawn as (faded) contours in the background. The 2 vertical grey lines show the left-mid-right segmentation of the SDF-L1 field. Therefore, part of the 1-2 edge falls into the mid section, where SDF-L1 have negative y- values, and the “lateral” steps are \(dx\). And part of the 1-2 edge falls into the right section, where SDF-L1 have positive x- values, and the “lateral” steps take \(dy\). The total area of these trapezoids can be approximated by 2 integrations \(\int SDF_x dy + \int SDF_y dx\).

In Figure 3, we covered part of the \(A_{extra}\) area by sub-dividing the 1-2 edge and integrating the SDF-L1 values. Now do the similar for the 2-3 edge, shown in Figure 4.

Again, part of the SDF-L1 heights are in the mid- section, where SDF values are positive and vertical, and part of the SDF-L1 heights are in the right- section, where SDF values are positive and horizontal. Then repeat for the 3-4 edge, shown in Figure 5.

Notice that SDF-L1 is not defined inside the box, so there are no SDF-L1 heights when the target box cuts into the prediction box, thereby we only account for the area outside of the target box.

Also notice that we are giving the SDF-L1 values signs (red color for positive, blue for negative). Why using signed values for the computation of area?

This is because sometimes there will be regions whose area get integrated more than once. For instance the lower-left triangle shown by blue SDF-L1 heights in Figure 5. This same triangle will also be covered when we take the integration process along the 4-1 edge, shown in Figure 6.

Therefore, we need some kind of mechanism to offset such double-counting. It turns out that by giving the trapezoid heights and bases signed values (effectively making them into vectors), and using cross product to compute signed areas, we can effectively get rid of the double-counted areas as we move around the box edge.

Now refer back to Figure 3–6 and check out the \(+\) or \(-\) signs shown next to the SDF-L1 and \(dx/dy\) vectors. Those are the signs of the cross products between those vector pairs, using the right-hand rule. There is only 1 cross product that has a minus sign, and that appears in the lower-left triangular region when we integrate the 3-4 edge. Notice that this same triangular area is offset by a positive areal integration, when we integrate along the 4-1 edge.

Figure 7 shows the entire integration process after finishing a whole loop around the prediction box, and the total integrated area gives an approximate to \(A_{extra}\). Lastly, the approximate IoU can be computed by substituting \(A_{extra}\) into Equation \eqref{orgfd579eb}, then Equation \eqref{org04b4f0b}.

Lastly, let’s call this area-finding method signed area function (SAF).

The function sdf_obox_l1() for computing SDF-L1 values to an oriented box is given in Signed distance function to oriented 2D boxes.

This is a helper function that converts an oriented box encoded in [x_center, y_center, width, height, angle] format into a 4-polygon:

def obbox2poly(xc, yc, w, h, theta, close=False):
    '''Convert oriented box to 4-polygon, single box
    Order: bottom-right, top-right, top-left, bottom-left'''

    cos, sin = np.cos(theta), np.sin(theta)
    rot = np.array([[cos, -sin], [sin, cos]])
    center = np.c_[xc, yc]
    corner_xy = np.array([
        [w/2, -h/2],
        [w/2,  h/2],
        [-w/2,  h/2],
        [-w/2, -h/2]])
    corner_xy = corner_xy.dot(rot.T)
    corner_xy = center + corner_xy

    if close:
        return np.r_[corner_xy, corner_xy[0:1]]
    else:
        return corner_xy

The following function saf_obox2obox() computes \(A_{extra}\):

def saf_obox2obox(pred_box, gt_box, n_samples=40):
    '''Signed area difference formed between a single pair of predition and target boxes

    Args:
        pred_box (ndarray): predition box in [xc, yc, w, h, angle].
        gt_box (ndarray): reference box in [xc, yc, w, h, angle].
    Keyword Args:
        n_samples (int): number of points to sample along box perimeter.
    Return:
        saf (float): mean sdf difference, i.e. the area of <gt_box> diff <pred_box>.
    '''

    pred_poly = obbox2poly(*pred_box, close=True)
    angle = gt_box[-1]
    rot = np.array([[np.cos(angle), -np.sin(angle)],
                    [np.sin(angle),  np.cos(angle)]])
    saf = 0
    for (x1, y1), (x2, y2) in zip(pred_poly[:-1], pred_poly[1:]):
        pii = np.linspace([x1, y1], [x2, y2], n_samples, True)
        sdfii, sdfxyii = sdf_obox_l1(pii,
            gt_box[0],
            gt_box[1],
            gt_box[2],
            gt_box[3],
            gt_box[4],
            True)
        pii = pii.dot(rot)
        dxyii = np.gradient(pii, axis=0)
        sdfii = np.cross(sdfxyii, dxyii)
        saf += sdfii.sum()

    return saf

A few points to note:

• We take a full loop around the 4 edges of the prediction box, for each edge, the starting point is (x1, y1) and the ending point (x2, y2).
• Then we evenly sample the edge into n_samples intervals.
• For those sampled points, compute their SDF-L1 values using sdf_obox_l1(), and return also the x- and y- components in sdfxyii.
• In the schematics in Figure 3–7, the target box (and therefore the SDF-L1 field) has no rotation. In the general case where the target box is oriented (as in Figure 1), we need to do an counter-rotation so the target box is not oriented. This is done in pii = pii.dot(rot).
• After the rotation, the \(dx\) and \(dy\) values are computed using np.gradient(pii, axis=0).
• Finally, we compute the signed area using np.cross(sdfxyii, dxyii), and integrate that to get the final saf.

Time for some test runs.

To validate the algorithm, I use the shapely module to compute the ground truth union area of 2 boxes, and compare that with the approximate solution. Code below:

from shapely.geometry import Polygon
import numpy as np
import matplotlib.pyplot as plt

figure = plt.figure(figsize=(10, 5))
# example 1
# target box
xc = 2
yc = 1
w = 4
h = 2.5
angle = 30/180*np.pi
gt_box = [xc, yc, w, h, angle]
gt_poly = obbox2poly(*gt_box, False)
gt_poly_sp = Polygon(gt_poly)

# prediction box
box = [4, 3, 4, 5, 145/180*np.pi]
box_poly = obbox2poly(*box, False)
box_poly_sp = Polygon(box_poly)

aa = saf_obox2obox(box, gt_box, 20)
union = gt_poly_sp.union(box_poly_sp).area
union_hat = gt_box[2] * gt_box[3] + aa
print('Union area:', union)
print('Union hat area:', union_hat)

ax = figure.add_subplot(1, 2, 1)
gt_poly = obbox2poly(*gt_box, True)
box_poly = obbox2poly(*box, True)
ax.plot(gt_poly[:, 0], gt_poly[:, 1], 'r-', label='Target')
ax.plot(box_poly[:, 0], box_poly[:, 1], 'b-', label='Prediction')
ax.set_aspect('equal')
ax.set_title('true union: %.3f, approx union: %.3f' % (union, union_hat))
ax.legend()

# example 2
# target box
xc = 2.2
yc = 1
w = 4
h = 2.5
angle = 30/180*np.pi
gt_box = [xc, yc, w, h, angle]
gt_poly = obbox2poly(*gt_box, False)
gt_poly_sp = Polygon(gt_poly)

# prediction box
box = [-2, 3, 4, 5, 45/180*np.pi]
box_poly = obbox2poly(*box, False)
box_poly_sp = Polygon(box_poly)

aa = saf_obox2obox(box, gt_box, 20)
union = gt_poly_sp.union(box_poly_sp).area
union_hat = gt_box[2] * gt_box[3] + aa
print('Union area:', union)
print('Union hat area:', union_hat)

ax = figure.add_subplot(1, 2, 2)
gt_poly = obbox2poly(*gt_box, True)
box_poly = obbox2poly(*box, True)
ax.plot(gt_poly[:, 0], gt_poly[:, 1], 'r-', label='Target')
ax.plot(box_poly[:, 0], box_poly[:, 1], 'b-', label='Prediction')
ax.set_aspect('equal')
ax.set_title('true union: %.3f, approx union: %.3f' % (union, union_hat))
ax.legend()

figure.show()

The output figure:

It is seen that the approximate solution using 20 sub-division steps for each edge is slightly bigger than the ground truth. This is because our integration of the trapezoid areas is not adjusting for the triangular top, but instead using rectangular approximations.

It is also noticed that when the 2 boxes have no intersection (Figure 8b), the method works as well.

In practice, one often needs to deal with multiple pairs of targets and predictions. It would be much more efficient to vectorize the computations.