寻找两个椭圆的交点

Question

我正在使用Python编写一个基本的2D形状库（主要用于操作SVG图形），但我不知道如何高效地计算两个椭圆的交点。

每个椭圆由以下变量定义（所有变量均为浮点数）：

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

忽略省略号相同的情况，可能存在0到4个交点（无交点、切点、部分重叠、内部切点和完全重叠）。

我找到了一些潜在的解决方案：

• SymPy几何模块 - 基本上只是将椭圆方程插入到 SymPy 的求解器中。如果没有已有的求解器，我不确定这是否有意义。（顺便说一下，我本来会用 SymPy 而不是自己写，但它在处理疯狂的浮点数时表现非常糟糕）
• 如何检测椭圆是否与圆相交？ - 这可能可以改为适用于两个椭圆，但我对如何将其转化为合理的代码有些模糊。
• 如何计算椭圆和椭圆的交点？ - 回答中提到的库（CADEMIA）可能有一个很好的算法，但我甚至无法确定它是否开源。
• 维基百科：两个圆锥曲线的交点 - 我对线性代数的掌握还不够，无法理解这个解决方案。

有没有关于如何计算相交的建议？速度（可能需要计算很多相交点）和优雅是主要考虑因素。提供代码是非常棒的，但甚至提供一个好的方向也会很有帮助。

- mrjogo

---

这可以简化为解决二元二次方程。 - thkang

大多数椭圆与圆的解决方案对您来说都不起作用，因为通常您需要通过参数化椭圆，然后找出其距离圆心的距离等于圆半径的位置。(如果您知道正确的相位以同步参数化两个椭圆，那么您可以这样做...但是我想，从我的经验来看，这并不比相交椭圆更容易...) - abarnert

此外，我认为这个问题应该发布在http://math.stackexchange.com上，而且实际上它是一个重复的问题，已经被迁移到那里了，http://math.stackexchange.com/questions/197982/calculate-intersection-of-two-ellipses。 - abarnert

这个可能不适合的唯一原因（也可能不是重复的）是因为这是专门用于编程的，所以类似于@HYRY在下面提供的折线解决方案可能不会被建议。 - mrjogo

Answer 1

在数学中，您需要将椭圆表示为二元二次方程，并解决它。我找到了一个文档。所有计算都在文档中，但在Python中实现可能需要一些时间。

因此，另一种方法是将椭圆近似为折线，并使用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])

并且输出结果:

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

- HYRY

---

我很担心我的两个选择是a)复杂或b)依赖项（特别是非Python）。为了这一个功能而需要整个库似乎太浪费了。 - mrjogo

Answer 2

看了一下 sympy，我认为它有你需要的一切。（我试图提供给你一些示例代码，但不幸的是，我无法使用无用的Python内置数学库安装gmpy2代替sympy）

你可以使用以下内容：

• 一个带有旋转方法的椭圆，可以与其他椭圆相交
• 或者是一个任意相交函数，它接受可变数量的...他们称之为“几何实体”。

从他们的示例中，我认为肯定可以相交两个椭圆：

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

编辑：由于Sympy不支持一般的椭圆（ax^2 + bx + cy^2 + dy + exy + f），您需要自己建立方程并进行转换，然后使用��们的求解器找到交点。

- thkang

---

注意：由于一般椭圆不受支持，因此椭圆的轴将不会被旋转。只有中心会旋转到新位置。因此，sympy无法找到一般椭圆的交点。 - HYRY

Sympy的问题不在于几何库，而在于系统建立在合理化具有许多小数点的浮点数时性能较差。 - mrjogo

Answer 3

您可以使用此处显示的方法：https://math.stackexchange.com/questions/864070/how-to-determine-if-two-ellipse-have-at-least-one-intersection-point/864186#864186 首先，您应该能够在一个方向上重新调整椭圆。为此，您需要计算椭圆的系数作为锥截面，重新调整大小，然后恢复椭圆的新几何参数：中心、轴、角度。

然后，您的问题就归结为查找椭圆到原点的距离。要解决这个问题，您需要进行一些迭代。以下是可能的自包含实现...

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

- Emanuele Paolini