3次元曲線の離散化について

ちょげん

曲線を離散化する際によく使われる方法は折れ線で近似する方法である．
もっと良さげな方法を思いついたのでメモする

0. フレネの公式

離散化はフレネの公式をベースにするためその簡単な解説から始める． フレネの公式は曲線に関する3つのベクトルの変化についての公式である． 3つのベクトルの1つ目は接線ベクトル\displaystyle{T}である．接線と同じ方向で大きさが1となるベクトルである．2つ目は主法線ベクトル\displaystyle{N}である．これは\displaystyle{T}の変化方向を向き大きさが1となるベクトルである．3つ目は従法線ベクトル\displaystyle{B}であり，これは\displaystyle{T}，\displaystyle{N}の両者に直交し大きさ1のベクトルである．この3つのベクトルは基本的に直交しており（直交を説明するにはより厳密な定義が必要），少なくとも独立している．フレネの公式は曲線\displaystyle{p(t)}上の各点で

\displaystyle{ \begin{align*} T' &= & \kappa vN & \\ N' &= -\kappa v T & & + \tau v B\\ B' &= & -\tau v N & \end{align*} }

が成り立つ法則である．ここでダッシュ記号は\displaystyle{t}についての微分を示し，\displaystyle{v, \kappa, \tau}はそれぞれtについての実数値関数で速度，曲率，捩じれ率と呼ばれる．

1.曲線のテイラー展開

曲線\displaystyle{p(t)}が\displaystyle{t=(a, b)}で定義されており，この区間を速度，曲率，捩じれ率がそれぞれ一定と見なせる程細かく区分する．その区間の幅を\displaystyle{l}とする．
この時，曲線を\displaystyle{p(t)}周りにテイラー展開すると，

\displaystyle{ p(t+l)=p(t)+p^{(1)}(t)l+\frac{1}{2!}p^{(2)}(t)l^2+\frac{1}{3!}p^{(3)}(t)l^3+... }

となる．ここで\displaystyle{p^{(k)}(t)}は\displaystyle{p(t)}のk回微分係数であり，これを求めることで\displaystyle{p(t+l)}が求まる．その導出にフレネの公式を用いる．

2.\displaystyle{p^{(k)}(t)}の導出

\displaystyle{p^{(1)}}は大きさが\displaystyle{v}で向きが\displaystyle{T}のベクトルなので

\displaystyle{ p^{(1)}=vT }

となる．検討している区間では\displaystyle{v, \kappa, \tau}がそれぞれ一定なので

\displaystyle{ p^{(2)}=vT' }

である．\displaystyle{T'}はフレネの公式より

\displaystyle{ p^{(2)}=vT'=v(\kappa vN)= \kappa v^2 N }

となる．同様にして

\displaystyle{ \begin{align*} p^{(1)} & & &=vT\\ p^{(2)} &=& vT'=v(\kappa vN) &= \kappa v^2 N\\ p^{(3)} &=& \kappa v^2 N' &= \kappa v^3 (-\kappa T + \tau B )\\ p^{(4)} &=& \kappa v^3 (-\kappa T' + \tau B' ) &= -\kappa v^4 (\kappa^2 + \tau^2) N\\ p^{(5)} &=& -\kappa v^4 (\kappa^2 + \tau^2) N' &= -\kappa v^5 (\kappa^2 + \tau^2) (-\kappa T + \tau B)\\ \vdots& & \end{align*} }

となる．これをそれぞれ代入して

\displaystyle{ \begin{align*} p(t+l) =&p(t) +vlT +\frac{1}{2!}\kappa v^2 Nl^2 +\frac{1}{3!}\kappa v^3 (-\kappa T + \tau B ) l^3 -\frac{1}{4!}\kappa v^4 (\kappa^2 + \tau^2) N l^4 +...\\ =&p(t)+vlT +\left\{ \frac{1}{2!}\kappa v^2 l^2 -\frac{1}{4!}\kappa v^4 (\kappa^2 + \tau^2) l^4 +\frac{1}{6!}\kappa v^6 (\kappa^2 + \tau^2)^2 l^6 - ... \right\}N\\ &+\left\{ \frac{1}{3!}\kappa v^3 l^3 -\frac{1}{5!}\kappa v^5 (\kappa^2 + \tau^2) l^5 +\frac{1}{7!}\kappa v^7 (\kappa^2 + \tau^2)^2 l^7 -... \right\}(-\kappa T +\tau B) \end{align*} }

となる．Nの係数は

\displaystyle{ \begin{align*} &\frac{1}{2!}\kappa v^2 l^2 -\frac{1}{4!}\kappa v^4 (\kappa^2 + \tau^2) l^4 +\frac{1}{6!}\kappa v^6 (\kappa^2 + \tau^2)^2 l^6 - ... +(-1)^{k-1}\frac{1}{2k!}\kappa v^{2k}(\kappa^2 + \tau^2)^{k-1}l^{2k}\\ &=\frac{\kappa}{(\kappa^2 + \tau^2)}\left\{ \frac{1}{2!} v^2 (\kappa^2 + \tau^2) l^2 -\frac{1}{4!} v^4 (\kappa^2 + \tau^2)^2 l^4+... +(-1)^{k-1}\frac{1}{2k!}v^{2k}(\kappa^2 + \tau^2)^{k}l^{2k} + ... \right\}\\ \end{align*} }

とでき，\displaystyle{d^ 2=\kappa^ 2+\tau^ 2}として，

\displaystyle{ \begin{align*} &=\frac{\kappa}{(\kappa^2 + \tau^2)}\left\{ \frac{1}{2!} v^2 d^2 l^2 -\frac{1}{4!} v^4 d^4 l^4+... +(-1)^{k-1}\frac{1}{2k!}v^{2k}d^{2k}l^{2k} + ... \right\}\\ &=\frac{\kappa}{d^2}\{1-\cos(vdl)\} \end{align*} }

である．\displaystyle{-\kappa T+ \tau B}の係数は

\displaystyle{ \begin{align*} &\left\{ \frac{1}{3!}\kappa v^3 l^3 -\frac{1}{5!}\kappa v^5 (\kappa^2 + \tau^2) l^5 +\frac{1}{7!}\kappa v^7 (\kappa^2 + \tau^2)^2 l^7 -... \right\}\\ &=\kappa \left\{ \frac{1}{3!} v^3 l^3 -\frac{1}{5!} v^5 d^2 l^5 +\frac{1}{7!} v^7 d^4 l^7 -... \right\}\\ &=\frac{\kappa}{d^3} \left\{ \frac{1}{3!} v^3 d^3 l^3 -\frac{1}{5!} v^5 d^5 l^5 +\frac{1}{7!} v^7 d^7 l^7 -... \right\}\\ &=\frac{\kappa}{d^3}\{vdl-\sin(vdl)\} \end{align*} }

である．まとめると，

\displaystyle{ \begin{align*} p(t+l)=&p(t)+vlT +\frac{\kappa}{d^2}\{1-\cos(vdl)\}N +\frac{\kappa}{d^3}\{vdl-\sin(vdl)\}(-\kappa T +\tau B)\\ =&p(t) -\frac{\kappa^2}{d^3}\{vdl-\frac{vd^3l}{\kappa^2}-\sin(vdl)\}T +\frac{\kappa}{d^2}\{1-\cos(vdl)\}N +\frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\}B\\ =&p(t) +\frac{\kappa^2}{d^3}\{vdl\frac{\tau^2}{\kappa^2}+\sin(vdl)\}T +\frac{\kappa}{d^2}\{1-\cos(vdl)\}N +\frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\}B\\ \end{align*} }

が得られる．

3.\displaystyle{l+t}でのフレーム

\displaystyle{p(t+l)}の計算では媒介係数\displaystyle{t}でのフレームが既知である前提で議論した．実際には求めた\displaystyle{p(t+l)}を\displaystyle{p(t)}に代入し\displaystyle{p(t+2l)}を求め，これを繰り返すことで目的の曲線を得る．そのためには媒介係数\displaystyle{t}の段階で\displaystyle{t+l}におけるフレーム\displaystyle{T(t+l),N(t+l),B(t+l)}が得られなければならない．本節ではそれを求める．
前節では以下の式を得た．

\displaystyle{ p(l+t)=p(t) +\frac{\kappa^2}{d^3}\{vdl\frac{\tau^2}{\kappa^2}+\sin(vdl)\}T +\frac{\kappa}{d^2}\{1-\cos(vdl)\}N +\frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\}B }

これは媒介係数を\displaystyle{l}とした曲線と見ることができる．また\displaystyle{T,N,B}の係数を取り出し，

\displaystyle{ c(l)=( \frac{\kappa^2}{d^3}\{vdl\frac{\tau^2}{\kappa^2}+\sin(vdl)\}, \frac{\kappa}{d^2}\{1-\cos(vdl)\}, \frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\}) }

は座標軸を\displaystyle{T,N,B}とした座標系における曲線と見ることができる.
\displaystyle{T(t+l)}は\displaystyle{c(l)}を通り，接線方向を向き，大きさが1のベクトルである． 接線方向は

\displaystyle{ \begin{align*} \frac{d}{dl}c(l)=& \frac{d}{dl}( \frac{\kappa^2}{d^3}\{vdl\frac{\tau^2}{\kappa^2}+\sin(vdl)\}, \frac{\kappa}{d^2}\{1-\cos(vdl)\}, \frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\} )\\ =&( \frac{\kappa^2}{d^3}\{vd\frac{\tau^2}{\kappa^2}+\cos(vdl)vd\}, \frac{\kappa}{d}\{\sin(vdl)v\}, \frac{\kappa\tau}{d^2}\{v-\cos(vdl)v\} )\\ =&\kappa v( \frac{\kappa}{d^2}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}, \frac{1}{d}\{\sin(vdl)\}, \frac{\tau}{d^2}\{1-\cos(vdl)\} )\\ \end{align*} }

