Check if an ellipse is within another ellipse

Question

I have an ellipse $E_1$ centered at $(h,k)$, with semi-major axis $r_x$, semi-minor axis $r_y$, both aligned with the Cartesian plane.

How do I determine if another ellipse $E_2$ is within this given ellipse $E_1$?

$E_2$ can be anywhere in the Cartesian plane. What is given, is the centerpoint at $(i,j)$, the semi-major axis $r_x$ and the semi-minor axis $r_y$ and a rotation angle $\alpha$ (can be $0$, so no rotation).

I need this for a computed algorithm. Given this computer science background, what i use right now is the formula from here, and choose a point on the ellipse $E_2$, check if its within $E_1$ and choose another point, $1°$ further and check that point again and so on, until i complete $360°$.

I was thinking, that there has to be a better solution, a more formal and complete one (i think it should be possible to get a wrong result with the current algorithm in some very special cases). Yet, i haven't found a better solution.

Answer 1

Hint:

Without loss of generality, one of the ellipses is the unit circle centered at the origin. (If not, you can translate its center to the origin, counter-rotate to bring its major axis horizontal and rescale non-isotropically by the axis lengths).

Then the problem reduces to checking if an ellipse is wholly contained in the unit circle.

Let the parametric equation of that ellipse be

$$\vec p=\vec p_c+\vec a\cos t+\vec b\sin t.$$ We need to guarantee the inequality

$$\vec p^2\le 1,$$ i.e.

$$\vec p_c^2+\vec a^2\cos^2t+\vec b^2\sin^2t+2\vec p_c\vec a\cos t+2\vec p_c\vec b\sin t\le 1$$ (we have $\vec a\vec b=0$).

The inequality is ensured by finding the maxima of the trigonometric polynomial and showing the they don't exceed $1$.

Unfortunately, this leads to a quartic equation, so that an analytical solution promises to be painful. The problem is very close to that of finding the shortest distance of a point to an ellipse, or to ellipse offsetting.

---

There is an analytical solution, but it is a little scary.

• by an affine transform, stretch the plane so that the contained ellipse becomes a circle, while the containing ellipse becomes centered at the origin and axis-aligned, hence of the form

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

