寻找椭圆和直线的交点答案

【问题标题】：Finding intersection between ellipse and a line寻找椭圆和直线的交点
【发布时间】：2013-08-21 14:45:16
【问题描述】：

我试图找到椭圆和直线之间的交点，但似乎无法做到。我尝试了两种方法，一种是通过尝试在下面的代码中找到 LineString 和 LinearRing 之间的交集，但没有从中得到任何可用的值。问题之一是椭圆总是偏离中心并且角度很小或很大

  # -*- coding: utf-8 -*-
"""
Created on Mon Aug 19 17:38:55 2013

@author: adudchenko
"""
from pylab import *
import numpy as np
from shapely.geometry.polygon import LinearRing
from shapely.geometry import LineString
def ellipse_polyline(ellipses, n=100):
    t = np.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, line):
    ea = LinearRing(a)
    eb = LinearRing(b)
    mp = ea.intersection(eb)
    print mp
    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)
line=LineString([[0,0],[4,4]])
x, y = intersections(a, line)
figure()
plot(x, y, "o")
plot(a[:,0], a[:,1])
plot(b[:,0], b[:,1])
show()

我也尝试使用 fsolve，如下例所示，但它找到了错误的交叉点（或者实际上是一个错误的点。

from pylab import *
from scipy.optimize import fsolve
import numpy as np

def ellipse_polyline(ellipses, n=100):
    t = np.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 * np.cos(t) - b * sa * np.sin(t)
        p[:, 1] = y0 + a * sa * np.cos(t) + b * ca * np.sin(t)
        result.append(p)
    return result
def ellipse_line(txy):
    t,x,y=txy
    x0, y0, a, b, angle,m,lb=1, 1, 2, 1, 45,0.5,0
    sa = np.sin(angle)
    ca = np.cos(angle)
    return (x0 + a * ca * np.cos(t) - b * sa * np.sin(t)-x,y0 + a * sa * np.cos(t) + b * ca * np.sin(t)-y,m*x+lb-y) 
a,b= ellipse_polyline([(1, 1, 2, 1, 45), (2, 0.5, 5, 1.5, -30)])
t,y,x=fsolve(ellipse_line,(0,0,0))
print t,y,x
#print a[:,0]
m=0.5
bl=0
xl,yl=[],[]
for i in range(10):
    xl.append(i)
    yl.append(m*i+bl)
figure()
plot(x, y, "o")
plot(a[:,0], a[:,1])
plot(xl,yl)

有什么帮助吗？

【问题讨论】：

似乎是一个非常适合math.stackexchange.com的问题
我不这么认为。假设这不能直接解析，因此您需要一种数值方法，在这种情况下，这是发布的正确位置。否则 math.stackexchange.com 是要走的路。
@usethedeathstar...好点！我没有这样想过……
An answer is here on math.stackexchange.com 获取基本解决方案或Wolfam Mathworld "Ellipse-Line" intersection。 Python 也不是，但快速编辑应该可以了解它的核心。我今天将在 Python 中解决这个商业问题并将其发布（因为这在 G**gle 搜索中很早就出现了）。

标签： python numpy shapely

【解决方案1】：

使用 Shapely 部分，可以修复 def intersections(a, line)。首先，有一个错误，它引用了全局b 而不是使用line 参数，该参数被忽略。因此，比较是在两个椭圆 a 和 b 之间进行，而不是在每个椭圆与 line 之间进行比较，正如我认为的那样。

此外，两条线之间的交点可以是以下几种结果之一：空集、点（1 个交点）、多点（超过 1 个交点）、线串（如果线串的任何部分重叠并且彼此平行），或点和线的集合。假设只有前三个结果：

def intersections(a, line):
    ea = LinearRing(a)
    mp = ea.intersection(line)
    if mp.is_empty:
        print('Geometries do not intersect')
        return [], []
    elif mp.geom_type == 'Point':
        return [mp.x], [mp.y]
    elif mp.geom_type == 'MultiPoint':
        return [p.x for p in mp], [p.y for p in mp]
    else:
        raise ValueError('something unexpected: ' + mp.geom_type)

所以现在这些看起来是正确的：

>>> intersections(a, line)
([2.414213562373095], [2.414213562373095])
>>> intersections(b, line)
([0.0006681263405436677, 2.135895843256409], [0.0006681263405436642, 2.135895843256409])

【讨论】：