である．大きさは

\displaystyle{ \begin{align*} ||\frac{1}{\kappa v}\frac{d}{dl}c(l)||^2 &= \frac{\kappa^2}{d^4}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}^2+ \frac{1}{d^2}\{\sin^2(vdl)\}+ \frac{\tau^2}{d^4}\{1-\cos(vdl)\}^2\\ & \text{↓第1,3項を展開}\\ &= \frac{\kappa^2}{d^4}\{ \frac{\tau^4}{\kappa^4} +\cos^2(vdl) +2\frac{\tau^2}{\kappa^2}\cos(vdl)\} +\frac{1}{d^2}\{\sin^2(vdl)\} +\frac{\tau^2}{d^4}\{1+\cos^2(vdl)-2\cos(vdl)\}\\ & \text{↓第3,7項を相殺，第2,6項を加算}\\ &= \frac{\tau^4}{\kappa^2 d^4} +\frac{1}{d^2}\cos^2(vdl) +\frac{1}{d^2}\{\sin^2(vdl)\} +\frac{\tau^2}{d^4}\\ & \text{↓第1,4項を加算}\\ &= \frac{\tau^4+\kappa^2\tau^2}{\kappa^2 d^4} +\frac{1}{d^2}\\ &= \frac{\tau^2}{\kappa^2 d^2} +\frac{1}{d^2}\\ &= \frac{\tau^2+\kappa^2}{\kappa^2 d^2}\\ ||\frac{1}{\kappa v}\frac{d}{dl}c(l)||^2 &= \frac{1}{\kappa^2}\\ \end{align*} }

