python - optimize - Ajustar una curva cerrada a un conjunto de puntos

optimize curve fit numpy

Question

Tengo un conjunto de puntos de puntos que forman un bucle y se ve así:

Esto es algo similar a 31243002 , pero en lugar de poner puntos entre pares de puntos, me gustaría ajustar una curva suave a través de los puntos (las coordenadas se dan al final de la pregunta), así que probé algo similar a la documentación scipy en 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)

pero me sale este error:

ValueError: error en datos de entrada

¿Hay alguna manera de encontrar tal ajuste?

Coordenadas de los puntos:

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

Answer 1

En realidad, no estabas lejos de la solución en tu pregunta.

Usar scipy.interpolate.splprep para la interpolación paramétrica B-spline sería el enfoque más simple. También admite de forma nativa curvas cerradas, si proporciona el parámetro 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()

Básicamente, este enfoque no es muy diferente al de la respuesta de @Joe Kington. Aunque, probablemente será un poco más robusto, porque el equivalente del vector i se elige, de forma predeterminada, en función de las distancias entre puntos y no simplemente su índice (consulte la scipy.interpolate.splprep para el parámetro u ).

Answer 2

Para ajustar una curva cerrada suave a través de N puntos, puede usar segmentos de línea con las siguientes restricciones:

• Cada segmento de línea tiene que tocar sus dos puntos finales (2 condiciones por segmento de línea)
• Para cada punto, el segmento de línea izquierdo y derecho debe tener la misma derivada (2 condiciones por punto == 2 condiciones por segmento de línea)

Para poder tener suficiente libertad para un total de 4 condiciones por segmento de línea, la ecuación de cada segmento de línea debe ser y = ax ^ 3 + bx ^ 2 + cx + d. (entonces la derivada es y ''= 3ax ^ 2 + 2bx + c)

Establecer las condiciones como se sugiere le dará N * 4 ecuaciones lineales para N * 4 incógnitas (a1..aN, b1..bN, c1..cN, d1..dN) solucionables por inversión de matriz (numpy).

Si los puntos están en la misma línea vertical, se requiere un manejo especial (pero simple) ya que la derivada será "infinita".

Answer 3

Su problema es porque está tratando de trabajar con x e y directamente. La función de interpolación que está llamando supone que los valores x están en orden y que cada valor x tendrá un valor y único.

En cambio, necesitará hacer un sistema de coordenadas parametrizado (por ejemplo, el índice de sus vértices) e interpolar x e y separadamente usándolo.

Para empezar, considere lo siguiente:

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

No cerré el polígono. Si lo desea, puede agregar el primer punto al final de la matriz (por ejemplo, pts = np.vstack([pts, pts[0]])

Si haces eso, notarás que hay una discontinuidad donde el polígono se cierra.

Esto se debe a que nuestra parametrización no tiene en cuenta el cierre del polígono. Una solución rápida es rellenar la matriz con los puntos "reflejados":

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

Alternativamente, puede usar un algoritmo especializado de suavizado de curvas como PEAK o un algoritmo de corte de esquinas.

Answer 4

Utilizando el marco ROOT y la interfaz pyroot pude generar la siguiente imagen

Con el siguiente código (convertí sus datos a un CSV llamado data.csv, por lo que leerlos en ROOT sería más fácil y le di a las columnas títulos 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'')