Ajustement d'une courbe fermée pour un ensemble de points

j'ai un ensemble de points pts qui forme une boucle et il ressemble à ceci:

enter image description here

c'est un peu similaire à 31243002 , mais au lieu de mettre des points entre les paires de points, je voudrais ajuster une courbe lisse à travers les points (les coordonnées sont données à la fin de la question), donc j'ai essayé quelque chose de similaire à scipy documentation sur 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)

mais je reçois cette erreur:

ValueError: Erreur sur les données d'entrée

Est-il possible de trouver une telle forme?

coordonnées des points:

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 ]])
23
demandé sur Community 2015-07-16 23:54:42

4 réponses

en Fait, vous n'étiez pas loin de la solution dans votre question.

utilisant scipy.interpolate.splprep pour l'interpolation B-spline paramétrique serait l'approche la plus simple. Il supporte aussi nativement les courbes fermées, si vous fournissez le paramètre 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

fondamentalement, cette approche n'est pas très différente de celle de la réponse de @Joe Kington. Bien que, il sera probablement un peu plus robuste, parce que l'équivalent du vecteur i est choisi, par défaut, basé sur les distances entre les points et pas simplement leur index (voir splprep documentation pour le paramètre u ).

21
répondu rth 2015-07-17 10:07:47

votre problème est que vous essayez de travailler avec x et y directement. La fonction d'interpolation que vous appelez suppose que les valeurs x sont en ordre et que chaque valeur x aura une valeur y unique.

au lieu de cela, vous aurez besoin de faire un système de coordonnées paramétré (par exemple l'index de vos sommets) et interpoler x et y séparément en l'utilisant.

pour commencer, considérez ce qui suit:

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

Je n'ai pas fermé le polygone. Si vous le souhaitez, vous pouvez ajouter le premier point à la fin du tableau (par exemple, pts = np.vstack([pts, pts[0]])

si vous faites ça, vous remarquerez qu'il y a une discontinuité là où le polygone se ferme.

enter image description here

c'est parce que notre paramétrage ne tient pas compte de la fermeture du polgyon. Rapide Fixer est de garnir le tableau avec les" reflected "points:

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

alternativement, vous pouvez utiliser un algorithme de lissage de courbe spécialisé tel que PEAK ou un algorithme de découpe d'angle.

18
répondu Joe Kington 2015-07-16 21:37:05

en utilisant le root Framework et l'interface pyroot j'ai pu générer l'image suivante Drawing Using pyroot

avec le code suivant(j'ai converti vos données en un CSV appelé data.csv donc la lecture dans la racine serait plus facile et a donné les titres de colonnes de 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')
8
répondu Matt 2015-07-16 22:24:47

pour ajuster une courbe fermée en douceur à travers N points, vous pouvez utiliser les segments de ligne avec les contraintes suivantes:

  • chaque segment de ligne doit toucher ses deux extrémités (2 conditions par segment de ligne)
  • pour chaque point, les segments de droite et de gauche doivent avoir la même dérivée (2 conditions par point == 2 conditions par segment de ligne)

pour être en mesure d'avoir assez de liberté pour au total 4 conditions par segment de ligne l'équation de chaque segment de ligne devrait être y = ax^3 + bx^2 + cx + D. (donc le dérivé est y' = 3ax^2 + 2bx + c)

la mise en place des conditions suggérées vous donnerait N * 4 équations linéaires pour N * 4 inconnues (a1 .. aN, b1..bN, c1 .. cN, d1..dN) résoluble par inversion de matrice (numpy).

si les points sont sur la même ligne verticale, une manipulation spéciale (mais simple) est nécessaire car la dérivée sera"infinie".

5
répondu Jacques de Hooge 2015-07-16 21:10:23