ellipse_intersection.rst

Ellipse Intersections

This section details the implementation of ellipse intersections in Sara. This section is extracted from the appendix of my thesis Ok:2013:phdthesis and has been retouched a bit.

Eberly provides a comprehensive study on the computation of ellipses intersection, namely the computation of its area and its intersection points. This is non a trivial geometric problem. We complement Eberly’s study with additional technical details about the area computation of two intersecting ellipses.

Origin-Centered Axis-Aligned Ellipses

Let ℰ be a ellipse with semi-major axis a and semi-minor axis b, i.e.,  a ≥ b > 0. Let us first suppose that ℰ is centered at the origin and is axis-aligned oriented, i.e., such that the axis a is along the x-axis and the axis b along the y-axis. Then,

Important

the equation of ellipse ℰ is

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$

Ellipse Area

Riemann sum approximating the upper quadrant area of the ellipse.

Using the symmetry in the ellipse, the area of ellipse ℰ is 4 times the upper quadrant area of the ellipse, i.e.,

$$\mathop{\mathrm{area}}(\mathcal{E}) = 4 \int_{0}^{a} y(x) \mathop{\mathrm{d}x} = 4 b \int_{0}^{a} \sqrt{1 - \frac{x^2}{a^2}} \mathop{\mathrm{d}x} = \pi a b$$

The integral is the limit of the Riemann sum as illustrated in Figure [figriemannsum].

Let us detail the computation. We use the 𝒞¹-diffeomorphism change of variable $\frac{x}{a} = \sin \theta$ which is valid for [0, a] → [0, π/2]. We recall that a 𝒞¹-diffeomorphism is an invertible differentiable function with continuous derivative.

The differential is dx  = acos (θ)dθ  and hence,

Important

.. math:

\begin{aligned}
\mathop{\mathrm{area}}(\mathcal{E})

&= 4ab int{0}^{pi/2} cos^2(theta) mathop{mathrm{d}theta} \ &= 4ab int{0}^{pi/2} frac{1 + cos(2theta)}{2} mathop{mathrm{d}theta} \ &= 4ab left[ frac{theta}{2} + frac{sin(2theta)}{4} right]_{0}^{pi/2} \ &= pi a b.

end{aligned}

Area of an Elliptical Sector

In this part, we review the computation of the area of an ellipse sector. We complement Eberly with a bit more details.

The ellipse sector delimited by the polar angles (θ₁, θ₂) is colored in blue

The elliptic sector area is delimited in polar coordinates by [θ₁, θ₂] (with θ₁ < θ₂) as illustrated in Figure [fig-ellsector]. Using polar coordinates, it equals to the following nonnegative integral

$$A(\theta_1, \theta_2) = \frac{1}{2} \int_{\theta_1}^{\theta_2} r^2 \mathop{\mathrm{d}\theta}.$$

The change of variable in polar coordinates is x = rcos θ and y = rsin θ and, thus with Equation eq-elleq, $\displaystyle\frac{r^2 \cos^2(\theta)}{a^2} + \frac{r^2 \sin^2(\theta)}{b^2} = 1$, therefore

$$\displaystyle r^2 = \frac{a^2 b^2}{b^2 \cos^2(\theta) + a^2 \sin^2(\theta)}.$$

Plugging the formula of r in the integral,

$$A(\theta_1, \theta_2) = \frac{a^2b^2}{2} \int_{\theta_0}^{\theta_1} \frac{\mathop{\mathrm{d}\theta}}{b^2 \cos^2(\theta) + a^2 \sin^2(\theta)}$$

Now the integrand $\frac{\mathop{\mathrm{d}\theta}}{b^2 \cos^2(\theta) + a^2 \sin^2(\theta)}$ is invariant by the transformation θ ↦ θ + π, i.e.,

$$\frac{\mathop{\mathrm{d}\theta}} {b^2 \cos^2(\theta) + a^2 \sin^2(\theta)} = \frac{\mathop{\mathrm{d}(\theta+\pi)}} {b^2 \cos^2(\theta+\pi) + a^2 \sin^2(\theta+\pi)}.$$

