8

私はPythonで基本的な2D形状ライブラリ（主にSVG図面を操作するため）を書いていますが、2つの楕円の交点を効率的に計算する方法に迷っています。

各楕円は、次の変数（すべて浮動小数点数）によって定義されます。

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の交点が存在する可能性があります（交点なし、接線、部分的に重なり合う、部分的に重なり合う、内部的に接する、完全に重なる）。

私はいくつかの潜在的な解決策を見つけました：

交差点の計算方法について何か提案はありますか？速度（多くの交差点を計算する必要があるかもしれません）と優雅さが主要な基準です。コードは素晴らしいでしょうが、進むべき良い方向性でさえも役に立ちます。

4

3 に答える 3

13

数学では、楕円を2変量2次方程式として表現し、それを解く必要があります。文書を見つけました。すべての計算はドキュメントに含まれていますが、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])

および出力：

ここに画像の説明を入力してください

私のPCでは約1.5msかかります。

于 2013-03-16T06:56:25.537 に答える
3

sympyを見ると、必要なものがすべて揃っていると思います。（私はあなたにサンプルコードを提供しようとしました;残念ながら、私は役に立たないpython組み込み数学の代わりにgmpy2でsympyをインストールすることに失敗しました）

あなたが持っている ：

  • 他の楕円と交差させることができる回転方式の楕円
  • または、 「幾何学的エンティティ」と呼ばれるものの可変個引数をとる任意の交差関数。

彼らの例から、2つの楕円を交差させることは間違いなく可能だと思います。

>>> 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ではサポートされていないため、

方程式を作成して自分で変換し、ソルバーを使用して交点を見つける必要があります。

于 2013-03-16T05:22:17.197 に答える
0

ここに示す方法を使用できます： 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()
于 2014-07-12T07:01:57.130 に答える