Julia Advent Calendar 2020/12/21

TL;DR

JuliaでBスプラインを扱うためのパッケージ BasicBSpline.jlを作りました。
B-splineは、コンピュータグラフィックス・CAD・数値解析など、様々な応用があります。
本記事の後半ではパッケージの使い方を書きます。

↓↓このパッケージのために作ったロゴです↓↓

B-spline基底関数の定義

B-spline基底関数とは, 与えられたノット列 $k=(k_1,\dots,k_l)$ に対して以下の漸化式(Cox-de Boorの公式)で定義される関数のことです。

\begin{aligned} {B}_{(i,p,k)}(t) &= \frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t) +\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\ {B}_{(i,0,k)}(t) &= \begin{cases} &1\quad (k_{i}\le t< k_{i+1})\\ &0\quad (\text{otherwise}) \end{cases} \end{aligned}

ノット列というのは、実軸上にならんだ(広義)単調増加有限列のことです。

例えば, $p=0$ に対してはこのような矩形の不連続関数になります。

$p=1$ に対してはこのような折れ線グラフになります。

$p=2$ だとこんな感じ。

Desmosで試すにはこちら↓

今回作ったパッケージBasicBSpline.jlはこれらの関数を扱うためのパッケージになります。
たったこれだけ、 「一つの漸化式(Cox-de Boorの公式)を扱うだけ」 なのですが、結構奥が深くて8ヶ月くらいかけてチマチマ実装しています。
そして先日ついに公式パッケージの方にv0.1.0が登録されました🎉

一例を見てみましょう。
JuliaでCox-de Boorの公式を素直にそのまま実装してみるとこのようになります。

function B(i,p,k,t)
    if p==0
        return float(k[i]≤t<k[i+1])
    else
        return B(i,p-1,k,t)*(t-k[i])/(k[i+p]-k[i]) + B(i+1,p-1,k,t)*(k[i+p+1]-t)/(k[i+p+1]-k[i+1])
    end
end

BasicBSpline.jlでは、この関数がすでにbsplinebasisとして用意されています。
両者の計算時間を比較してみましょう。

using BenchmarkTools

knotvector = sort(rand(10))
p = 3
i = 1
t = 0.2
@benchmark B(i,p,knotvector,t)

k = Knots(knotvector)
P = FastBSplineSpace(p,k)
@benchmark bsplinebasis(P,i,t)

素直な実装: 118.757 ns
BasicBSpline.jl: 46.722 ns

ダンチでBasicBSpline.jlの勝ち！
速度面もそうですが、機能面でも色々と充実してます。

そもそもB-splineって何に使うの？

まずはBézier曲線から

Bézier曲線をご存知の方はたぶん多いと思います。
IllustratorやInkscapeなど等のドロー系ソフトで使われているアレですね。
このようなイラスト等では、複数のBézier曲線を組み合わせて形状が決定されています。

Bézier曲線には次のような性質があります。

曲線の形状は制御点から決定される
制御点を平行移動すると、それに追従して曲線も平行移動する

じつは、このBézier曲線は $p=3$ 次の多項式曲線で、以下のように表すことができます。

\begin{aligned} \bm{p}(t) &= \sum_{i=1}^n B_{(i,p)}(t)\bm{a}_i & (n=p+1) \\ B_{(i,p)}(t) &= \binom{p}{i}(1-t)^{p-i-1}t^{i} \end{aligned}

ここで $B_{(i,p)}$ はBernstein多項式と呼ばれるもので、 $n$ は制御点の数です。

さて、Bézier曲線を複数繋げたとして、その滑らかさは一般には保証されません。
(滑らかさとは、数学的には「どれくらい微分することができるか」という情報のことです)

接ベクトルが滑らかに繋がらない例↓

曲線の平面曲率が不連続関数になる例↓

では、ある程度の滑らかさを保証してBézier曲線を繋げることはできないのでしょうか?
この問題に対する一つの答えがB-spline曲線です。
(ここでは $p$ 次多項式を繋げて区分多項式を作ることを考えています。全体での $C^p$ 級を仮定すれば区分化する意味が無いので、 $C^{p-1}$ 級以下の滑らかさを要求することになります。)

B-spline曲線

上記のBézier曲線でのBernstein基底関数 $B_{(i,p)}$ をB-spline基底関数 $B_{(i,p,k)}$ に取り替えてみましょう！

\begin{aligned} \bm{p}(t) &= \sum_{i=1}^n B_{(i,p,k)}(t)\bm{a}_i & (k=(k_1,\dots, k_l), n = l-p-1) \end{aligned}

はい、これがB-spline曲線です。
等間隔にノット列を取って、 $p=2$ とした場合は以下のようなB-spline曲線が得られます。

B-spline基底関数は区分多項式なので、B-spline曲線は結局、複数のBézier曲線をつなげたものということになります。

さらに、その繋ぎ目(ノット)である程度滑らかになっています。これがかなり嬉しい性質です。

その一例として、iPhoneの外周を紹介しましょう。
iPhoneの外周は線分と円弧で

