BasicBSpline.jlで始めるB-spline

1.

BasicBSpline.jlで始める B‑spline 堀川由人 (@Hyrodium) 2023 年 7 月 11 日

2.

スライド公開先 https://www.docswell.com/s/hyrodium/5Q89MJ‑B‑spline 1

https://www.docswell.com/s/hyrodium/5Q89MJ-B-spline

3.

目次 1. 自己紹介 2. B‑spline 最速入門 3. B‑spline もう少し入門 4. B‑spline の応用例 5. まとめ 2

4.

目次 1. 自己紹介 2. B‑spline 最速入門 3. B‑spline もう少し入門 4. B‑spline の応用例 5. まとめ 3

5.

経歴・所属名前 • 堀川由人 • @Hyrodium 経歴 • 大阪府立工業高等専門学校 • 大阪大学 (学部・修士) • DMG 森精機現所属 • Rist • 大阪大学大学院基礎工学研究科 (招聘研究員) 4

6.

Julia と私 Zenn.dev • Julia の行列・ベクトルを完全に理解すっぞ！！ • FastGaussQuadrature.jl で数値積分しましょう • Quaternions.jl をメンテナンスしてる話 • Julia 言語における中置演算子の扱い GitHub • JuliaGeometry/Rotations.jl • JuliaApproximation/FastGaussQuadrature.jl • JuliaArrays/StaticArrays.jl • JuliaGeometry/Quaternions.jl • JuliaMath/IntervalSets.jl • hyrodium/ImageClipboard.jl • hyrodium/BasicBSpline.jl 5

7.

B‑spline と私これまでに執筆した記事など • B‑spline 入門（線形代数がすこし分かる人向け） • NURBS 多様体による形状表現 • BasicBSpline.jl のドキュメント • BasicBSpline.jl を作ったので宣伝です！ • 新しくなった BasicBSpline.jl(v0.8.3) の紹介 • なめらかな曲線を Julia で SVG 出力する B‑spline のモチベーションを伝えるための動画証明を詳しく書いた PDF 資料パッケージの使用方法と簡単な数学的解説 BasicBSpline.jl を作ったアナウンス on Zenn.dev バージョンが上がった件について解説 BasicBSpline.jl の応用例の紹介記事 • BasicBSpline.jl ‑ Basic mathematical operations for B‑spline • Plotting smooth graphs with Julia Discourse でのパッケージ告知「なめらかな曲線を Julia で SVG 出力する」の英訳版 6

8.

目次 1. 自己紹介 2. B‑spline 最速入門 B‑spline とは？関数近似 B‑spline 基底関数 3. B‑spline もう少し入門 4. B‑spline の応用例 5. まとめ 7

9.

B‑spline とは？区分多項式を使って関数近似するための道具 • 多項式閉区間上の任意の連続関数は多項式で近似可能 (Weierstrass の多項式近似定理) • 区分関数定義域を適当な区間に分割し、各区間で関数を近似 B‑spline で実現できること： • 高すぎない多項式次数 • 多すぎない区間分割数 • 区分点で滑らかに繋げる 8

10.

関数近似 ➀ (折れ線) 適当な点を折れ線で繋げば近似になる 9

11.

関数近似 ➀ (折れ線) 適当な点を折れ線で繋げば近似になる多項式の繋ぎ目 (区分点) を「ノット (knot)」と呼ぶ 10

12.

関数近似 ➁ (放物線) 各区間で 3 点取って放物線 (2 次多項式) を描けば近似になる 11

13.

関数近似 ➁ (放物線) 各区間で 3 点取って放物線 (2 次多項式) を描けば近似になる 12

14.

関数近似 ➁ (放物線) 各区間で 3 点取って放物線 (2 次多項式) を描けば近似になる区分点で C 0 級になる 13

15.

関数近似 ➂ (放物線、ただし C 1 級) 放物線で近似し、区分点で C 1 級も要求 14

16.

関数近似 ➂ (放物線、ただし C 1 級) 放物線で近似し、区分点で C 1 級も要求 15