より\displaystyle{||\frac{d}{dl}c(l)||=v}である．これより\displaystyle{c(l)}が前提である「領域内での速度は\displaystyle{v}で一定」と矛盾しない事が確認された．また\displaystyle{T(t+l)}を大きさ1にするため，\displaystyle{\frac{d}{dl}c(l)}を\displaystyle{v}で割って，\displaystyle{T(t+l)}は

\displaystyle{ \begin{align*} T(t+l)&= \kappa( \frac{\kappa}{d^2}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}, \frac{1}{d}\sin(vdl), \frac{\tau}{d^2}\{1-\cos(vdl)\} )\\ &= ( \frac{\tau^2}{d^2}+\frac{\kappa^2}{d^2}\cos(vdl), \frac{\kappa}{d}\sin(vdl), \frac{\kappa\tau}{d^2}\{1-\cos(vdl)\} ) \end{align*} }

である．\displaystyle{N(t+l)}は\displaystyle{T(t+l)}の変化方向で大きさが1のベクトルであるため，\displaystyle{T(t+l)}の微分は

\displaystyle{ \begin{align*} \frac{d}{dl}T(t+l) &= \kappa ( -\frac{\kappa}{d^2}\sin(vdl)vd, \frac{1}{d}\cos(vdl)vd, \frac{\tau}{d^2}\sin(vdl)vd )\\ & \text{↓係数をまとめる}\\ &= \kappa v ( -\frac{\kappa}{d}\sin(vdl), \cos(vdl), \frac{\tau}{d}\sin(vdl) )\\ \end{align*} }

