BasicBSpline.jlを作ったので宣伝です!
Julia Advent Calendar 2020/12/21
TL;DR
- JuliaでBスプラインを扱うためのパッケージ BasicBSpline.jlを作りました。
- B-splineは、コンピュータグラフィックス・CAD・数値解析など、様々な応用があります。
- 本記事の後半ではパッケージの使い方を書きます。
↓↓このパッケージのために作ったロゴです↓↓
B-spline基底関数の定義
B-spline基底関数とは, 与えられたノット列
ノット列というのは、実軸上にならんだ(広義)単調増加有限列のことです。
例えば,
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曲線は
ここで
さて、Bézier曲線を複数繋げたとして、その滑らかさは一般には保証されません。
(滑らかさとは、数学的には「どれくらい微分することができるか」という情報のことです)
接ベクトルが滑らかに繋がらない例↓
曲線の平面曲率が不連続関数になる例↓
では、ある程度の滑らかさを保証してBézier曲線を繋げることはできないのでしょうか?
この問題に対する一つの答えがB-spline曲線です。
(ここでは
B-spline曲線
上記のBézier曲線でのBernstein基底関数
はい、これが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曲面というものもあります。
定義式は
です。
B-spline曲面
同様に基底関数をBernstein基底関数からB-spline基底関数に取り替えましょう!
曲線のときと同様に、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)
ここでのフィッティングの計算には最小二乗法を使って以下の値を最小化するような制御点の位置
ノット列の位置を変更して、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-細分と呼ばれます。(ノット列の間隔を表す記号に
一方で、次数を増やす細分の操作は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