Approximate area of overlap of two rotated rectangles

Question

I need to estimate the overlap ratio of two rectangles, each one with arbitrary size and orientation.

I know how to perform the exact computation, using the Sutherland-Hodgman algorithm, which can be optimized for this case.

Anyway as I need to use that function intensively and perfect accuracy isn't required (say 10% error can be tolerated), I was wondering if it cannot be evaluated in a faster way.

If that helps, one can assume the same aspect ratio for both rectangles, and ratio of the areas not exceeding $4$.

Answer 1

Here is an untested idea:

Take the thinnest rectangle (in red). Intersect one of its long edges with the other rectangle (in blue). Estimate the overlap using the subrectangle (in green) of the first rectangle defined by the two interaction points. If there is only one intersection point, use the vertex that is inside the other rectangle as the second point.

Answer 2

This is not an answer, but may be of use to others as a starting point for fast exact solutions or an approximation.

Edited to add: Using this approach the exact solution is pretty cheap, computationally: a few dozen multiplications and additions, and a handful of divisions, plus a dozen or so comparisons.

---

Let's say each rectangle is defined by one vertex and two edge vectors; $\vec{C}$, $\vec{U}$, and $\vec{V}$ for one, and $\vec{c}$, $\vec{u}$, and $\vec{v}$ for the other.

If the two are rectangles, then the edge vectors are perpendicular, $\vec{U}\cdot\vec{V} = 0$, and $\vec{u} \cdot \vec{v} = 0$. Note that the rectangles areas are $\vec{U}\times\vec{V}$ and $\vec{u}\times\vec{v}$, respectively, where $\times$ is the two-dimensional analog of vector cross product ($x_1 y_2 - x_2 y_1$).

If we transform the first rectangle to the unit square at $(0,0)-(1,1)$, the second rectangle becomes a parallelogram, with the vertex at $$\vec{c}' = \left( \frac{\vec{V} \times \left( \vec{C} - \vec{c} \right) }{\vec{U} \times \vec{V}} ,\; \frac{ \vec{U} \times \left( \vec{c} - \vec{C}\right)}{\vec{U} \times \vec{V} } \right)$$ and the two edge vectors becoming $$\vec{u}' = \left( \frac{\vec{u} \times \vec{V}}{\vec{U} \times \vec{V}} ,\; \frac{\vec{U} \times \vec{u}}{\vec{U} \times \vec{V}} \right)$$ and $$\vec{v}' = \left( \frac{\vec{v} \times \vec{V}}{\vec{U} \times \vec{V}} ,\; \frac{\vec{U} \times \vec{v}}{\vec{U} \times \vec{V}} \right)$$ A straightforward implementation using scalars (and Cartesian coordinates) uses 14 multiplications, 6 divisions, and 9 subtractions (or, alternatively, 20 multiplications, one reciprocal (or division), and 9 subtractions, if you replace the divisions by the slightly less precise multiplication by reciprocal).

Although the other rectangle is now a parallelogram, we can reflect the parallelogram around $x = 0.5$, $y = 0.5$, and $x = y$, in any combination, due to the symmetries.

If we find the fraction $f$ of the unit square $(0,0)-(1,0)$ covered by the parallelogram, the area of the intersection is the area of the first rectangle multiplied by that fraction, i.e. $\left(\vec{U}\times\vec{V}\right) f$.

(The area of this parallelogram, $\left({\vec{u}'} \times {\vec{v}'}\right)$, is not interesting at all.)

Note that if the parallelogram vertices have all $x \le 0$, or all $y \le 0$, or all $x \ge 1$, or all $y \ge 1$, the parallelogram does not intersect with the unit square $(0,0)-(1,1)$ at all. This test is cheap, and rules out most (but not all) non-intersecting parallelograms.

---

