Interpolación B-spline con Python

Estoy tratando de reproducir un ejemplo de Mathematica para un B-spline con Python.

El código del ejemplo matemático dice:

pts = {{0, 0}, {0, 2}, {2, 3}, {4, 0}, {6, 3}, {8, 2}, {8, 0}}; Graphics[{BSplineCurve[pts, SplineKnots -> {0, 0, 0, 0, 2, 3, 4, 6, 6, 6, 6}], Green, Line[pts], Red, Point[pts]}] 

y produce lo que espero. Ahora trato de hacer lo mismo con Python / scipy:

 import numpy as np import matplotlib.pyplot as plt import scipy.interpolate as si points = np.array([[0, 0], [0, 2], [2, 3], [4, 0], [6, 3], [8, 2], [8, 0]]) x = points[:,0] y = points[:,1] t = range(len(x)) knots = [2, 3, 4] ipl_t = np.linspace(0.0, len(points) - 1, 100) x_tup = si.splrep(t, x, k=3, t=knots) y_tup = si.splrep(t, y, k=3, t=knots) x_i = si.splev(ipl_t, x_tup) y_i = si.splev(ipl_t, y_tup) print 'knots:', x_tup fig = plt.figure() ax = fig.add_subplot(111) plt.plot(x, y, label='original') plt.plot(x_i, y_i, label='spline') plt.xlim([min(x) - 1.0, max(x) + 1.0]) plt.ylim([min(y) - 1.0, max(y) + 1.0]) plt.legend() plt.show() 

