Finding intersection points of two ellipses

Question

I'm writing a basic 2D shape library in Python (primarily for manipulating SVG drawings), and I'm at a loss for how to efficiently calculate the intersection points of two ellipses.

Each ellipse is defined by the following variables (all floats):

c: center point (x, y)
hradius: "horizontal" radius
vradius: "vertical" radius
phi: rotation from coordinate system's x-axis to ellipse's horizontal axis

Ignoring when the ellipses are identical, there could be 0 through 4 intersection points (no intersection, tangent, partially overlapping, partially overlapping and internally tangent, and fully overlapping).

I've found a few potential solutions:

SymPy geometry module - This basically just plugs the ellipse equations into SymPy's solver. I'm not sure whether this makes sense without already having the solver. (Incidentally, I would have used SymPy instead of rolling my own, but it performs horribly when dealing with crazy floats)
How to detect if an ellipse intersects(collides with) a circle - This could probably be adapted for two ellipses, but I'm a little fuzzy on how to turn it into sensible code.
How Ellipse to Ellipse intersection? - The library the answer references (CADEMIA) might have a good algorithm, but I can't even figure out if it's open source.
Wikipedia: Intersecting Two Conics - I don't have enough of a grasp of linear algebra to understand this solution.

Any suggestions on how I should go about calculating the intersections? Speed (it might have to calculate a lot of intersections) and elegance are the primary criteria. Code would be fantastic, but even a good direction to go in would be helpful.

this can be reduced into solving two-variable quadratic equations — thkang, Mar 16 '13 at 04:36
Most ellipse-vs.-circle solutions won't work for you, because generally you do that by parameterizing the ellipse and then just finding out where its distance from the circle's center equals the circle's radius. (If you knew the right phase to parameterize the two ellipses in lockstep, you could do them… but off the top of my head, I suspect that's no easier than just intersecting the ellipses…) — abarnert, Mar 16 '13 at 06:23
Also, I think this belongs on http://math.stackexchange.com instead of here—and in fact it's a dup of a question that was already migrated there, http://math.stackexchange.com/questions/197982/calculate-intersection-of-two-ellipses — abarnert, Mar 16 '13 at 06:30
The only reason it might not belong there (and might not be a dup) is since this is explicitly for programming, polyline solutions such as @HYRY's below probably would not be suggested — mrjogo, Mar 23 '13 at 22:36

score 14 · Accepted Answer · answered Mar 16 '13 at 06:56

In math, you need to express the ellipses as bivariate quadratic equations, and solve it. I found a doucument. All the calculations are in the document, but it may take a while to implement it in Python.

So another method is to approximate the ellipses as polylines, and use shapely to find the intersections, here is the code:

import numpy as np
from shapely.geometry.polygon import LinearRing

def ellipse_polyline(ellipses, n=100):
    t = linspace(0, 2*np.pi, n, endpoint=False)
    st = np.sin(t)
    ct = np.cos(t)
    result = []
    for x0, y0, a, b, angle in ellipses:
        angle = np.deg2rad(angle)
        sa = np.sin(angle)
        ca = np.cos(angle)
        p = np.empty((n, 2))
        p[:, 0] = x0 + a * ca * ct - b * sa * st
        p[:, 1] = y0 + a * sa * ct + b * ca * st
        result.append(p)
    return result

def intersections(a, b):
    ea = LinearRing(a)
    eb = LinearRing(b)
    mp = ea.intersection(eb)

    x = [p.x for p in mp]
    y = [p.y for p in mp]
    return x, y

ellipses = [(1, 1, 2, 1, 45), (2, 0.5, 5, 1.5, -30)]
a, b = ellipse_polyline(ellipses)
x, y = intersections(a, b)
plot(x, y, "o")
plot(a[:,0], a[:,1])
plot(b[:,0], b[:,1])

and the output:

enter image description here

It takes about 1.5ms on my PC.

I was worried my two options were a) complicated or b) dependencies (particularly non-Python). It seems such a shame to have to require an entire library for this one feature. — mrjogo, Mar 23 '13 at 22:40

thkang · Answer 2 · 2013-03-16T06:05:50.110

looking at sympy I thinks it has everything you need. (I tried to provide you with example codes; unfortunately, I failed at installing sympy with gmpy2 instead of useless python built-in mathematics)

you have :

