BasicBSpline.jlで始めるB-spline

13.7K Views

July 11, 23

スライド概要

2023-07-11に開催される「数学と物理におけるJuliaの活用」で発表するスライドです。
https://akio-tomiya.github.io/julia_imi_workshop2023/

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

関連スライド

各ページのテキスト
1.

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

2.

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

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