である．大きさは

\displaystyle{ \begin{align*} ||\frac{d}{dl}T(t+l)||^2 &= \kappa^2v^2 ( \frac{\kappa^2}{d^2}\sin^2(vdl) +\cos^2(vdl) +\frac{\tau^2}{d^2}\sin^2(vdl) )\\ &=\kappa^2v^2 (\cos^2(vdl)+\sin^2(vdl))=\kappa^2v^2\\ \end{align*} }

より\displaystyle{||\frac{d}{dl}T(t+l)||=\kappa v}であり，\displaystyle{N=1/\kappa v\frac{dT}{dl}}として，

\displaystyle{ \begin{align*} N(t+l)&=( -\frac{\kappa}{d}\sin(vdl), \cos(vdl), \frac{\tau}{d}\sin(vdl) )\\ \end{align*} }

である． ちなみに大きさは

\displaystyle{ \begin{align*} ||N(t+l)||^2 &= \cos^2(vdl) +\frac{\kappa^2}{d^2}\sin^2(vdl) +\frac{\tau^2}{d^2}\sin^2(vdl))=1\\ \end{align*} }

と確かに1である．また，\displaystyle{T(t+l)}と\displaystyle{N(t+l)}の内積は

\displaystyle{ \begin{align*} &T(t+l)\cdot N(t+l)\\ &=\kappa( -\frac{\kappa}{d^2}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}\frac{\kappa}{d}\sin(vdl) +\frac{1}{d}\sin(vdl)\cos(vdl) +\frac{\tau}{d^2}\{1-\cos(vdl)\}\frac{\tau}{d}\sin(vdl) )\\ & \text{↓第1,3項を展開}\\ &=\kappa( -\frac{\tau^2}{d^3}\sin(vdl) -\frac{\kappa^2}{d^3}\sin(vdl)\cos(vdl) +\frac{1}{d}\sin(vdl)\cos(vdl) +\frac{\tau^2}{d^3}\sin(vdl) -\frac{\tau^2}{d^3}\sin(vdl)\cos(vdl) )\\ & \text{↓第1,4項を相殺，2,5項を加算}\\ &=\kappa( -\frac{1}{d}\sin(vdl)\cos(vdl) +\frac{1}{d}\sin(vdl)\cos(vdl) )=0\\ \end{align*} }

より直交している事が確認できる． \displaystyle{B(t+l)}は\displaystyle{T(t+l)}と\displaystyle{N(t+l)}の両方に直交しているベクトルなので外積を使い，

