Bezier and B-Spline Techniques (Mathematics and Visualization)
目次
はじめに
様々なデータを補間したり、スムージングしたりするのに、
三次スプライン曲線はよく使われますが、
より複雑な曲線を表現するのに、BSplineという曲線があり、
Pythonの科学技術計算用ライブラリであるSciPyでは、
このBSpline用のクラスや便利関数を提供しています。
今回の記事では、SciPyが提供するBSpline関連のツールの使い方を
紹介したいと思います。
BSplineとは?
BSplineは区間多項式曲線(piecewise polynomial curve)の一つで、
下記の式で表されます。
ここで、
S(x)はxにおけるスプライン上の補間点、
ciは制御点と呼ばれるスプラインを生成する代表点、
p+1はBSplineの次元
kはノットの数
B_ip(x)は基底関数(Basis Function)と呼ばれる関数
です。
(BSplineのBはこのBasis Functionから来ているそうです)
上記の基底関数は、
下記のDe Boorの再帰式で表現することができます。
ここで、
- ti はノットベクトルの各要素
です。
上記の各iに対するBSplineの基底関数をグラフ化すると、
下記のようなグラフになります。
このグラフはこちらのコードで作成できます。
scipyのBSplineは基底関数単体を計算する関数を提供していないため、
指定したi以外の制御点を0とすることで、
擬似的にある基底関数単体を計算できる関数Bを定義しています。
関数B_origは前述の定義式をベースとした基底関数の計算関数ですが、
一般的に再帰を使うため、計算が遅いとされています。
from scipy.interpolate import BSpline def B_orig(x, k, i, t): if k == 0: return 1.0 if t[i] <= x < t[i + 1] else 0.0 if t[i + k] == t[i]: c1 = 0.0 else: c1 = (x - t[i]) / (t[i + k] - t[i]) * B(x, k - 1, i, t) if t[i + k + 1] == t[i + 1]: c2 = 0.0 else: c2 = (t[i + k + 1] - x) / (t[i + k + 1] - t[i + 1]) * B(x, k - 1, i + 1, t) return c1 + c2 def B(x, k, i, t): c = np.zeros_like(t) c[i] = 1 return BSpline(t, c, k)(x) def main(): k = 3 # degree of the spline t = [0, 1, 2, 3, 4, 5] # knots vector x = np.linspace(0, 5, 1000, endpoint=False) t = np.r_[[np.min(t)]*k, t, [np.max(t)]*k] n = len(t) - k - 1 for i in range(n): y = np.array([B(ix, k, i, t) for ix in x]) plt.plot(x, y, label=f'i = {i}') plt.title(f'Basis functions (k = {k}, knots = {t})') plt.show()
ここで重要なのは、すべての基底関数を足すと、
どのxの値でも1になることです。
これにより、各制御点をこの基底関数で重み付け足し算した時に、
滑らかな曲線を作ることができます。
SciPyにおける1次元のBスプライン曲線を生成するルーチン
SciPyで1次元のBスプライン曲線を作るルーチンとして、
下記のスプライン補間とスプライン近似(スムージング)の2種類が提供されています。
また、SciPyのBSpline関連のコードでは、歴史的な理由から、
1: 1980年代に開発されたFORTRANのライブラリであるFITPACKをPythonでラップしたAPI
https://netlib.org/dierckx/readme
と
2: 近年新しく内部ロジックをPythonで実装した新しいAPI
があります。
1は長年使われていて、安定していますが、
拡張性やメンテナンスの問題、
int32を使っているため大規模なデータでの処理に問題があることから、
最近は、2のAPIの方がアクティブに開発されていますので、
本記事では2のAPIをメインに説明します。
両方の種類のAPIともに、補間と近似の両方の機能を提供していますが、
1はAPIがオブジェクト指向であり、それぞれの補間クラスが補間を実施しています。
2はAPIがジェネレータ設計になっており、生成された曲線は、
すべてBSplineという単一のクラスのオブジェクトとして生成されます。
1. スプライン補間
与えられたx,yの点列を必ず通るスプライン補間としては、
scipy.interpolate.make_interp_splineという関数が提供されています。
下記は、あるx,yデータを
1, 2, 3次元のBSplineで補間するコードです。
import numpy as np from scipy.interpolate import make_interp_spline, BSpline rng = np.random.default_rng() x = np.linspace(-3, 3, 10) y = np.exp(-x**2) + 0.1 * rng.standard_normal(10) xs = np.linspace(-3, 3, 100) plt.plot(x, y, 'ro') for k in [1, 2, 3]: spl_i = make_interp_spline(x, y, k=k) plt.plot(xs, spl_i(xs), label=f'Interpolation spline (k={k}') plt.legend() plt.show()
上記のコードを実行すると、
下記のように各指定した点を通る
各次元のBSpline曲線が引けることがわかります。
また、生成されたBSplineオブジェクトには、
derivativeというメンバー関数が存在しており、
これを使うことで、任意の点の任意の次数の微分を計算できます。
下記は、上記のスプライン補間の各xにおける微分値を計算するコードです。
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import make_interp_spline, BSpline rng = np.random.default_rng(1234) x = np.linspace(-3, 3, 10) y = np.exp(-x**2) + 0.1 * rng.standard_normal(10) xs = np.linspace(-3, 3, 100) for k in [1, 2, 3]: spl_i = make_interp_spline(x, y, k=k) dspl_i = spl_i.derivative(1) plt.plot(xs, dspl_i(xs), label=f'Interpolation spline derivative (k={k})') plt.legend() plt.show()
1次のスプライン関数は折れ線グラフなので、
微分値が同じ区間では定数であることがわかります。
また、make_interp_splineは、端点の条件をbc_typeという
オプション引数で指定することができます。
デフォルトは"not-a-knot"(None)で端点の制約無し、
"clampled"は端点の微分が０、
"natural"は端点の微分の変化量が０
になる設定です。
("periodicという設定もありますが、
これは始点と終点が同じで閉じるような
スプラインを作りたい時に、始点と終点を滑らかにつなぐ設定です。
下記は、各bc_typeの時のスプラインを比較したものになります。
端点の制約が変更されているので、始点と終点周りを見ると
"clamped"の場合は、補間が水平になる(微分値が0)になり、
"natural"の場合は、微分値の変化量が0なので、直線の補間になっていることがわかります。
FITPACKのAPIによるスプライン補間
FITPACKのAPIによるスプライン補間としては、
scipy.interpolate.InterpolatedUnivariateSplineが提供されています。
こちらのAPIもmake_interp_splineと同じように利用することが可能です。
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import InterpolatedUnivariateSpline rng = np.random.default_rng() x = np.linspace(-3, 3, 10) y = np.exp(-x**2) + 0.1 * rng.standard_normal(10) xs = np.linspace(-3, 3, 100) plt.plot(x, y, 'ro') for k in [1, 2, 3]: spl_i = InterpolatedUnivariateSpline(x, y, k=k) plt.plot(xs, spl_i(xs), label=f'Interpolation spline (k={k})') plt.legend() plt.show()
2. スプライン近似(スムージング)
続いて、与えられてたx-yの点を必ず通るスプラインを作るのではなく、
そのx-yの点からの偏差が小さくなるようなスプライン曲線を生成するのには
scipy.interpolate.make_lsq_splineを利用できます。
このようなスプライン曲線の生成をスムージングやスプライン近似(Approximation)と呼びます。
下記がスプライン近似のサンプルコードです。
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import make_lsq_spline, BSpline rng = np.random.default_rng() x = np.linspace(-3, 3, 50) y = np.exp(-x**2) + 0.1 * rng.standard_normal(50) t = [-1, 0, 1] k = 3 t = np.r_[(x[0],)*(k+1), t, (x[-1],)*(k+1)] spl = make_lsq_spline(x, y, t, k=k) xs = np.linspace(-3, 3, 100) plt.subplots() plt.plot(x, y, 'ro', ms=5, label='data') plt.plot(xs, spl(xs), label='LSQ spline') plt.legend(loc='best') plt.show()
上記のコードを実行すると、下記のような結果が得られます。
入力の重みのWを調整することで、
一部の点にだけはできるだけ近づけるなどが可能になります。
また、scipy 3.10.0からは、scipy.interpolate.make_smoothing_spline
という新しいルーチンが追加されました。
こちらは、先程のscipy.interpolate.make_lsq_splineと異なり、
必要な入力は、x,yのデータのみで、
ノットベクトルtを評価関数に基づいて
自動設定してくれるため、使いやすいと思います。
先程の例で、scipy.interpolate.make_smoothing_splineを利用すると、
下記のようになります。
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import make_smoothing_spline rng = np.random.default_rng() x = np.linspace(-3, 3, 50) y = np.exp(-x**2) + 0.1 * rng.standard_normal(50) spl = make_smoothing_spline(x, y) xs = np.linspace(-3, 3, 100) plt.subplots() plt.plot(x, y, 'ro', ms=5, label='data') plt.plot(xs, spl(xs), label='Smoothing spline') plt.legend(loc='best') plt.show()
このルーチンはオプショナル引数である
lamの値を大きな値に設定すると、
下記のように、よりなめらかな曲線になり、
入力データに追従しないようになります。
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import make_smoothing_spline rng = np.random.default_rng() x = np.linspace(-3, 3, 50) y = np.exp(-x**2) + 0.1 * rng.standard_normal(50) spl = make_smoothing_spline(x, y) xs = np.linspace(-3, 3, 100) plt.subplots() plt.plot(x, y, 'ro', ms=5, label='data') for l in np.arange(-5.0, 5.0): spl = make_smoothing_spline(x, y, lam=10**l) plt.plot(xs, spl(xs), label=f'lam={10**l:.2f}') plt.legend(loc='best') plt.show()
FITPACKのAPIによるスプライン近似
FITPACKのAPIによるスプライン近似としては、
scipy.interpolate.LSQUnivariateSplineが提供されています。
こちらのAPIもmake_lsq_splineと同じように利用することが可能です。
import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import LSQUnivariateSpline rng = np.random.default_rng() x = np.linspace(-3, 3, 50) y = np.exp(-x**2) + 0.1 * rng.standard_normal(50) t = [-1, 0, 1] k = 3 t = np.r_[(x[0],)*(k+1), t, (x[-1],)*(k+1)] spl = LSQUnivariateSpline(x, y, t, k=k) xs = np.linspace(-3, 3, 100) plt.subplots() plt.plot(x, y, 'ro', ms=5, label='data') plt.plot(xs, spl(xs), label='LSQ spline') plt.legend(loc='best') plt.show()
参考資料
Bezier and B-Spline Techniques (Mathematics and Visualization)
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。