17.

関数近似 ➂ (放物線、ただし C 1 級) 放物線で近似し、区分点で C 1 級も要求 → どうすれば実現できる⋯？ 16

18.

B‑spline 基底関数 ➀ (近似の方針) 「滑らかさを保証した区分多項式」をたくさん用意して、線型結合すれば OK！これが B‑spline 基底関数！ 17

19.

B‑spline 基底関数 ➁ (定義) 定義 (B‑spline 基底関数) ノット列 k = (k1 , . . . , kl ) に対する p 次 B‑spline 基底関数 B(i,p,k) は以下で定義される t − ki ki+p+1 − t B(i,p−1,k) (t) + B(i+1,p−1,k) (t) ki+p − ki ki+p+1 − ki+1   1 (k ≤ t < k ) i i+1 B(i,0,k) (t) =  0 (otherwise) B(i,p,k) (t) = 18

20.

B‑spline 基底関数 ➁ (定義) 19

21.

B‑spline 基底関数 ➂ (B‑spline 基底関数の張る空間) 定理 (B‑spline 基底関数の張る空間) 狭義単調増加なノット列 k = (k1 , . . . , kl ) で構成される B‑spline 基底関数について以下が成立 ( ) ∀i, f |(ki ,ki+1 ) は p 次多項式 p−1 span B(i,p,k) = f ∈ C (R) i supp(f ) ⊆ [k1 , kl ] この線形空間を P[p, k] で表す 20

22.

B‑spline 基底関数 ➃ (疑問点いろいろ) • 区分多項式空間 P[p, k] の次元は何で決まる？ • 2 つの線形空間 P[p1 , k 1 ], P[p2 , k 2 ] の包含関係は何で決まる？ • 特定のノット ki での滑らかさの制約を変更できる？ • B‑spline 基底関数 B(i,p,k) はどんな良い性質がある？ • B‑spline 基底関数 B(i,p,k) の導関数は？ • 21

23.

B‑spline 基底関数 ➃ (疑問点いろいろ) • 区分多項式空間 P[p, k] の次元は何で決まる？ • 2 つの線形空間 P[p1 , k 1 ], P[p2 , k 2 ] の包含関係は何で決まる？ • 特定のノット ki での滑らかさの制約を変更できる？ • B‑spline 基底関数 B(i,p,k) はどんな良い性質がある？ • B‑spline 基底関数 B(i,p,k) の導関数は？ • BasicBSpline.jl で何ができるのか？ 22

24.

[beta]

B‑spline 基底関数 ➃ (疑問点いろいろ)
• 区分多項式空間 P[p, k] の次元は何で決まる？
• 2 つの線形空間 P[p1 , k 1 ], P[p2 , k 2 ] の包含関係は何で決まる？
• 特定のノット ki での滑らかさの制約を変更できる？
• B‑spline 基底関数 B(i,p,k) はどんな良い性質がある？
• B‑spline 基底関数 B(i,p,k) の導関数は？
• BasicBSpline.jl で何ができるのか？
using BasicBSpline, Plots
P = BSplineSpace{2}(KnotVector([0,1,2,3,4,5,6]))
n = dim(P)
P ⊆ BSplineSpace{2}(KnotVector(0:8))
plot(P)
plot!(BSplineSpace{2}(KnotVector([0,1,2,3,3,4,5,6])))
plot!(t -> sum(bsplinebasis(P,i,t) for i in 1:n))
plot!(BSplineDerivativeSpace{1}(P))

↑ 完全理解しなくて OK。雰囲気だけ
23

25.

目次 1. 自己紹介 2. B‑spline 最速入門 3. B‑spline もう少し入門 B‑spline 基底関数の性質 B‑spline 多様体リファインメント 4. B‑spline の応用例 5. まとめ 24

26.

[beta]