Esto da como resultado algo que también está interpolado pero no se ve bien. Puedo parametrizar y dividir los componentes x e y por separado, utilizando los mismos nudos que matemática. Sin embargo, me sobran y me faltan los cordones, lo que hace que mi curva interpolada se arquee fuera del casco convexo de los puntos de control. ¿Cuál es la forma correcta de hacer esto / cómo las matemáticas lo hacen?

    Pude recrear el ejemplo de Mathematica que pregunté en el post anterior usando Python / scipy. Aquí está el resultado:

    B-Spline, Aperiodic

    Spline a través de una curva 2D.

    El truco consistía en interceptar los coeficientes, es decir, el elemento 1 de la tupla devuelta por scipy.interpolate.splrep , y reemplazarlos con los valores del punto de control antes de entregarlos a scipy.interpolate.splev , o, si está bien con la creación Los nudos usted mismo, también puede prescindir de splrep y crear la tupla completa usted mismo.

    Sin embargo, lo extraño de todo esto es que, según el manual, splrep devuelve (y splev espera) una tupla que contiene, entre otros, un vector de coeficientes de spline con un coeficiente por nudo. Sin embargo, según todas las fonts que encontré, una spline se define como la sum ponderada de las splines basadas en N_control_points , por lo que esperaría que el vector de coeficientes tenga tantos elementos como puntos de control, no posiciones de nudo.

    De hecho, cuando se suministra la tupla de resultados de splrep con el vector de coeficientes modificado como se describe anteriormente a scipy.interpolate.splev , resulta que los primeros N_control_points de ese vector en realidad son los coeficientes esperados para las splines de base de N_control_points . El último grado + 1 elementos de ese vector parecen no tener efecto. Estoy sorprendido de por qué se hace de esta manera. Si alguien puede aclarar eso, sería genial. Aquí está la fuente que genera los gráficos anteriores:

     import numpy as np import matplotlib.pyplot as plt import scipy.interpolate as si points = [[0, 0], [0, 2], [2, 3], [4, 0], [6, 3], [8, 2], [8, 0]]; points = np.array(points) x = points[:,0] y = points[:,1] t = range(len(points)) ipl_t = np.linspace(0.0, len(points) - 1, 100) x_tup = si.splrep(t, x, k=3) y_tup = si.splrep(t, y, k=3) x_list = list(x_tup) xl = x.tolist() x_list[1] = xl + [0.0, 0.0, 0.0, 0.0] y_list = list(y_tup) yl = y.tolist() y_list[1] = yl + [0.0, 0.0, 0.0, 0.0] x_i = si.splev(ipl_t, x_list) y_i = si.splev(ipl_t, y_list) #============================================================================== # Plot #============================================================================== fig = plt.figure() ax = fig.add_subplot(231) plt.plot(t, x, '-og') plt.plot(ipl_t, x_i, 'r') plt.xlim([0.0, max(t)]) plt.title('Splined x(t)') ax = fig.add_subplot(232) plt.plot(t, y, '-og') plt.plot(ipl_t, y_i, 'r') plt.xlim([0.0, max(t)]) plt.title('Splined y(t)') ax = fig.add_subplot(233) plt.plot(x, y, '-og') plt.plot(x_i, y_i, 'r') plt.xlim([min(x) - 0.3, max(x) + 0.3]) plt.ylim([min(y) - 0.3, max(y) + 0.3]) plt.title('Splined f(x(t), y(t))') ax = fig.add_subplot(234) for i in range(7): vec = np.zeros(11) vec[i] = 1.0 x_list = list(x_tup) x_list[1] = vec.tolist() x_i = si.splev(ipl_t, x_list) plt.plot(ipl_t, x_i) plt.xlim([0.0, max(t)]) plt.title('Basis splines') plt.show() 

    B-Spline, Periódico

    Ahora para crear una curva cerrada como la siguiente, que es otro ejemplo de Mathematica que se puede encontrar en la web, Curva b-spline cerrada

    es necesario establecer el parámetro per en la llamada de splrep , si utiliza eso. Después de rellenar la lista de puntos de control con valores de grado + 1 al final, esto parece funcionar lo suficientemente bien, como muestran las imágenes.

    La siguiente peculiaridad aquí, sin embargo, es que los elementos de primer y último grado en el vector de coeficientes no tienen ningún efecto, lo que significa que los puntos de control deben colocarse en el vector a partir de la segunda posición, es decir, la posición 1. Sólo entonces los resultados De acuerdo. Para grados k = 4 y k = 5, esa posición incluso cambia a la posición 2.

    Aquí está la fuente para generar la curva cerrada:

     import numpy as np import matplotlib.pyplot as plt import scipy.interpolate as si points = [[-2, 2], [0, 1], [-2, 0], [0, -1], [-2, -2], [-4, -4], [2, -4], [4, 0], [2, 4], [-4, 4]] degree = 3 points = points + points[0:degree + 1] points = np.array(points) n_points = len(points) x = points[:,0] y = points[:,1] t = range(len(x)) ipl_t = np.linspace(1.0, len(points) - degree, 1000) x_tup = si.splrep(t, x, k=degree, per=1) y_tup = si.splrep(t, y, k=degree, per=1) x_list = list(x_tup) xl = x.tolist() x_list[1] = [0.0] + xl + [0.0, 0.0, 0.0, 0.0] y_list = list(y_tup) yl = y.tolist() y_list[1] = [0.0] + yl + [0.0, 0.0, 0.0, 0.0] x_i = si.splev(ipl_t, x_list) y_i = si.splev(ipl_t, y_list) #============================================================================== # Plot #============================================================================== fig = plt.figure() ax = fig.add_subplot(231) plt.plot(t, x, '-og') plt.plot(ipl_t, x_i, 'r') plt.xlim([0.0, max(t)]) plt.title('Splined x(t)') ax = fig.add_subplot(232) plt.plot(t, y, '-og') plt.plot(ipl_t, y_i, 'r') plt.xlim([0.0, max(t)]) plt.title('Splined y(t)') ax = fig.add_subplot(233) plt.plot(x, y, '-og') plt.plot(x_i, y_i, 'r') plt.xlim([min(x) - 0.3, max(x) + 0.3]) plt.ylim([min(y) - 0.3, max(y) + 0.3]) plt.title('Splined f(x(t), y(t))') ax = fig.add_subplot(234) for i in range(n_points - degree - 1): vec = np.zeros(11) vec[i] = 1.0 x_list = list(x_tup) x_list[1] = vec.tolist() x_i = si.splev(ipl_t, x_list) plt.plot(ipl_t, x_i) plt.xlim([0.0, 9.0]) plt.title('Periodic basis splines') plt.show() 

    B-Spline, Periódico, Grado Superior

    Por último, hay un efecto que tampoco puedo explicar, y esto es cuando va al grado 5, hay una pequeña discontinuidad que aparece en la curva dividida, vea el panel superior derecho, que es un primer plano de esa mitad. -mone-with-nose-shape ‘. El código fuente que produce esto se enumera a continuación.

    Discontinuidad.

     import numpy as np import matplotlib.pyplot as plt import scipy.interpolate as si points = [[-2, 2], [0, 1], [-2, 0], [0, -1], [-2, -2], [-4, -4], [2, -4], [4, 0], [2, 4], [-4, 4]] degree = 5 points = points + points[0:degree + 1] points = np.array(points) n_points = len(points) x = points[:,0] y = points[:,1] t = range(len(x)) ipl_t = np.linspace(1.0, len(points) - degree, 1000) knots = np.linspace(-degree, len(points), len(points) + degree + 1).tolist() xl = x.tolist() coeffs_x = [0.0, 0.0] + xl + [0.0, 0.0, 0.0] yl = y.tolist() coeffs_y = [0.0, 0.0] + yl + [0.0, 0.0, 0.0] x_i = si.splev(ipl_t, (knots, coeffs_x, degree)) y_i = si.splev(ipl_t, (knots, coeffs_y, degree)) #============================================================================== # Plot #============================================================================== fig = plt.figure() ax = fig.add_subplot(231) plt.plot(t, x, '-og') plt.plot(ipl_t, x_i, 'r') plt.xlim([0.0, max(t)]) plt.title('Splined x(t)') ax = fig.add_subplot(232) plt.plot(t, y, '-og') plt.plot(ipl_t, y_i, 'r') plt.xlim([0.0, max(t)]) plt.title('Splined y(t)') ax = fig.add_subplot(233) plt.plot(x, y, '-og') plt.plot(x_i, y_i, 'r') plt.xlim([min(x) - 0.3, max(x) + 0.3]) plt.ylim([min(y) - 0.3, max(y) + 0.3]) plt.title('Splined f(x(t), y(t))') ax = fig.add_subplot(234) for i in range(n_points - degree - 1): vec = np.zeros(11) vec[i] = 1.0 x_i = si.splev(ipl_t, (knots, vec, degree)) plt.plot(ipl_t, x_i) plt.xlim([0.0, 9.0]) plt.title('Periodic basis splines') plt.show() 

    Dado que las b-splines son omnipresentes en la comunidad científica, y que Scipy es una caja de herramientas tan completa, y que no he podido encontrar mucho sobre lo que estoy pidiendo aquí en la web, me hace creer que estoy en El camino equivocado o pasar por alto algo. Cualquier ayuda sería apreciada.

    Usa esta función que escribí para otra pregunta que hice aquí.

    En mi pregunta, estaba buscando formas de calcular bsplines con scipy (así es como me topé con tu pregunta).

    Después de mucha obsesión, se me ocurrió la siguiente función. Evaluará cualquier curva hasta el grado 20 (mucho más de lo que necesitamos). Y en cuanto a velocidad, lo probé para 100,000 muestras y tomé 0.017s.

     import numpy as np import scipy.interpolate as si def bspline(cv, n=100, degree=3, periodic=False): """ Calculate n samples on a bspline cv : Array ov control vertices n : Number of samples to return degree: Curve degree periodic: True - Curve is closed False - Curve is open """ # If periodic, extend the point array by count+degree+1 cv = np.asarray(cv) count = len(cv) if periodic: factor, fraction = divmod(count+degree+1, count) cv = np.concatenate((cv,) * factor + (cv[:fraction],)) count = len(cv) degree = np.clip(degree,1,degree) # If opened, prevent degree from exceeding count-1 else: degree = np.clip(degree,1,count-1) # Calculate knot vector kv = None if periodic: kv = np.arange(0-degree,count+degree+degree-1,dtype='int') else: kv = np.concatenate(([0]*degree, np.arange(count-degree+1), [count-degree]*degree)) # Calculate query range u = np.linspace(periodic,(count-degree),n) # Calculate result return np.array(si.splev(u, (kv,cv.T,degree))).T 

    Resultados tanto para curvas abiertas como periódicas:

     cv = np.array([[ 50., 25.], [ 59., 12.], [ 50., 10.], [ 57., 2.], [ 40., 4.], [ 40., 14.]]) 

    Curva periódica (cerrada) Curva abierta

    Creo que la biblioteca fitpack de scipy está haciendo algo más complicado que lo que hace Mathematica. Estaba confundido en cuanto a lo que estaba pasando también.

    Existe el parámetro de suavizado en estas funciones, y el comportamiento de interpolación predeterminado es tratar de hacer que los puntos pasen por las líneas. Eso es lo que hace este software de Fitpack, así que supongo que Scipy simplemente lo heredó. ( http://www.netlib.org/fitpack/all – No estoy seguro de que este sea el Fitpack correcto)

    Tomé algunas ideas de http://research.microsoft.com/en-us/um/people/ablake/contours/ y codifiqué su ejemplo con las B-splines allí.

    Ajuste de spline

    funciones básicas

     import numpy import matplotlib.pyplot as plt # This is the basis function described in eq 3.6 in http://research.microsoft.com/en-us/um/people/ablake/contours/ def func(x, offset): out = numpy.ndarray((len(x))) for i, v in enumerate(x): s = v - offset if s >= 0 and s < 1: out[i] = s * s / 2.0 elif s >= 1 and s < 2: out[i] = 3.0 / 4.0 - (s - 3.0 / 2.0) * (s - 3.0 / 2.0) elif s >= 2 and s < 3: out[i] = (s - 3.0) * (s - 3.0) / 2.0 else: out[i] = 0.0 return out # We have 7 things to fit, so let's do 7 basis functions? y = numpy.array([0, 2, 3, 0, 3, 2, 0]) # We need enough x points for all the basis functions... That's why the weird linspace max here x = numpy.linspace(0, len(y) + 2, 100) B = numpy.ndarray((len(x), len(y))) for k in range(len(y)): B[:, k] = func(x, k) plt.plot(x, B.dot(y)) # The x values in the next statement are the maximums of each basis function. I'm not sure at all this is right plt.plot(numpy.array(range(len(y))) + 1.5, y, '-o') plt.legend('B-spline', 'Control points') plt.show() for k in range(len(y)): plt.plot(x, B[:, k]) plt.title('Basis functions') plt.show() 

    De todos modos, creo que otras personas tienen los mismos problemas, eche un vistazo a: Comportamiento de splrep de scipy