A closed form solution for the intersections of two ellipses

ellipseIntersectionThere intersections of two ellipses are determined using the following steps

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

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

  1. 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

(A2/b+B2/d)x2 + (-2AB/b+2AB/d)xy + (B2/b+A2/d)y2 + (-2aA/b-2cB/d)x + (2aB/b-2cA/d)y + a2/b+c2/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 = A2/b + B2/d
BB = -2*A*B/b + 2*A*B/d
CC = B2/b + A2/d
FF = -1

in the equation AA”x2 + 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

AAx2 + BBxy + CCy2 -(2AA a + BB c)x – (BB a + 2CC c)y + AA a2 + BB ac + CC c2 – 1 = 0

We have now standardized the equations defining the individual ellipses.

  1. 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 * y2 + z3 * y3 + z4 * y4  = 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.

  1. 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(bb2 – 4 * a * cc)) / (2 * a) or x = (-bb – sqrt(bb2 – 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.

Advertisements
Comments
9 Responses to “A closed form solution for the intersections of two ellipses”
  1. Bob Lyon says:

    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

    • Elliot Noma says:

      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.

  2. Bob Lyon says:

    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

  3. Alfredo Lamadrid says:

    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).

  4. Stephan says:

    @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?

  5. Stephan says:

    @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.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: