用numpy找到xy数据点图的所有交点？

Question

我正在分析循环拉伸试验的数据。 巨大的 x 和 y 值列表作为输入。 为了描述材料是变硬还是变软，我需要得到每个循环循环的蓝色斜率。

下坡是儿童节，上坡就是挑战。

到目前为止，我已经采用了这种方法，将每个循环的局部最大值以下几个点的循环切掉，并使红线从 hardnumbered 点计数。通过poly1d(polyfit(x1,x2,1)) 近似红线，然后使用fsolve 得到交点。但是它并不总是正确地工作，因为点的分布并不总是相同的。

问题是如何正确定义两条（红色）相交线的间隔。上图中是 3 个实验以及平均斜率。我花了几天时间试图为每个循环找到 4 个最近的点，决定这不是最好的方法。最后，我在 stackowerflow 结束了。

所需的输出是带有交叉点近似坐标的列表 - 如果你想玩，here 是曲线的数据 (0,[[xvals],[yvals]])。使用

可以轻松阅读这些内容

import csv
import sys
csv. field_size_limit(sys.maxsize)     

csvfile = 'data.csv'
tc_data = {}
for key, val in csv.reader(open(csvfile, "r")):
    tc_data[key] = val
for key in tc_data:
  tc = eval(tc_data[key])

x = tc[0]
y = tc[1]

Answer 1

这可能有点矫枉过正，但是找到交点的正确方法是，一旦你将曲线分成块，看看第一个块的任何部分是否与第二个块的任何部分相交。

我要给自己做一些简单的数据，一个prolate cycloid的一部分，我要找到y坐标从增加到减少的地方，类似于here：

a, b = 1, 2
phi = np.linspace(3, 10, 100)
x = a*phi - b*np.sin(phi)
y = a - b*np.cos(phi)
y_growth_flips = np.where(np.diff(np.diff(y) > 0))[0] + 1

plt.plot(x, y, 'rx')
plt.plot(x[y_growth_flips], y[y_growth_flips], 'bo')
plt.axis([2, 12, -1.5, 3.5])
plt.show()

如果您有两条线段，一条从点P0 到P1，另一条从点Q0 到Q1，您可以通过求解向量方程P0 + s*(P1-P0) = Q0 + t*(Q1-Q0) 找到它们的交点，并且如果s 和t 都在[0, 1] 中，则这两个段实际上会相交。尝试所有细分市场：

x_down = x[y_growth_flips[0]:y_growth_flips[1]+1]
y_down = y[y_growth_flips[0]:y_growth_flips[1]+1]
x_up = x[y_growth_flips[1]:y_growth_flips[2]+1]
y_up = y[y_growth_flips[1]:y_growth_flips[2]+1]

def find_intersect(x_down, y_down, x_up, y_up):
    for j in xrange(len(x_down)-1):
        p0 = np.array([x_down[j], y_down[j]])
        p1 = np.array([x_down[j+1], y_down[j+1]])
        for k in xrange(len(x_up)-1):
            q0 = np.array([x_up[k], y_up[k]])
            q1 = np.array([x_up[k+1], y_up[k+1]])
            params = np.linalg.solve(np.column_stack((p1-p0, q0-q1)),
                                     q0-p0)
            if np.all((params >= 0) & (params <= 1)):
                return p0 + params[0]*(p1 - p0)

>>> find_intersect(x_down, y_down, x_up, y_up)
array([ 6.28302264,  1.63658676])

crossing_point = find_intersect(x_down, y_down, x_up, y_up)
plt.plot(crossing_point[0], crossing_point[1], 'ro')
plt.show()

在我的系统上，它每秒可以处理大约 20 个交叉点，这不是超快的，但可能足以不时地分析图形。您可以通过向量化 2x2 线性系统的解决方案来加快速度：

def find_intersect_vec(x_down, y_down, x_up, y_up):
    p = np.column_stack((x_down, y_down))
    q = np.column_stack((x_up, y_up))
    p0, p1, q0, q1 = p[:-1], p[1:], q[:-1], q[1:]
    rhs = q0 - p0[:, np.newaxis, :]
    mat = np.empty((len(p0), len(q0), 2, 2))
    mat[..., 0] = (p1 - p0)[:, np.newaxis]
    mat[..., 1] = q0 - q1
    mat_inv = -mat.copy()
    mat_inv[..., 0, 0] = mat[..., 1, 1]
    mat_inv[..., 1, 1] = mat[..., 0, 0]
    det = mat[..., 0, 0] * mat[..., 1, 1] - mat[..., 0, 1] * mat[..., 1, 0]
    mat_inv /= det[..., np.newaxis, np.newaxis]
    import numpy.core.umath_tests as ut
    params = ut.matrix_multiply(mat_inv, rhs[..., np.newaxis])
    intersection = np.all((params >= 0) & (params <= 1), axis=(-1, -2))
    p0_s = params[intersection, 0, :] * mat[intersection, :, 0]
    return p0_s + p0[np.where(intersection)[0]]

是的，它很乱，但它可以工作，而且速度快 100 倍：

find_intersect(x_down, y_down, x_up, y_up)
Out[67]: array([ 6.28302264,  1.63658676])

find_intersect_vec(x_down, y_down, x_up, y_up)
Out[68]: array([[ 6.28302264,  1.63658676]])

%timeit find_intersect(x_down, y_down, x_up, y_up)
10 loops, best of 3: 66.1 ms per loop

%timeit find_intersect_vec(x_down, y_down, x_up, y_up)
1000 loops, best of 3: 375 us per loop

Answer 2

您可以非常简单地做到这一点，只需使用 scipy 中的 interp1d 函数以相同的 x 值重新采样所有行的 y 值。

http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html