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

Question

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

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

到目前为止，我已经采用了这种方法，将每个循环的局部最大值以下几个点的循环切片，并从硬编号的点数中制作红线。红线的近似是由poly1d(polyfit(x1,x2,1))然后fsolve用于获得交点。然而，它并不总是正确地工作，因为点的分布并不总是相同的。

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

所需的输出是具有交叉点近似坐标的列表 - 如果你想玩，这里是曲线的数据 (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

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

我将为自己制作一些简单的数据，一条长摆线，并找到 y 坐标从增加到减少的地方，类似于这里：

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