• "deflate" the circle while you "inflate" the ellipse; that means that you reduce the radius of the circle until it reduces to a point, while you compute the corresponding offset curve of the ellipse (https://en.wikipedia.org/wiki/Parallel_curve).

• the offset curve has a known implicit equation, given in "Brief Atlas of Offset Curves, Juana Sandra", available online.

• the point is contained in the offset curve, hence the circle in the ellipse, when the implicit expression has the desired sign.

Answer 2

Hint :

In the general case (any position of center, any length of axes, any angle of orientation) for both ellipses :

Express the equations of the ellipses on the general form : $$y^2+A_1xy+B_1x^2+C_1y+D_1x+F_1=0 \qquad\qquad \text{Ellipse } E_1$$ $$y^2+A_2xy+B_2x^2+C_2y+D_2x+F_2=0 \qquad\qquad \text{Ellipse } E_1$$

For example if $E_1$ is centered at $(h,k)$, with semi-major axis $r_x$, semi-minor axis $r_y$, both aligned with the Cartesian plane : $$\left(\frac{y-k}{r_y}\right)^2+\left(\frac{x-h}{r_x}\right)^2=1$$ $A_1=0\quad;\quad B_1=\left(\frac{r_y}{r_x}\right)^2\quad ;\quad C_1=2k \quad;\quad D_1=-2h\left(\frac{r_y}{r_x}\right)^2 \quad;\quad F_1=k^2+h^2\left(\frac{r_y}{r_x}\right)^2-1$

Do the same for the parameters of $E_2$

Of course, it is a bit more complicated if the axes of the ellipses are not parallel to the Cartesian axes.

So, at this point, the numerical values of the ten coefficients are known.

Solve for $(x,y)$ (with numerical calculus) the system of equations : $$\begin{cases} y^2+A_1xy+B_1x^2+C_1y+D_1x+F_1=0 \\ y^2+A_2xy+B_2x^2+C_2y+D_2x+F_2=0 \end{cases}$$

Four roots $(x,y)$ are obtained (real and\or complex).

The result is : $\begin{cases} \text{Case 1 : If one or more root is real, the ellipses intersect.} \\ \text{Case 2 : If all roots are complex, see case below:} \\ \end{cases}$

In case 2, determine if the center of $E_2$ is outside $E_1$. And determine if the center of $E_2$ is outside $E_1$. If both are outside, the two ellipses are outside each other. If not, one is inside the other.

NOTE :

To determine if a point $(x,y)$ is outside $E_1$, compute $$\quad G_1=y^2+A_1xy+B_1x^2+C_1y+D_1x+F_1$$ If $G_1>0$ the point is outside.

Similarly to determine if a point is outside $E_2$ .

Answer 3

I'll suppose too the enclosing ellipse scaled as a circle of radius $1$ and the figure rotated so that the (possibly) enclosed ellipse will have a vertical semi-major axis $b$ and an horizontal semi-minor axis $a$. Since an ellipse stretched in any direction remains an ellipse this can be done even if the ellipses axis are not aligned.

I will consider only the case $\;0<a\le b<1\;$ (for $\,b=1$ we need $\,x_0=y_0=0$) and start with an illustration for $\;a=\frac 1{10},\ b=\frac 3{10}$ :

Let's first search the orbit of the center $(x_0,y_0)$ of the enclosed ellipse when remaining tangent to the circle in the first quadrant and obtain the system (with $(x,y)$ the common tangent point and $\,R=1$) :

\begin{align} \tag{1}1&=\left(\frac{x-x_0}a\right)^2+\left(\frac{y-y_0}b\right)^2\\ \tag{2}R^2&=x^2+y^2\\ \tag{3}0&=2\frac{x-x_0}{a^2}dx+2\frac{y-y_0}{b^2}dy\\ \end{align} The third one was obtained by differentiation of $(1)$. Differentiation of $(2)$ leaving to $\;\displaystyle \frac {dy}{dx}=-\frac xy\;$ gives us following relation at the tangent point (for $y\neq 0$) :

$$-\frac {dy}{dx}=\frac xy=\frac{x-x_0}{y-y_0}\frac{b^2}{a^2}$$ that I will rewrite as : $$\tag{4}by\frac{x-x_0}{a}=ax\frac{y-y_0}{b}$$ Its LHS may be substituted in $(1)$ multiplied by $\,b^2y^2\,$ giving :

\begin{align} b^2y^2 &=\left(ax\frac{y-y_0}{b}\right)^2+b^2y^2\left(\frac{y-y_0}b\right)^2\\ \left(b^2y\right)^2 &=\left(a^2x^2+b^2y^2\right)\left(y-y_0\right)^2\\ \tag{5}|y-y_0|&=\frac{b^2\,|y|}{\sqrt{a^2x^2+b^2y^2}}\\ \text{while from }&\ (4)\ \text{(or by symmetry) :}\\ \tag{6}|x-x_0|&=\frac{a^2\,|x|}{\sqrt{a^2x^2+b^2y^2}}\\ \end{align} The absolute values are not really necessary and we obtain the parametric answer (for $R=1$ and $\;x=\cos(t),\,y= \sin(t)\;$) : $$\tag{7}x_0(t)=\cos(t)-\frac{a^2\;\cos(t)}{\sqrt{a^2\cos(t)^2+b^2\sin(t)^2}},\ y_0(t)=\sin(t)-\frac{b^2\;\sin(t)}{\sqrt{a^2\cos(t)^2+b^2\sin(t)^2}}$$ leading to the black curve in the first picture and this one for $\;a=\frac 3{10},\ b=\frac 3{4}$ :

This picture illustrates that in the case $\,a<b^2\,$ the value $\,x_0\,$ has to be restrained (because the enclosed ellipse will get tangent at two points instead of one when reaching $\,y_0=0$ as shown).
Setting $y_0=0$ in the previous solution we find that we must impose (with $\,R=1$) : $$\tag{8}|x_0|\le \sqrt{\left(\frac {R^2}{b^2}-1\right)\left(b^2-a^2\right)},\quad\text{for}\ a\le b^2\quad$$ We found the external bounds for $(x_0,\ y_0)$ wih $(7)$ and $(8)$ but the initial problem remains :
what conditions are required for the center of an ellipse to remain inside the bounds $(7)$ ?

A mathematical solution is to use Cauchy's integral formula with the (parametric) contour given by $\;z_0:=x_0(t)+i\,y_0(t)$, using $(7)$, to obtain : $$\tag{9}f(z):=\frac 1{2\pi i}\int_0^{2\pi} \frac{dz_0}{z_0-z}=\frac 1{2\pi}\int_0^{2\pi}\frac {e^{it}\left(1-\large{\frac{a^2b^2}{\sqrt{a^2\cos(t)^2+b^2\sin(t)^2}^3}}\right)}{e^{it}-\large{\frac{a^2\cos(t)+ib^2\sin(t)}{\sqrt{a^2\cos(t)^2+b^2\sin(t)^2}}}-z_0}dt$$ which returns indeed $f(z)=1$ if the center $z$ of the ellipse is inside the black curve in the complex plane and $0$ if outside (supposing that $(8)$ applies for the real part of $z$).
Of course a simplification of this integral would be welcome! $$-$$ To find if the ellipse is 'inside' it may seem easier to use $(4)$ to find the farthest point of the ellipse $(x,y)$ (the tangent one verifying $\,x^2+y^2=R^2\,$) and see if $R^2>1$.
Let's try this using $\,\displaystyle u:=\frac{x-x_0}a,\ v:=\frac{y-y_0}b\,$ and $\,u^2+v^2=1\,$ from $(1)$ : \begin{align} by\frac{x-x_0}{a}&=ax\frac{y-y_0}{b}\\ b\,(y_0+b\,v)u&=a\,(x_0+a\,u)v\\ b\,y_0 u&=((b^2-a^2)u+a\,x_0)v\\ (b\,y_0 u)^2&=((b^2-a^2)u+a\,x_0)^2(1-u^2)\\ 0&=\Delta^2 u^4+2ax_0 \Delta\, u^3+((ax_0)^2+(by_0)^2-\Delta^2)u^2-2ax_0\Delta \,u-(a x_0)^2\\ &\quad\text{for }\ \Delta=b^2-a^2 \end{align}

and that's where you obtain a quartic in $u$ with an awful looking explicit solution but you may apply the algorithm (numerically!) directly or find the solution in $(0,1)$ by iterations (say using Newton's method). Once $u$ and then $v=\sqrt{1-u^2}$ obtained you may deduce $x$ and $y$ and conclude with $R^2=x^2+y^2$ (sorry if this conclusion doesn't add much to the other answers, I hope it may help anyway...).

Answer 4

First, express $E_2$ in cartesian coordinates $(x,y)$. Then substitute $X=\dfrac{x+h}{r_x},Y=\dfrac{y+k}{r_y}$. This substitution transforms $E_1$ to the unit circle centred at the origin, so the problem reduces to deciding whether the transformed ellipse $E_2'$ is contained in this unit circle.

To do this, you can use Lagrange multipliers to maximise $X^2+Y^2$ subject to the transformed equation of $E_2'$. I don't have time for this at the moment, but I don't foresee any problems.

Edited to add: I have tried to work through the Lagrange multiplier solution, and unfortunately it doesn't look any easier than the solutions proposed in other answers.

Answer 5

I came up with a method for a $n$-ellipse. The 2D case is a particular one and the same method can be applied.

It supposes that you know the center of the ellipses and the axis.

• Refer to this question if you know the coeffs

$$ax^2 + bxy + cy^2 +dx + ey + f = 0$$

• To compute the matrix $A$ by using major axis $a$, minor axis $b$ and angle $\varphi$:

$$A = \begin{bmatrix}a \cdot \cos \varphi & a \cdot \sin \varphi \\ -b \cdot \sin \varphi & b \cdot \cos \varphi\end{bmatrix}$$

General method

Let $E_1$ and $E_2$ be ellipses of center $\alpha$ and $\beta$ and axis the matrices $\left[A\right]$ and $\left[B\right]$:

$$E_1: \left\{\left(\alpha + A \cdot p\right) \in \mathbb{R}^{n} : \|p\| \le 1\right\}$$ $$E_2: \left\{\left(\beta + B \cdot q\right) \in \mathbb{R}^{n} : \|q\| \le 1\right\}$$

I use the notation $a_i = \left[A\right]_i$, $b_i = \left[B\right]_{i}$, and $\langle u, \ v\rangle = u^{T} \cdot v$ refers to the inner product between $u$ and $v$.

To know if a point $P \in E_1$, it must satisfy

$$\sum_{i=1}^{n} \left(\dfrac{\langle P - \alpha, \ a_i\rangle}{\langle a_i, \ a_i \rangle}\right)^2 \le 1 \tag{3}\label{3}$$

So, to $E_2 \subset E_1$, then \eqref{3} must be satisfied by every point of $E_2$, which leads to \eqref{4}

$$\sum_{i=1}^{n} \left(\dfrac{\left\langle \left(\beta - \alpha + \sum_{j=1}^{n} q_j \cdot b_j\right) , \ a_i\right\rangle}{\langle a_i, \ a_i \rangle}\right)^2 \le 1 \ \ \ \forall \ \|q\| \le 1 \tag{4}\label{4}$$

Expanding \eqref{4} and setting new variables to get \eqref{5}

$$ \sum_{i=1}^{n} \left(C_i + \sum_{j=1}^{n} M_{ij} \cdot q_j \right)^2 \le 1 \ \ \ \ \forall \ \|q\| \le 1 \tag{5}\label{5} $$ $$C_{i} = \dfrac{\left\langle a_i, \ \beta - \alpha \right\rangle}{\left\langle a_i, \ a_i \right\rangle} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ M_{ij} = \dfrac{\left\langle a_i, \ b_i \right\rangle}{\left\langle a_i, \ a_i \right\rangle}$$

Or also using matrix notation

$$C = \left(A \cdot A^{T} \right)^{-1} \cdot A \cdot \left(\beta - \alpha\right)$$ $$M = \left(A \cdot A^{T} \right)^{-1} \cdot A \cdot B^{T}$$

\eqref{5} can be rewrote as \eqref{6} by using matrix notation

$$\left(C + Mq\right)^{T} \cdot \left(C + Mq\right) \le 1 \ \ \ \ \ \ \ \forall \|q\| \le 1 $$ $$q^T \cdot X \cdot q + 2 \cdot Y \cdot q + Z \le 1 \ \ \ \ \ \ \ \forall \|q\| \le 1 \tag{6}\label{6}$$ $$X = M^{T} \cdot M \ \ \ \ \ \ \ \ \ \ Y = C^{T} \cdot M \ \ \ \ \ \ \ \ \ \ Z = C^{T} \cdot C $$

I tried to find $q$ such maximizes the left side of \eqref{6} by using lagrange multiplier (to apply the constraint). Which involves solving

$$\begin{cases}X\cdot q + Y + \mu q = 0 \\ q^T \cdot q - 1 =0\end{cases}$$

Unfortunatelly it was harder than I thought to find a solution. So I used a different approach which is true for most cases, for 'almost' centered ellipses.

Completing squares to isolate $q$ we get \eqref{7} $$w^T \cdot X \cdot w \le W \ \ \ \ \ \ \ \forall \|w - N\| \le 1 \tag{7}\label{7}$$

$$w = q + N \ \ \ \ \ \ \ \ \ \ \ N = X^{-1} Y \ \ \ \ \ \ \ \ \ \ \ W = 1 + Y^{T} \cdot X^{-1} \cdot Y - Z$$

Since $X$ is a positive defined matrix, the left side is always positive. Finding the maximum of the left side happens when the vector $w$ is parallel to the eigenvector $v$ that gives the maximum eigenvalue $\lambda$.

This hypothesis works only when $N^T \cdot N < 1 + \left(V^T N\right)^2$. This doesn't happen when the centers $\alpha$ and $\beta$ are distant.

$$q + N = \mu \cdot v$$ $$\underbrace{q^T \cdot q}_{1} = \mu^2 \cdot \underbrace{v^T \cdot v}_1 - 2 v^T \cdot N + N^T \cdot N$$ $$\mu = v^T \cdot N \pm \sqrt{1 + \left(v^T N\right)^2 - N^T \cdot N}\tag{8}\label{8}$$

Applying \eqref{8} in \eqref{7} we obtain the final verification

$$\boxed{\mu^2 \cdot \lambda \le W}$$

Algorithm

I made a python algorithm that works with the hypothesis made above.

import numpy as np
from matplotlib import pyplot as plt

A = [[ 2.74630871, -3.17694060], [-2.22172112, -1.92056851]]
B = [[ 1.62742095, -2.34414119], [-1.97924175, -1.37408937]]
alpha = (0, 0)
beta = (0, 0)

AAT = np.dot(A, np.transpose(A))
ABT = np.dot(A, np.transpose(B))
AdC = np.dot(A, np.array(alpha)-beta)
C = np.linalg.solve(AAT, AdC)
M = np.linalg.solve(AAT, ABT)
X = np.dot(np.transpose(M), M)
Y = np.dot(C, M)
Z = np.dot(C, C)
N = np.linalg.solve(X, Y)
W = 1 + np.dot(Y, N) - Z
eigvals, eigvecs = np.linalg.eigh(X)
index = tuple(eigvals).index(max(eigvals))
lamda, v = eigvals[index], eigvecs[:, index]
b = np.dot(v, N)
c = np.dot(N, N) - 1
mu1 = b + np.sqrt(b**2 - c)
mu2 = b - np.sqrt(b**2 - c)

print("X = ")
print(X)
print(f"Y = {Y}")
print(f"Z = {Z}")
print(f"N = {N}")
print(f"W = {W}")
print(f"mu1 = {mu1}")
print(f"mu2 = {mu2}")
print(mu1**2 * lamda <= W)
print(mu2**2 * lamda <= W)

theta = np.linspace(0, 2*np.pi, 129)
points = np.tensordot(np.cos(theta), A[0], axes=0)
points += np.tensordot(np.sin(theta), A[1], axes=0)
xvals, yvals = np.transpose(points)
plt.plot(xvals, yvals, label="A")
points = np.tensordot(np.cos(theta), B[0], axes=0)
points += np.tensordot(np.sin(theta), B[1], axes=0)
xvals, yvals = np.transpose(points)
plt.plot(xvals, yvals, label="B")
plt.legend()
plt.show()
```