According to Bioche’s rule, a relevant change of variable is the 𝒞¹-diffeomorphism change of variable t = tan (θ) which is valid for ] − π/2, π/2[ → ] − ∞, ∞[. Let us first rewrite

$$\begin{aligned} \begin{aligned} A(\theta_1, \theta_2) &= \frac{a^2b^2}{2} \int_{\theta_1}^{\theta_2} \frac{\mathop{\mathrm{d}\theta}}{b^2 \cos^2(\theta) + a^2 \sin^2(\theta)}\\\ &= \frac{a^2b^2}{2} \int_{\theta_1}^{\theta_2} \frac{\frac{\mathop{\mathrm{d}\theta}}{\cos^2(\theta)}}{b^2 + a^2 \tan^2(\theta)}\\\ &= \frac{\cancel{a^2}b^2}{2} \int_{\theta_1}^{\theta_2} \frac{\frac{\mathop{\mathrm{d}\theta}}{\cos^2(\theta)}}{\cancel{a^2} (b/a)^2 + \tan^2(\theta))}\\\ \end{aligned} \end{aligned}$$

Differentiating t = tan θ, $\mathop{\mathrm{d}t} = \frac{\mathop{\mathrm{d}\theta}}{\cos^2(\theta)}$, thus

$$\begin{aligned} \begin{aligned} A(\theta_1, \theta_2) &= \frac{b^2}{2} \int_{\tan\theta_1}^{\tan\theta_2} \frac{\mathop{\mathrm{d}t}}{(b/a)^2 + t^2}\\\ &= \frac{b^{\cancel{2}}}{2} \left[ \frac{a}{\cancel{b}} \arctan\left(\frac{a}{b} t\right) \right]_{\tan\theta_1}^{\tan\theta_2}\\\ &= \frac{ab}{2} \left[ \arctan\left(\frac{a}{b} t\right) \right]_{\tan\theta_1}^{\tan\theta_2} \\\ &= \frac{ab}{2} \left( \arctan\left(\frac{a}{b} \tan\theta_2\right) - \arctan\left(\frac{a}{b} \tan\theta_1\right) \right)\end{aligned} \end{aligned}$$

Hence,

$$A(\theta_1, \theta_2) = \frac{ab}{2} \left( \arctan\left(\frac{a}{b} \tan\theta_2\right) - \arctan\left(\frac{a}{b} \tan\theta_1\right) \right)$$

Warning

The integral is properly defined for (θ₁, θ₂) ∈ ] − π/2, π/2[. But, using symmetry properties of the ellipse, we can easily retrieve the elliptical sector for any (θ₁, θ₂) ∈ ] − π, π[.

Alternatively, Eberly provides a more convenient antiderivative because it is defined in ] − π, π] as follows

$$F(\theta) = \frac{ab}{2} \left[ \theta - \arctan \left( \frac{(b-a) \sin 2\theta}{(b+a) + (b-a)\cos 2 \theta} \right) \right].$$

Hence, the elliptic sector area equals to the following nonnegative quantity

Important

.. math:

\forall (\theta_1, \theta_2) \in ]-\pi, \pi], \ A(\theta_1, \theta_2) =
\left| F(\theta_2) - F(\theta_1) \right|.

Area Bounded by a Line Segment and an Elliptical Arc

The ellipse sector bounded by a line segment and the elliptical arc (θ₁, θ₂) is colored in blue.

We are interested in computing the elliptic portion by a line segment and the elliptical arc (θ₁, θ₂) such that

|θ₂ − θ₁| ≤ π

This condition is important as a such elliptic portion always corresponds to the blue elliptic portion in Figure [figellsector2]. Let us denote the area of such portion by B(θ₁, θ₂). Geometrically, we see that, if |θ₂ − θ₁| ≤ π, then

