2015-07-16 11 views
19

Mam zbiór punktów pts które tworzą pętlę i wygląda to tak:Montaż zamkniętą krzywą do zbioru punktów

enter image description here

ten jest nieco podobny do 31243002, ale zamiast umieszczać punkty pomiędzy parami punktów, chciałbym aby dopasować gładką krzywą przez punkty (współrzędne podane są na końcu tego pytania), więc próbowałem coś podobnego do scipy dokumentacji na Interpolation:

values = pts 
tck = interpolate.splrep(values[:,0], values[:,1], s=1) 
xnew = np.arange(2,7,0.01) 
ynew = interpolate.splev(xnew, tck, der=0) 

ale otrzymuję ten błąd:

ValueError: Error on input data

Czy istnieje jakiś sposób, aby znaleźć taki atak?

Współrzędne punktów:

pts = array([[ 6.55525 , 3.05472 ], 
    [ 6.17284 , 2.802609], 
    [ 5.53946 , 2.649209], 
    [ 4.93053 , 2.444444], 
    [ 4.32544 , 2.318749], 
    [ 3.90982 , 2.2875 ], 
    [ 3.51294 , 2.221875], 
    [ 3.09107 , 2.29375 ], 
    [ 2.64013 , 2.4375 ], 
    [ 2.275444, 2.653124], 
    [ 2.137945, 3.26562 ], 
    [ 2.15982 , 3.84375 ], 
    [ 2.20982 , 4.31562 ], 
    [ 2.334704, 4.87873 ], 
    [ 2.314264, 5.5047 ], 
    [ 2.311709, 5.9135 ], 
    [ 2.29638 , 6.42961 ], 
    [ 2.619374, 6.75021 ], 
    [ 3.32448 , 6.66353 ], 
    [ 3.31582 , 5.68866 ], 
    [ 3.35159 , 5.17255 ], 
    [ 3.48482 , 4.73125 ], 
    [ 3.70669 , 4.51875 ], 
    [ 4.23639 , 4.58968 ], 
    [ 4.39592 , 4.94615 ], 
    [ 4.33527 , 5.33862 ], 
    [ 3.95968 , 5.61967 ], 
    [ 3.56366 , 5.73976 ], 
    [ 3.78818 , 6.55292 ], 
    [ 4.27712 , 6.8283 ], 
    [ 4.89532 , 6.78615 ], 
    [ 5.35334 , 6.72433 ], 
    [ 5.71583 , 6.54449 ], 
    [ 6.13452 , 6.46019 ], 
    [ 6.54478 , 6.26068 ], 
    [ 6.7873 , 5.74615 ], 
    [ 6.64086 , 5.25269 ], 
    [ 6.45649 , 4.86206 ], 
    [ 6.41586 , 4.46519 ], 
    [ 5.44711 , 4.26519 ], 
    [ 5.04087 , 4.10581 ], 
    [ 4.70013 , 3.67405 ], 
    [ 4.83482 , 3.4375 ], 
    [ 5.34086 , 3.43394 ], 
    [ 5.76392 , 3.55156 ], 
    [ 6.37056 , 3.8778 ], 
    [ 6.53116 , 3.47228 ]]) 
+1

