找到两个椭圆的交点(Python)

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.

我正在用Python编写一个基本的2D形状库(主要是用于操作SVG绘图)，并且我对于如何有效地计算两个省略号的交集点感到茫然。

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

忽略椭圆是相同的，可以有0到4个相交点(没有相交、切线、部分重叠、部分重叠和内部相切，以及完全重叠)。

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)
• 这基本上就是把椭圆方程插入到y的求解中。我不确定这是否有意义，而不是已经有了求解器。(顺便说一句，我会用“好”而不是“我自己”，但在处理疯狂的花车时它表现得很糟糕)
• 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.
• 椭圆与椭圆相交如何?-答案参考文献(CADEMIA)可能有一个很好的算法，但我甚至不知道它是否是开源的。
• 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.

关于如何计算交叉口，有什么建议吗?速度(可能需要计算很多交集)和优雅是主要的标准。代码会很棒，但即使是一个好的方向也会有帮助。

3 个解决方案

#1

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.

在数学中，你需要把椭圆表示成双变量二次方程，并求解。我找到了一个doucument。所有的计算都在文档中，但是在Python中实现它可能需要一段时间。

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

另一种方法是将椭圆近似为折线，用shapely来找到交集，这里是代码:

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:

和输出:

It takes about 1.5ms on my PC.

我的电脑大约需要1.5毫秒。

#2

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)

看着很好，我认为它有你需要的一切。(我试着给你提供示例代码;不幸的是，我没能通过使用gmpy2而不是无用的python内置数学来安装“症状”。

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,

编辑:由于一般椭圆(ax ^ 2 + bx + cy ^ 2 + dy + exy + f)在sympy不支持,

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

你应该自己建立方程，并把它们转化，用它们的求解器找到交点。

#3

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

您可以使用这里所示的方法:https://math.stackexchange.com/questions/864070/ howto - default -if-two- - - - - -至少是一个intersectionpoint/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()