B‑spline 基底関数の性質 ➀
• 台
supp(B(i,p,k) ) = [ki , ki+p+1 ]
• 1 の分割
P
i B(i,p,k) (t) = 1 (t ∈ [k1+p , kl−p ])
k = KnotVector(1:8)
P = BSplineSpace{2}(k)
n = dim(P)
plot(P)
scatter!(
collect(k),
zeros(length(k))
)
bsplinesupport(P, 2)
plot!(t->sum(
bsplinebasis(P, i, t)
for i in 1:n)
)
25

27.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 26

28.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 27

29.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 28

30.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 29

31.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 30

32.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 31

33.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 32

34.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 33

35.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ 34

36.

B‑spline 基底関数の性質 ➁ ノット列が重なるとどうなるか？ • 各点収束している • 収束先で滑らかさが落ちている 35

37.

B‑spline 基底関数の定義 (ふたたび) ノット列の重複があるとき： • 各点収束している • 収束先で滑らかさが落ちているしかし、B‑spline 基底関数の定義はノット列の重複を許さないはず ki+p+1 − t t − ki B(i,p−1,k) (t) + B(i+1,p−1,k) (t) B(i,p,k) (t) = ki+p − ki ki+p+1 − ki+1   1 (k ≤ t < k ) i i+1 B(i,0,k) (t) =  0 (otherwise) 36

38.

B‑spline 基底関数の定義 (ふたたび) ノット列の重複があるとき： • 各点収束している • 収束先で滑らかさが落ちているしかし、B‑spline 基底関数の定義はノット列の重複を許さないはず ki+p+1 − t t − ki B(i,p−1,k) (t) + B(i+1,p−1,k) (t) B(i,p,k) (t) = ki+p − ki ki+p+1 − ki+1   1 (k ≤ t < k ) i i+1 B(i,0,k) (t) =  0 (otherwise) ゼロ除算をゼロとすれば、収束先の関数が得られる！区分多項式空間 P[p, k] の定義も拡張される   ∀i, f |(ki ,ki+1 ) は p 次多項式       supp(f ) ⊆ [k1 , kl ] P[p, k] = span B(i,p,k) = f : R → R   i     f は t の近傍で C p−nk (t) 級 nk (t) = # { i | t = ki } 37

39.

一旦復習これまでのあらすじ： (1) 与えられた関数 f : I → R を区分多項式で近似したい (I: 有界閉区間) (2) B‑spline 基底関数 B(i,p,k) の線型結合で近似するのが良さそう X f (t) ≈ ai B(i,p,k) (t) i (3) ノット列 k に重複があると B(i,p,k) の滑らかさが落ちる (4) B(i,p,k) の張る関数空間が P[p, k] 次に気になること： • 関数だけでなく、曲線 f : I → V を近似できるか？係数 ai ∈ R を線形空間の元 ai ∈ V に置き換え X f (t) ≈ ai B(i,p,k) (t) i • 1 変数関数だけでなく、2 変数関数 f : I × J → R を近似できるか？ B‑spline 基底関数のテンソル積 X f (u1 , u2 ) ≈ aij B(i,p1 ,k1 ) (u1 )B(j,p1 ,k2 ) (u2 ) i,j 38

40.

一旦復習これまでのあらすじ： (1) 与えられた関数 f : I → R を区分多項式で近似したい (I: 有界閉区間) (2) B‑spline 基底関数 B(i,p,k) の線型結合で近似するのが良さそう X f (t) ≈ ai B(i,p,k) (t) i (3) ノット列 k に重複があると B(i,p,k) の滑らかさが落ちる (4) B(i,p,k) の張る関数空間が P[p, k] 次に気になること： • 関数だけでなく、曲線 f : I → V を近似できるか？係数 ai ∈ R を線形空間の元 ai ∈ V に置き換え X f (t) ≈ ai B(i,p,k) (t) i • 1 変数関数だけでなく、2 変数関数 f : I × J → R を近似できるか？ B‑spline 基底関数のテンソル積 X f (u1 , u2 ) ≈ aij B(i,p1 ,k1 ) (u1 )B(j,p1 ,k2 ) (u2 ) i,j 39

41.