のようには作られていません。「円弧と線分」だとつなぎ目で曲率が不連続になるため、光の反射が不自然になってしまいます。

実際は以下のようになっています↓

画像はAccessory Design Guidelines for Apple Devices (Release R13)より

このような複雑な形状は、B-splineを用いて設計されており、その形状をもとに図面に落とされています。(たぶん)

Bézier曲面

Bézier曲線に関連して、Bézier曲面というものもあります。

定義式は

\begin{aligned} \bm{p}(t_1,t_2) &= \sum_{i=1}^{n_1}\sum_{j=1}^{n_2} B_{(i,p_1)}(t_1)B_{(j,p_2)}(t_2)\bm{a}_{ij} & (n_1=p_1+1, n_2=p_2+1) \end{aligned}

です。

B-spline曲面

同様に基底関数をBernstein基底関数からB-spline基底関数に取り替えましょう！

\begin{aligned} \bm{p}(t_1,t_2) &= \sum_{i=1}^{n_1}\sum_{j=1}^{n_2} B_{(i,p_1,k_1)}(t_1)B_{(j,p_2,k_2)}(t_2)\bm{a}_{ij} & (n_1=l_1-p_1-1, n_2=l_2-p_2-1) \end{aligned}

曲線のときと同様に、B-spline曲線はBézier曲面を繋げたものなのですが、こちらの画像からも「滑らか」に繋がっていることが分かると思います。

また、埋め込み先の空間が2次元の場合でも"曲面"と呼ぶこととします。
(後のコード例ではこの曲面を例にします。)

結局、B-splineとは

数学的には区分多項式
Bézier曲線を「滑らかに繋げて」一般化したもの
曲線(curve)・曲面(surface)だけではなく、立体(中身が詰まったもの, solid)も表現可能
なのでコンピュータ上のモデリングなどによく使われる
さらにそのような形状を解析するための数値計算にも使われます
(発展)より一般化したNURBS、T-splineとかもあります

BasicBSpline.jlのできること

ようやく本題です。
以下はあくまで一例です。

🄋インストール

公式パッケージに入ったので！以下のコマンドで入ります。

pkg> add BasicBSpline

①B-spline基底関数の計算

Plots.jlでプロットまで試してみましょう

using BasicBSpline
using Plots
gr()

k = Knots([0.00,1.50,2.50,5.50,8.00,9.00,9.50,10.0])
P0 = BSplineSpace(0,k) # 0th degree piecewise polynomial space
P1 = BSplineSpace(1,k) # 1st degree piecewise polynomial space
P2 = BSplineSpace(2,k) # 2nd degree piecewise polynomial space
P3 = BSplineSpace(3,k) # 3rd degree piecewise polynomial space
plot(
    plot([t->bsplinebasis(P0,i,t) for i in 1:dim(P0)], 0, 10, ylims=(0,1), legend=false),
    plot([t->bsplinebasis(P1,i,t) for i in 1:dim(P1)], 0, 10, ylims=(0,1), legend=false),
    plot([t->bsplinebasis(P2,i,t) for i in 1:dim(P2)], 0, 10, ylims=(0,1), legend=false),
    plot([t->bsplinebasis(P3,i,t) for i in 1:dim(P3)], 0, 10, ylims=(0,1), legend=false),
    layout=(2,2),
)

②B-spline曲線の描画

描画にはBasicBSpline.jlとは別のパッケージExportNURBS.jlを使います。
こちらはまだ公式パッケージに登録してないので、以下のコマンドでインストールして下さい。

pkg> add https://github.com/hyrodium/ExportNURBS.jl

BasicBSpline.jlで出力するには以下のコードを実行します。

## 1-dim B-spline manifold
p = 2 # degree of polynomial
k = Knots(1:12) # knot vector
P = FastBSplineSpace(p, k) # B-spline space
a = [[i-5, 3*sin(i^2)] for i in 1:dim(P)] # control points
M = BSplineCurve([P], a) # Define B-spline manifold
save_png("docs/src/img/1dim.png", M, unitlength = 50)

③B-spline曲線のフィッティング

ここではsin関数を例にフィッテングさせてみます。
BasicBSpline.jlでは以下のようにしてフィッテングして制御点が決定できます。

p = 3
k = Knots(range(-2π,2π,length=8))+p*Knots(-2π,2π)
P = FastBSplineSpace(p, k)

f(t) = [t,sin(t)]

a = fittingcontrolpoints(f, P)
M = BSplineManifold([P],a)
save_svg("sine_curve.svg", M, unitlength=50, up=2, down=-2, left=-8, right=8)

ここでのフィッティングの計算には最小二乗法を使って以下の値を最小化するような制御点の位置 $\bm{a}_i$ を決定しています。

\begin{aligned} &\int_{I}\left\|\bm{f}(t) - \sum_i {B_{(i,p,k)}(t) \bm{a}_i}\right\|^2 dt & \bigl(I = [k_{1+p},k_{l-p}], k &= (k_1,\dots,k_l)\bigr) \end{aligned}

ノット列の位置を変更して、sin関数を近似する区間を変更すると以下のようになります。

この動画を生成するコードはは以下のgistにあります。

④B-spline曲面の描画

以下のコードで実行できます。

using BasicBSpline
using ExportNURBS

p = 2 # degree of polynomial
k = Knots(1:8) # knot vector
P = BSplineSpace(p,k) # B-spline space
rand_a = [rand(2) for i in 1:dim(P), j in 1:dim(P)]
a = [[2*i-6.5,2*j-6.5] for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
M = BSplineManifold([P,P],a) # Define B-spline manifold
save_png("2dim.png", M) # save image

上記のコードでやっていることは

多項式次数の定義 (p)
ノット列の定義 (k)
区分多項式空間の定義 (P)
B-spline曲面の定義 (BSplineManifold)
png出力 (save_png)

です。
出力結果(png)はこちら↓

⑤B-spline曲面のフィッテング

以下のコードで実行できます。

p1 = 2
p2 = 2
k1 = Knots(-10:10)+p1*Knots(-10,10)
k2 = Knots(-10:10)+p2*Knots(-10,10)
P1 = FastBSplineSpace(p1, k1)
P2 = FastBSplineSpace(p2, k2)

f(u1, u2) = [2u1 + sin(u1) + cos(u2) + u2 / 2, 3u2 + sin(u2) + sin(u1) / 2 + u1^2 / 6] / 5

a = fittingcontrolpoints(f, P1, P2)
M = BSplineManifold([P1,P2],a)
save_png("fitting.png", M, unitlength=50, up=10, down=-10, left=-10, right=10)

上記のコードでやっていることは

多項式次数の定義 (p1, p2)
ノット列の定義 (k1, k2)
区分多項式空間の定義 (P1, P2)
フィッテング対象の関数の定義 (f)
フィッテングの計算 (fittingcontolpoints)
B-spline曲面の定義 (BSplineManifold)
png出力 (save_png)

です。
出力結果(png)はこちら↓

Desmosでの確認はこちら↓

⑥B-spline多様体の細分

「形状を保ったまま、次数やノットの数を増やす操作」のことを細分(refinement)と呼びます。
理論はともかくとして、BasicBSpline.jlで計算してみましょう。

using BasicBSpline
using ExportNURBS

p = 2 # degree of polynomial
k = Knots(1:8) # knot vector
P = BSplineSpace(p,k) # B-spline space
rand_a = [rand(2) for i in 1:dim(P), j in 1:dim(P)]
a = [[2*i-6.5,2*j-6.5] for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
M = BSplineManifold([P,P],a) # Define B-spline manifold
save_png("2dim.png", M) # save image

これは上の「④B-spline曲面の描画」での例と同じです。
これの細分を計算しましょう。

k₊=[Knots(3.3,4.2),Knots(3.8,3.2,5.3)] # additional knots
M′ = refinement(M,k₊=k₊) # refinement of B-spline manifold
save_png("2dim_refinement.png", M′) # save image

これら2つの画像での「赤色の形状」は全く同じですが、制御点の数は細分後に増えています。
ここの例のような、ノット列の挿入による細分はh-細分と呼ばれます。(ノット列の間隔を表す記号に $h$ が使われるためのようです。)
一方で、次数を増やす細分の操作はp-細分と呼ばれます。(こちらの例も入れたかったのですが、バグが入ってて微妙に間に合わず..あとで修正します)

このように制御点を増やすことはいくつかの利点があります:

制御点の細かい調整でより精度の高いモデリングができる
B-spline多様体上の数値計算の精度が向上できる

今後の方針

B-spline関連の実装は一通り終わった気がしますが..以下のような変更予定があります
- Knots型をパラメトリック型に変更してFloat64以外も対応できるようにする
- 等間隔のノット列のための型を用意して、一部の計算を高速化させる
NURBS、T-splineへの一般化
- T-splineはまだ私が完全には理解できてないです。今後の課題
- そもそもT-splineは特許関連があるのでややこしそう。
- 2003年に特許申請されたらしいのでそろそろ切れそう…？
- 実装するにしても別パッケージになる予定
形状出力用のパッケージ(ExportNURBS.jl)の充実化
- NURBSに対応してから充実させるつもりで命名していました
- しかし、BasicBSpline.jl用にパッケージを分けた方が良さそう
- 現在はPNG, SVG出力のみに対応していますが、将来的にはPOV-Ray, Blender, STL, Paraview向けの出力などにも対応したい
数値計算の具体例を書く
- 数値計算AdventCalendar2018に書いた記事を書き直したい…。

もっと学びたい人へ

BasicBSpline.jlのドキュメント
- https://hyrodium.github.io/BasicBSpline.jl/dev/
- Mathematical properties of B-splineの節とかどうかな
Geometric Modeling with Splines
- https://www.routledge.com/p/book/9780367447243
- 英語ですが数学的に厳密に書かれてて分かりやすいです。
- 日本語の良い書籍って無いんですか？ → 知らないです
NURBSによる形状表現
- https://hyrodium.github.io/pdf/
- 私の書いた資料です。(全40ページ)
- マジで読んでください！
- マジで読んでください！！

Discussion