给定两个从计算中获得的数据数组,可以分别用连续可微函数来建模。这两个“线段”的数据点在(至少)一个点上相交,现在的问题是这两个数据集背后的函数实际上是交叉还是错位。
2、解决方案: 一种潜在的解决方案是:
查看一条曲线的一条线段与另一条曲线的一条线段的每种组合。 边界框有重叠的线段组合有可能包含交点。 求解一个线性方程组来计算两个直线是否相交以及相交位置。 如果方程组无解,则直线平行但不重叠,忽略这种情况。 如果有一个解,检查它是否真正落在线段内,如果是,则记录此交点。 如果存在无限多个交点,则直线相同。这也不是真正的交叉,可以忽略。 对所有线段组合执行此操作,并消除孪生情况,即当两条曲线在某一段的起点或终点相交时。
以下是基于提供的数据的代码示例:
import math
import matplotlib.pyplot as plt
import numpy as np
def intersect_curves(x1, y1, x2, y2):
N1 = x1.shape[0]
N2 = x2.shape[0]
xs1 = np.vstack((x1[:-1], x1[1:]))
ys1 = np.vstack((y1[:-1], y1[1:]))
xs2 = np.vstack((x2[:-1], x2[1:]))
ys2 = np.vstack((y2[:-1], y2[1:]))
mix1 = np.tile(np.amin(xs1, axis=0), (N2-1,1))
max1 = np.tile(np.amax(xs1, axis=0), (N2-1,1))
miy1 = np.tile(np.amin(ys1, axis=0), (N2-1,1))
may1 = np.tile(np.amax(ys1, axis=0), (N2-1,1))
mix2 = np.transpose(np.tile(np.amin(xs2, axis=0), (N1-1,1)))
max2 = np.transpose(np.tile(np.amax(xs2, axis=0), (N1-1,1)))
miy2 = np.transpose(np.tile(np.amin(ys2, axis=0), (N1-1,1)))
may2 = np.transpose(np.tile(np.amax(ys2, axis=0), (N1-1,1)))
idx = np.where((mix2 <= max1) & (max2 >= mix1) & (miy2 <= may1) & (may2 >= miy1))
x0 = []
y0 = []
for (i, j) in zip(idx[0], idx[1]):
xa = xs1[:, j]
ya = ys1[:, j]
xb = xs2[:, i]
yb = ys2[:, i]
a = np.array([[xa[1] - xa[0], xb[0] - xb[1]], [ya[1] - ya[0], yb[0]- yb[1]]])
b = np.array([xb[0] - xa[0], yb[0] - ya[0]])
r, residuals, rank, s = np.linalg.lstsq(a, b)
if rank == 2 and not residuals and r[0] >= 0 and r[0] < 1 and r[1] >= 0 and r[1] < 1:
if r[0] == 0 and r[1] == 0 and i > 0 and j > 0:
angle_a1 = math.atan2(ya[1] - ya[0], xa[1] - xa[0])
angle_b1 = math.atan2(yb[1] - yb[0], xb[1] - xb[0])
xa2 = xs1[:, j-1]
ya2 = ys1[:, j-1]
xb2 = xs2[:, i-1]
yb2 = ys2[:, i-1]
angle_a2 = math.atan2(ya2[0] - ya2[1], xa2[0] - xa2[1])
angle_b2 = math.atan2(yb2[0] - yb2[1], xb2[0] - xb2[1])
if angle_a2 < angle_a1:
h = angle_a1
angle_a1 = angle_a2
angle_a2 = h
if (angle_b1 > angle_a1 and angle_b1 < angle_a2 and (angle_b2 < angle_a1 or angle_b2 > angle_a2)) or \
((angle_b1 < angle_a1 or angle_b1 > angle_a2) and angle_b2 > angle_a1 and angle_b2 < angle_a2):
x0.append(xa[0])
y0.append(ya[0])
else:
x0.append(xa[0] + r[0] * (xa[1] - xa[0]))
y0.append(ya[0] + r[0] * (ya[1] - ya[0]))
return (x0, y0)
x1 = np.arange(-10, 10, .5)
x2 = x1
y1 = [np.absolute(x**3)+100*np.absolute(x) for x in x1]
y2 = [-np.absolute(x**3)-100*np.absolute(x) for x in x2][::-1]
plt.plot(x1, y1, 'b-o')
plt.plot(x2, y2, 'r-o')
x0, y0 = intersect_curves(x1, y1, x2, y2)
print('{} intersection points'.format(len(x0)))
plt.plot(x0, y0, 'go')
plt.show()
这种方法具有以下优点:
它不限于形式为 y=f(x) 的曲线,而是适用于二维中的任意曲线。 它还适用于自相交曲线。