Numpy和线交叉点

Question

我如何使用numpy来计算两个线段之间的交点？

在代码中我有segment1 =((x1,y1),(x2,y2))和segment2 =((x1,y1),(x2,y2)).注意段1不等于segment2.所以在我的代码中我也一直在计算斜率和y截距,如果可以避免这种情况会很好但我不知道如何.

我一直在使用Cramer的规则和我在Python中编写的函数,但我想找到一种更快的方法.

Answer 1

直接从http://www.cs.mun.ca/~rod/2500/notes/numpy-arrays/numpy-arrays.html窃取

#
# line segment intersection using vectors
# see Computer Graphics by F.S. Hill
#
from numpy import *
def perp( a ) :
    b = empty_like(a)
    b[0] = -a[1]
    b[1] = a[0]
    return b

# line segment a given by endpoints a1, a2
# line segment b given by endpoints b1, b2
# return 
def seg_intersect(a1,a2, b1,b2) :
    da = a2-a1
    db = b2-b1
    dp = a1-b1
    dap = perp(da)
    denom = dot( dap, db)
    num = dot( dap, dp )
    return (num / denom.astype(float))*db + b1

p1 = array( [0.0, 0.0] )
p2 = array( [1.0, 0.0] )

p3 = array( [4.0, -5.0] )
p4 = array( [4.0, 2.0] )

print seg_intersect( p1,p2, p3,p4)

p1 = array( [2.0, 2.0] )
p2 = array( [4.0, 3.0] )

p3 = array( [6.0, 0.0] )
p4 = array( [6.0, 3.0] )

print seg_intersect( p1,p2, p3,p4)

Answer 2

import numpy as np

def get_intersect(a1, a2, b1, b2):
    """ 
    Returns the point of intersection of the lines passing through a2,a1 and b2,b1.
    a1: [x, y] a point on the first line
    a2: [x, y] another point on the first line
    b1: [x, y] a point on the second line
    b2: [x, y] another point on the second line
    """
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)

if __name__ == "__main__":
    print get_intersect((0, 1), (0, 2), (1, 10), (1, 9))  # parallel  lines
    print get_intersect((0, 1), (0, 2), (1, 10), (2, 10)) # vertical and horizontal lines
    print get_intersect((0, 1), (1, 2), (0, 10), (1, 9))  # another line for fun

说明

注意,一条线的方程是ax+by+c=0.所以,如果一个点在这一行,那么它是一个解决方案(a,b,c).(x,y,1)=0(.是点积)

让l1=(a1,b1,c1),l2=(a2,b2,c2)是两行p1=(x1,y1,1),p2=(x2,y2,1)是两个点.

找到通过两点的线:

let t=p1xp2(两点的叉积)是表示直线的向量.

我们知道这p1是t因为t.p1 = (p1xp2).p1=0.我们也知道这p2是t因为t.p2 = (p1xp2).p2=0.所以t必须是通过p1和p2.

这意味着我们可以通过获取该线上两点的叉积来获得线的矢量表示.

找到交叉点:

现在让r=l1xl2(两条线的叉积)成为表示点的向量

我们知道r谎言是l1因为r.l1=(l1xl2).l1=0.我们也知道r谎言l2因为r.l2=(l1xl2).l2=0.所以,r必须在线路的交叉点l1和l2.

有趣的是,我们可以通过取两条线的叉积来找到交点.

Answer 3

这可能是一个迟到的反应,但这是我用谷歌搜索'numpy line intersection'时的第一个响应.在我的情况下,我在一个平面上有两条线,我想快速得到它们之间的任何交叉点,而Hamish的解决方案会很慢 - 需要在所有线段上嵌套for循环.

以下是没有for循环的方法(它非常快):

from numpy import where, dstack, diff, meshgrid

def find_intersections(A, B):

    # min, max and all for arrays
    amin = lambda x1, x2: where(x1<x2, x1, x2)
    amax = lambda x1, x2: where(x1>x2, x1, x2)
    aall = lambda abools: dstack(abools).all(axis=2)
    slope = lambda line: (lambda d: d[:,1]/d[:,0])(diff(line, axis=0))

    x11, x21 = meshgrid(A[:-1, 0], B[:-1, 0])
    x12, x22 = meshgrid(A[1:, 0], B[1:, 0])
    y11, y21 = meshgrid(A[:-1, 1], B[:-1, 1])
    y12, y22 = meshgrid(A[1:, 1], B[1:, 1])

    m1, m2 = meshgrid(slope(A), slope(B))
    m1inv, m2inv = 1/m1, 1/m2

    yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
    xi = (yi - y21)*m2inv + x21

    xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12), 
              amin(x21, x22) < xi, xi <= amax(x21, x22) )
    yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
              amin(y21, y22) < yi, yi <= amax(y21, y22) )

    return xi[aall(xconds)], yi[aall(yconds)]

然后使用它,提供两行作为参数,其中arg是一个2列矩阵,每行对应一个(x,y)点:

# example from matplotlib contour plots
Acs = contour(...)
Bsc = contour(...)

# A and B are the two lines, each is a 
# two column matrix
A = Acs.collections[0].get_paths()[0].vertices
B = Bcs.collections[0].get_paths()[0].vertices

# do it
x, y = find_intersections(A, B)

玩得开心

Answer 4

这是 @Hamish Grubijan 答案的一个版本，它也适用于每个输入参数中的多个点，即 、、 、a1可以a2是2D 点的 Nx2 行数组。该函数被点积取代。b1b2perp

T = np.array([[0, -1], [1, 0]])
def line_intersect(a1, a2, b1, b2):
    da = np.atleast_2d(a2 - a1)
    db = np.atleast_2d(b2 - b1)
    dp = np.atleast_2d(a1 - b1)
    dap = np.dot(da, T)
    denom = np.sum(dap * db, axis=1)
    num = np.sum(dap * dp, axis=1)
    return np.atleast_2d(num / denom).T * db + b1