Сплайны с Python (с использованием узлов управления и конечных точек)

Я пытаюсь сделать что-то вроде следующего (изображение извлечено из википедии)

#!/usr/bin/env python
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

# sampling
x = np.linspace(0, 10, 10)
y = np.sin(x)

# spline trough all the sampled points
tck = interpolate.splrep(x, y)
x2 = np.linspace(0, 10, 200)
y2 = interpolate.splev(x2, tck)

# spline with all the middle points as knots (not working yet)
# knots = x[1:-1]  # it should be something like this
knots = np.array([x[1]])  # not working with above line and just seeing what this line does
weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y, t=knots, w=weights)
x3 = np.linspace(0, 10, 200)
y3 = interpolate.splev(x2, tck)

# plot
plt.plot(x, y, 'go', x2, y2, 'b', x3, y3,'r')
plt.show()

Первая часть кода - это код, извлеченный из основной ссылки, но он не объяснил, как использовать точки в качестве узлов управления.

Результатом этого кода является следующее изображение.

Точками являются образцы, синяя линия - сплайн с учетом всех точек. И красная линия - это та, которая не работает для меня. Я стараюсь учитывать все промежуточные точки в качестве узлов управления, но я просто не могу. Если я пытаюсь использовать knots=x[1:-1], это просто не работает. Я был бы признателен за любую помощь.

Короче говоря: как использовать все промежуточные точки в качестве узлов управления в функции сплайна?

Примечание: это последнее изображение именно то, что мне нужно, и это разница между тем, что у меня есть (сплайн, передающий все точки) и тем, что мне нужно (сплайн с контрольными узлами). Есть идеи?

Ответ 1

Я только что нашел что-то действительно интересное с ответом, который мне нужен с bézier в этой ссылке. Затем я использовал код, чтобы попробовать самостоятельно. Очевидно, он работает нормально. Это моя реализация:

#! /usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import binom

def Bernstein(n, k):
    """Bernstein polynomial.

    """
    coeff = binom(n, k)

    def _bpoly(x):
        return coeff * x ** k * (1 - x) ** (n - k)

    return _bpoly


def Bezier(points, num=200):
    """Build Bézier curve from points.

    """
    N = len(points)
    t = np.linspace(0, 1, num=num)
    curve = np.zeros((num, 2))
    for ii in range(N):
        curve += np.outer(Bernstein(N - 1, ii)(t), points[ii])
    return curve
xp = np.array([2,3,4,5])
yp = np.array([2,1,4,0])
x, y = Bezier(list(zip(xp, yp))).T

plt.plot(x,y)
plt.plot(xp,yp,"ro")
plt.plot(xp,yp,"b--")

plt.show()

И образ для примера.

Красные точки представляют контрольные точки. Что это =)

Ответ 2

В этом IPython Notebook http://nbviewer.ipython.org/github/empet/geom_modeling/blob/master/FP-Bezier-Bspline.ipynb вы можете найти подробное описание данных, связанных с созданием кривой B-сплайна, а также Реализация Python алгоритма де Бура.

Ответ 3

Я думаю, что проблема заключается в том, чтобы сделать с вами узел узлов. Кажется, что это вызывает проблемы, если вы выбираете слишком много узлов, у него должна быть точка данных между узлами. Этот вопрос решает проблему Ошибка (?) При выборе узлов в функции splrep scipy.insterpolate

#!/usr/bin/env python
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

# sampling
x = np.linspace(0, 10, 10)
y = np.sin(x)

# spline trough all the sampled points
tck = interpolate.splrep(x, y)
print tck
x2 = np.linspace(0, 10, 200)
y2 = interpolate.splev(x2, tck)

# spline with all the middle points as knots (not working yet)
knots = np.asarray(x[1:-1])  # it should be something like this
#knots = np.array([x[1]])  # not working with above line and just seeing what this line does
nknots = 5
idx_knots = (np.arange(1,len(x)-1,(len(x)-2)/np.double(nknots))).astype('int')
knots = x[idx_knots]
print knots

weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y,  t=knots, w=weights)
x3 = np.linspace(0, 10, 200)
y3 = interpolate.splev(x2, tck)

# plot
plt.plot(x, y, 'go', x2, y2, 'b', x3, y3,'r')
plt.show()

кажется, что он выбирает 5 узлов, 6 дает более неожиданные результаты, а также дает ошибки.

Ответ 4

Если вы хотите оценить bspline, вам нужно выяснить соответствующий вектор узла для вашего сплайна, а затем вручную перестроить tck в соответствии с вашими потребностями.

tck обозначает узлы t + коэффициенты c + степень кривой k. splrep вычисляет tck для кубической кривой, проходящей через заданные контрольные точки. Поэтому вы не можете использовать его для чего хотите.

Функция ниже покажет вам мое решение для аналогичного вопроса, который я задал некоторое время назад., адаптированный для того, что вы хотите.

Интересный факт: код работает для кривых любого измерения (1D, 2D, 3D,..., nD)

import numpy as np
import scipy.interpolate as si


def bspline(cv, n=100, degree=3):
    """ Calculate n samples on a bspline

        cv :      Array ov control vertices
        n  :      Number of samples to return
        degree:   Curve degree
    """
    cv = np.asarray(cv)
    count = cv.shape[0]

    # Prevent degree from exceeding count-1, otherwise splev will crash
    degree = np.clip(degree,1,count-1)

    # Calculate knot vector
    kv = np.array([0]*degree + range(count-degree+1) + [count-degree]*degree,dtype='int')

    # Calculate query range
    u = np.linspace(0,(count-degree),n)

    # Calculate result
    points = np.zeros((len(u),cv.shape[1]))
    for i in xrange(cv.shape[1]):
        points[:,i] = si.splev(u, (kv,cv[:,i],degree))

    return points

Проверьте это:

import matplotlib.pyplot as plt
colors = ('b', 'g', 'r', 'c', 'm', 'y', 'k')

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

plt.plot(cv[:,0],cv[:,1], 'o-', label='Control Points')

for d in range(1,5):
    p = bspline(cv,n=100,degree=d,periodic=True)
    x,y = p.T
    plt.plot(x,y,'k-',label='Degree %s'%d,color=colors[d%len(colors)])

plt.minorticks_on()
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(35, 70)
plt.ylim(0, 30)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

Результат:

Ответ 5

функция вашего примера периодична, и вам нужно добавить параметр per=True к методу interpolate.splrep.

knots = x[1:-1]
weights = np.concatenate(([1],np.ones(x.shape[0]-2)*.01,[1]))
tck = interpolate.splrep(x, y, t=knots, w=weights, per=True)

Это дает мне следующее:

Изменить: Это также объясняет, почему он работал с knots = x[-2:2], который является непериодическим подмножеством полного диапазона.