If we label the transformed vertices of the parallelogram, $$\begin{array}{ll} \vec{p}_0 = \vec{c}' & = ( x_0 , y_0 )\\ \vec{p}_1 = \vec{c}' + \vec{u}' & = ( x_1 , y_1 )\\ \vec{p}_2 = \vec{c}' + \vec{u}' + \vec{v}' & = ( x_2 , y_2 ) \\ \vec{p}_3 = \vec{c}' + \vec{v}' & = ( x_3 , y_3 ) \end{array}$$ with $$\begin{array}{l} \vec{c}' = ( x_c , y_c ) \\ \vec{u}' = ( x_u , y_u ) \\ \vec{v}' = ( x_v , y_v ) \end{array}$$ it turns out that we can find the intersections between the parallelogram and the unit square $(0,0)-(1,1)$ very cheaply: $$\begin{array}{l|l|l} \text{edge} & \text{intersection} & \text{solution} \\ \hline \vec{p}_0 - \vec{p}_1 & x = 0 & y = y_c - x_c \frac{y_u}{x_u} \\ \; & x = 1 & y = y_c - (x_c - 1) \frac{y_u}{x_u} \\ \; & y = 0 & x = x_c - y_c \frac{x_u}{y_u} \\ \; & y = 1 & x = x_c - (y_c - 1) \frac{x_u}{y_u} \\ \hline \vec{p}_1 - \vec{p}_2 & x = 0 & y = y_c + y_u - (x_c + x_u) \frac{y_v}{x_v} \\ \; & x = 1 & y = y_c + y_u - (x_c + x_u - 1) \frac{y_v}{x_v} \\ \; & y = 0 & x = x_c + x_u - (y_c + y_u) \frac{x_v}{y_v} \\ \; & y = 1 & x = x_c + x_u - (y_c + y_u - 1) \frac{x_v}{y_v} \\ \hline \vec{p}_2 - \vec{p}_3 & x = 0 & y = y_c + y_v - (x_c + x_v)\frac{y_u}{x_u} \\ \; & x = 1 & y = y_c + y_v - (x_c + x_v - 1)\frac{y_u}{x_u} \\ \; & y = 0 & x = x_c + x_v - (y_c + y_v)\frac{x_u}{y_u} \\ \; & y = 1 & x = x_c + x_v - (y_c + y_v - 1)\frac{x_u}{y_u} \\ \hline \vec{p}_3 - \vec{p}_0 & x = 0 & y = y_c - x_c \frac{y_v}{x_v} \\ \; & x = 1 & y = y_c - (x_c - 1) \frac{y_v}{x_v} \\ \; & y = 0 & x = x_c - y_c \frac{x_v}{y_v} \\ \; & y = 1 & x = x_c - (y_c - 1) \frac{x_v}{y_v} \\ \end{array}$$

In other words, there are only four unique multipliers: $\frac{x_u}{y_u}$, $\frac{y_u}{x_u}$, $\frac{x_v}{y_v}$, and $\frac{y_v}{x_v}$; the $1$-intersections are found by adding the multiplier to the $0$-intersections.

Thus, only four divisions, eight multiplications, and twelve additions or subtractions are needed to find the intersection coordinates. (You'll need to avoid division by zero, so add eight comparisons: $\ge -\epsilon$ and $\le \epsilon$ per divisor, where $\epsilon$ depends on the numerical precision.)

Of course, this does not yield a fast approximation; but this does yield a sufficiently fast "exact" intersection area -- about 22 multiplications, 5 divisions, and 21 additions or subtractions; and about a dozen or so comparisons or checks -- that I do not believe an approximation is worth pursuing.

(Of course, calculating the area of the clipped convex polygon will then require a few 2D analogs of vector cross products, two multiplications and a subtraction each. I don't see that as a much of an additional cost. In practice, one should further consider the comparisons/checks, to see if some of them can be avoided (using eg. masking), as conditional jumps tend to be costly on current architectures.)

---

If additions and comparisons are very cheap (relative to multiplications or divisions), then one possibility to estimate the area of the intersection, is to transform the larger rectangle to unit square $(0,0)-(1,1)$, and use a simple nested loop to count the number of points inside the unit square, using a regular grid over the parallelogram. In pseudocode:

Function unit_fraction(o.x, o.y, u.x, u.y, v.x, v.y, nu, nv):
    du.x = u.x / (nu-1)
    du.y = u.y / (nu-1)
    dv.x = u.x / (nv-1)
    dv.y = u.y / (nv-1)
    p.x = o.x
    p.y = o.y
    h = 0
    For i = 1..nv:
       c.x = p.x
       c.y = p.y
       For j = 1..nu:
           If (c.x >= 0) And (c.x <= 1) And
              (c.y >= 0) And (c.y <= 1):
               h = h + 1
           End If
       End For
       c.x = c.x + dv.x
       c.y = c.y + dv.y
    End For
    Return h / (nu * nv)
End Function

The absolute error in the intersection area can exceed the area of the smaller rectangle divided by nu and nv, if the sample points fall just outside or just inside the intersection, but I'm too lazy to try to calculate the bounds for the absolute error based on nu and nv.

On typical architectures, this is not that feasible, because multiplications and divisions are not that expensive compared to additions and subtractions.