$$\begin{aligned} \begin{aligned} B(\theta_1, \theta_2) &= \mathop{\mathrm{area}}(\mathrm{sector(\theta_1, \theta_2)}) - \mathop{\mathrm{area}}(\mathrm{triangle(\theta_1, \theta_2)})\\\ &= A(\theta_1, \theta_2) - \frac{1}{2} |x_2y_1 - x_1y_2|\end{aligned} \end{aligned}$$

where (x_i, y_i) = (r_icos θ_i, r_isin θ_i) and $\displaystyle r_i = \frac{ab}{\sqrt{b^2 \cos^2(\theta_i)+a^2 \sin^2(\theta_i)}}$ for i ∈ {1, 2}.

Note that the other portion corresponding to the red one in Figure 3 has an area which equals to πab − B(θ₁, θ₂) ≥ B(θ₁, θ₂) if |θ₂ − θ₁| ≤ π.

To summarize, our portion of interest, illustrated by the blue elliptic portion in Figure 3, has an area which equals to

Important

For any (θ₁, θ₂) ∈ ] − π, π],

$$\begin{aligned} \ B(\theta_1, \theta_2) = \left\{ \begin{array}{cl} \displaystyle A(\theta_1, \theta_2) - \frac{1}{2} |x_2y_1 - x_1y_2| & \textrm{if} \ |\theta_2 - \theta_1| \leq \pi \\\ \displaystyle \pi a b - A(\theta_1, \theta_2) + \frac{1}{2} |x_2y_1 - x_1y_2| & \textrm{otherwise} \end{array} \right. . \end{aligned}$$

General Ellipse Parameterization

The previous sections has provided the basis for area of intersecting ellipses. However, ellipses are neither centered at the origin nor aligned with the axes of the reference frame in general. Therefore, an ellipse ℰ is entirely defined by the following geometric information

• a center x_ℰ,
• axis radii (a_ℰ, b_ℰ),
• an orientation θ_ℰ, i.e., the oriented angle between the x-axis and the axis of radius a_ℰ.

or more concisely by the pair (x_ℰ, Σ_ℰ) where the positive definite matrix Σ_ℰ ∈ 𝒮₂^+ + is such that

Σ_ℰ = R_ℰD_ℰR_ℰ^T

where R_ℰ is a rotation matrix defined as

$$\begin{aligned} \mathbf{R}_{\mathcal{E}} \overset{\textrm{def}}{=} \begin{bmatrix} \cos\theta_{\mathcal{E}} & -\sin\theta_{\mathcal{E}}\\\ \sin\theta_{\mathcal{E}} & \cos\theta_{\mathcal{E}} \end{bmatrix} \end{aligned}$$

and D_ℰ is the diagonal matrix defined as

$$\begin{aligned} \mathbf{D}_{\mathcal{E}} \overset{\textrm{def}}{=} \begin{bmatrix} 1/b_{\mathcal{E}}^2 & 0\\\ 0 & 1/a_{\mathcal{E}}^2 & \\\ \end{bmatrix} \end{aligned}$$

Note that Equation eq-sigma_eps is the singular value decomposition of Σ_ℰ if the axis radii satisfy a_ℰ ≥ b_ℰ. Thus more generally,

Important

The ellipse ℰ is characterized by the equation

(x − x_ℰ)^TΣ_ℰ(x − x_ℰ) = 1

Or

x^TA_ℰx + b_ℰ^Tx + c_ℰ = 0

where A_ℰ = Σ_ℰ, b_ℰ = 2Σ_ℰx_ℰ and c_ℰ = x_ℰ^TΣ_ℰx_ℰ − 1. Denoting x^T = [x, y], ellipse ℰ can be defined algebraically as

E(x, y) = e₁x² + e₂xy + e₃y² + e₄x + e₅y + e₆ = 0,

where $\mathbf{A}_{\mathcal{E}} = \begin{bmatrix} e_1 & e_2/2 \\ e_2/2 & e_3 \end{bmatrix}$, b_ℰ^T = [e₄, e₅] and c_ℰ = e₆. This algebraic form is the convenient one that we will use in order to compute the intersection points of two intersecting ellipses.

Intersection Points of Two Ellipses

We explain how we can retrieve the intersection points of two ellipses. Our presentation complements Eberly.

First let (ℰ_i)_{1 ≤ i ≤ 2} be two ellipses defined as

(x, y) ∈ ℰ_i ⇔ E_i(x, y) = e_i1x² + e_i2xy + e_i3y² + e_i4x + e_i5y + e_i6 = 0

The intersection points of ellipses (ℰ_i)_{1 ≤ i ≤ 2} satisfy Equation eq-twoellipses for i ∈ {1, 2}, i.e., the following equation system holds for intersection points

$$\begin{aligned} \left\{ \begin{matrix} E_1(x,y) = 0 \\ E_2(x,y) = 0 \end{matrix} \right. \end{aligned}$$

Now let us rewrite E_i(x, y) as a quadratic polynomial in x, i.e.,

E_i(x, y) = e_i1x² + (e_i2y + e_i4)x + (e_i3y² + e_i5y + e_i6) = 0

Conveniently we define auxiliary polynomials in y

$$\begin{aligned} \begin{aligned} p_0(y) &= e_{13} y^2 + e_{15} y + e_{16} & q_0(y) &= e_{23} y^2 + e_{25} y + e_{26} \\\ p_1(y) &= e_{12} y + e_{14} & q_1(y) &= e_{22} y + e_{24} \\\ p_2(y) &= e_{11} & q_2(y) &= e_{21} \end{aligned} \end{aligned}$$

Introducing the polynomials above, Equation eq-twoellipses is rewritten as

$$\begin{aligned} \left\{ \begin{matrix} p_2(y) x^2 + p_1(y) x + p_0(y) = 0 \\\ q_2(y) x^2 + q_1(y) x + q_0(y) = 0 \end{matrix} \right. \end{aligned}$$

Suppose we know the y-coordinate of an intersection point, we can calculate the x-coordinate of this intersection point.

Indeed we multiply the first equation by q₂(y) and the second equation by p₂(y).

$$\begin{aligned} \left\{ \begin{matrix} q_2(y) \times \left( p_2(y) x^2 + p_1(y) x + p_0(y) \right)= 0\times q_2(y)\\\ p_2(y) \times \left( q_2(y) x^2 + q_1(y) x + q_0(y) \right)= 0\times p_2(y) \end{matrix} \right. \end{aligned}$$

Then subtracting the first equation from the second equation, the monomial x² disappears. Thus:

Important

$$x = \frac{p_0(y)q_2(y) - p_2(y)q_0(y)}{p_1(y)q_2(y) - p_2(y)q_1(y)}.$$

Furthermore, Equation eq-system is equivalent to the following augmented equation system

$$\begin{aligned} \left\{ \begin{array}{rl} E_1(x,y) &= 0 \\\ x\times E_1(x,y) &= 0 \\\ E_2(x,y) &= 0 \\\ x\times E_2(x,y) &= 0 \\\ \end{array} \right., \end{aligned}$$

And we see more clearly in matrix notation that

Important

[1, x, x², x³]^T is in the nullspace of B(y), where B(y) is defined as

$$\begin{aligned} \underbrace{ \begin{bmatrix} p_{0}(y) & p_{1}(y) & p_{2}(y) & 0 \\\ 0 & p_{0}(y) & p_{1}(y) & p_{2}(y) \\\ q_{0}(y) & q_{1}(y) & q_{2}(y) & 0 \\\ 0 & q_{0}(y) & q_{1}(y) & q_{2}(y) \end{bmatrix} }_{\mathbf{B}(y)} \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \end{aligned}$$

We observe that the vector [1, x, x², x³]^T is never zero for any real value x. Thus necessarily the nullspace Null(B(y)) is always nontrivial and that means the determinant of B(y) has to be zero.

Important

Let the polynomial R be defined as

$$R \overset{\textrm{def}}{=} \left( p_{0}q_{2} - p_{2}q_{0} \right)^2 - \left( p_{0}q_{1} - p_{1}q_{0} \right) \left( p_{1}q_{2} - p_{2}q_{1} \right),$$

Equation eq-system is equivalent to the following quartic equation in y.

det (B(y)) = R(y) = 0,

Using any polynomial solver, we get the 4 roots (y_i)_{1 ≤ i ≤ 4} of the quartic polynomial R and only keep those that are real. Finally (x_i)_{1 ≤ i ≤ 4} are deduced from Equation eq:xinter.

Implementation Notes

In Sara, we can use several solvers to retrieve the roots of polynomial R.

1. Companion matrix approach: since Sara depends on Eigen, Eigen has an unsupported Polynomial solver using this simple approach.
2. Jenkins-Traub iterative but very accurate approach also available in Sara.
3. Ferrari’s method available in Sara.

The implementation in Sara uses Ferrari's method. While more tedious to implement, the method has the advantage of being direct. Also, we experimentally observe Ferrari’s method can sometimes be numerically inaccurate in particular situations where for example one of the ellipse is quasi-degenerate.

In the future, depending on the use case, we can polish the roots to refine the root values.

Intersection Area of Two Ellipses

Our presentation complements Eberly. In the rest of the section, we consider two ellipses (ℰ_i)_{1 ≤ i ≤ 2} and we respectively denote

• the axes of ellipse ℰ_i by (a_i, b_i), the ellipse center by x_i, the orientation by θ_i, and the direction vectors of axis a_i and b_i by

  
  $$\begin{aligned} \begin{aligned} \mathbf{u}_i &\overset{\textrm{def}}{=}\begin{bmatrix} \cos(\theta_i) \\ \sin(\theta_i) \end{bmatrix} & \mathbf{v}_i &\overset{\textrm{def}}{=}\begin{bmatrix} -\sin(\theta_i) \\ \cos(\theta_i) \end{bmatrix}\end{aligned} \end{aligned}$$
  

• the area of the elliptic portion bounded a line segment and an arc for ellipse ℰ_i by B_i,
• the number of intersection points by L,

• the intersection points by p_i for i ∈ ⟦1, L⟧, sorted in a counter-clockwise order, i.e.,

  
  ∀i ∈ ⟦1, L − 1⟧, ∠([1,0]^T,p_i) < ∠([1,0]^T,p_i + 1)
  

  where ∠(., .) denotes the angle between two vectors in the plane ℝ².

• the polar angles of points (p_i)_{1 ≤ i ≤ L} with respect to ellipses ℰ₁ and ℰ₂ by (ϕ_i)_{1 ≤ i ≤ 2} and (ψ_i)_{1 ≤ i ≤ 2}, i.e.,

  
  $$\begin{aligned} \begin{gathered} \forall i \in \llbracket 1, L\rrbracket, \phi_i \overset{\textrm{def}}{=}\angle\left(\mathbf{u}_1, \mathbf{p}_i - \mathbf{x}_1\right) \\ \forall i \in \llbracket 1, L\rrbracket, \psi_i \overset{\textrm{def}}{=}\angle\left(\mathbf{u}_2, \mathbf{p}_i - \mathbf{x}_2\right)\end{gathered} \end{aligned}$$
  

Retrieving the polar angles

To retrieve the polar angles, we need to place ourselves in the coordinate system (x_i, u_i, v_i). Using the convenient function atan2 giving values in ] − π, π], we have

$$\begin{aligned} \begin{aligned} \phi_i &= \mathrm{atan2} \left( \langle \mathbf{p}_i-\mathbf{x}_1, \mathbf{v}_1 \rangle, \langle \mathbf{p}_i-\mathbf{x}_1, \mathbf{u}_1 \rangle \right)\\\ \psi_i &= \mathrm{atan2} \left( \langle \mathbf{p}_i-\mathbf{x}_2, \mathbf{v}_2 \rangle, \langle \mathbf{p}_i-\mathbf{x}_2, \mathbf{u}_2 \rangle \right)\end{aligned} \end{aligned}$$

0 or 1 intersection point

Either one ellipse is contained in the other or there are separated as illustrated in Figure [figinter01].

An ellipse, say ℰ₁, is contained in the other ℰ₂ if and only if its center satisfies E₂(x₁) < 0. In that case, the area of the intersection is just the area of ellipse ℰ₁. Otherwise, if there is no containment, the intersection area is zero. In summary,

$$\begin{aligned} \mathop{\mathrm{area}}(\mathcal{E}_1 \cap \mathcal{E}_2) = \left\{ \begin{array}{ll} \pi a_1 b_1 & \textrm{if}\ E_2(\mathbf{x}_1) < 0\\\ \pi a_2 b_2 & \textrm{if}\ E_1(\mathbf{x}_2) < 0\\\ 0 & \textrm{otherwise} \end{array} \right. \end{aligned}$$

2 intersection points

We will not detail the case when Polynomial eq-detBy have 2 roots with multiplicity 2. This still corresponds to the case where there are two intersection points. But because of the root multiplicities, one ellipse is contained in the other one and then Equation eq:eq-area01 gives the correct intersection area.

Otherwise, we have to consider two cases as illustrated in Figure [figinter2], which Eberly apparently forgot to consider. Namely, the cases correspond to whether the center of ellipses ℰ₁ and ℰ₂ are on the same side or on opposite side with respect to the line (p₁, p₂).

Denoting a unit normal of the line going across the intersection points (p₁, p₂) by n (cf. Figure 1.9). If the ellipse centers x₁ and x₂ are on opposite side with respect to the line (p₁, p₂), i.e.,

⟨n, x₁ − p₁⟩⟨n, x₂ − p₁⟩ < 0,

then

area (ℰ₁ ∩ ℰ₂) = B₁(ϕ₁, ϕ₂) + B₂(ψ₁, ψ₂)

If they are on the same side with respect to the line (p₁, p₂), i.e.,

⟨n, x₁ − p₁⟩⟨n, x₂ − p₁⟩ > 0,

then

$$\begin{aligned} \mathop{\mathrm{area}}(\mathcal{E}_1 \cap \mathcal{E}_2) = \left\{ \begin{array}{ll} \displaystyle \left( \pi a_1 b_1 - B_1(\phi_1, \phi_2) \right) + B_2(\psi_1, \psi_2) & \textrm{if} |\langle\mathbf{n},\mathbf{x}_1-\mathbf{p}_1\rangle| \leq |\langle\mathbf{n},\mathbf{x}_2-\mathbf{p}_1\rangle| \\\ \\\ \displaystyle B_1(\phi_1, \phi_2) + \left( \pi a_2 b_2 - B_2(\psi_1, \psi_2) \right) & \textrm{otherwise}. \end{array} \right. \end{aligned}$$

Note that the condition |⟨n, x₁ − p₁⟩| ≤ |⟨n, x₂ − p₁⟩| in Equation eqinter2b just expresses the fact that the distance of ellipse center x₁ to the line (p₁, p₂) is smaller than the distance of ellipse center x₂ to the line (p₁, p₂).

3 and 4 intersection points

These cases are rather easy to handle. Indeed, we see geometrically from Figure [fig-inter34],

$$\mathop{\mathrm{area}}(\mathcal{E}_1 \cap \mathcal{E}_2) = \sum_{i=1}^{L} \underbrace{\min \left( B_1(\phi_i, \phi_{i+1}), B_2(\psi_i, \psi_{i+1}) \right)}_{\textrm{smallest of elliptic portion area}} + \underbrace{\frac{1}{2} \sum_{i=1}^{L} \left| \det\left(\mathbf{p}_i, \mathbf{p}_{i+1}\right) \right|}_{\textrm{area of polygon}\ (\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_L)}$$

with ϕ_L + 1 = ϕ₁, ψ_L + 1 = ψ₁ and p_L + 1 = p₁.