B-сплайн-интерполяция с Python

Ответ 1

Мне удалось воссоздать пример Mathematica, о котором я спросил в предыдущем сообщении, используя Python/scipy. Здесь результат:

B-сплайн, апериодический

Хитрость заключалась в том, чтобы либо перехватить коэффициенты, т.е. элемент 1 кортежа, возвращаемый scipy.interpolate.splrep, и заменить их значениями контрольной точки, прежде чем передать их в scipy.interpolate.splev, или, если вы в порядке с созданием узлы, вы также можете обойтись без splrep и создать весь кортеж самостоятельно.

Однако это странно, так как согласно руководству splrep возвращает (и splev ожидает) кортеж, содержащий, среди прочего, вектор сплайн-коэффициентов с одним коэффициентом за узел. Однако, по всем найденным мною источникам, сплайн определяется как взвешенная сумма сплайнов базиса N_control_points, поэтому я ожидаю, что вектор коэффициентов будет иметь столько элементов, сколько контрольных точек, а не узлов.

Фактически, при поставке splrep результата кортеж с вектором коэффициентов, измененным, как описано выше, на scipy.interpolate.splev, оказывается, что первые N_control_points этого вектора фактически являются ожидаемыми коэффициентами для сплайнов базиса N_control_points. Последняя степень + 1 элемент этого вектора, по-видимому, не имеет никакого эффекта. Я в тупике, почему это так. Если кто-нибудь сможет это прояснить, это будет здорово. Здесь источник, который генерирует приведенные выше графики:

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-сплайн, периодический

Теперь, чтобы создать замкнутую кривую, подобную следующей, что является другим примером Mathematica, который можно найти в Интернете,

необходимо установить параметр per в вызове splrep, если вы его используете. После заполнения списка контрольных точек со степенью + 1 значения в конце, это, кажется, работает достаточно хорошо, как показывают изображения.

Следующая особенность здесь, однако, состоит в том, что элементы первой и последней степени в векторе коэффициентов не влияют, что означает, что контрольные точки должны быть помещены в вектор, начиная со второй позиции, то есть положению 1. Только тогда это результаты ok. Для степеней k = 4 и k = 5 эта позиция даже изменяется в положение 2.

Здесь источник генерации замкнутой кривой:

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-сплайн, периодический, высшая степень

Наконец, есть эффект, который я тоже не могу объяснить, и это при переходе к степени 5, есть небольшой разрыв, который появляется в сплайсированной кривой, см. верхнюю правую панель, которая является крупным планом это "полумесяц с носовой формой". Исходный код, который производит это, приведен ниже.

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

Учитывая, что b-сплайны являются повсеместными в научном сообществе, и этот scipy является таким всеобъемлющим набором инструментов и что я не смог много узнать о том, что я прошу здесь в Интернете, заставляет меня поверить, что я на неправильном пути или что-то не замечаешь. Любая помощь будет оценена.

Ответ 2

Используйте эту функцию, которую я написал для еще один вопрос, который я задал здесь.

В моем вопросе я искал способы вычисления bsplines с scipy (вот как я на самом деле наткнулся на ваш вопрос).

После большой одержимости, я придумал функцию ниже. Он будет оценивать любую кривую до 20-й степени (путь больше, чем нам нужно). И скорость мудрая я протестировал ее на 100 000 образцов, и это заняло 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.array([0]*degree + range(count-degree+1) + [count-degree]*degree,dtype='int')

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


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

    return points

Результаты как для открытых, так и для периодических кривых:

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

Ответ 3

Я считаю, что scipy fitpack Library делает что-то более сложное, чем то, что делает Mathematica. Я был в замешательстве относительно того, что происходило.

В этих функциях есть параметр сглаживания, а поведение интерполяции по умолчанию заключается в попытке сделать точки пройденными линиями. Это то, что делает это программное обеспечение fitpack, поэтому я думаю, что scipy просто унаследовал его? (http://www.netlib.org/fitpack/all - я не уверен, что это подходящий пакет)

Я взял некоторые идеи http://research.microsoft.com/en-us/um/people/ablake/contours/ и закодировал ваш пример с B-сплайнами там.

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

В любом случае, я думаю, что у других людей такие же проблемы, посмотрите: Поведение screy's sprep