16

Python で B スプラインの Mathematica の例を再現しようとしています。

Mathematica の例のコードは

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

そして私が期待するものを生み出します。今、私は 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()

これにより、補間も行われますが、見た目が正しくありません。Mathematica と同じノットを使用して、x 成分と y 成分を別々にパラメーター化し、スプライン化します。ただし、オーバーシュートとアンダーシュートが発生し、補間された曲線が制御点の凸包の外側で曲がってしまいます。これを行う正しい方法は何ですか/Mathematica はどのようにそれを行いますか?

4

3 に答える 3

13

私がここで尋ねた別の質問のために私が書いたこの関数を使用してください。

私の質問では、scipy で bsplines を計算する方法を探していました (これが実際にあなたの質問に出くわした方法です)。

多くの執着の後、私は以下の機能を思いつきました。20次までの曲線を評価します（必要以上に）。速度に関しては、100,000 サンプルでテストしたところ、0.017 秒かかりました

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

開曲線と周期曲線の両方の結果:

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

周期的 (閉じた) 曲線 開いた曲線

于 2016-01-15T09:01:59.570 に答える