一旦復習これまでのあらすじ： (1) 与えられた関数 f : I → R を区分多項式で近似したい (I: 有界閉区間) (2) B‑spline 基底関数 B(i,p,k) の線型結合で近似するのが良さそう X f (t) ≈ ai B(i,p,k) (t) i (3) ノット列 k に重複があると B(i,p,k) の滑らかさが落ちる (4) B(i,p,k) の張る関数空間が P[p, k] 次に気になること： • 関数だけでなく、曲線 f : I → V を近似できるか？係数 ai ∈ R を線形空間の元 ai ∈ V に置き換え X f (t) ≈ ai B(i,p,k) (t) i • 1 変数関数だけでなく、2 変数関数 f : I × J → R を近似できるか？ B‑spline 基底関数のテンソル積 X f (u1 , u2 ) ≈ aij B(i,p1 ,k1 ) (u1 )B(j,p1 ,k2 ) (u2 ) i,j 40

42.

B‑spline 曲線 ➀ 定義 (B‑spline 曲線) V を R‑線形空間、ai を V 上の点とする。このとき X p(t) = ai B(i,p,k) (t) (t ∈ [k1+p , kl−p ]) i で表される曲線を B‑spline 曲線と呼ぶ。ai を制御点と呼ぶ。 using StaticArrays # ノット列の定義 k = KnotVector([0,1,2,3,3,4,5,6,7]) # 区分多項式空間の定義 P = BSplineSpace{2}(k) n = dim(P) # 制御点の定義 a = [SVector(cos(i),sin(i)) for i in 0:n-1] # B-spline 曲線の定義 C = BSplineManifold(a, P) plot(C) 41

43.

B‑spline 曲線 ➁ 定義 (B‑spline 曲線) V を R‑線形空間、ai を V 上の点とする。このとき X p(t) = ai B(i,p,k) (t) (t ∈ [k1+p , kl−p ]) i で表される曲線を B‑spline 曲線と呼ぶ。ai を制御点と呼ぶ。端点を通るとは限らない！ 42

44.

[beta]

B‑spline 曲線 ➂
制御点の平行移動で B‑spline 曲線も平行移動する
a2 = [p+SVector(1,2) for p in a]
C2 = BSplineManifold(a2, P)
plot!(C2)

より一般に次の命題が成立
命題 (Affine 可換性)
f : x 7→ Ax + b を Affine 変換とし、制御点 {ai } から構成される B‑spline
曲線を p({ai }; t) で表すとする。このとき次が成立。
p({f (ai )}; t) = f (p({ai }; t))
証明. 1 の分割の性質を使う
X
X
X
p({f (ai )}; t) =
(Aai + b)B(i,p,k) (t) = A
ai B(i,p,k) (t) + b
B(i,p,k) (t)
i

=A

X
i

i

ai B(i,p,k) (t) + b = f (p({ai }; t))

i

□
43

45.

B‑spline 曲面定義 (B‑spline 曲面) V を R‑線形空間、aij を V 上の点とする。このとき X p(u1 , u2 ) = aij B(i,p1 ,k1 ) (u1 )B(j,p2 ,k2 ) (u2 ) i,j 1 1 2 2 ((u , u ) ∈ [k1+p 1 , kl1 −p1 ] × [k1+p2 , kl2 −p2 ]) 1 2 で表される曲面を B‑spline 曲面と呼ぶ。aij を制御点と呼ぶ。 # P # a 区分多項式空間の定義 = BSplineSpace{3}(KnotVector(0:7)) 制御点の定義 = [SVector(i-1.5, j-1.5, 1≤i≤2&&1≤j≤2) for i in 0:3, j in 0:3] # B-spline 曲線の定義 M = BSplineManifold(a, P, P) plotly() plot(M) gr() 44

46.

B‑spline 曲面定義 (B‑spline 曲面) V を R‑線形空間、aij を V 上の点とする。このとき X p(u1 , u2 ) = aij B(i,p1 ,k1 ) (u1 )B(j,p2 ,k2 ) (u2 ) i,j 1 1 2 2 ((u , u ) ∈ [k1+p 1 , kl1 −p1 ] × [k1+p2 , kl2 −p2 ]) 1 2 で表される曲面を B‑spline 曲面と呼ぶ。aij を制御点と呼ぶ。 3 次元以上も同様に定義可能 45

