### A closed form solution for the intersections of two ellipses

There intersections of two ellipses are determined using the following steps

- Rewrite each ellipse as a second degree equations in x and y coordinates
- Convert the intersection of the two second degree equations into a fourth degree polynomial which can be solved using closed form methods
- Compute the x values for each y.

In theory all these steps can be computed using closed form solutions.

- The ellipse centered at (a,c) with axes determined by coefficients b and d is specified:

(x’ – a)^{2}/b + (y’ – c)^{2}/d = 1

In addition, the xy coordinates can be rotated by

x’ = x * cos(q) – y * sin(q)

y’ = x * sin(q) + y * cos(q)

To simplify the notation, set A = cos(q) and B = sin(q). Substituting for x’ and y’ produces the general formula for an ellipse

(A * x –B * y – a)^{2}/b + (B * x + A * y – c)^{2}/d = 1

This can be rewritten

(A^{2}/b+B^{2}/d)x^{2} + (-2AB/b+2AB/d)xy + (B^{2}/b+A^{2}/d)y^{2} + (-2aA/b-2cB/d)x + (2aB/b-2cA/d)y + a^{2}/b+c^{2}/d-1 = 0

This method first translates and then rotates the ellipse.

An alternative method first rotates the ellipse then recenters it. With A = cos(q) and B = sin(q), the polynomial coefficients for the initial rotation are

AA = A^{2}/b + B^{2}/d

BB = -2*A*B/b + 2*A*B/d

CC = B^{2}/b + A^{2}/d

FF = -1

in the equation AA”x^{2} + BBx”y” + CCy”^{2} + F = 0.

By substituting x” = x – a and y” = y – c to shift the center of the ellipse and using the rotates variables AA, BB, CC, we get the polynomial describing the ellipse

AAx^{2 }+ BBxy + CCy^{2 }-(2AA a + BB c)x – (BB a + 2CC c)y + AA a^{2} + BB ac + CC c^{2} – 1 = 0

We have now standardized the equations defining the individual ellipses.

- The two ellipses are now written as quadratic conic sections. Their xy coordinates can be solved using a solution from http://www.math.niu.edu/~rusin/known-math/99/2ellipses. This formula solves for y values using the coefficients for each ellipse determined in the previous step. The coefficients for each ellipse, f1 and f2, are identified in the following equalities

f1 : a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0

f2 : a1*x^2 + b1*x*y + c1*y^2 + d1*x + e1*y + fq = 0

The solution is the zeros of a fourth degree polynomial z0 + z1 * y + z2 * y^{2} + z3 * y^{3} + z4 * y^{4} = 0

with the following coefficients derived from the expressions of f1 and f2,

z0 = f*a*d1^2+a^2*fq^2-d*a*d1*fq+a1^2*f^2-2*a*fq*a1*f-d*d1*a1*f+a1*d^2*fq

z1 = e1*d^2*a1-fq*d1*a*b-2*a*fq*a1*e-f*a1*b1*d+2*d1*b1*a*f+2*e1*fq*a^2+d1^2*a*e-e1*d1*a*d-

2*a*e1*a1*f-f*a1*d1*b+2*f*e*a1^2-fq*b1*a*d-e*a1*d1*d+2*fq*b*a1*d

z2 = e1^2*a^2+2*c1*fq*a^2-e*a1*d1*b+fq*a1*b^2-e*a1*b1*d-fq*b1*a*b-2*a*e1*a1*e+

2*d1*b1*a*e-c1*d1*a*d-2*a*c1*a1*f+b1^2*a*f+2*e1*b*a1*d+e^2*a1^2-c*a1*d1*d-

e1*b1*a*d+2*f*c*a1^2-f*a1*b1*b+c1*d^2*a1+d1^2*a*c-e1*d1*a*b-2*a*fq*a1*c

z3 = -2*a*a1*c*e1+e1*a1*b^2+2*c1*b*a1*d-c*a1*b1*d+b1^2*a*e-e1*b1*a*b-2*a*c1*a1*e-

e*a1*b1*b-c1*b1*a*d+2*e1*c1*a^2+2*e*c*a1^2-c*a1*d1*b+2*d1*b1*a*c-c1*d1*a*b

z4 = a^2*c1^2-2*a*c1*a1*c+a1^2*c^2-b*a*b1*c1-b*b1*a1*c+b^2*a1*c1+c*a*b1^2

There are analytic solutions to this equation which produce four solutions in the complex plane. Only the solutions with zero imaginary components are of interest when computing the intersection points.

- The corresponding x values are computed by

x = -(a*fq+a*c1*y^2-a1*c*y^2+a*e1*y-a1*e*y-a1*f)/(a*b1*y+a*d1-a1*b*y-a1*d)

There are, however, cases where the denominator is zero when the main axes of the ellipses are horizontal and vertical. In this case we can solve for x using the formula

bb = b * y + d

cc = c * y^2 + e * y + f

x = (-bb + sqrt(bb^{2} – 4 * a * cc)) / (2 * a) or x = (-bb – sqrt(bb^{2} – 4 * a * cc)) / (2 * a)

There are two intersection points with identical y values and positive and negative solutions for x.

By solving for the intersections of two ellipses, we have also solved for the intersections of any two conic sections showing there are a maximum of four intersections for any two distinct conic figures.

##### Comments

**9 Responses to “A closed form solution for the intersections of two ellipses”**