an ellipse with rotate method, which can be intersected with other ellipses
or an arbitrary intersection function that takes variadic number of.. what they call as 'Geometric Entities'.

from their examples, I think it is definitely possible to intersect two ellipses:

>>> from sympy import Ellipse, Point, Line, sqrt
>>> e = Ellipse(Point(0, 0), 5, 7)
...
>>> e.intersection(Ellipse(Point(1, 0), 4, 3))
[Point(0, -3*sqrt(15)/4), Point(0, 3*sqrt(15)/4)]
>>> e.intersection(Ellipse(Point(5, 0), 4, 3))
[Point(2, -3*sqrt(7)/4), Point(2, 3*sqrt(7)/4)]
>>> e.intersection(Ellipse(Point(100500, 0), 4, 3))
[]
>>> e.intersection(Ellipse(Point(0, 0), 3, 4))
[Point(-363/175, -48*sqrt(111)/175), Point(-363/175, 48*sqrt(111)/175),
Point(3, 0)]
>>> e.intersection(Ellipse(Point(-1, 0), 3, 4))
[Point(-17/5, -12/5), Point(-17/5, 12/5), Point(7/5, -12/5),
Point(7/5, 12/5)]

edit : since general ellipse (ax^2 + bx + cy^2 + dy + exy + f) is not supported in sympy,

you should build equations and transform them yourself, and use their solver to find intersection points.

Note: since the general ellipse is not supported, the axes of the ellipse will not be rotated. Only the center is rotated to a new position. So, general ellipse intersection can't be found by sympy. — HYRY, Mar 16 '13 at 05:54
The problem with sympy wasn't in the geometry libraries, it was in the poor performance rationalizing floats with a lot of decimal points which the system was built on — mrjogo, Mar 23 '13 at 22:32

score 3 · Answer 3 · edited Apr 13 '17 at 12:19

You can use the method shown here: https://math.stackexchange.com/questions/864070/how-to-determine-if-two-ellipse-have-at-least-one-intersection-point/864186#864186

First you should be able to rescale an ellipse in one direction. To do this you need to compute the coefficients of the ellipse as a conic section, rescale, and then recover the new geometric parameters of the ellipse: center, axes, angle.

Then your problem reduces to that of finding the distance from an ellipse to the origin. To solve this last problem you need some iteration. Here is a possible self contained implementation...

from math import *

def bisect(f,t_0,t_1,err=0.0001,max_iter=100):
    iter = 0
    ft_0 = f(t_0)
    ft_1 = f(t_1)
    assert ft_0*ft_1 <= 0.0
    while True:
        t = 0.5*(t_0+t_1)
        ft = f(t)
        if iter>=max_iter or ft<err:
            return t
        if ft * ft_0 <= 0.0:
            t_1 = t
            ft_1 = ft
        else:
            t_0 = t
            ft_0 = ft
        iter += 1

class Ellipse(object):
    def __init__(self,center,radius,angle=0.0):
        assert len(center) == 2
        assert len(radius) == 2
        self.center = center
        self.radius = radius
        self.angle = angle

    def distance_from_origin(self):
        """                                                                           
        Ellipse equation:                                                             
        (x-center_x)^2/radius_x^2 + (y-center_y)^2/radius_y^2 = 1                     
        x = center_x + radius_x * cos(t)                                              
        y = center_y + radius_y * sin(t)                                              
        """
        center = self.center
        radius = self.radius

        # rotate ellipse of -angle to become axis aligned                             

        c,s = cos(self.angle),sin(self.angle)
        center = (c * center[0] + s * center[1],
                  -s* center[0] + c * center[1])

        f = lambda t: (radius[1]*(center[1] + radius[1]*sin(t))*cos(t) -
                       radius[0]*(center[0] + radius[0]*cos(t))*sin(t))

        if center[0] > 0.0:
            if center[1] > 0.0:
                t_0, t_1 = -pi, -pi/2
            else:
                t_0, t_1 = pi/2, pi
        else:
            if center[1] > 0.0:
                t_0, t_1 = -pi/2, 0
            else:
                t_0, t_1 = 0, pi/2

        t = bisect(f,t_0,t_1)
        x = center[0] + radius[0]*cos(t)
        y = center[1] + radius[1]*sin(t)
        return sqrt(x**2 + y**2)

print Ellipse((1.0,-1.0),(2.0,0.5)).distance_from_origin()

Finding intersection points of two ellipses

3 Answers3

Linked