\displaystyle{ \begin{align*} T(t+l)\times N(t+l) =&\kappa( \frac{1}{d}\{\sin(vdl)\}\frac{\tau}{d}\sin(vdl) -\frac{\tau}{d^2}\{1-\cos(vdl)\}\cos(vdl),\\ &-\frac{\kappa}{d^2}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}\frac{\tau}{d}\sin(vdl) -\frac{\tau}{d^2}\{1-\cos(vdl)\}\frac{\kappa}{d}\sin(vdl),\\ &\frac{\kappa}{d^2}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}\cos(vdl) +\frac{1}{d}\{\sin(vdl)\}\frac{\kappa}{d}\sin(vdl) )\\ =&\kappa( \frac{\tau}{d^2}\sin^2(vdl) -\frac{\tau}{d^2}\{1-\cos(vdl)\}\cos(vdl),\\ &-\frac{\kappa\tau}{d^3}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}\sin(vdl) -\frac{\kappa\tau}{d^3}\{1-\cos(vdl)\}\sin(vdl),\\ &\frac{\kappa}{d^2}\{\frac{\tau^2}{\kappa^2}+\cos(vdl)\}\cos(vdl) +\frac{\kappa}{d^2}\{\sin(vdl)\}\sin(vdl) )\\ =&\kappa( \frac{\tau}{d^2}\{1-\cos(vdl)\},\\ &-\frac{\tau^3}{\kappa d^3}\sin(vdl) -\frac{\kappa\tau}{d^3}\sin(vdl)\cos(vdl) -\frac{\kappa\tau}{d^3}\sin(vdl) +\frac{\kappa\tau}{d^3}\sin(vdl)\cos(vdl),\\ &\frac{\kappa}{d^2} +\frac{\kappa\tau^2}{d^2}\cos(vdl)) )\\ B(t+l)=&( \frac{\kappa\tau}{d^2}\{1-\cos(vdl)\}, -\frac{\tau}{d}\sin(vdl), \frac{\kappa^2}{d^2} +\frac{\tau^2}{d^2}\cos(vdl))\\ \end{align*} }

である．ちなみに大きさは

\displaystyle{ \begin{align*} ||B(t+l)||^2 =& \frac{\kappa^2\tau^2}{d^4}\{1-\cos(vdl)\}^2 +\frac{\tau^2}{d^2}\sin^2(vdl) +\left\{\frac{\kappa^2}{d^2}+\frac{\tau^2}{d^2}\cos(vdl)\right \}^2\\ & \text{↓第1項と第3項を展開}\\ =& \frac{\kappa^2\tau^2}{d^4}\{1-2\cos(vdl)+\cos^2(vdl)\} +\frac{\tau^2}{d^2}\sin^2(vdl) +\left\{\frac{\kappa^4}{d^4}+\frac{\tau^4}{d^4}\cos^2(vdl)+2\frac{\kappa^2\tau^2}{d^4}\cos(vdl) \right\}\\ & \text{↓$\cos(vdl)$の項を消去}\\ =& \frac{\tau^2}{d^2}\sin^2(vdl) +\frac{\kappa^2\tau^2+\kappa^4}{d^4} +\frac{\kappa^2\tau^2+\tau^4}{d^4}\cos^2(vdl)\\ & \text{↓$d^2=\kappa^2+\tau^2$を利用}\\ =& \frac{\tau^2}{d^2}\sin^2(vdl) +\frac{\kappa^2}{d^2} +\frac{\tau^2}{d^2}\cos^2(vdl)=1 \end{align*} }

より確かに1である．

4.まとめ

以上の議論をまとめると，\displaystyle{p(t+l)}は

\displaystyle{ \begin{align} p(l+t)=p(t) +\frac{\kappa^2}{d^3}\{vdl\frac{\tau^2}{\kappa^2}+\sin(vdl)\}T +\frac{\kappa}{d^2}\{1-\cos(vdl)\}N +\frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\}B\\ \end{align} }

であり，\displaystyle{T(t+l),N(t+l),B(t+l)}は，

\displaystyle{ \begin{align} T(t+l)&= \left[\frac{\tau^2}{d^2}+\frac{\kappa^2}{d^2}\cos(vdl)\right ]T(t) & +\left[\frac{\kappa}{d}\sin(vdl)\right ]N(t) & +\left[\frac{\kappa\tau}{d^2}\{1-\cos(vdl)\}\right ]B(t) \\ N(t+l)&= \left[-\frac{\kappa}{d}\sin(vdl)\right ]T(t) & +\left[\cos(vdl)\right ]N(t) & +\left[\frac{\tau}{d}\sin(vdl)\right ]B(t) \\ B(t+l)&= \left[\frac{\kappa\tau}{d^2}\{1-\cos(vdl)\}\right ]T(t) & +\left[-\frac{\tau}{d}\sin(vdl)\right ]N(t) & +\left[\frac{\kappa^2}{d^2} +\frac{\tau^2}{d^2}\cos(vdl))\right ]B(t)\\ \end{align} }