There sure seems to be a lot of Step 1’s here.

Wrt

f1 : a*x^2 + b*x*y + c*y^2 + d*x + e*y + f = 0

f2 : a1*x^2 + b1*x*y + c1*y^2 + d1*x + e1*y + fq = 0

There must be some subtlety that I don’t get. Why is the last term called “fq” and not “f1”.

This also feeds into another problem I’m having. Do we “shift” both ellipses by their a and c coefficients? I suspect not, but then I do not see how to get both of them into their standard function forms of f1 and f2.

Cluelessly,

/Bob

Since the equations were numbered f1 and f2, the constant coefficient for the second equations needs to be identified using another label. fq is used in this case.

The equations are the most general bivariate quadratic form, so they are a general solution for the intersection of two conic sections.

Well, that clears it all up!

Here’s a living, breathing demo where you can look at the (Javascript) code that implements the technique discussed in this page: https://www.khanacademy.org/cs/ellipse-collision-detector/5514890244521984

This is an old post but maybe you’ll read it.

First of all thanks for writing this down. I got this to work in code.

But this fails for me when intersecting two congruent ellipses (same radiusses (what’s the word?)) and same rotation, different origins), as all quartics are basically zero (I check against an epsilon of 10^-7). Is this expected?

The good thing is, that this special case can be reduced to an intersection between a line and one of the ellipses.

Let the two ellipses have the axes determined by the coefficients b and d. Let the origin of the first ellipse be at and the origin of the second ellipse be at . The rotation of both is q.

The line in hesse normal form x * n = d is easily calculated.

n is the vector . Normalize that to length 1.

d = Z – W

U = -2 * AA * a + BB * c

V = BB * a + 2 * CC * c

W = AA * a^2 + BB * a * c + CC * c^2

X = -2 * AA * o + BB * p

Y = BB * o + 2 * CC * p

Z = AA * o^2 + BB * o * p + CC * p^2

This also applies when abs(q1 – q2) = PI/2 and the axis coefficients b1 = d1 and b2 = d2.

These could be further simplified by ignoring the rotation altogether and later rotate the resulting points. You’d need to rotate one origin around the other.

Hello. Thanks for all the work that is in this page. It has been very helpful in my own work.

I’m trying to implement it as a DEM, but I get errors that don’t know where they come from. Maybe your solution serves me, but I can not understand it well.

Does the vector defines itself as tangent to one of the ellipses?

What is p?

An example of error is as follows:

a: 0.006195215690411591189

b: 0.0038067195103558847981

c: .0065136915639502088529

d: -5.1654958580412948521

e: -3.0844415930179378549

f: 1170.2827541311000914

a1: .0046922628185795929937

b1: -0.0018817617213015217902

c1: .0080166444357822070482

d1: -3.5625973379117654716

e1: -1.0678944804070475971

fq: 776.67002592325184196

z0: 0.12093560247970655541

z1: -0.0053551065513971329163

z2: 8.7342719462466861984e-005

z3: -6.212609087657343584e-007

z4: 1.6275442949293521071e-009

Drawing these ellipses, is that they are contacting in only two points, but the result of the quartic throws me 4 real roots. The quartic solver is fine…, i’ve contrasted it with this numbers and have the same result.

The difference in the angle of these ellipses is not pi / 2 and in i’ve not experimenting any problem to find right solutions for z coefficients smalls as this are.

I’ve also checked with this part of code from: https://www.khanacademy.org/computer-programming/dr-daves-ellipse-collision-detector/5776508379463680.

* Then get the differnce of the two equations.

var deltaB = (b1 / a1) – (b / a);

var deltaC = (c1 / a1) – (c / a);

var deltaD = (d1 / a1) – (d / a);

var deltaE = (e1 / a1) – (e / a);

var deltaF = (f1 / a1) – (f / a);

* Special case for b’s and d’s being equal */

if (deltaB === 0 && deltaD === 0) {

return realRoot(0, 0, deltaC, deltaE, deltaF);

I don’t understand this, but is not the solution for the case i’ve already posted, where i get deltaB and delta D not even close to 0 (deltaB : -1.015496196358965042 and deltaD : 74.538579851653636865).

@Alfredo: (Can’t seem to reply to your post directly).

ATTENTION: Be careful to only read my second post as the first got destroyed a bit by wordpress.

p is the y-Coordinate of one of the ellipses.

What “vector” do you mean? The vector n is the normal of the line that goes through the two intersection points of the ellipses.

n.x = (U-X)

n.y = (Y-V)

d is the distance of that line from the origin (0, 0)

d = Z – W

U, V, W, X, Y, Z are calculated as follows. (AA, BB, CC) are directly taken from the original blog post:

U = -2 * AA * a + BB * c

V = BB * a + 2 * CC * c

W = AA * a^2 + BB * a * c + CC * c^2

X = -2 * AA * o + BB * p

Y = BB * o + 2 * CC * p

Z = AA * o^2 + BB * o * p + CC * p^2

What is a DEM?

@Alfredo: Here’s my Typescript/Javascript-code for intersecting two ellipses: https://gist.github.com/drawable/92792f59b6ff8869d8b1

I do some things to make this more robust:

I check zeros against an epsilon of 10e-7. That is what the eq and is0 -function do.

I check if one of the ellipses has an axis parallel to the x- or y-axis. If that is the case I rotate the two ellipses by some angle. That saves me the checking of some special cases.