47.

GUI アプリケーションでの例 Blender Inkscape B‑spline 曲面を編集可能 Bézier 曲線を編集可能 46

48.

リファインメント制御点の数を増やす操作のこと元の B‑spline 曲線 (P[p, k]) リファインメント後 (P[p′ , k ′ ]) (1) 新しい次数 p′ とノット列 k′ を決める (2) 次数とノット列に合わせて制御点 a′j を配置 p(t) = n X ′ ai B(i,p,k) (t) = i=1 B(i,p,k) = n′ X j=1 n X a′j B(j,p′ ,k′ ) (t) j=1 Aij B(j,p′ ,k′ ) a′j = n X i=1 ai Aij 47

49.

[beta]

区分多項式空間の包含関係
P
B(i,p,k) = j Aij B(j,p′ ,k′ ) を満たす行列 A が存在するには P[p, k] ⊆ P[p′ , k ′ ]
でなければならない
定理 (区分多項式空間の包含関係)
P[p, k] ⊆ P[p, k ′ ] ⇔ k ⊆ k′
p < p′ の場合はもう少し複雑な条件になる
k = KnotVector(0:10)
k′= k + KnotVector([3.4, 5.0, 6.3])
P = BSplineSpace{3}(k)
P′= BSplineSpace{3}(k′)
k ⊆ k′
P ⊆ P′
plot(P)
plot!(P′)
A = changebasis(P, P′)
48

50.

目次 1. 自己紹介 2. B‑spline 最速入門 3. B‑spline もう少し入門 4. B‑spline の応用例関数近似滑らかなグラフを出力する内挿補間数値計算の例：編み紙 5. まとめ 49

51.

関数を近似するには？(伏線回収) 最小 2 乗法で次式を最小化して係数 {ai } を決めれば良い !2 Z X f (t) − ai B(i,p,k) (t) dt I i using BasicBSplineFitting f(x) = (1+x^2/10)*sin(exp(-(x-12)/4))+x/2+3 a = fittingcontrolpoints_I(f, P) 50

52.

滑らかなグラフを出力するには？➀ Plots.jl の出力グラフは折れ線 plot(sin,-8,8) の出力拡大滑らかなグラフが欲しい！ 51

53.

[beta]

滑らかなグラフを出力するには？➁
using BasicBSpline
using BasicBSplineFitting
using BasicBSplineExporter
using StaticArrays
f(t) = SVector(t,sin(t))
a, b = -8, 8
p = 3
k = KnotVector(a:b)+p*KnotVector([a,b])
P = BSplineSpace{p}(k)
a = fittingcontrolpoints_I(f,(P,))
C = BSplineManifold(a,(P,))
save_svg("sin.svg", C, xlims=(-10,10), ylims=(-2,2))

BasicBSplineExporter.jl パッケージで画像保存が可能！

52

54.

点列を内挿補間するには？➀ 以下の条件の補間を考える • p 次区分多項式で内挿補間 • 区分点とデータ点が一致 • 区分点の近傍で C p−1 級 1 次区分多項式では自明 53

55.

点列を内挿補間するには？➁ 2 次区分多項式では一意的ではないデータ点を通る区分多項式は無数に存在 (1 自由度) 54

56.

点列を内挿補間するには？➂ 3 次区分多項式でも一意的ではないデータ点を通る区分多項式は無数に存在 (2 自由度) 55

57.

点列を内挿補間するには？➃ 3 次区分多項式でも一意的ではないデータ点を通る区分多項式は無数に存在 (2 自由度) 両端での 2 階微分を 0 とすれば一意的！ 56

58.

[beta]

