查找兩個橢圓的交點（Python）

Question

我正在用Python編寫一個基本的二維形狀庫（主要用於處理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 

當橢圓是相同的忽略，有可能通過4個交叉點（沒有交點，正切，部分重疊的，部分重疊的是0和內部相切，並完全重疊）。

我已經發現了一些可能的解決方案：

• SymPy geometry module - 這基本上只是插頭橢圓方程爲SymPy的求解。如果沒有求解者，我不確定這是否合理。 （順便說一句，我會用SymPy代替滾動我自己的，但它在處理瘋狂的浮動時表現可怕）
• How to detect if an ellipse intersects(collides with) a circle - 這可能適用於兩個省略號，但我對如何打開它有點模糊變成明智的代碼。
• How Ellipse to Ellipse intersection? - 答案參考庫（CADEMIA）可能有一個很好的算法，但我甚至不知道它是否是開源的。
• Wikipedia: Intersecting Two Conics - 我沒有足夠的線性代數來理解這個解決方案。

任何關於如何計算交叉點的建議？速度（可能需要計算很多交叉點），優雅是主要標準。代碼會很棒，但即使是一個好的方向，也會有所幫助。

Answer 1

在數學中，您需要將橢圓表示爲二元二次方程，並解決它。我發現了一個doucument。所有的計算都在文檔中，但是它可能需要一段時間才能在Python中實現。

所以另一種方法是近似橢圓的折線，並使用勻稱找到路口，這裏是代碼：

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.5ms。

Answer 2

看着sympy我認爲它有你需要的一切。 （我想向您提供例如代碼;不幸的是，我沒有在與gmpy2安裝sympy而不是內置的數學無用蟒）

您有：

• 的ellipse與rotate方法，它可以是intersected與其他橢圓
• 或arbitrary intersection函數，它採用可變數字..他們稱之爲'幾何實體'。

從它們的實例，我認爲這是絕對有可能相交的兩個橢圓：

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

編輯：因爲一般橢圓（AX^2 + BX + CY^2 + DY + EXY + F）在sympy中不支持，

您應該自己構建方程並進行轉換，並使用它們的求解器找到相交點。

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