はじめに
筆者は数学を勉強しているわけではないので、間違いがあれば指摘してください。ちゃんと勉強したい人はSEGAさんの基礎線形代数講座を読むことをお勧めします。
3dcgで回転を表現する方法として、回転行列、オイラー角、クォータニオンが主に挙げられると思います。それぞれの特徴をまとめました。
回転行列
3×3行列の直交行列で回転を表し、9つのパラメータを使います。直行行列である拘束条件(\(\mathrm{R R^T}=\mathrm{R^T R}=\mathrm{E}\))により、自由度は3になります。また、任意軸周りの回転を表すロドリゲスの回転行列に対応させることができ、以下のように\(n_x, n_y, n_z, θ\) の4つのパラメータで回転を表せます。回転軸nのノルムが1であるという拘束条件により、自由度が3つになることが分かります。
ロドリゲスの回転公式
\(
R_n(θ)= \begin{pmatrix}
\cosθ+n_x^2(1-\cosθ) & n_x n_y(1-\cosθ)-n_z\sinθ & n_z n_x(1-\cosθ)+n_y\sinθ \\
n_x n_y(1-\cosθ)+n_z\sinθ & \cosθ+n_y^2(1-\cosθ) & n_y n_z(1-\cosθ)-n_x\sinθ \\
n_z n_x(1-\cosθ)-n_y\sinθ & n_y n_z(1-\cosθ)-n_x\sinθ & \cosθ+n_z^2(1-\cosθ)
\end{pmatrix}
\)
※今回は右手系、内因性(intrinsic)の場合を考えます。
オイラー角では、3つの角度の組み合わせで回転を表現します。 剛体に固定された座標系で、「X軸周りにα、Y軸周りにβ、Z軸周りにγ回転させる」といった感じで表せ、直感的に分かりやすいです。
型が多い
まず、オイラー角の回転には親子関係があります。具体的にXYZ順の場合、Y軸を回転すればX軸が回転し、Z軸を回転すればY軸が回転し、それに起因してX軸も回転します。こちらにXYZ順の回転を実装したものを置いておきました。
親子関係によって軸が回転するということは、XYZの順とYXZの順では結果が変わってきます。この時点で6通りの表し方があります。また、最初の回転と最後の回転で同じ座標軸を使うことができます。1番目と2番目で同じ座標軸を使っても2番目は冗長になるので、1つ飛ばして同じ回転軸を使用するということです。1番目と3番目で同じ軸を使うものは狭義のオイラー角といい、違う軸を使う場合はタイトブライアン角とも言います。結局これらを考慮すると、以下のように12通りの回転の表し方が存在することなります。
ABC順 |
XYZ |
XZY |
YXZ |
YZX |
ZXY |
ZYX |
ABA順 |
XYX |
XZX |
YXY |
YZY |
ZXZ |
ZYZ |
※UnityではZXY順、blenderでは選択できるようになっています。
blenderでの回転順の選択箇所
XYZ順の場合、2番目の軸を±π/2 回転させた時にX軸とZ軸が平行になり、自由度が1つ失われるという現象が起きます。XYX順の場合は2番目の軸を0またはπ回転させた時に同様の現象が起きます。これをジンバルロックと言います。こちらのXYZ順でもβの値を90度か270度にしたとき、αとγによる回転がどちらもz軸周りの回転になることが確認できます。3軸の回転を表現したかったのに、2軸の回転に成り下がってしまうというのは計算上かなり問題になってきます。実際、アポロ計画ではあらかじめジンバルロックにならないようなコースを選び、手動で調節していたらしいです。
※航空学ではピッチ角(β)が90度になる可能性が非常に低いので今でもよく用いられるようです。
XYZ順の回転で、ジンバルロックが起きたときの回転行列を考えてみます。β=π/2 であるとき、
\begin{eqnarray}
R&=& \begin{pmatrix}
\cosγ & -\sinγ & 0 \\
\sinγ & \cosγ & 0 \\
0 & 0 & 1\end{pmatrix}\begin{pmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
-1 & 0 & 0
\end{pmatrix}\begin{pmatrix}
1 & 0 & 0 \\
0 & \cosα & -\sinα \\
0 & \sinα & \cosα
\end{pmatrix}\\[15px]&=& \begin{pmatrix}
0 & -\sinγ & \cosγ \\
0 & \cosγ & \sinγ \\
-1 & 0 & 0
\end{pmatrix}\begin{pmatrix}
1 & 0 & 0 \\
0 & \cosα & -\sinα \\
0 & \sinα & \cosα
\end{pmatrix}\\[15px]&=& \begin{pmatrix}
0 & \sin(α-γ) & \cos(α-γ) \\
0 & \cos(α-γ) & -\sin(α-γ) \\
-1 & 0 & 0
\end{pmatrix}
\end{eqnarray}
上記より、αとγの変化が回転行列Rに与える影響は同じであることが分かります。γの符号が負なので、αとγによる回転の向きは反対になります。x軸の回転はz軸周りの回転と等しくなることが分かっているので、αによる回転は余分になります。
せっかくなので回転行列からオイラー角(XYZ順)を求めてみましょう。
\begin{eqnarray}
R&=& \begin{pmatrix}
\cosγ & -\sinγ & 0 \\
\sinγ & \cosγ & 0 \\
0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
\cosβ & 0 & \sinβ \\
0 & 1 & 0 \\
-\sinβ & 0 & \cosβ
\end{pmatrix}\begin{pmatrix}
1 & 0 & 0 \\
0 & \cosα & -\sinα \\
0 & \sinα & \cosα
\end{pmatrix}\\[15px]&=& \begin{pmatrix}
\cosβ\cosγ & \sinα\sinβ\sinγ-\cosα\cosγ & \cosα\sinβ\cosγ+\sinα\sinγ\\
\cosβ\sinγ & \sinα\sinβ\sinγ+\cosα\cosγ & \cosα\sinβ\sinγ-\sinα\cosγ \\
-\sinβ & \sinα\cosβ & \cosα\cosβ
\end{pmatrix}
\end{eqnarray}
よって、
\(
m_{20}=-\sinβ
\)
\(\displaystyle β= \mathrm{arc} \sin (-m_{20}) \)
また、β≠±π/2のとき、
\(
m_{21}= \sinα\cosβ, m_{22}= \cosα\cosβ
\)
\(\displaystyle \sin α=\frac{m_{21}}{\cosβ}\)
\(\displaystyle \cos α = \frac{m_{22}}{\cosβ}\)
\(\displaystyle \tan α=\frac{\sinα}{\cosα} = \frac{m_{21}}{m_{22}}\)
\( \displaystyle α= \mathrm{arc} \tan ( \frac{m_{21}}{m_{22}} ) \)
同様に、
\(m_{10}= \cosβ\sinγ, m_{00}= \cosβ\cosγ\)
\(\displaystyle \sin γ=\frac{m_{10}}{\cosβ}\)
\(\displaystyle \cos γ = \frac{m_{00}}{\cosβ}\)
\(\displaystyle \tan γ=\frac{\sinγ}{\cosγ} = \frac{m_{10}}{m_{00}}\)
\( \displaystyle γ= \mathrm{arc} \tan ( \frac{m_{10}}{m_{00}} ) \)
まとめ
β≠±π/2のとき、
\(\text{}
\begin{cases}
\displaystyle α= \mathrm{arc} \tan ( \frac{m_{21}}{m_{22}} ) \\
\displaystyle β= \mathrm{arc} \sin (-m_{20}) \\
\displaystyle γ= \mathrm{arc} \tan ( \frac{m_{10}}{m_{00}} )
\end{cases}\)
β=±π/2のとき、
\(\text{}
\begin{cases}
\displaystyle α= 0 \\
\displaystyle γ= \mathrm{arc} \tan ( -\frac{m_{01}}{m_{11}} )
\end{cases}\)
先ほど取り上げた回転行列ではパラメータが9つあり、余分にメモリが取られます。またオイラー角には、ジンバルロックが存在します。こういった問題に遭遇しない、ちょうど良い回転表現がクォータニオンになります。
クォータニオンは四元数とも呼ばれます。複素数は虚数成分iを用いて\(x+iy\)と表すことができますが、これに虚数成分jとkを加えたものがクォータニオンになります。複素数は2次元の回転を扱うことができるので、それを応用して3次元の回転を表現しようとしたわけです。本来なら虚数成分はiとjの2つで済ませたいところですが、それでは上手くいかなかったためkを追加しました。それで四元数となっています。
2次元の複素数平面では1-i平面の1平面のみで回転させていました。クォータニオンは4次元なので、6平面で回転させる必要があります。
平面 |
乗算 |
計算結果 |
条件 |
1-i |
×i |
1→i→-1→-i→1 |
|
1-j |
×j |
1→j→-1→-j→1 |
|
1-k |
×k |
1→k→-1→-k→1 |
|
j-k |
×i |
j→k→-j→-k→j |
ij=k, ik=-j |
k-i |
×j |
k→i→-k→-i→k |
jk=i, ji=-k |
i-j |
×k |
i→j→-i→-j→i |
ki=j, kj=-i |
※i→j→-i→-j→iとなるのが理想でしたが、j→ij→-j→-ij→jとなって上手く行かなかったので、kを導入して4次元にすることで辻褄を合わせたという感じです。
しかし、iをかける(\(\frac{π}{2}\) 回転させる)と、1-i平面とj-k平面で回転してしまうことが分かります。同様に、jをかけると1-j平面とk-i平面で回転し、kをかけると1-k平面とi-j平面で回転してしまいます。2つの平面で回転するのは困るので、何とかして1つの平面だけを回転させたいです。ここで、右から-iをかけたときを考えてみます。
平面 |
乗算 |
計算結果 |
1-i |
×(-i) |
1→-i→-1→i→1 |
j-k |
×(-i) |
j→k→-j→-k→j |
左からiをかけたときと比較して、1-i平面では反対向き、j-k平面では同じ向きに回っていることが分かります。このことから、回転を表現するクォータニオンを、\( q= a\cos \frac{θ}{2}+i b\sin \frac{θ}{2}+ j c\sin \frac{θ}{2}+k d\sin \frac{θ}{2}\)と定義し、左からクォータニオン、右から共役クォータニオンを乗算することで回転を表現するようにしました。\(\frac{θ}{2}\)となっているのは、左右からクォータニオンを乗算したときに回したい角度の2倍回転してしまうからです。点p=(x,y,z)を回転させたものをp'と置いたとき、
\begin{equation}\boldsymbol{p'= qp\bar{q}} \tag{1}\end{equation}
p'のi, j, k成分が回転後の座標になっています。
回転行列による表現
地道に直接計算しても良いです。
(1)式に対して、\(q=q_{0}+ \boldsymbol{q}\)と置くと(\(q_{0}\)はスカラー部)
\begin{eqnarray}
\boldsymbol{p'} &=& q\boldsymbol{p}\bar{q} \\
&=& (q_{0}+\boldsymbol{q})\boldsymbol{p}(q_{0}-\boldsymbol{q}) \\
&=& q_{0}^2 \boldsymbol{p} + q_{0}\boldsymbol{(qp-pq)} -\boldsymbol{qpq} \\&=& q_{0}^2 \boldsymbol{p} + 2q_{0}(\boldsymbol{q×p})-\{-\boldsymbol{(p・q)q+(q・q)p-(q・p)q-q・(p×q)} \tag{2}\}\\&=& (q_{0}^2-\boldsymbol{q・q}) \boldsymbol{p} + 2\boldsymbol{(q・p)q}+2q_{0}(\boldsymbol{q×p}) \tag{3}
\end{eqnarray}
※(2)式は\(\boldsymbol{pq}=-\boldsymbol{p・q}+\boldsymbol{p×q}\)、(3)式はq同士の外積が0になることを利用
\( \boldsymbol{p'}=\begin{pmatrix}p_{1}' \\ p_{2}' \\ p_{3}' \end{pmatrix}\) \(\boldsymbol{p}=\begin{pmatrix} p_{1} \\ p_{2} \\ p_{3} \end{pmatrix}\) \(\boldsymbol{q}=\begin{pmatrix} q_{1} \\ q_{2} \\ q_{3} \end{pmatrix}\)と置くと、以下のように計算できます。
\begin{eqnarray}
p_{i}'&=&(q_{0}^2-\boldsymbol{q・q}) p_{i}+ 2\sum_{j=1}^{3} q_{j}p_{j}q_{i}+2q_{0}\sum_{k,j=1}^{3}ε_{ikj}q_{k}p_{j} \\ &=& \sum_{j=1}^{3} \{ (q_{0}^2-\boldsymbol{q・q}) δ_{ij}+2q_{i}q_{j}-2q_{0}\sum_{k=1}^{3}ε_{ijk}q_{k} \}p_{j}\end{eqnarray}
\begin{eqnarray}R_{ij}=(q_{0}^2-\boldsymbol{q・q})δ_{ij}+2q_{i}q_{j}-2q_{0}\sum_{k=1}^{3}ε_{ijk}q_{k}\end{eqnarray}
※\(ε_{ijk}\)は3階のエディントンのイプシロン
ここから、以下の行列式が求まります。
\begin{eqnarray}\begin{pmatrix}p_{1}' \\p_{2}' \\p_{3}' \end{pmatrix}=\begin{pmatrix}
q_{0}^2+q_{1}^2-q_{2}^2-q_{3}^2&2(q_{1}q_{2}-q_{0}q_{3})&2(q_{1}q_{3}+q_{0}q_{2}) \\
2(q_{2}q_{1}+q_{0}q_{3})&q_{0}^2-q_{1}^2+q_{2}^2-q_{3}^2&2(q_{2}q_{3}-q_{0}q_{1}) \\
2(q_{3}q_{1}-q_{0}q_{2})&2(q_{3}q_{2}+q_{0}q_{1})&q_{0}^2-q_{1}^2-q_{2}^2+q_{3}^2
\end{pmatrix}\begin{pmatrix}p_{1} \\p_{2} \\p_{3} \end{pmatrix}\end{eqnarray}
これでクォータニオンから回転行列への変換ができました。
ロドリゲスの回転公式
(3)式で、\(q_{0}=\cos \frac{θ}{2}、\boldsymbol{q}=\sin \frac{θ}{2}\boldsymbol{u}\)と置くと、
\begin{eqnarray}
\boldsymbol{p'} &=& (q_{0}^2-\boldsymbol{q・q}) \boldsymbol{p} + 2\boldsymbol{(q・p)q}+2q_{0}(\boldsymbol{q×p})\\
&=& (\cos^2 \frac{θ}{2}-\sin^2 \frac{θ}{2}) \boldsymbol{p} + 2\sin^2 \frac{θ}{2}\boldsymbol{(u・p)u}+2\cos \frac{θ}{2} \sin \frac{θ}{2}(\boldsymbol{u×p})\\&=& \boldsymbol{p}\cosθ +(1-\cos θ)\boldsymbol{(p・u)u}+(\boldsymbol{u×p})\sin θ
\end{eqnarray}
この式はロドリゲスの回転公式のベクトル表記になります。このことから、回転行列とクォータニオンは本質的には同じものであることが分かります。
球面線形補間
クォータニオンのメリットとしてよく挙げられるのが球面線形補間(slerp)ができることです。考え方についてはCGのための数学のクォータニオンの部分で分かりやすく説明されています。2つのクォータニオン\(q_{1}、q_{2}\)に対して、補間されたクォータニオン\(q_{t}\)が以下のように表現されます。
\begin{eqnarray}\boldsymbol{q_{t}}=slerp(q_{1},q_{2},t)=\frac{\sin(θ(1-t))}{\sinθ}\boldsymbol{q_{1}}+\frac{\sin(θt)}{\sinθ}\boldsymbol{q_{2}}
\end{eqnarray}
しかし2つ以上のキーがある回転の間を移動する場合、slerpでは補間の間にガタツキが生じます。ガタツキの例はこちらで分かりやすく紹介されています。これを避ける方法がEberlyさんにより提供されています。この方法では以下のように、クォータニオン\(a_{i}、a_{i+1}\)を\(q_{i}、q_{i+1}\)の間に導入します。
\begin{eqnarray}\boldsymbol{a_{i}}=\boldsymbol{q_{i}} exp[-\frac{log(\boldsymbol{q_{i}^{-1} q_{i-1}})+log(\boldsymbol{{q_{i}^{-1} q_{i+1}}}) } {4} ] \end{eqnarray}
3次スプラインによるクォータニオンの球面線形補間であるsquad(Shoemakeのクォータニオン曲線)を以下のように求めることで、上手く球面線形補間できます。
\begin{eqnarray}squad(q_{i},q_{i+1},a_{i},a_{i+1},t)=slerp(slerp(q_{i},q_{i+1},t),slerp(a_{i},a_{
i+1},t),2t(1-t)) \end{eqnarray}
逆行列のように、あるクォータニオンに乗算したときに1になるクォータニオンのことを逆クォータニオンと言います。\(q=q_{0}+ \boldsymbol{q}\)と置くと、共役クォータニオンと掛け合わせたとき、
\begin{eqnarray}q\bar{q}=(q_{0}+\boldsymbol{q})(q_{0}-\boldsymbol{q})=(q_{0}^2+\boldsymbol{q・q},-q_{0}\boldsymbol{q}+q_{0}\boldsymbol{q}-\boldsymbol{q×q})=|q|^2 \end{eqnarray}
※\(\boldsymbol{pq}=-\boldsymbol{p・q}+\boldsymbol{p×q}\)を利用
以上より、逆クォータニオンはクォータニオンのノルムと共役を用いて以下のように表せます。
\begin{eqnarray}q^{-1}=\frac{\bar{q}}{|q|^2} \end{eqnarray}
デュアルクォータニオン
デュアルクォータニオンとは
先ほどはクォータニオンを使って回転を表現しましたが、回転に加え平行移動の表現も可能にするのがデュアルクォータニオン(クォータニオンの二次元ベクトル)です。3dcgのスキニングでは、従来のリニアブレンディングスキンニング(LBS)に代わる方法として、クォータニオンを用いたデュアルクォータニオンスキニング(DQS)があります。LBSにはcandy wrapper artifactと呼ばれる、曲げたときに体積が失われる問題が存在します。これは、LBSでは回転行列と移動ベクトルを並べた剛体変換を重み付けし、その線形和を各々の頂点に適用するのですが、この変換が剛体変換であることは保証されていないために発生する問題です。DQSはcandy wrapper artifactを避けてくれます。ちなみに計算コストはLBSの1.5倍未満程度のようです。
※DQSにも、joint-bulging artefactと呼ばれる、曲げたときに関節の周りに膨らみが生じるという問題があります。これを避けるために、LBSとDQSを組み合わせる方法や、scale付きのデュアルクォータニオンを使う手法、元の距離を維持するために修正を加える手法があります。最近の手法では、Disney Researchの回転中心スキニングがあります。
デュアルクォータニオンの定義
元の状態を\(s_{0}\)、変化後の状態を\(s_{tr}\)とします。どちらもクォータニオンの二次元ベクトルです。回転を表すクォータニオンを\(r\)、移動を表すクォータニオンを\(t\)とすると、以下のように表すことができます。tとrの乗算する順が気になるところですが、回転と移動を作用させる順序に限らずこの形になります。これは嬉しいですね。
\begin{eqnarray}s_{tr}=\{\begin{pmatrix}1&0 \\0&1\end{pmatrix}r+\begin{pmatrix}0&0 \\1&0 \end{pmatrix}tr\}s_{0}\end{eqnarray}
ここで、\(\begin{pmatrix}0&0 \\1&0 \end{pmatrix}\)と同じ性質を持つようなε(\(ε^2=0\))を定義します。そうしてデュアルクォータニオンを以下のように表すことにしました。
\begin{eqnarray}d=r+εtr\end{eqnarray}
εの固有ベクトルを変更しても問題ないので、ε→\(\frac{ε}{2}\)として再定義します。\(d=r+\frac{ε}{2}tr\)としたとき、
\begin{eqnarray}|d|^2=(r+\frac{ε}{2}tr)(\bar{r}+\frac{ε}{2}\bar{t}\bar{r})=r\bar{r}+\frac{ε}{2}(tr\bar{r}+\bar{t}\bar{r}r) =1+(t+\bar{t})ε=1\end{eqnarray}
※回転クォータニオンrのノルムは1、tは移動量よりi、j、k成分しかないので\(\bar{t}=-t\)が成り立つ
デュアルクォータニオンとその共役の積が1となり、\(d\)は単位デュアルクォータニオンとなることが分かります。
参考
qiita.com
qiita.com
en.wikipedia.org
www.sky-engin.jp
techblog.sega.jp
en.wikipedia.org
kamino.hatenablog.com
www.f-sp.com
rodolphe-vaillant.fr