点列を内挿補間するには？➄
BasicBSpline.jl では内挿補間の API は無いが、実装は容易
function interpolate(xs, fs::AbstractVector{T}) where T
k = KnotVector(xs) + KnotVector([xs[1],xs[end]]) * 3
P = BSplineSpace{3}(k)
ddP = BSplineDerivativeSpace{2}(P)
dda = [bsplinebasis(ddP,j,xs[1]) for j in 1:dim(P)]
ddb = [bsplinebasis(ddP,j,xs[end]) for j in 1:dim(P)]
M = [bsplinebasis(P,j,x) for x in xs, j in 1:dim(P)]
M = vcat(dda', M, ddb')
y = vcat(zero(T), fs, zero(T))
return BSplineManifold(M\y, P)
end
ts = 0:8
F = interpolate(ts, sin.(ts))
scatter(ts, sin.(ts), label="data points")
plot!(t->F(t), 0, 8, label="interpolation")
57

59.

編み紙 ➀ 紙を編んで滑らかな曲面を作る話埋め込み形状 (上図中央) を B‑spline 多様体として数値計算 • 曲面片の弾性変形エネルギーを最小化 (偏微分方程式) • 制御点を未知変数として Galerkin 法が使える • B‑spline によって滑らかな (C 2 級の) 解が得られる • 3 次 B‑spline 曲線は SVG 画像として出力可能 58

60.

編み紙 ➁ 本日、曲面模型を持ってきています！詳細は JuliaCon2023 にて！ 59

61.

目次 1. 自己紹介 2. B‑spline 最速入門 3. B‑spline もう少し入門 4. B‑spline の応用例 5. まとめ 60

62.

まとめ • B‑spline は滑らかな区分多項式を扱うための道具「区分関数」と「多項式」両方の性質を持つ ♣ • B‑spline は数学的に奥深い区分多項式空間 P[p, k] の性質、など • B‑spline は応用先が広い幾何形状表現、内挿補間、Galerkin 法など • B‑spline を Julia で扱うには BasicBSpline.jl が便利！あなたの Star をお待ちしてます！ 61

63.

Julia 言語の良いところ計算が速い、以外にも色々ある • Unicode コーディング k1 ⊆ k2 , k1 + k2 , P1 ⊆ P2 • 抽象型・具象型 KnotVector <: AbstractKnotVector , UniformKnotVector <: AbstractKnotVector • パラメトリック型 BSplineSpace{p} • メタプログラミング P = BSplineSpace{3}(KnotVector(1:8)) @less bsplinebasisall(P, 1, 4.2) 62

64.

他パッケージとの比較 • sostock/BSplines.jl 22 • B‑spline の基本的な操作を実装 • 多項式次数 p が型パラメータに含まれないので遅い • jipolanco/BSplineKit.jl 40 • 応用を重視した B‑spline パッケージ • 悪い歴史的経緯を踏襲している (e.g. BSplineOrder ) • HoBeZwe/NURBS.jl 2 • 幾何形状表現への応用を意図したパッケージ • 実装途中の箇所が多い • kbarbary/Dierckx.jl 141 • Fortran ライブラリの DIERCKX のラッパー • Pure Julia ではない (e.g. 有理数型が使えない) • hyrodium/BasicBSpline.jl 76 • 数学志向で高速な B‑spline の実装 (公平な比較になっていないことに注意) 63

65.

参考文献 • Elaine Cohen, Richard F. Riesenfeld, Gershon Elber (2001) “Geometric Modeling with Splines: An Introduction” • Larry Schumaker (1989) “Spline Functions: Basic Theory” • 市田浩三, 吉本富士市 (1979) 『スプライン関数とその応用』 64

BasicBSpline.jlで始めるB-spline

Yuto Horikawa

関連スライド

逆数の作図からCayley変換まで[第一回日曜数学会]

Riemann球面に内接する直方体[第四回日曜数学会]

Creating smooth surface models with ElasticSurfaceEmbedding.jl

学振特別研究員になるために～2025年度申請版

研究に使える便利なフリーソフト ImageJ

大規模言語モデルに追加学習で専門知識を教える試み (2023, arXiv:2312.03360)

各ページのテキスト