I have seen a lot of papers on how to find points of intersection between two ellipses for 2D case, but i only need to check if two ellipses are in collision. I don't need to know points of intersection if there are any. Is there simplified algorithm for this. Thanks.
I know center and two radii for every ellipse. Both ellipses can be rotated.
Checking whether two ellipses intersect, or more generally two n-dimenensional ellipsoids, can be done by solving a 1D smooth convex optimization problem. This is very easy and fast to do computationally.
First, the data format that ellipsoids are given in should be a a vector which defines the center of the ellipse, and a symmetric positive definite matrix which defines the shape of the ellipse. In particular, let $a,b \in \mathbb{R}^n$, let $A,B \in \mathbb{R}^{n \times n}$ be symmetric positive definite matrices, and define the ellipsoids as follows: \begin{align*} E_A &= \{x : (x-a)^T A (x-a) \le 1\} \\ E_B &= \{x : (x-b)^T B (x-b) \le 1\}. \end{align*} The vector $a$ is the centerpoint of the ellipsoid $E_A$. The eigenvectors of the matrix $A$ are the unit vectors that point in the directions of the primary axes of $E_A$. The inverse square roots of the corresponding eigenvalues ($1/\sqrt{\lambda_i}$) of $A$ are the lengths of the primary axes of $E_A$. This is shown in the image below. The vector $b$ and matrix $B$ characterize $E_B$ in the same manner.
The vectors $a$ and $b$ and matrices $A$ and $B$ uniquely define the ellipsoids $E_A$ and $E_B$, respectively. Also, for every ellipsoid, there is a corresponding vector and matrix (like $a$ and $A$) that can put it in this format. For more on this, you can see the following wikipedia article: https://en.wikipedia.org/wiki/Ellipsoid
If your ellipses (2D) or ellipsoids (nD) are given in another format, such as axis lengths and angles, the first step is to convert them to this format by computing the vectors $a$ and $b$, and matrices $A$ and $B$.
Now define the following convex scalar function $K:(0,1)\rightarrow \mathbb{R}$, $$K(s) := 1 - (b-a)^T\left(\frac{1}{1-s}A^{-1} + \frac{1}{s}B^{-1}\right)^{-1}(b-a).$$
The following result holds: $$E_A \cap E_B = \{\} \quad \text{if and only if}\quad K(s)<0 \quad \text{for some}\quad s \in (0,1).$$
This is proven in Proposition 2 of the following paper:
Gilitschenski, Igor, and Uwe D. Hanebeck. "A robust computational test for overlap of two arbitrary-dimensional ellipsoids in fault-detection of kalman filters." 2012 15th International Conference on Information Fusion. IEEE, 2012.
So you can check whether the two ellipsoids intersect by finding the minimum of $K$ on $(0,1)$ using any 1D convex minimization routine, then check if $K(s)$ is greater than 0 at the minima.
I like to use Brent's algorithm to find the minimum, because it is fast and implementations are available in many software libraries such as scipy. https://en.wikipedia.org/wiki/Brent%27s_method
The minimization algorithm will evaluate $K(s)$ many times for different values of $s$. You can avoid performing matrix inverses/linear solves at each evaluation of $K(s)$ by precomputing the generalized eigenvalue decomposition of $A$ with respect to $B$, then working in the basis of generalized eigenvectors, because both matrices become diagonal in this basis. http://fourier.eng.hmc.edu/e161/lectures/algebra/node7.html
Here is some python code I wrote to do this. In this code, $\Sigma_A := A^{-1}$, $\Sigma_B := B^{-1}$, $\mu_A = a$, and $\mu_B = b$, and $\tau$ is an optional scalar that makes the ellipses uniformly bigger or smaller.
import numpy as np
from scipy.linalg import eigh
from scipy.optimize import minimize_scalar
def ellipsoid_intersection_test(Sigma_A, Sigma_B, mu_A, mu_B, tau=1.0):
lambdas, Phi = eigh(Sigma_A, b=Sigma_B)
v_squared = np.dot(Phi.T, mu_A - mu_B) ** 2
res = minimize_scalar(K_function,
bracket=[0.0, 0.5, 1.0],
args=(lambdas, v_squared, tau))
return (res.fun[0] >= 0)
def K_function(s, lambdas, v_squared, tau):
return 1.-(1./tau*2)np.sum(v_squared((s(1.-s))/(1.+s*(lambdas-1.))))
Here are a sequence of pictures showing how the function $K(s)$ changes as one ellipsoid slides past another:
The complete code I used to generate these pictures is in the following jupyter notebook: https://github.com/NickAlger/nalger_helper_functions/blob/master/tutorial_notebooks/ellipsoid_intersection_test_tutorial.ipynb
By a suitable stretching of the plane in the direction of the axis of one of the ellipses, you can turn this ellipse to a circle, while the other remains an ellipse.
Now checking if the circle and the ellipse have a nonempty intersection is the same as checking if the center of the circle is inside the outward offset curve of the ellipse, at distance $r$.
Unfortunately, such an offset curve is known to have an octic ($8^{th}$) degree implicit equation, which is just horrible. Check "Brief Atlas of Offset Curves", example 4.
This tends to show that there is no easy exact analytical solution. If you can do with an approximate solution, just replace the offset curve by another ellipse of axis $a+r$ and $b+r$.
Let we suppose that $E_1$ is an ellipse with equation $f(x,y)=\frac{x^2}{a}+\frac{y^2}{b}-1=0$ and $E_2$ is another ellipse. To check if $E_1$ and $E_2$ intersect, it is sufficient to check if $f(x,y)$ takes only positive values on $\partial E_2$. So we can take a parametrization of $\partial E_2$ and compute the stationary points for the quadratic function $f(x,y)$ on $\partial E_2$. If we values of $f$ in such points are positive, $E_1$ and $E_2$ do not intersect, otherwise they intersect.
Here I assumed that the ellipses lie on the euclidean plane, but the same argument can be extended also to check if two ellipses in $\mathbb{R}^3$ are "linked" or not.
Here is a C++ class for iteratively testing collision between two ellipses on a two-dimensional plane, given their center points, half-width vectors (each corresponding to a major or minor radius), and height/width ratios. The algorithm transforms the coordinates of the two-ellipse system so that the other ellipse becomes a circle, and sandwiches the other between two stretched and sheared (due to stretching of the system along other than one the axes of the ellipse) regular polygons that touch on the ellipse. This encloses the curve of the ellipse inside triangles. The algorithm iteratively splits the triangle in which collision might take place into two smaller triangles until it is known whether a collision takes place, or until an iteration limit is reached. No division or square roots are needed in the calculation. If one of the ellipses is inside the other, that is also detected as a collision.
To make one of the ellipses a circle, its coordinates and those of the polygons enclosing the other ellipse are transformed by:
$$\left[\begin{array}{c}x\\ y\end{array}\right] \longmapsto \left[\begin{array}{cc}hw\ wx^2 + wy^2 & hw\ wx\ wy - wx\ wy\\ hw\ wx\ wy - wx\ wy& hw\ wy^2 + wx^2\end{array}\right] \left[\begin{array}{c}x\\ y\end{array}\right] = \left[\begin{array}{c}hw\ wx\ (wy\ y + wx\ x) - wy\ (wx\ y - wy\ x)\\ hw\ wy\ (wy\ y + wx\ x) + wx\ (wx\ y - wy\ x)\end{array}\right],$$
where $hw$ is the height/width ratio of the ellipse and $(wx, wy)$ is its half-width vector. The squared radius of the circle resulting from the transformation is $hw^2\ (wx^2 + wy^2)^3$. At every iteration of the algorithm, to create a $2n$-gon from an $n$-gon in the transformed coordinates, neighboring center–vertex vectors are linearly combined using constants from a small precomputed table.
Source: http://yehar.com/blog/?p=2926
#include <math.h>
class EllipseCollisionTest {
private:
double *innerPolygonCoef;
double *outerPolygonCoef;
int maxIterations;
bool iterate(double x, double y, double c0x, double c0y, double c2x, double c2y, double rr) const {
for (int t = 1; t <= maxIterations; t++) {
double c1x = (c0x + c2x)*innerPolygonCoef[t];
double c1y = (c0y + c2y)*innerPolygonCoef[t];
double tx = x - c1x;
double ty = y - c1y;
if (tx*tx + ty*ty <= rr) {
return true;
}
double t2x = c2x - c1x;
double t2y = c2y - c1y;
if (tx*t2x + ty*t2y >= 0 && tx*t2x + ty*t2y <= t2x*t2x + t2y*t2y &&
(ty*t2x - tx*t2y >= 0 || rr*(t2x*t2x + t2y*t2y) >= (ty*t2x - tx*t2y)*(ty*t2x - tx*t2y))) {
return true;
}
double t0x = c0x - c1x;
double t0y = c0y - c1y;
if (tx*t0x + ty*t0y >= 0 && tx*t0x + ty*t0y <= t0x*t0x + t0y*t0y &&
(ty*t0x - tx*t0y <= 0 || rr*(t0x*t0x + t0y*t0y) >= (ty*t0x - tx*t0y)*(ty*t0x - tx*t0y))) {
return true;
}
double c3x = (c0x + c1x)*outerPolygonCoef[t];
double c3y = (c0y + c1y)*outerPolygonCoef[t];
if ((c3x-x)*(c3x-x) + (c3y-y)*(c3y-y) < rr) {
c2x = c1x;
c2y = c1y;
continue;
}
double c4x = c1x - c3x + c1x;
double c4y = c1y - c3y + c1y;
if ((c4x-x)*(c4x-x) + (c4y-y)*(c4y-y) < rr) {
c0x = c1x;
c0y = c1y;
continue;
}
double t3x = c3x - c1x;
double t3y = c3y - c1y;
if (ty*t3x - tx*t3y <= 0 || rr*(t3x*t3x + t3y*t3y) > (ty*t3x - tx*t3y)*(ty*t3x - tx*t3y)) {
if (tx*t3x + ty*t3y > 0) {
if (fabs(tx*t3x + ty*t3y) <= t3x*t3x + t3y*t3y || (x-c3x)*(c0x-c3x) + (y-c3y)*(c0y-c3y) >= 0) {
c2x = c1x;
c2y = c1y;
continue;
}
} else if (-(tx*t3x + ty*t3y) <= t3x*t3x + t3y*t3y || (x-c4x)*(c2x-c4x) + (y-c4y)*(c2y-c4y) >= 0) {
c0x = c1x;
c0y = c1y;
continue;
}
}
return false;
}
return false; // Out of iterations so it is unsure if there was a collision. But have to return something.
}
public:
// Test for collision between two ellipses, "0" and "1". Ellipse is at (x, y) with major or minor radius
// vector (wx, wy) and the other major or minor radius perpendicular to that vector and hw times as long.
bool collide(double x0, double y0, double wx0, double wy0, double hw0,
double x1, double y1, double wx1, double wy1, double hw1) const {
float rr = hw1*hw1*(wx1*wx1 + wy1*wy1)*(wx1*wx1 + wy1*wy1)*(wx1*wx1 + wy1*wy1);
float x = hw1*wx1*(wy1*(y1 - y0) + wx1*(x1 - x0)) - wy1*(wx1*(y1 - y0) - wy1*(x1 - x0));
float y = hw1*wy1*(wy1*(y1 - y0) + wx1*(x1 - x0)) + wx1*(wx1*(y1 - y0) - wy1*(x1 - x0));
float temp = wx0;
wx0 = hw1*wx1*(wy1*wy0 + wx1*wx0) - wy1*(wx1*wy0 - wy1*wx0);
float temp2 = wy0;
wy0 = hw1*wy1*(wy1*wy0 + wx1*temp) + wx1*(wx1*wy0 - wy1*temp);
float hx0 = hw1*wx1*(wy1*(temp*hw0)-wx1*temp2*hw0)-wy1*(wx1*(temp*hw0)+wy1*temp2*hw0);
float hy0 = hw1*wy1*(wy1*(temp*hw0)-wx1*temp2*hw0)+wx1*(wx1*(temp*hw0)+wy1*temp2*hw0);
if (wx0*y - wy0*x < 0) {
x = -x;
y = -y;
}
if ((wx0 - x)*(wx0 - x) + (wy0 - y)*(wy0 - y) <= rr) {
return true;
} else if ((wx0 + x)*(wx0 + x) + (wy0 + y)*(wy0 + y) <= rr) {
return true;
} else if ((hx0 - x)*(hx0 - x) + (hy0 - y)*(hy0 - y) <= rr) {
return true;
} else if ((hx0 + x)*(hx0 + x) + (hy0 + y)*(hy0 + y) <= rr) {
return true;
} else if (x*(hy0 - wy0) + y*(wx0 - hx0) <= hy0*wx0 - hx0*wy0 &&
y*(wx0 + hx0) - x*(wy0 + hy0) <= hy0*wx0 - hx0*wy0) {
return true;
} else if (x*(wx0-hx0) - y*(hy0-wy0) > hx0*(wx0-hx0) - hy0*(hy0-wy0)
&& x*(wx0-hx0) - y*(hy0-wy0) < wx0*(wx0-hx0) - wy0*(hy0-wy0)
&& (x*(hy0-wy0) + y*(wx0-hx0) - hy0*wx0 + hx0*wy0)*(x*(hy0-wy0) + y*(wx0-hx0) - hy0*wx0 + hx0*wy0)
<= rr*((wx0-hx0)*(wx0-hx0) + (wy0-hy0)*(wy0-hy0))) {
return true;
} else if (x*(wx0+hx0) + y*(wy0+hy0) > -wx0*(wx0+hx0) - wy0*(wy0+hy0)
&& x*(wx0+hx0) + y*(wy0+hy0) < hx0*(wx0+hx0) + hy0*(wy0+hy0)
&& (y*(wx0+hx0) - x*(wy0+hy0) - hy0*wx0 + hx0*wy0)*(y*(wx0+hx0) - x*(wy0+hy0) - hy0*wx0 + hx0*wy0)
<= rr*((wx0+hx0)*(wx0+hx0) + (wy0+hy0)*(wy0+hy0))) {
return true;
} else {
if ((hx0-wx0 - x)*(hx0-wx0 - x) + (hy0-wy0 - y)*(hy0-wy0 - y) <= rr) {
return iterate(x, y, hx0, hy0, -wx0, -wy0, rr);
} else if ((hx0+wx0 - x)*(hx0+wx0 - x) + (hy0+wy0 - y)*(hy0+wy0 - y) <= rr) {
return iterate(x, y, wx0, wy0, hx0, hy0, rr);
} else if ((wx0-hx0 - x)*(wx0-hx0 - x) + (wy0-hy0 - y)*(wy0-hy0 - y) <= rr) {
return iterate(x, y, -hx0, -hy0, wx0, wy0, rr);
} else if ((-wx0-hx0 - x)*(-wx0-hx0 - x) + (-wy0-hy0 - y)*(-wy0-hy0 - y) <= rr) {
return iterate(x, y, -wx0, -wy0, -hx0, -hy0, rr);
} else if (wx0*y - wy0*x < wx0*hy0 - wy0*hx0 && fabs(hx0*y - hy0*x) < hy0*wx0 - hx0*wy0) {
if (hx0*y - hy0*x > 0) {
return iterate(x, y, hx0, hy0, -wx0, -wy0, rr);
}
return iterate(x, y, wx0, wy0, hx0, hy0, rr);
} else if (wx0*x + wy0*y > wx0*(hx0-wx0) + wy0*(hy0-wy0) && wx0*x + wy0*y < wx0*(hx0+wx0) + wy0*(hy0+wy0)
&& (wx0*y - wy0*x - hy0*wx0 + hx0*wy0)*(wx0*y - wy0*x - hy0*wx0 + hx0*wy0) < rr*(wx0*wx0 + wy0*wy0)) {
if (wx0*x + wy0*y > wx0*hx0 + wy0*hy0) {
return iterate(x, y, wx0, wy0, hx0, hy0, rr);
}
return iterate(x, y, hx0, hy0, -wx0, -wy0, rr);
} else {
if (hx0*y - hy0*x < 0) {
x = -x;
y = -y;
}
if (hx0*x + hy0*y > -hx0*(wx0+hx0) - hy0*(wy0+hy0) && hx0*x + hy0*y < hx0*(hx0-wx0) + hy0*(hy0-wy0)
&& (hx0*y - hy0*x - hy0*wx0 + hx0*wy0)*(hx0*y - hy0*x - hy0*wx0 + hx0*wy0) < rr*(hx0*hx0 + hy0*hy0)) {
if (hx0*x + hy0*y > -hx0*wx0 - hy0*wy0) {
return iterate(x, y, hx0, hy0, -wx0, -wy0, rr);
}
return iterate(x, y, -wx0, -wy0, -hx0, -hy0, rr);
}
return false;
}
}
}
// Test for collision between an ellipse of horizontal radius w0 and vertical radius h0 at (x0, y0) and
// an ellipse of horizontal radius w1 and vertical radius h1 at (x1, y1)
bool collide(double x0, double y0, double w0, double h0, double x1, double y1, double w1, double h1) const {
double x = fabs(x1 - x0)*h1;
double y = fabs(y1 - y0)*w1;
w0 *= h1;
h0 *= w1;
double r = w1*h1;
if (x*x + (h0 - y)*(h0 - y) <= r*r || (w0 - x)*(w0 - x) + y*y <= r*r || x*h0 + y*w0 <= w0*h0
|| ((x*h0 + y*w0 - w0*h0)*(x*h0 + y*w0 - w0*h0) <= r*r*(w0*w0 + h0*h0) && x*w0 - y*h0 >= -h0*h0 && x*w0 - y*h0 <= w0*w0)) {
return true;
} else {
if ((x-w0)*(x-w0) + (y-h0)*(y-h0) <= r*r || (x <= w0 && y - r <= h0) || (y <= h0 && x - r <= w0)) {
return iterate(x, y, w0, 0, 0, h0, r*r);
}
return false;
}
}
// Test for collision between an ellipse of horizontal radius w and vertical radius h at (x0, y0) and
// a circle of radius r at (x1, y1)
bool collide(double x0, double y0, double w, double h, double x1, double y1, double r) const {
double x = fabs(x1 - x0);
double y = fabs(y1 - y0);
if (x*x + (h - y)*(h - y) <= r*r || (w - x)*(w - x) + y*y <= r*r || x*h + y*w <= w*h
|| ((x*h + y*w - w*h)*(x*h + y*w - w*h) <= r*r*(w*w + h*h) && x*w - y*h >= -h*h && x*w - y*h <= w*w)) {
return true;
} else {
if ((x-w)*(x-w) + (y-h)*(y-h) <= r*r || (x <= w && y - r <= h) || (y <= h && x - r <= w)) {
return iterate(x, y, w, 0, 0, h, r*r);
}
return false;
}
}
EllipseCollisionTest(int maxIterations) {
this->maxIterations = maxIterations;
innerPolygonCoef = new double[maxIterations+1];
outerPolygonCoef = new double[maxIterations+1];
for (int t = 0; t <= maxIterations; t++) {
int numNodes = 4 << t;
innerPolygonCoef[t] = 0.5/cos(4*acos(0.0)/numNodes);
outerPolygonCoef[t] = 0.5/(cos(2*acos(0.0)/numNodes)*cos(2*acos(0.0)/numNodes));
}
}
~EllipseCollisionTest() {
delete[] innerPolygonCoef;
delete[] outerPolygonCoef;
}
};
Thought I'd give another option in case anyone stumbles across this post in the future. I think my implementation should be faster than the other methods posted as it contains no iterative calculations.
"""
Based on the publication:
Fernando Etayo, Laureano Gonzalez-Vega, Natalia del Rio.
A new approach to characterizing the relative position of two
ellipses depending on one parameter. Computer Aided Geometric Design. 2006.
Jason Conradt, 2023
"""
import numpy as np
def detectOverlap(e1,e2,p1,p2):
#e1 is the tuple containing the axes a and b of ellipse 1
#e2 is the tuple containing the axes a and b of ellipse 2
#p1 is the tuple containing the position and orientation of ellipse 1
#p2 is the tuple containing the position and orientation of ellipse 2
a1, b1 = e1
a2, b2 = e2
x1, y1, theta1 = p1
x2, y2, theta2 = p2
#Represent the ellipse in matricial form
A1 = np.array([[1/a1**2,0,0],[0,1/b1**2,0],[0,0,-1]])
A2 = np.array([[1/a2**2,0,0],[0,1/b2**2,0],[0,0,-1]])
#Define relevant matrices for a rigid body transformation
M1 = np.array([[np.cos(theta1),-np.sin(theta1),x1],[np.sin(theta1),np.cos(theta1),y1],[0,0,1]])
M2 = np.array([[np.cos(theta2),-np.sin(theta2),x2],[np.sin(theta2),np.cos(theta2),y2],[0,0,1]])
M1_inv = np.linalg.inv(M1)
M2_inv = np.linalg.inv(M2)
#Both ellipses have been transformed to the same coordinate system
A = np.matmul(np.matmul(np.transpose(M1_inv),A1),M1_inv)
B = np.matmul(np.matmul(np.transpose(M2_inv),A2),M2_inv)
#Explicate the indices for convenience
a11 = A[0,0]
a12 = A[0,1]
a13 = A[0,2]
a21 = A[1,0]
a22 = A[1,1]
a23 = A[1,2]
a31 = A[2,0]
a32 = A[2,1]
a33 = A[2,2]
b11 = B[0,0]
b12 = B[0,1]
b13 = B[0,2]
b21 = B[1,0]
b22 = B[1,1]
b23 = B[1,2]
b31 = B[2,0]
b32 = B[2,1]
b33 = B[2,2]
#Now calculate the coffecients of the characteristic equation: det(L*A+B)
#
#Coefficients calculated using SymPy
#
#pen = sp.Matrix(L*A+B).evalf()
#char = pen.det()
#simp = collect(simplify(char),[L**3,L**2,L])
#
#Which yields the coefficients of L**3, L**2, L, and L**0
factor = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31)
a = (a11*a22*b33 - a11*a23*b32 - a11*a32*b23 + a11*a33*b22 - a12*a21*b33 + a12*a23*b31 + a12*a31*b23 - a12*a33*b21 + a13*a21*b32 - a13*a22*b31 - a13*a31*b22 + a13*a32*b21 + a21*a32*b13 - a21*a33*b12 - a22*a31*b13 + a22*a33*b11 + a23*a31*b12 - a23*a32*b11)/factor
b = (a11*b22*b33 - a11*b23*b32 - a12*b21*b33 + a12*b23*b31 + a13*b21*b32 - a13*b22*b31 - a21*b12*b33 + a21*b13*b32 + a22*b11*b33 - a22*b13*b31 - a23*b11*b32 + a23*b12*b31 + a31*b12*b23 - a31*b13*b22 - a32*b11*b23 + a32*b13*b21 + a33*b11*b22 - a33*b12*b21)/factor
c = (b11*b22*b33 - b11*b23*b32 - b12*b21*b33 + b12*b23*b31 + b13*b21*b32 - b13*b22*b31)/factor
#Now check for the overlap conditions
if a >= 0:
cond1 = -3*b+a**2
cond2 = 3*a*c+b*a**2-4*b**2
cond3 = -27*c**2 + 18*c*a*b + a**2*b**2 - 4*a**3*c - 4*b**3
if cond1 > 0 and cond2 < 0 and cond3 > 0:
return('Separated')
else:
return('Overlapping')
else:
cond1 = -3*b+a**2
cond2 = -27*c**2 + 18*c*a*b + a**2*b**2 - 4*a**3*c - 4*b**3
if cond1 > 0 and cond2 > 0:
return('Separated')
else:
return('Overlapping')
#Example parameters
#Let's make both ellipses identical
a = 3 #Ellipse axis (x)
b = 1 #ELlipse axis (y)
#Position and orientation of ellipse 1
x1 = 3
y1 = 0
theta1 = 0
#Position and orientation of ellipse 2
x2 = 1
y2 = 2
theta2 = np.pi/8
#Determine whether they overlap
out = detectOverlap((a,b),(a,b),(x1,y1,theta1),(x2,y2,theta2))
print(out)
```
For axis-aligned ellipses, a good approximation can be to check that the center point of the first ellipse is outside of the augmented ellipse centered on the second ellipse. The augmented ellipse being the ellipse with added width/height.
For first ellipse being centered on $(x_1, y_1)$ with width $a_1$ and height $b_1$ and second ellipse centered on $(x_2, y_2)$ with width $a_2$ and height $b_2$:
$$ \frac{(x_2 - x_1)^2}{(a_1 + a_2)^2} + \frac{(y_2 - y_1)^2}{(b_1 + b_2)^2} \leq 1 $$
For the first ellipse being on the edge of the augmented ellipse, see the following error estimation:
