Пусть 0 <= x <= 1. У меня два столбца f и g длины 5000 соответственно. Теперь я рисую:
plt.plot(x, f, '-')
plt.plot(x, g, '*')
Я хочу найти точку "x", где пересекает кривая. Я не хочу найти пересечение f и g. Я могу сделать это просто:
set(f) & set(g)
Вы можете использовать np.sign в сочетании с np.diff и np.argwhere для получения индексов точек, где пересекаются линии (в этом случае точки [ 0, 149, 331, 448, 664, 743]):
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0, 1000)
f = np.arange(0, 1000)
g = np.sin(np.arange(0, 10, 0.01) * 2) * 1000
plt.plot(x, f, '-')
plt.plot(x, g, '-')
idx = np.argwhere(np.diff(np.sign(f - g))).flatten()
plt.plot(x[idx], f[idx], 'ro')
plt.show()
Сначала он вычисляет f - g и соответствующие знаки, используя np.sign. Применение np.diff показывает все позиции, где знак изменяется (например, пересекаются линии). Используя np.argwhere дает нам точные индексы.
Ну, я искал matplotlib для двух кривых, которые отличались по размеру и не имели одинаковых значений x. Вот что я придумал:
import numpy as np
import matplotlib.pyplot as plt
import sys
fig = plt.figure()
ax = fig.add_subplot(111)
# x1 = [1,2,3,4,5,6,7,8]
# y1 = [20,100,50,120,55,240,50,25]
# x2 = [3,4,5,6,7,8,9]
# y2 = [25,200,14,67,88,44,120]
x1=[1.4,2.1,3,5.9,8,9,12,15]
y1=[2.3,3.1,1,3.9,8,9,11,9]
x2=[1,2,3,4,6,8,9,12,14]
y2=[4,12,7,1,6.3,7,5,6,11]
ax.plot(x1, y1, color='lightblue',linewidth=3, marker='s')
ax.plot(x2, y2, color='darkgreen', marker='^')
y_lists = y1[:]
y_lists.extend(y2)
y_dist = max(y_lists)/200.0
x_lists = x1[:]
x_lists.extend(x2)
x_dist = max(x_lists)/900.0
division = 1000
x_begin = min(x1[0], x2[0]) # 3
x_end = max(x1[-1], x2[-1]) # 8
points1 = [t for t in zip(x1, y1) if x_begin<=t[0]<=x_end] # [(3, 50), (4, 120), (5, 55), (6, 240), (7, 50), (8, 25)]
points2 = [t for t in zip(x2, y2) if x_begin<=t[0]<=x_end] # [(3, 25), (4, 35), (5, 14), (6, 67), (7, 88), (8, 44)]
# print points1
# print points2
x_axis = np.linspace(x_begin, x_end, division)
idx = 0
id_px1 = 0
id_px2 = 0
x1_line = []
y1_line = []
x2_line = []
y2_line = []
xpoints = len(x_axis)
intersection = []
while idx < xpoints:
# Iterate over two line segments
x = x_axis[idx]
if id_px1>-1:
if x >= points1[id_px1][0] and id_px1<len(points1)-1:
y1_line = np.linspace(points1[id_px1][1], points1[id_px1+1][1], 1000) # 1.4 1.401 1.402 etc. bis 2.1
x1_line = np.linspace(points1[id_px1][0], points1[id_px1+1][0], 1000)
id_px1 = id_px1 + 1
if id_px1 == len(points1):
x1_line = []
y1_line = []
id_px1 = -1
if id_px2>-1:
if x >= points2[id_px2][0] and id_px2<len(points2)-1:
y2_line = np.linspace(points2[id_px2][1], points2[id_px2+1][1], 1000)
x2_line = np.linspace(points2[id_px2][0], points2[id_px2+1][0], 1000)
id_px2 = id_px2 + 1
if id_px2 == len(points2):
x2_line = []
y2_line = []
id_px2 = -1
if x1_line!=[] and y1_line!=[] and x2_line!=[] and y2_line!=[]:
i = 0
while abs(x-x1_line[i])>x_dist and i < len(x1_line)-1:
i = i + 1
y1_current = y1_line[i]
j = 0
while abs(x-x2_line[j])>x_dist and j < len(x2_line)-1:
j = j + 1
y2_current = y2_line[j]
if abs(y2_current-y1_current)<y_dist and i != len(x1_line) and j != len(x2_line):
ymax = max(y1_current, y2_current)
ymin = min(y1_current, y2_current)
xmax = max(x1_line[i], x2_line[j])
xmin = min(x1_line[i], x2_line[j])
intersection.append((x, ymin+(ymax-ymin)/2))
ax.plot(x, y1_current, 'ro') # Plot the cross point
idx += 1
print "intersection points", intersection
plt.show()
Для массивов f и g мы могли бы просто сделать следующее:
np.pad(np.diff(np.array(f > g).astype(int)), (1,0), 'constant', constant_values = (0,))
Это даст массив всех точек кроссовера. Каждый 1 является кроссовером снизу вверх и каждым -1 кроссовером сверху вниз.
Там может быть несколько пересечений, вы можете найти точку (x,y) на каждом пересечении следующим пониманием списка
intersections = [(x[i], f[i]) for i,_ in enumerate(zip(f,g)) if f[i] == g[i]]
В качестве простого примера
>>> x = [1,2,3,4,5]
>>> f = [2,4,6,8,10]
>>> g = [10,8,6,4,2]
>>> [(x[i], f[i]) for i,_ in enumerate(zip(f,g)) if f[i] == g[i]]
[(3, 6)]
Таким образом, найдена точка пересечения в точке x = 3, y = 6. Обратите внимание: если вы используете float, два значения могут быть не равными, поэтому вы можете использовать некоторый допуск вместо ==.
Даже если f и g пересекаются, вы не можете быть уверены, что f [i] == g [i] для целого я (пересечение, вероятно, происходит между точками).
Вместо этого вы должны проверить, как
# detect intersection by change in sign of difference
d = f - g
for i in range(len(d) - 1):
if d[i] == 0. or d[i] * d[i + 1] < 0.:
# crossover at i
x_ = x[i]
Пересечение, вероятно, происходит между точками. Давайте рассмотрим пример ниже.
import numpy as np
import matplotlib.pyplot as plt
xs=np.arange(0, 20)
y1=np.arange(0, 20)*2
y2=np.array([1, 1.5, 3, 8, 9, 20, 23, 21, 13, 23, 18, 20, 23, 24, 31, 28, 30, 33, 37, 36])
построив 2 кривые выше вместе с их пересечениями, используя в качестве пересечения средние координаты до и после предложенного из пересечения idx, все точки ближе к первой кривой.
idx=np.argwhere(np.diff(np.sign(y1 - y2 )) != 0).reshape(-1) + 0
plt.plot(xs, y1)
plt.plot(xs, y2)
for i in range(len(idx)):
plt.plot((xs[idx[i]]+xs[idx[i]+1])/2.,(y1[idx[i]]+y1[idx[i]+1])/2., 'ro')
plt.legend(['Y1', 'Y2'])
plt.show()
используя в качестве пересечения средние координаты до и после, но для обеих кривых y1 и y2 обычно ближе к истинному пересечению
plt.plot(xs, y1)
plt.plot(xs, y2)
for i in range(len(idx)):
plt.plot((xs[idx[i]]+xs[idx[i]+1])/2.,(y1[idx[i]]+y1[idx[i]+1]+y2[idx[i]]+y2[idx[i]+1])/4., 'ro')
plt.legend(['Y1', 'Y2'])
plt.show()
Для еще более точной оценки пересечения мы могли бы использовать интерполяцию.
Здесь есть решение, которое:
distance_upper_bound в запросе KD-дерева, чтобы определить, насколько близко достаточно близко. 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import cKDTree
from scipy import interpolate
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], projection='3d')
ax.axis('off')
def upsample_coords(coord_list):
# s is smoothness, set to zero
# k is degree of the spline. setting to 1 for linear spline
tck, u = interpolate.splprep(coord_list, k=1, s=0.0)
upsampled_coords = interpolate.splev(np.linspace(0, 1, 100), tck)
return upsampled_coords
# target line
x_targ = [1, 2, 3, 4, 5, 6, 7, 8]
y_targ = [20, 100, 50, 120, 55, 240, 50, 25]
z_targ = [20, 100, 50, 120, 55, 240, 50, 25]
targ_upsampled = upsample_coords([x_targ, y_targ, z_targ])
targ_coords = np.column_stack(targ_upsampled)
# KD-tree for nearest neighbor search
targ_kdtree = cKDTree(targ_coords)
# line two
x2 = [3,4,5,6,7,8,9]
y2 = [25,35,14,67,88,44,120]
z2 = [25,35,14,67,88,44,120]
l2_upsampled = upsample_coords([x2, y2, z2])
l2_coords = np.column_stack(l2_upsampled)
# plot both lines
ax.plot(x_targ, y_targ, z_targ, color='black', linewidth=0.5)
ax.plot(x2, y2, z2, color='darkgreen', linewidth=0.5)
# find intersections
for i in range(len(l2_coords)):
if i == 0: # skip first, there is no previous point
continue
distance, close_index = targ_kdtree.query(l2_coords[i], distance_upper_bound=.5)
# strangely, points infinitely far away are somehow within the upper bound
if np.isinf(distance):
continue
# plot ground truth that was activated
_x, _y, _z = targ_kdtree.data[close_index]
ax.scatter(_x, _y, _z, 'gx')
_x2, _y2, _z2 = l2_coords[i]
ax.scatter(_x2, _y2, _z2, 'rx') # Plot the cross point
plt.show()