xyデータポイントグラフのすべての交差点をnumpyで区切りますか？

Question

私は周期的引張試験からデータを分析しています。入力として使用される巨大な xとy値のリスト。 材料が硬化するか軟化するかを記述するために、各サイクルループの青い勾配を取得する必要があります。

斜面の下のポイントを取得

は、子祭りですが、上側の一つは、それが課題です。

Iは、各ループの極大以下、いくつかのポイントでループをスライスし、点のhardnumberedカウントから赤い線を作り、これまでのところ、このアプローチを行いました。赤線の近似は、poly1d(polyfit(x1,x2,1))によって行われ、fsolveは交点を得るために使用されます。ただし、ポイントの分布が常に同じではないため、常に正しく機能しません。

問題が正しく線と交差間隔2つ（赤）のを定義する方法あります。上の写真では、3回の実験と平均傾斜を示しています。私は数日間、ループごとに4つの最も近い点を見つけようとしましたが、これは最良のアプローチではないと判断しました。そして最後に、私はここでstackowerflowで終わった。

希望の出力は、交点のおおよその座標を持つリストです。再生する場合は、hereはカーブ（0、[xvals]、[yvals]）のデータです。 Theeseは簡単にあなたが同じxの値で、すべてのラインのy値をリサンプリングするためにscipyのダウンロードからinterp1d機能を使って、非常に簡単にこれを行うことができます

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

これは少しやり過ぎが、あなたの交点を見つけるための適切な方法、あなたはチャンクにあなたの曲線を分割したら、最初のチャンクから任意のセグメントは、任意のセグメントと交差するかどうかであること2番目のチャンクから。

私は、自分自身いくつかの簡単なデータにするためにprolate cycloidの一部をつもりです、とhereと同様に減少に増加からフリップを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のセグメントの2つのセグメントがある場合は、ベクトル式を解くことによって交点を見つけることができます両方のセグメントが実際に交差するのは、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の線形システムの解をベクトル化し物事をsppedすることができる場合があります

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]] 

はい、それは厄介だが、それが動作し、より高速なので、X100回を行います。

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

で読み取ることができます。

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