である．
これを曲線\displaystyle{p(t):(a, b)\rightarrow R^ 3}を求める計算に利用する事を考える．そのためにフレームの式の係数をまとめて，

\displaystyle{ R(vl, \kappa, \tau)=\left[ \begin{array}{ccc} \frac{\tau^2}{d^2}+\frac{\kappa^2}{d^2}\cos(vdl) & -\frac{\kappa}{d}\sin(vdl) & \frac{\kappa\tau}{d^2}\{1-\cos(vdl)\} \\ \frac{\kappa}{d}\sin(vdl) & \cos(vdl) & -\frac{\tau}{d}\sin(vdl) \\ \frac{\kappa\tau}{d^2}\{1-\cos(vdl)\} & \frac{\tau}{d}\sin(vdl) & \frac{\kappa^2}{d^2} +\frac{\tau^2}{d^2}\cos(vdl)) \\ \end{array} \right] }

とする．ちなみにこれは特殊直交行列となる．これよりフレームの式は

\displaystyle{ F(t+l)=F(t)R(vl, \kappa, \tau) }

となる．ここでフレームをまとめて\displaystyle{F(t)=[T(t),N(t),B(t)]}とした．また，

\displaystyle{ C(vl,\kappa,\tau)=\left[ \begin{array}{c} \frac{\kappa^2}{d^3}\{vdl\frac{\tau^2}{\kappa^2}+\sin(vdl)\} \\ \frac{\kappa}{d^2}\{1-\cos(vdl)\} \\ \frac{\kappa\tau}{d^3}\{vdl-\sin(vdl)\} \\ \end{array} \right] }

とすると，\displaystyle{p(t+l)}は

\displaystyle{ p(t+l)=p(t)+F(t)C(vl, \kappa, \tau) }

と表現できる． 区間\displaystyle{(a, b)}をそれぞれ速度，曲率，捩じれ率がそれぞれ一定と見なせるような領域になるように媒介係数を\displaystyle{n}個に等分割する．1つの領域の幅\displaystyle{l}は

\displaystyle{ l=L/n=(b-a)/n }

である．ここで\displaystyle{L=b-a}である．領域の境界点は

\displaystyle{ t_i=a+li\quad(i=0,1,...,n) }

である．この時微小領域は

\displaystyle{ I_i=(t_{i-1}, t_{i})\quad(i=1,2,...,n) }

内 と表記できる．\displaystyle{I_i}内でのパラメータを\displaystyle{v_i},\displaystyle{\kappa_i},\displaystyle{\tau_i}で与えられるとする．また，初期位置\displaystyle{p(a)}，初期フレーム\displaystyle{F(a)}は与えられるとする． i番目境界点でのフレームは

\displaystyle{ F_i=F(t_i)=F(a)\prod_{j=1}^iR(v_jl, \kappa_j, \tau_j) }

と書ける． また曲線は

\displaystyle{ p_i=p(t_i)=p(a)+\sum_{j=1}^iF_jC(v_jl, \kappa_j, \tau_j) =p(a)+F(a)\sum_{j=1}^i\prod_{k=1}^jR(v_kl, \kappa_k, \tau_k)C(v_jl, \kappa_j, \tau_j) }

となる．これをi=1からnまで繰り返し計算することで離散化された曲線を得られる．

謝辞

このブログページを作製するにあたり，自前の環境での数式からはてなブログの数式へ変換するためano3様のこちらのツールを使用させて頂きました．非常に有用なツールをありがとうございます．

参考文献

微分幾何の勉強はこれでしました．
Neill, Barrett. Elementary differential geometry. Burlington, MA San Diego, CA London: Academic Press, an imprint of Elsevier, 2006. Print.
https://www.amazon.co.jp/Elementary-Differential-Geometry-Revised-Second/dp/0120887355
英語で長く比較的丁寧に書いているのでフレネの公式だけについて知りたい方はwikipedia等のページを参照することをお勧めします． 一方で微分幾何について応用するのに十分に理解できるので時間がある方には一読をお勧めします．