Czy chcesz zainstalować nowy pakiet/framework? Jeśli jesteś w rodzaju dopasowania, o którym mówisz, jest dostępne przez [ROOT-Framework] (https://root.cern.ch), a także wiele innych opcji dopasowania. Powinno być całkiem łatwo dostosować przykład [Histogram 2D] (https://root.cern.ch/root/htmldoc/tutorials/fit/fit2dHist.C.html), aby dopasować dane do PyROOT (interfejs Pythona do ROOT który używa składni python zamiast interpretera C++). Jeśli jest to coś, czemu nie jesteś przeciwny, mogę wysłać prawidłową odpowiedź i przykład. – Matt

+1

@Matt: Dziękuję za komentarz. Nie mam nic przeciwko instalowaniu nowego pakietu, chociaż martwię się, że dane wyjściowe mogą być użyte w 'matplotlib' (mam kilka zdjęć i chciałbym zachować ten sam styl w tym). – Mahdi

+1

To najwyraźniej dotyczyło kogoś innego, ponieważ był tam [post na temat używania matplotlib w/ROOT] (http://www.rootpy.org/auto_examples/plotting/plot_matplotlib_hist.html). ROOT to bardzo potężne narzędzie i polecam wypróbowanie go, istnieje wiele wspaniałych funkcji do analizy danych i wizualizacji. – Matt

Odpowiedz

18

W rzeczywistości nie byłeś zbytnio od rozwiązania w swoim pytaniu.

Przy użyciu parametru scipy.interpolate.splprep w przypadku parametrycznej interpolacji B-splajnu można zastosować najprostsze podejście. To również natywnie obsługuje zamknięte krzywe, jeśli podasz parametr per=1,

import numpy as np 
from scipy.interpolate import splprep, splev 
import matplotlib.pyplot as plt 

# define pts from the question 

tck, u = splprep(pts.T, u=None, s=0.0, per=1) 
u_new = np.linspace(u.min(), u.max(), 1000) 
x_new, y_new = splev(u_new, tck, der=0) 

plt.plot(pts[:,0], pts[:,1], 'ro') 
plt.plot(x_new, y_new, 'b--') 
plt.show() 

enter image description here

Zasadniczo podejście to nie bardzo różni się od tej w odpowiedzi @Joe kington użytkownika. Chociaż prawdopodobnie będzie to nieco bardziej stabilne, ponieważ odpowiednik wektora i jest wybierany domyślnie na podstawie odległości między punktami, a nie tylko ich indeksu (patrz splprep documentation dla parametru u).

+0

nice! thx do wysyłania. –

+1

Jest zwięzły i działa świetnie, czego jeszcze chcę? – Mahdi

+0

do czego odnosi się ten magiczny "pts.T"? –

5

Aby dopasować gładką krzywą zamkniętą przez N punktów można wykorzystać segmenty linii z następującymi ograniczeniami:

Segment
  • Każda linia ma dotknąć jej dwa końcowe punktów (2 warunki na odcinek linii)
  • Dla każdego punktu, lewy i prawy odcinek linii muszą mieć tę samą pochodną (2 warunki na punkt == 2 warunki na odcinek linii)

Aby mieć wystarczającą swobodę w sumie 4 warunki na odcinek linii, równanie każdego odcinka linii powinno być y = ax^3 + bx^2 + cx + d. (więc pochodną jest y = 3ax^2 + 2bx + c)

Ustawienie warunków zgodnie z sugestiami dałoby N * 4 równania liniowego dla N * 4 niewiadomych (a1..aN, b1..bN, c1 ..cN, d1..dN) rozwiązywalne poprzez inwersję macierzy (numpy).

Jeśli punkty znajdują się na tej samej linii pionowej, wymagana jest specjalna (ale prosta) obsługa, ponieważ pochodna będzie "nieskończona".

16

Twój problem polega na tym, że próbujesz bezpośrednio pracować z X i Y. Wywoływana funkcja interpolacji zakłada, że ​​wartości x są posortowane według kolejności i że każda wartość x będzie miała unikalną wartość y.

Zamiast tego należy wykonać sparametryzowany układ współrzędnych (na przykład indeks wierzchołków) i interpolować xiy oddzielnie.

Aby rozpocząć, należy rozważyć następujące kwestie:

import numpy as np 
from scipy.interpolate import interp1d # Different interface to the same function 
import matplotlib.pyplot as plt 

#pts = np.array([...]) # Your points 

x, y = pts.T 
i = np.arange(len(pts)) 

# 5x the original number of points 
interp_i = np.linspace(0, i.max(), 5 * i.max()) 

xi = interp1d(i, x, kind='cubic')(interp_i) 
yi = interp1d(i, y, kind='cubic')(interp_i) 

fig, ax = plt.subplots() 
ax.plot(xi, yi) 
ax.plot(x, y, 'ko') 
plt.show() 

enter image description here

nie zamknąć wielokąt. Jeśli chcesz, możesz dodać pierwszy punkt do końca tablicy (na przykład: pts = np.vstack([pts, pts[0]])

Jeśli to zrobisz, zauważysz, że istnieje nieciągłość, w której wielokąt się zamyka.

enter image description here

To dlatego nasza parametryzacja nie bierze pod uwagę zamknięcie polgyon. Szybkim rozwiązaniem jest pad tablica z „odbitych” punktów:

import numpy as np 
from scipy.interpolate import interp1d 
import matplotlib.pyplot as plt 

#pts = np.array([...]) # Your points 

pad = 3 
pts = np.pad(pts, [(pad,pad), (0,0)], mode='wrap') 
x, y = pts.T 
i = np.arange(0, len(pts)) 

interp_i = np.linspace(pad, i.max() - pad + 1, 5 * (i.size - 2*pad)) 

xi = interp1d(i, x, kind='cubic')(interp_i) 
yi = interp1d(i, y, kind='cubic')(interp_i) 

fig, ax = plt.subplots() 
ax.plot(xi, yi) 
ax.plot(x, y, 'ko') 
plt.show() 

enter image description here

Alternatywnie, można użyć specjalistycznego algorytmu wygładzania krzywej takich jak PEAK lub algorytm narożną cięcia.

+0

Dziękuję za wspaniałą odpowiedź i szczegółowe wyjaśnienie. – Mahdi

8

Korzystanie z interfejsu ROOT Framework i pyroot udało mi się wygenerować następujący obraz Drawing Using pyroot

z następującego kodu (I konwertowane danych do CSV nazywa Data.csv więc czytając go do ROOT byłoby łatwiejsze i nadał kolumnom tytuły xp, yp)

from ROOT import TTree, TGraph, TCanvas, TH2F 

c1 = TCanvas('c1', 'Drawing Example', 200, 10, 700, 500) 
t=TTree('TP','Data Points') 
t.ReadFile('./data.csv') 
t.SetMarkerStyle(8) 
t.Draw("yp:xp","","ACP") 
c1.Print('pydraw.png')