我将如何使用 numpy 计算两条线段之间的交点?
在我的代码中segment1 = ((x1,y1),(x2,y2))
和segment2 = ((x1,y1),(x2,y2))
. 注意segment1
不等于segment2
。因此,在我的代码中,我也一直在计算斜率和 y 截距,如果可以避免这种情况会很好,但我不知道如何做。
我一直在将 Cramer 规则与我用 Python 编写的函数一起使用,但我想找到一种更快的方法。
我将如何使用 numpy 计算两条线段之间的交点?
在我的代码中segment1 = ((x1,y1),(x2,y2))
和segment2 = ((x1,y1),(x2,y2))
. 注意segment1
不等于segment2
。因此,在我的代码中,我也一直在计算斜率和 y 截距,如果可以避免这种情况会很好,但我不知道如何做。
我一直在将 Cramer 规则与我用 Python 编写的函数一起使用,但我想找到一种更快的方法。
#
# 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)
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)
是两点。
让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
有趣的是,我们可以通过两条线的叉积找到交点。
也许这是一个迟到的回应,但当我在谷歌上搜索“numpy line intersections”时,这是第一个命中。就我而言,我在一个平面上有两条线,我想快速获得它们之间的任何交点,而 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)
玩得开心
这是@Hamish Grubijan 的答案的一个版本,它也适用于每个输入参数中的多个点,即 , ,a1
可以a2
是二维点的 Nx2 行数组。该函数由点积代替。b1
b2
perp
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
这是一个(有点强迫的)单行:
import numpy as np
from scipy.interpolate import interp1d
x = np.array([0, 1])
segment1 = np.array([0, 1])
segment2 = np.array([-1, 2])
x_intersection = interp1d(segment1 - segment2, x)(0)
# if you need it:
y_intersection = interp1d(x, segment1)(x_intersection)
对差值进行插值(默认为线性),并找到一个 0 的倒数。
干杯!
我想在这里添加一些小东西。最初的问题是关于线段的。我来到这里,因为我正在寻找线段交点,这在我的情况下意味着我需要过滤那些不存在线段交点的情况。这是一些执行此操作的代码:
def line_intersection(x1, y1, x2, y2, x3, y3, x4, y4):
"""find the intersection of line segments A=(x1,y1)/(x2,y2) and
B=(x3,y3)/(x4,y4). Returns a point or None"""
denom = ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))
if denom==0: return None
px = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) / denom
py = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) / denom
if (px - x1) * (px - x2) < 0 and (py - y1) * (py - y2) < 0 \
and (px - x3) * (px - x4) < 0 and (py - y3) * (py - y4) < 0:
return [px, py]
else:
return None
如果您正在寻找我们可以排除垂直线段的矢量化版本。
def intersect(a):
# a numpy array with dimension [n, 2, 2, 2]
# axis 0: line-pair, axis 1: two lines, axis 2: line delimiters axis 3: x and y coords
# for each of the n line pairs a boolean is returned stating of the two lines intersect
# Note: the edge case of a vertical line is not handled.
m = (a[:, :, 1, 1] - a[:, :, 0, 1]) / (a[:, :, 1, 0] - a[:, :, 0, 0])
t = a[:, :, 0, 1] - m[:, :] * a[:, :, 0, 0]
x = (t[:, 0] - t[:, 1]) / (m[:, 1] - m[:, 0])
y = m[:, 0] * x + t[:, 0]
r = a.min(axis=2).max(axis=1), a.max(axis=2).min(axis=1)
return (x >= r[0][:, 0]) & (x <= r[1][:, 0]) & (y >= r[0][:, 1]) & (y <= r[1][:, 1])
示例调用将是:
intersect(np.array([
[[[1, 2], [2, 2]],
[[1, 2], [1, 1]]], # I
[[[3, 4], [4, 4]],
[[4, 4], [5, 6]]], # II
[[[2, 0], [3, 1]],
[[3, 0], [4, 1]]], # III
[[[0, 5], [2, 5]],
[[2, 4], [1, 3]]], # IV
]))
# returns [False, True, False, False]
可视化(我需要更多的声誉才能在这里发布图片)。
这就是我用来查找线交点的方法,它可以使用每条线的 2 个点,或者只有一个点及其斜率。我基本上解决了线性方程组。
def line_intersect(p0, p1, m0=None, m1=None, q0=None, q1=None):
''' intersect 2 lines given 2 points and (either associated slopes or one extra point)
Inputs:
p0 - first point of first line [x,y]
p1 - fist point of second line [x,y]
m0 - slope of first line
m1 - slope of second line
q0 - second point of first line [x,y]
q1 - second point of second line [x,y]
'''
if m0 is None:
if q0 is None:
raise ValueError('either m0 or q0 is needed')
dy = q0[1] - p0[1]
dx = q0[0] - p0[0]
lhs0 = [-dy, dx]
rhs0 = p0[1] * dx - dy * p0[0]
else:
lhs0 = [-m0, 1]
rhs0 = p0[1] - m0 * p0[0]
if m1 is None:
if q1 is None:
raise ValueError('either m1 or q1 is needed')
dy = q1[1] - p1[1]
dx = q1[0] - p1[0]
lhs1 = [-dy, dx]
rhs1 = p1[1] * dx - dy * p1[0]
else:
lhs1 = [-m1, 1]
rhs1 = p1[1] - m1 * p1[0]
a = np.array([lhs0,
lhs1])
b = np.array([rhs0,
rhs1])
try:
px = np.linalg.solve(a, b)
except:
px = np.array([np.nan, np.nan])
return px
为了解决这个问题,我们必须将我们的行转换为以下形式:ax+by=c。在哪里
a = y1 - y2
b = x1 - x2
c = ax1 + by1
如果我们对每条线应用这个方程,我们将得到两条线方程。a1x+b1y=c1和a2x+b2y=c2。
现在,当我们得到两行的表达式时。
首先,我们必须检查线是否平行。为了检验这一点,我们想找到行列式。如果行列式等于零,则这些线是平行的。
我们通过求解以下表达式找到行列式:
det = a1 * b2 - a2 * b1
如果行列式等于 0,那么这些线是平行的并且永远不会相交。如果线不平行,它们必须在某个点相交。
使用以下公式找到线的相交点:
class Point:
def __init__(self, x, y):
self.x = x
self.y = y
'''
finding intersect point of line AB and CD
where A is the first point of line AB
and B is the second point of line AB
and C is the first point of line CD
and D is the second point of line CD
'''
def get_intersect(A, B, C, D):
# a1x + b1y = c1
a1 = B.y - A.y
b1 = A.x - B.x
c1 = a1 * (A.x) + b1 * (A.y)
# a2x + b2y = c2
a2 = D.y - C.y
b2 = C.x - D.x
c2 = a2 * (C.x) + b2 * (C.y)
# determinant
det = a1 * b2 - a2 * b1
# parallel line
if det == 0:
return (float('inf'), float('inf'))
# intersect point(x,y)
x = ((b2 * c1) - (b1 * c2)) / det
y = ((a1 * c2) - (a2 * c1)) / det
return (x, y)
我为 line 编写了一个模块来计算这个和其他一些简单的 line 操作。它是用 C++ 实现的,所以运行速度非常快。您可以通过 pip 安装FastLine,然后以这种方式使用它:
from FastLine import Line
# define a line by two points
l1 = Line(p1=(0,0), p2=(10,10))
# or define a line by slope and intercept
l2 = Line(m=0.5, b=-1)
# compute intersection
p = l1.intersection(l2)
# returns (-2.0, -2.0)