Rozwinąć wszystkie punkty przecięcia wykresu punktów danych Xy za pomocą numpy?

Question

Analizuję dane z cyklicznych prób rozciągania. Jako dane wejściowe są używane ogromnych list wartości x i y. Aby opisać, czy materiał twardnieje lub zmiękcza, muszę uzyskać niebieskie nachylenie każdej pętli cyklu.

Getting dolny punkt nachylenia jest festiwal dziecko, ale górna, która jest wyzwaniem.

Mam tego podejścia do tej pory, krojenie na pętle z kilku punktów poniżej lokalnego maksimum każdej pętli i podejmowania czerwone linie z hardnumbered liczbę punktów. Aproximation czerwonych linii jest przez poly1d(polyfit(x1,x2,1)) i fsolve jest używany następnie, aby uzyskać punkt przecięcia. Jednak nie działa zawsze poprawnie, ponieważ rozkład punktów nie zawsze jest taki sam.

Problemem jest to, jak prawidłowo określić interwał dwóch (czerwony) przecinających się linii. Na powyższym obrazku znajdują się 3 eksperymenty z uśrednionym nachyleniem. Spędziłem kilka dni próbując znaleźć 4 najbliższe punkty dla każdej pętli, decydując, że nie jest to najlepsze podejście. I w końcu skończyłem tutaj na stosie upustu.

Żądanymi danymi wyjściowymi są lista z przybliżonymi współrzędnymi punktów przecięcia - jeśli chcesz odtworzyć, here są danymi dla krzywej (0, [[xvals], [yvals]]). Theese można łatwo odczytać z

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

To może być trochę przesada, ale właściwa droga, aby znaleźć punkt przecięcia, gdy już podzielone krzywą na kawałki, aby sprawdzić, czy każdy odcinek od pierwszego kawałka przecina się z każdym segmencie z drugiej porcji.

mam zamiar zrobić sobie kilka prostych danych, kawałek o prolate cycloid i zamierzam znaleźć miejsca gdzie współrzędna y koziołki ze zwiększenia się do zmniejszenia podobnie do 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() 

Jeśli masz dwa segmenty, jeden przechodzący od punktu P0 do P1, a drugi biegnący od punktu Q0 do Q1, możesz znaleźć punkt przecięcia, rozwiązując równanie wektorowe P0 + s*(P1-P0) = Q0 + t*(Q1-Q0), a dwa segmenty faktycznie przecinają się, jeśli oba s i t są w [0, 1].Starając się to dla wszystkich segmentów:

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() 

W moim systemie, to może obsłużyć około 20 skrzyżowań na sekundę, co nie jest super, ale na tyle pewnie, aby analizować wykresy co jakiś czas. Możesz być w stanie spped rzeczy przez Wektoryzacja rozwiązania systemów 2x2 liniowej:

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

Tak, to bałagan, ale to działa, i robi tak x100 razy szybciej:

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

Można to zrobić w bardzo prosty sposób za pomocą funkcji interp1d z scipy do resample Y wartości wszystkich liniach w tych samych wartości x.

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