数学的に難しい話は別の方にお任せして，なんとな～くクォータニオンの片鱗に触れられたらいいなと思います．

📌 はじめに

物体を回転させるときに，よくXYZの各軸に対して回転値を設定すると思います．これは直感的にわかりやすい方法であり，このような表現を オイラー角表現 といいます．では実際にオイラー角で指定した回転を処理する場合，通常は行列やクォータニオンに変換して計算します．これは行列やクォータニオンが便利な性質を持っているからです．特にクォータニオンは回転を表すのに行列よりも少ない情報量と計算負荷が軽いといった特徴があります．また，クォータニオンを拡張して回転と平行移動を表すデュアルクォータニオンがあります．これはスキニングで使われています．ここでは，はじめにクォータニオン，その後にデュアルクォータニオンについて扱います．

📌 特性

クォータニオンの特性としてよく挙げられるのは次の２つです．

• ジンバルロックが発生しない
• 球面線形補間ができる

この辺りの詳細は割愛します．

📌 複素数

クォータニオンに入る前に，少し複素数について復習しておきます．2乗して-1になるものを 虚数 といい，iで表します．

i^2 = -1

複素数は実部 a と虚部 bi で構成された数です．

z = a+bi

これは，実部をX軸，虚部をY軸にとった平面で表現することができます．

これを複素平面またはガウス平面といいます．ここで，原点から z までの距離は

|z|=\sqrt{a^2+b^2}

で得られます．これを複素数の 絶対値 といいます．実際に複素数の絶対値を求める場合， a^2+b^2 を得るためには a+bi に虚数部の符号が反転した a-bi をかけます．これら複素数を 共役 といいます． a+bi の共役な複素数は a-bi であり， a-bi の共役な複素数は a+bi です．共役な複素数は \overline{z} または z* で表します．複素数 z=a+bi ，その共役な複素数 \overline{z}=a-bi を使うと

|z| = \sqrt{z\overline{z}} = \sqrt{(a+bi)(a-bi)} = \sqrt{a^2-b^2i^2} = \sqrt{a^2+b^2}

となります．

📌 クォータニオン

最初に複素数を考えると，これは実数の世界に虚数 i という概念を導入して構築されたものです．これと同様に虚数成分 i,j,k を導入して構築したものが クォータニオン(Quaternion) または ４元数 と呼ばれるものです．クォータニオン q は虚数成分を使って

q = w+xi+yj+zk

と表現します． i,j,k は虚数の性質

i^2 = j^2 = k^2 = -1

を持っています．また，

ijk = -1

という性質もあります．さらに

ij = -ji = k \qquad jk = -kj = i \qquad ki = -ik = j

という関係があります．これは例えば ijk=-1 の両辺に k を掛けると

ijk^2 = -k \qquad ij(-1) = -k \qquad -ij = -k \\[1ex] \therefore ij=k

となります．

加減乗算

まず加減算ですが，クォータニオンをベクトルとして考えるとわかりやすいです．ここで

q_1 = w_1+x_1i+y_1j+z_1k, \qquad q_2 = w_2+x_2i+y_2j+z_2k

とすると

\begin{aligned} q_1+q_2 &= (w_1+w_2)+(x_1+x_2)i+(y_1+y_2)j+(z_1+z_2)k \\[1ex] q_1-q_2 &= (w_1-w_2)+(x_1-x_2)i+(y_1-y_2)j+(z_1-z_2)k \end{aligned}

となります．次にスカラー倍は

sq = sw+sxi+syj+szk

となります．ベクトルと変わりませんね．次は乗算です．ベクトルではベクトル同士の積だと内積や外積がありましたが，クォータニオンでは素直に展開します．

\begin{aligned} q_1 q_2 &= (w_1+x_1i+y_1j+z_1k)(w_2+x_2i+y_2j+z_2k) \\[1ex] &= w_1 w_2 + w_1 x_2i + w_1 y_2j + w_1 z_2k \\[1ex] &\quad+ x_1i w_2 + x_1i x_2i + x_1i y_2j + x_1i z_2k \\[1ex] &\quad+ y_1i w_2 + y_1i x_2i + y_1i y_2j + y_1i z_2k \\[1ex] &\quad+ z_1i w_2 + z_1i x_2i + z_1i y_2j + z_1i z_2k \end{aligned}

i,j,k で囲むと

\begin{aligned} q_1 q_2 &= (w_1 w_2 - x_1 x_2 - y_1 y_2 - z_1 z_2) \\[1ex] &\quad+ (w_1 x_2 + w_2 x_1 + y_1 z_2 - z_1 y_2)i \\[1ex] &\quad+ (w_1 y_2 + w_2 y_1 + z_1 x_2 - x_1 z_2)j \\[1ex] &\quad+ (w_1 z_2 + w_2 z_1 + x_1 y_2 - y_1 x_2)k \end{aligned}

となります．

ベクトルを用いた表記

加減算とスカラー倍ではベクトルの演算と同じような方法でした．乗算については実はベクトルを使って表すことができます．そこで，クォータニオンをベクトルを用いて表現してみます． i,j,k それぞれに掛かる係数を要素とした3次ベクトル \vec{v} を使って

\vec{v} = <x,y,z>, \quad q = (w,\vec{v})

と表記します．このとき w の部分をスカラー部， \vec{v} の部分をベクトル部と呼ぶこともあります．これを用いると加減算は

\begin{aligned} q_1+q_2 &= (w_1+w_2,\vec{v_1}+\vec{v_2}) \\[1ex] q_1-q_2 &= (w_1-w_2,\vec{v_1}-\vec{v_2}) \end{aligned}

そして，乗算は

q_1 q_2 = (w_1 w_2 - \vec{v}_1\cdot\vec{v}_2, w_1 \vec{v}_2 + w_2 \vec{v}_1 + \vec{v}_1\times\vec{v}_2)

となって，スカラー部のところに内積，ベクトル部のところに外積を使ってすっきりと表記することができます．ここで w=0 の場合を考えると

q_1 q_2 = (- \vec{v}_1\cdot\vec{v}_2, \vec{v}_1\times\vec{v}_2)

となって，内積と外積だけになりました．クォータニオンとベクトル，内積と外積の関係が浮かび上がってきた感じがします．ちなみに， w が 0 のクォータニオンのことを 純クォータニオン または 純四元数 といいます．また，クォータニオンの乗算は交換法則が成り立ちません．これは外積が含まれていることを考えれば予想できますね．

逆クォータニオン

次に除算を考えます．行列には逆数と同じ振舞いをする逆行列があったように，クォータニオンには同じ振舞いをする 逆クォータニオン というのがあります．除算は逆クォータニオンを乗算で行います．それでは逆クォータニオンを求めるためにいくつか準備していきます．ここで複素数のときに触れたことを思い出してみましょう．まずは 絶対値 です．クォータニオンの絶対値は複素数のように

|q| = \sqrt{w^2+x^2+y^2+z^2}

で得られます．絶対値いわば長さが求められるということは，正規化することができます．スカラー倍が可能なら，スカラー除算も可能なので

normalize(q) = \frac{q}{|q|}

とすることで長さが1のクォータニオンを作ることができます．これを 単位クォータニオン といいます．
次に，複素数には共役があったように，クォータニオンにも共役が存在します．クォータニオンの共役はベクトル部にマイナスを掛けたクォータニオンのことです．

\overline{q} = (w, -\vec{v})

この共役の性質を見ていきます．まず，クォータニオンとその共役を掛け合わせると

\begin{aligned} q\overline{q} &= (w,\vec{v})(w,-\vec{v}) \\[1ex] &= (w^2+\vec{v}\cdot\vec{v},-w\vec{v}+w\vec{v}+\vec{v}\times(-\vec{v})) \\[1ex] &= w^2 + \vec{v}\cdot\vec{v} \\[1ex] &= |q|^2 \end{aligned}

となって，スカラーになりました．しかも，絶対値（長さ）の2乗です．この式を変形すると

\frac{q\overline{q}}{|q|^2} = 1

となります．同じものを割っているので1ですね．これは

q\frac{\overline{q}}{|q|^2} = 1

と見ることができます． q にあるものをかけると 1 になる…それはまさに逆クォータニオンのことです．よって

q^{-1} = \frac{\overline{q}}{|q|^2}

となります．また，クォータニオンが正規化されていた場合，絶対値は1なので，その逆クォータニオンは共役そのものになります．

📌 回転軸表現

クォータニオンによる回転について触れていく前に，回転軸表現(Axis-angle representation)について見ていきます．オイラー角表現ではXYZ軸における回転角度によって表現していましたが，回転軸表現では回転軸となるベクトルと，回転角度によって表現します．このような回転を 任意軸回転 と呼ぶこともあります．この回転を下図を見ながら考えてみます．

回転軸となる正規化ベクトルを \vec{n} として角度 \theta 回転を行うとします．回転の対象となる点を指すベクトル \vec{r} ，回転後のベクトルを \vec{r}' とします．次にベクトル \vec{r} を２つのベクトルに分解します．回転軸 \vec{n} に沿った平行成分ベクトル \vec{r}\_{\parallel} と， \vec{n} に対して垂直な回転面に沿った垂直成分ベクトル \vec{r}\_{\perp} に分解すると，

\begin{aligned} \vec{r}_{\parallel} &= (\vec{n}\cdot\vec{r})\vec{n} \\[1ex] \vec{r}_{\perp} &= \vec{r}-(\vec{n}\cdot\vec{r})\vec{n} \end{aligned}

となります．ベクトル \vec{r} を回転しても，平行成分ベクトル \vec{r}\_{\parallel} は変化しませんので，回転面のみの計算になります．右図は回転面を上から見たものです．ここで，回転軸 \vec{n} と垂直成分ベクトル \vec{r}\_{\perp} に直交していて回転面上にあるベクトル \vec{v} を導入します．このようなベクトルは外積で求まります．

\vec{v} = \vec{n}\times\vec{r}_{\perp} = \vec{n}\times\vec{r}

このベクトル \vec{v} を使うと回転後のベクトルの垂直成分 \vec{r}_{\perp}' は次のようになります．

\vec{r}_{\perp}' = (\cos\theta)\vec{r}_{\perp} + (\sin\theta)\vec{v}

あとは，垂直成分ベクトル \vec{r}_{\parallel} と合成して，展開すると

\begin{aligned} \vec{r}' &= \vec{r}_{\parallel} + \vec{r}_{\perp}' \\[1ex] &= \vec{r}_{\parallel} + (\cos\theta)\vec{r}_{\perp} + (\sin\theta)\vec{v} \\[1ex] &= (\vec{n}\cdot\vec{r})\vec{n} + (\cos\theta)\left\{\vec{r}-(\vec{n}\cdot\vec{r})\vec{n}\right\}+(\sin\theta)(\vec{n}\times\vec{r}) \\[1ex] &= (\cos\theta)\vec{r} + (\vec{n}\cdot\vec{r})\vec{n} - (\cos\theta)(\vec{n}\cdot\vec{r})\vec{n}+(\sin\theta)(\vec{n}\times\vec{r}) \\[1ex] &= (\cos\theta)\vec{r} + (1-\cos\theta)(\vec{n}\cdot\vec{r})\vec{n}+(\sin\theta)(\vec{n}\times\vec{r}) \qquad \qquad \text{\color{blue}(1)} \end{aligned}

となります．これは回転軸 \vec{n} と回転角度 \theta によるベクトルの回転を表しています．

クォータニオンでの回転

それではクォータニオンによる回転を見ていきます．まず，回転の対象となるベクトル \vec{r} をクォータニオン p で表すと

p = (0,\vec{r})

となります．そして，ある正規化されたクォータニオン q を使って，次のように考えます．

p' = qpq^{-1}

この式が一体なにをしているのか見ていきます．これを展開してみると

\begin{aligned} qpq^{-1} &= (w,\vec{v})(0,\vec{r})(w,-\vec{v}) \\[1ex] &= (-\vec{v}\cdot\vec{r},w\vec{r}+\vec{v}\times\vec{r})(w,-\vec{v}) \\[1ex] &= ((-\vec{v}\cdot\vec{r})w-(w\vec{r}+\vec{v}\times\vec{r})\cdot(-\vec{v}), \quad (-\vec{v}\cdot\vec{r})(-\vec{v})+ w(w\vec{r}+\vec{v}\times\vec{r})+ (w\vec{r}+\vec{v}\times\vec{r})\times(-\vec{v})) \qquad \qquad \text{\color{blue}(2)} \end{aligned}

となります．この式を変形するために内積や外積の性質，スカラー3重積，ベクトル3重積を使います．これまで扱っていない用語もありますがとりあえずそういうものだと思ってください．説明のためにベクトルはアルファベットの大文字 A,B,C で表し，スカラーは小文字 k を使います．

---

内積の性質

A\cdot B = B\cdot A \\[1ex] A\cdot(B+C) = A\cdot B + A\cdot C, \qquad (A+B)\cdot C = A\cdot C + B\cdot C \\[1ex] (kA)\cdot B = A\cdot (kB) = k(A\cdot B)

---

外積の性質

A\times B = -B\times A \\[1ex] A\times(B+C) = A\times B + A\times C, \qquad (A+B)\times C = A\times C + B\times C \\[1ex] (kA)\times B = A\times (kB) = k(A\times B)

---

スカラー3重積

A\cdot(B\times C) = B\cdot(C\times A) = C\cdot(A\times B)

---

ベクトル3重積

\begin{aligned} A\times(B\times C) &= (A\cdot C)B - (A\cdot B)C \\[1ex] (B\times C)\times A &= -(A\cdot C)B + (A\cdot B)C \\[1ex] (A\times B)\times C &= -(C\cdot B)A + (C\cdot A)B = (A\cdot C)B - (B\cdot C)A \end{aligned}

---

式(2)のスカラー部を見てみると

\begin{aligned} (-\vec{v}\cdot\vec{r})w-(w\vec{r}+\vec{v}\times\vec{r})\cdot(-\vec{v}) &= (-\vec{v}\cdot\vec{r})w-\left\{(w\vec{r})\cdot(-\vec{v})+(\vec{v}\times\vec{r})\cdot(-\vec{v})\right\} \\[1ex] &= (-\vec{v}\cdot\vec{r})w-(-\vec{v}\cdot\vec{r})w-(\vec{v}\times\vec{r})\cdot(-\vec{v}) \\[1ex] &= (\vec{v}\times\vec{r})\cdot(\vec{v}) \\[1ex] &= (\vec{v}\times\vec{v})\cdot\vec{r} \\[1ex] &= 0 \end{aligned}

となります．最後の \vec{v}\times\vec{v} は同じベクトルの外積なので 0 になりますね．
次にベクトル部です．

\begin{aligned} (-\vec{v}\cdot\vec{r})&(-\vec{v})+ w(w\vec{r}+\vec{v}\times\vec{r})+ (w\vec{r}+\vec{v}\times\vec{r})\times(-\vec{v}) \\[1ex] &= (\vec{v}\cdot\vec{r})\vec{v}+ w^2\vec{r}+ w(\vec{v}\times\vec{r})+ w(\vec{r}\times(-\vec{v}))+ (\vec{v}\times\vec{r})\times(-\vec{v}) \\[1ex] &= (\vec{v}\cdot\vec{r})\vec{v}+ w^2\vec{r}+ 2w(\vec{v}\times\vec{r})+ (\vec{v}\times\vec{r})\times(-\vec{v}) \\[1ex] &= (\vec{v}\cdot\vec{r})\vec{v}+ w^2\vec{r}+ 2w(\vec{v}\times\vec{r})+ (\vec{v}\cdot(-\vec{v}))\vec{r}-(\vec{r}\cdot(-\vec{v}))\vec{v} \\[1ex] &= (\vec{v}\cdot\vec{r})\vec{v}+ w^2\vec{r}+ 2w(\vec{v}\times\vec{r})+ -(\vec{v}\cdot\vec{v})\vec{r}+(\vec{r}\cdot\vec{v})\vec{v} \\[1ex] &= w^2\vec{r}-(\vec{v}\cdot\vec{v})\vec{r}+ 2(\vec{v}\cdot\vec{r})\vec{v}+ 2w(\vec{v}\times\vec{r}) \\[1ex] &= (w^2-\vec{v}\cdot\vec{v})\vec{r}+ 2(\vec{v}\cdot\vec{r})\vec{v}+ 2w(\vec{v}\times\vec{r}) \end{aligned}

よって

\begin{aligned} qpq^{-1} &= (w,\vec{v})(0,\vec{r})(w,-\vec{v}) \\[1ex] &= (0, \quad (w^2-\vec{v}\cdot\vec{v})\vec{r} + 2\vec{v}(\vec{v}\cdot\vec{r})+2w(\vec{v}\times\vec{r})) \end{aligned}

となります．ここで， q を

q = (\cos\phi, \sin\phi \vec{n})

とすると（ \vec{n} は正規化ベクトル）

\begin{aligned} p' &= qpq^{-1} \\[1ex] &= (w,\vec{v})(0,\vec{r})(w,-\vec{v}) \\[1ex] &= (0, \quad (w^2-\vec{v}\cdot\vec{v})\vec{r} + 2\vec{v}(\vec{v}\cdot\vec{r})+2w(\vec{v}\times\vec{r})) \\[1ex] &= (0, \quad (\cos^2\phi-\sin^2\phi(\vec{n}\cdot\vec{n}))\vec{r}+ 2\sin^2\phi(\vec{n}\cdot\vec{r})\vec{n}+ 2\cos\phi\sin\phi(\vec{n}\times\vec{r})) \\[1ex] &= (0, \quad \cos2\phi\vec{r}+ (1-\cos2\phi)(\vec{n}\cdot\vec{r})\vec{n}+ \sin2\phi(\vec{n}\times\vec{r})) \qquad \qquad \text{\color{blue}(3)} \end{aligned}

になります．ここでは，以下の公式（2倍角の公式など）を使いました．

\sin 2A = 2\sin A\cos A, \quad \cos 2A = \cos^2 A-\sin^2 A, \quad 1-\cos 2A = 2\sin^2

この式(3)と任意軸回転における式(1)

\vec{r}' = (\cos\theta)\vec{r} + (1-\cos\theta)(\vec{n}\cdot\vec{r})\vec{n}+(\sin\theta)(\vec{n}\times\vec{r})

を比べてみましょう．ほとんど同じですね．気になるところは純クォータニオンになっていること，そして，回転角度が 2\phi になっていることです．これを考慮して

q = \left(\cos\frac{\theta}{2}, \sin\frac{\theta}{2}\vec{n}\right)

とします．このクォータニオンを使って qpq^{-1} の演算を行うとベクトルを回転させることができるわけです．

📌 行列表記

クォータニオンによるベクトルの回転を見てきましたが，ベクトルの回転は行列を使って行うこともできます．行列では回転だけでなく平行移動や拡大縮小といった変換を適用することができます．場合によっては回転を行列にして扱い時もあるでしょう．そこでクォータニオンから行列に変換します．まずは，クォータニオン同士の乗算を考えてみます．クォータニオン同士の乗算は

\begin{aligned} q_1 q_2 &= (w_1 w_2 - x_1 x_2 - y_1 y_2 - z_1 z_2) \\[1ex] &\quad+ (w_1 x_2 + w_2 x_1 + y_1 z_2 - z_1 y_2)i \\[1ex] &\quad+ (w_1 y_2 + w_2 y_1 + z_1 x_2 - x_1 z_2)j \\[1ex] &\quad+ (w_1 z_2 + w_2 z_1 + x_1 y_2 - y_1 x_2)k \end{aligned}

でした．クォータニオンの4つの要素を i 軸, j 軸, k 軸，実軸の順番に並べたベクトルとして行列を作ると

q_1 q_2 = \begin{pmatrix}w_1 x_2 + w_2 x_1 + y_1 z_2 - z_1 y_2 \\ w_1 y_2 + w_2 y_1 + z_1 x_2 - x_1 z_2 \\ w_1 z_2 + w_2 z_1 + x_1 y_2 - y_1 x_2 \\ w_1 w_2 - x_1 x_2 - y_1 y_2 - z_1 z_2 \end{pmatrix}

となります．ここで，この乗算を q_1 の行列 L と q_2 のベクトルの積で表すと

q_1 q_2 = \begin{pmatrix}w_1 x_2 - z_1 y_2 + y_1 z_2 + x_1 w_2 \\ z_1 x_2 + w_1 y_2 - x_1 z_2 + y_1 w_2 \\ - y_1 x_2 + x_1 y_2 + w_1 z_2 + z_1 w_2 \\ - x_1 x_2 - y_1 y_2 - z_1 z_2 + w_1 w_2\end{pmatrix} = L(q_1)\begin{pmatrix}x_2 \\ y_2 \\ z_2 \\ w_2\end{pmatrix}

となって，行列 L は

L(q) = \begin{pmatrix}w & -z & y & x \\ z & w & -x & y \\ -y & x & w & z \\ -x & -y & -z & w\end{pmatrix}

となります．同じように今度は右側から掛けたことと同じになる行列 R を求めると

q_1 q_2 = \begin{pmatrix}w_2 x_1 + z_2 y_1 - y_2 z_1 + x_2 w_1 \\ -z_2 x_1 + w_2 y_1 + x_2 z_1 + y_2 w_1 \\ y_2 x_1 - x_2 y_1 + w_2 z_1 + z_2 w_1 \\ -x_2 x_1 - y_2 y_1 - z_2 z_1 + w_2 w_1\end{pmatrix} = R(q_2)\begin{pmatrix}x_1 \\ y_1 \\ z_1 \\ w_1\end{pmatrix} = \begin{pmatrix}w & z & -y & x \\ -z & w & x & y \\ y & -x & w & z \\ -x & -y & -z & w\end{pmatrix} \begin{pmatrix}x_1 \\ y_1 \\ z_1 \\ w_1\end{pmatrix}

となります．これでクォータニオンを左からベクトルに乗算する行列 L と，右からベクトルに乗算する行列 R が得られました．この2つの行列を使ってベクトル p を回転させる場合

\begin{aligned} qpq^{-1} &= qp\overline{q}\\[1ex] & = L(q)R(\overline{q})p = \begin{pmatrix}w^2+x^2-y^2-z^2 & 2(xy-wz) & 2(xz+wy) & 0 \\ 2(xy+wz) & w^2-x^2+y^2-z^2 & 2(yz-wx) & 0 \\ 2(xz-wy) & 2(yz+wx) & w^2-x^2-y^2+z^2 & 0 \\ 0 & 0 & 0 & 1\end{pmatrix}\begin{pmatrix}x_r \\ y_r \\ z_r \\ 0\end{pmatrix} \end{aligned}

と計算することができます．この行列がクォータニオン q の回転行列です．

📌 回転の合成

クォータニオンは行列と似たように，回転の合成を行うことができます．最初に q_1 回転して，続けて q_2 回転を行う場合

q_2\left(q_1 p q_1^{-1}\right)q_2^{-1} = \left(q_2 q_1\right)p\left(q_1^{-1} q_2^{-1}\right) = \left(q_2 q_1\right)p\left(q_2 q_1\right)^{-1}

となって， q_2 q_1 というクォータニオンで一度に回転を行うことができます．例えば，次のようなXYZ軸回転のクォータニオンを作成することができます．

\begin{aligned} q_x &= \left(\cos\frac{\phi}{2},<\sin\frac{\phi}{2},0,0>\right) \\[1ex] q_y &= \left(\cos\frac{\phi}{2},<0,\sin\frac{\phi}{2},0>\right) \\[1ex] q_z &= \left(\cos\frac{\phi}{2},<0,0,\sin\frac{\phi}{2}>\right) \\[1ex] q_{rot} &= q_z q_y q_x \end{aligned}

📌 球面線形補間

クォータニオンの良いところにこの 球面線形補間(Slerp) が使えることがあります．次の図を見てください．

2次元ベクトル p_1 と p_2 があり，それぞれ正規化されたベクトルとします． t は 0 から 1 の任意の数です．このとき，中間ベクトル p_t は半径 1 の円弧上を移動します．単純に線形補間だと円弧上を線形に移動しません．ですが，球面線形補間を使うと正しく線形補間することができます．それでは球面線形補間について幾何学的に考えてみます．下図を見てください．

ベクトル A からベクトル B への補間を考えます．2つとも正規化されているとします． t は 0 から 1 までの値を取るパラメータで， t のときの補間ベクトルを P とします．この補間ベクトルを求めるには次の図のように考えることができます．

\vec{v}_x は \overline{OA} に対して平行で， \vec{v}_y は \overline{OB} に対して平行です．この2つのベクトル \vec{v}_x と \vec{v}_y を求めて足し合わせれば補間ベクトル P になります．まずは \vec{v}_y を求めてみます．次の図を見てください．

B から \overline{OA} に対して垂直に下ろした位置を C ， P から \overline{OA} に対して垂直に下ろした位置を E とします． \overline{OB} の長さは 1 ですから，三角形COBにおいて \overline{BC} の長さは \sin\theta です．そして，三角形EOPでは \overline{OP} の長さが 1 なので， \overline{PE} の長さが \sin(t\theta) となります．青い三角形は相似の関係になっているので， \overline{OB}:\overline{BC} と \overline{CP}:\overline{PE} の比は同じです．よって， y の長さは y:1=\sin(t\theta):\sin\theta から

y = \frac{\sin(t\theta)}{\sin\theta} \qquad \therefore \vec{v}_y = \frac{\sin(t\theta)}{\sin\theta}\vec{B}

となります．次に \vec{v}_x を求めます．次の図を見てください．

A から \overline{OB} に対して垂直に下ろした位置を F ， P から \overline{OB} に対して垂直に下ろした位置を H ， P から \overline{OB} に向かって， \overline{OA} に平行な線を引き， \overline{OB} との交点を G とします． \overline{OA} の長さは 1 ですから，三角形AOFにおいて \overline{AF} の長さは \sin\theta です．そして，三角形HOPでは \overline{OP} の長さが 1 なので， \overline{PH} の長さが \sin(1-t)\theta となります．青い三角形は相似の関係になっているので， \overline{OA}:\overline{AF} と \overline{GP}:\overline{PH} の比は同じです．よって， x の長さは x:1=\sin((1-t)\theta):\sin\theta から

x = \frac{\sin((1-t)\theta)}{\sin\theta} \qquad \therefore \vec{v}_x = \frac{\sin((1-t)\theta)}{\sin\theta}\vec{A}

となります．よって

において，補間ベクトル \vec{p}_t は

\vec{p}_t = \vec{p}_1 \frac{\sin((1-t)\theta)}{\sin\theta} + \vec{p}_2 \frac{\sin(t\theta)}{\sin\theta}

となります．実際には2つのクォータニオン q_1, q_2 を指定して補間ベクトルを求めるといった使い方をするので，その場合の \theta は内積で求められます．

\theta = \cos^{-1}(q_1 \cdot q_2)

また， q_1 \cdot q_2 < 0 の場合，外側に向かって補間されてしまいます．クォータニオン q と -q は同じ回転を表すことを利用して， q_1 \cdot q_2 \geq 0 を満たすようにすれば，補間ベクトルは最短の経路で行われます．

📌 デュアルクォータニオン

クォータニオンは複素数を拡張して回転を表したものでした．このクォータニオンをさらに拡張して，回転と平行移動を表すデュアルクォータニオンというものがあります．このデュアルクォータニオンは行列と比べてパラメータが少ないことや回転と平行移動を回転軸に沿って補間する ScLerp(Screw Linear Interpolation) が使えます．このScLerpはスキニングアニメーションで使われることがあり，例えばねじるような回転をしたときに関節の部分が縮まってしまって，その部分の体積が無くなってしまう問題に効果的です．ただし，計算コストが若干増えるのと，デュアルクォータニオンによって別の問題が発生してしまい，それを解決するためにまた計算コストが高くなるという問題や，回転と移動だけなので，スケール情報がないといった問題もありますので，中々使いどころが難しい気がします．スキニングアニメーション以外ですと，スケールを考慮しない剛体の並進・回転計算に使われることもあるようです．

Mayaでのデュアルクォータニオン設定

二元数

クォータニオンは4つの元（実数，虚数i,j,k）つまり四元数です．これは複素数の拡張と考えることができます．では，複素数はというと実数と虚数の２つの元で表します．このように2つの元で表す数を二元数といいます．実数と次のような新しい元

\varepsilon^2 = 0

で構成された二元数があり，

z = a+\varepsilon b \qquad \varepsilon \neq 0

と表します．これを 二重数 といいます．複素数のように a を実部(Real Part)， \varepsilon b を虚部(Dual Part/Non Real Part)， \varepsilon は虚数単位です．二重数は複素数と同じように計算します．

\begin{aligned} z_1 + z_2 &= (a_1 + a_1) + \varepsilon(b_1+b_2) \\[2ex] z_1 z_2 &= a_1 a_2 + \varepsilon a_1 b_2 + \varepsilon a_2 b_1 + \varepsilon^2 a_2 b_2 \\[2ex] &= a_1 a_2 + \varepsilon(a_1 b_2 + a_2 b_1) \end{aligned}

デュアルクォータニオン

デュアルクォータニオンを二元数で次のように表します．

\sigma = p + \varepsilon q

ここで p と q はクォータニオンで

\begin{aligned} p &= p_w + \vec{p} = p_w + p_x i + p_y j + p_z k \\[2ex] q &= q_w + \vec{q} = q_w + q_x i + q_y j + q_z k \end{aligned}

となっています．加減乗算は次のようになります．

\begin{aligned} \sigma_1 + \sigma_2 &= (p_1 + p_2) + \varepsilon(q_1 + q_2) \\[2ex] \sigma_1 \sigma_2 &= p_1 p_2 + \varepsilon(p_1 q_2 + p_2 q_1) \end{aligned}

また，デュアルクォータニオンの実部と虚部のスカラー成分の和を

\hat{d} = p_w + \varepsilon q_w

ベクトル成分の和を

\vec{d} = \vec{p_v} + \varepsilon \vec{q_v}

とおくと，次のようにクォータニオンと同じ形に表せます．

\sigma = (\hat{d}, \ \vec{d}).

これを使うと乗算はクォータニオンと同じように計算できます．

\begin{aligned} \sigma_1 \sigma_2 &= (\hat{d_1}, \ \vec{d_1})(\hat{d_2}, \ \vec{d_2}) \\[2ex] &= (\hat{d_1}\hat{d_2} - \vec{d_1}\cdot\vec{d_2}, \ \hat{d_1}\vec{d_2} + \hat{d_2}\vec{d_1}+\vec{d_1}\times\vec{d_2}) \end{aligned}

共役

デュアルクォータニオンでは二重数の共役，クォータニオンの共役，その２つを合わせた共役の３種類あります．まず，二重数の共役は複素数と同じように虚部の符号が反転したものです

\sigma^{\bullet} = p - \epsilon q.

デュアルクォータニオンとその二重数の共役を乗算すると

\begin{aligned} \sigma \sigma^{\bullet} &= (p+\varepsilon q)(p-\varepsilon q) \\[2ex] &= pp + \varepsilon(qp - pq). \end{aligned}

となります．次に，クォータニオンの共役の場合は

\sigma^{\ast} = p^{\ast} + \epsilon q^{\ast}.

となって，乗算すると

\begin{aligned} \sigma \sigma^{\ast} &= (p+\varepsilon q)(p^{\ast} + \varepsilon q^{\ast}) \\[2ex] &= pp^{\ast} + \varepsilon (pq^{\ast} + qp^{\ast}) \\[2ex] &= (p_w^2+p_x^2+p_y^2+p_z^2) + 2\varepsilon(p_w q_w + p_x q_x + p_y q_y + p_z q_z). \\[2ex] &= ||p||^2 + 2\varepsilon(p_w q_w + \vec{p}_v \vec{q}_v). \end{aligned}

となります．三つ目は二重数とクォータニオンの共役を合わせたもので

\sigma^{\diamond} = p^{\ast} - \varepsilon q^{\ast}.

この乗算は

\begin{aligned} \sigma \sigma^{\diamond} &= (p + \varepsilon q)(p^{\ast} - \varepsilon q^{\ast}) \\[2ex] &= pp^{\ast} + \varepsilon(qp^{\ast}-pq^{\ast}). \end{aligned}

となります．
次に，デュアルクォータニオン同士を乗算したものの各共役を求めます．その前にクォータニオンの次の性質を確認します．

\begin{aligned} (pq)^{\ast} &= \left\{(p_w,\vec{p}_v)(q_w,\vec{q}_v)\right\}^{\ast} \\[2ex] &= (p_w q_w - \vec{p}_v\cdot\vec{q}_v,\ p_w \vec{q}_v + q_w \vec{p}_v + \vec{p}_v \times \vec{q}_v)^{\ast} \\[2ex] &= (p_w q_w - \vec{p}_v\cdot\vec{q}_v,\ - p_w \vec{q}_v - q_w \vec{p}_v - \vec{p}_v \times \vec{q}_v) \\[2ex] q^{\ast}p^{\ast} &= (q_w,\vec{q}_v)^{\ast}(p_w,\vec{p}_v)^{\ast} \\[2ex] &= (q_w,-\vec{q}_v)(p_w,-\vec{p}_v) \\[2ex] &= (q_w p_w - (-\vec{q}_v)\cdot(-\vec{p}_v),\ q_w (-\vec{p}_v) + p_w (-\vec{q}_v) + (-\vec{q}_v) \times (-\vec{p}_v) \\[2ex] &= (q_w p_w - \vec{q}_v\cdot\vec{p}_v,\ - q_w \vec{p}_v - p_w \vec{q}_v + \vec{q}_v \times \vec{p}_v) \end{aligned}

よって

\therefore (pq)^{\ast} = q^{\ast}p^{\ast}

となります．この関係を利用して

\begin{aligned} (\sigma_1 \sigma_2)^{\bullet} &= \left\{p_1 p_2 + \varepsilon(p_1 q_2 + q_1 p_2)\right\}^{\bullet} \\[2ex] &= p_1 p_2 + \varepsilon(p_1 q_2 + q_1 p_2) \\[2ex] &= (p_1 - \varepsilon q_1)(p_2 - \varepsilon q_2) \\[2ex] &= \sigma_1^{\bullet} \sigma_2^{\bullet} \end{aligned}

\begin{aligned} (\sigma_1 \sigma_2)^{\ast} &= \left\{p_1 p_2 + \varepsilon(p_1 q_2 + q_1 p_2)\right\}^{\ast} \\[2ex] &= (p_1 p_2)^{\ast} + \varepsilon(p_1 q_2 + q_1 p_2)^{\ast} \\[2ex] &= (p_1 p_2)^{\ast} + \varepsilon\left\{(p_1 q_2)^{\ast} + (q_1 p_2)^{\ast}\right\} \\[2ex] &= p_2^{\ast} p_1^{\ast} + \varepsilon(q_2^{\ast} p_1^{\ast} + p_2^{\ast} q_1^{\ast}) \\[2ex] &= (p_2^{\ast} + \varepsilon q_2^{\ast})(p_1^{\ast} + \varepsilon q_1^{\ast}) \\[2ex] &= \sigma_2^{\ast} \sigma_1^{\ast} \end{aligned}

\begin{aligned} (\sigma_1 \sigma_2)^{\diamond} &= \left\{p_1 p_2 + \varepsilon(p_1 q_2 + q_1 p_2)\right\}^{\diamond} \\[2ex] &= (p_1 p_2)^{\ast} - \varepsilon(p_1 q_2 + q_1 p_2)^{\ast} \\[2ex] &= (p_1 p_2)^{\ast} - \varepsilon\left\{(p_1 q_2)^{\ast} + (q_1 p_2)^{\ast}\right\} \\[2ex] &= p_2^{\ast} p_1^{\ast} - \varepsilon(q_2^{\ast} p_1^{\ast} + p_2^{\ast} q_1^{\ast}) \\[2ex] &= (p_2^{\ast} - \varepsilon q_2^{\ast})(p_1^{\ast} - \varepsilon q_1^{\ast}) \\[2ex] &= \sigma_2^{\diamond} \sigma_1^{\diamond} \end{aligned}

となります．まとめると次のようになっています．

(\sigma_1 \sigma_2)^{\bullet} = \sigma_1^{\bullet} \sigma_2^{\bullet}, \quad (\sigma_1 \sigma_2)^{\ast} = \sigma_2^{\ast} \sigma_1^{\ast}, \quad (\sigma_1 \sigma_2)^{\diamond} = \sigma_2^{\diamond} \sigma_1^{\diamond}

単位デュアルクォータニオン

クォータニオンの場合と同様に，デュアルクォータニオンとその共役を掛けたら 1 になるとき

\sigma \sigma^{\ast} = 1

このデュアルクォータニオンを 単位デュアルクォータニオン といいます．ノルムが 1 なので逆デュアルクォータニオンはその共役と同じになるからです．
単位デュアルクォータニオンは

p_w^2 + p_x^2 + p_y^2 + p_z^2 = 1, \qquad p_w q_w + p_x q_x + p_y q_y + p_z q_z = 0

を満たします．これは実部の大きさが 1 で，虚部が 0 であることを意味しています．

回転と平行移動

それではデュアルクォータニオンを使って回転と平行移動を表してみます．まず，回転軸の単位ベクトル \vec{n} ，回転角度 \theta のクォータニオン r は

r = \cos\frac{\theta}{2} + \vec{n}\frac{\theta}{2}

です．このクォータニオンの共役は

r^{\ast} = \cos\frac{\theta}{2} - \vec{n}\frac{\theta}{2}

で， rr^{\ast} = r^{\ast}r = 1 の関係になります．このクォータニオンの回転行列を R とします．次に，平行移動 \vec{t} は

\vec{t} = (t_1,t_2,t_3)

とします．あるベクトル \vec{v} を回転行列 R と平行移動 \vec{t} で座標変換すると

\vec{v}' = R\vec{v} + \vec{t}

と表せます．回転のクォータニオン r ，平行移動 \vec{t} の純クォータニオンからデュアルクォータニオン \sigma は

\begin{aligned} \sigma &= r + \frac{\varepsilon}{2}\vec{t}r \\[2ex] &= \cos\frac{\theta}{2} + \vec{n}\sin\frac{\theta}{2} + \frac{\varepsilon}{2}\left( -\sin\frac{\theta}{2}(\vec{t}\cdot\vec{n}) + \cos\frac{\theta}{2}\vec{t} + \sin\frac{\theta}{2}\vec{t}\times\vec{n} \right) \end{aligned}

となります．ここで \vec{t}=(0,0,0) なら

\sigma = r

また， \theta=0 なら

\sigma = 1 + \frac{\varepsilon}{2}\vec{t}

となります．回転と平行移動を表したクォータニオンとその共役を掛けてみると

\begin{aligned} \sigma \sigma^{\ast} &= \left(r + \frac{\varepsilon}{2}tr\right)\left(r^{\ast}+\frac{\varepsilon}{2}(tr)^{\ast}\right) \\[2ex] &= \left(r+\frac{\varepsilon}{2}tr\right)\left(r^{\ast}+\frac{\varepsilon}{2}r^{\ast}t^{\ast}\right) \\[2ex] &= rr^{\ast} + \frac{\varepsilon}{2}(rr^{\ast}t + trr^{\ast}) \\[2ex] &= 1 + \frac{\varepsilon}{2}(t^{\ast}+t) \\[2ex] &= 1 \end{aligned}

となって，単位デュアルクォータニオンということがわかります．途中式の t^{\ast}+t のところでは， t が純クォータニオンですから，その共役は -t となって結局その項は 0 になります．ここで単位デュアルクォータニオン \sigma から回転・平行移動を取り出すには

\begin{aligned} \sigma &= p + \varepsilon q \\[2ex] r &= p \qquad \qquad \text{\color{blue}(23)} \\[2ex] t &= 2qp^{\ast} \qquad \qquad \text{\color{blue}(24)} \end{aligned}

となります．式(24)は， q=\frac{1}{2}tr に式(23)を代入して変形したものです．この時， r で割る必要がありますので，その共役 r^{\ast} を掛けています．これまでで回転と平行移動を表すデュアルクォータニオンが得られました．次にこのデュアルクォータニオンを使って座標変換してみます．まず，座標変換したいベクトル \vec{v} を純クォータニオン (0,\vec{v}) とし，平行移動だけのデュアルクォータニオン \sigma = 1+\varepsilon \vec{v} とします．そして，クォータニオンで回転させるときのように，変換したいベクトルにデュアルクォータニオンとその共役を掛けます．

\begin{aligned} \sigma (1+\varepsilon \vec{v}) \sigma^{\diamond} &= \left(r+\frac{\varepsilon}{2}tr\right)(1+\varepsilon \vec{v})\left(r^{\ast}-\frac{\varepsilon}{2}(tr)^{\ast}\right) \\[2ex] &= \left(r + \frac{\varepsilon}{2}tr + \varepsilon r\vec{v}\right)\left(r^{\ast}-\frac{\varepsilon}{2}r^{\ast}t^{\ast}\right) \\[2ex] &= rr^{\ast} + \varepsilon\left(\frac{1}{2}(t-t^{\ast})+r\vec{v}r^{\ast}\right) \\[2ex] &= 1 + \varepsilon(r\vec{v}r^{\ast} + t) \qquad \qquad \text{\color{blue}(25)} \end{aligned}

式(25)にある r\vec{v}r^{\ast} はクォータニオンによる回転そのものです．それに平行移動 t を足していて，結局この式は平行移動のみのデュアルクォータニオンの形になっています．この演算は R\vec{v}+t と同じであることがわかります．

ついでに，回転する前に平行移動する場合のデュアルクォータニオン \xi を考えると

\xi = r + \frac{\varepsilon}{2}rt.

となって，次のような演算になります．

\begin{aligned} \xi (1+\varepsilon \vec{v}) \xi^{\diamond} &= \left(r+\frac{\varepsilon}{2}rt\right)(1+\varepsilon \vec{v})\left(r^{\ast}-\frac{\varepsilon}{2}(rt)^{\ast}\right) \\[2ex] &= \left(r + \frac{\varepsilon}{2}rt + \varepsilon r\vec{v}\right)\left(r^{\ast}-\frac{\varepsilon}{2}t^{\ast}r^{\ast}\right) \\[2ex] &= rr^{\ast} + \varepsilon\left(\frac{1}{2}(rtr^{\ast}-rt^{\ast}r^{\ast})+r\vec{v}r^{\ast}\right) \\[2ex] &= 1 + \varepsilon(rtr^{\ast} + r\vec{v}r^{\ast}) \\[2ex] &= 1 + \varepsilon r(\vec{v}+t)r^{\ast} \qquad \qquad \text{\color{blue}(27)} \end{aligned}

式(27)では r(\vec{v}+t)r^{\ast} となって，回転をする前に平行移動する形になっているのがわかります．最後に回転と平行移動を表すデュアルクォータニオン

\sigma = r + \frac{\varepsilon}{2}tr

の右辺を r でくくると

\sigma = \left(1+\frac{\varepsilon}{2}t\right)r

となって，回転のクォータニオン r と，平行移動 1+\frac{\varepsilon}{2}t で構成されていることがわかりやすくなります．

行列表記

デュアルクォータニオンから行列 M に変換する式です．

\sigma = p + \varepsilon q

M = \begin{pmatrix} 1 - 2p_y^2 - 2p_z^2 & 2 p_x p_y - 2 p_w p_z & 2 p_x p_z + 2 p_w p_y & t_x \\[2ex] 2 p_x p_y + 2 p_w p_z & 1 - 2 p_x^2 - 2p_z^2 & 2 p_y p_z - 2 p_w p_x & t_y \\[2ex] 2 p_x p_z - 2 p_w p_y & 2 p_y p_z + 2 p_w p_x & 1 - 2 p_x^2 - 2 p_y ^ 2 & t_z \\[2ex] 0 & 0 & 0 & 1 \end{pmatrix}

\begin{aligned} t_x &= 2(-q_w p_x + q_x p_w - q_y p_z + q_z p_y) \\[2ex] t_y &= 2(-q_w p_y + q_x p_z + q_y p_w - q_z p_x) \\[2ex] t_z &= 2(-q_w p_z - q_x p_y + q_y p_x + q_z p_w) \end{aligned}

ScLerp(Screw Linear Interpolation)

デュアルクォータニオンによる補間について少し触れておきます．下図を見ながら説明していきます．

A から B に変換するデュアルクォータニオンを考えます．回転軸は \vec{n} ，回転角度は \theta ，平行移動 \vec{t} です．この変換の補間を考えると，パラメータ \alpha を使って単純な線形補間の場合， A から回転軸 \vec{n} に対して \alpha\theta 分回転した後， \alpha\vec{t} 移動することになります．クォータニオンによる回転になっていますが，平行移動は A と B を結ぶ直線上で線形に補間されます．ScLerpではこの平行移動の補間を回転軸 \vec{n} に沿って行うことで，ねじるような軌跡になります．このようなデュアルクォータニオンを作成するために，まず d を計算します．

d = \vec{t}\cdot\vec{n} = (2qp^{\ast})\cdot\vec{n}

で求まります．そして，移動ベクトル \vec{m} を求めると

\begin{aligned} m &= \frac{1}{2}\left(\vec{t}\times\vec{n} + \vec{n}\times(\vec{t}\times\vec{n})\cot\frac{\theta}{2}\right) \\[2ex] &= \frac{1}{2}\left(\vec{t}\times\vec{n} + (\vec{t}-d\vec{n})\cot\frac{\theta}{2}\right) \end{aligned}

となります． d , \vec{m} , \vec{n} , \theta をscrew parameterと言います．この式の両辺に \sin\frac{\theta}{2} を掛けて変形すると

\sin\frac{\theta}{2}\vec{m} + \frac{d}{2}\cos\frac{\theta}{2}\vec{n} = \frac{1}{2}\left(\vec{t}\times\vec{n}\sin\frac{\theta}{2}+\cos\frac{\theta}{2}\vec{t}\right)

となります．この式を使って虚部を求めると

\begin{aligned} \frac{1}{2}\vec{t}r &= \frac{1}{2}\vec{t}\left(\cos\frac{\theta}{2}+\vec{n}\sin\frac{\theta}{2}\right) \\[2ex] &= -\frac{1}{2}\vec{t}\cdot\vec{n}\frac{\theta}{2} + \frac{1}{2}\vec{t}\cos\frac{\theta}{2} + \frac{1}{2}\vec{t}\times\vec{n}\sin\frac{\theta}{2} \\[2ex] &= -\frac{d}{2}\sin\frac{\theta}{2} + \frac{1}{2}\left(\vec{t}\cos\frac{\theta}{2} + \vec{t}\times\vec{n}\sin\frac{\theta}{2}\right) \\[2ex] &= -\frac{d}{2}\sin\frac{\theta}{2} + \sin\frac{\theta}{2}\vec{m} + \frac{d}{2}\cos\frac{\theta}{2}\vec{n} \end{aligned}

よって，デュアルクォータニオンは

\begin{aligned} \sigma &= r + \frac{\varepsilon}{2}\vec{t}r \\[2ex] &= \left(\cos\frac{\theta}{2}+\vec{n}\sin\frac{\theta}{2}\right) + \varepsilon\left(-\frac{d}{2}\sin\frac{\theta}{2} + \sin\frac{\theta}{2}\vec{m} + \frac{d}{2}\cos\frac{\theta}{2}\vec{n}\right) \\[2ex] &= \left(\cos\frac{\theta}{2}-\varepsilon\frac{d}{2}\sin\frac{d}{2}\right)+\left( \vec{n}\sin\frac{\theta}{2} + \varepsilon\left(\sin\frac{\theta}{2}\vec{m}+\frac{d}{2}\cos\frac{\theta}{2}\vec{n}\right)\right) \end{aligned}

ここで三角関数の公式を使うと

\begin{aligned} \cos\frac{\theta+\varepsilon d}{2} &= \cos\frac{\theta}{2}-\varepsilon\frac{d}{2}\sin\frac{\theta}{2} \\[2ex] \sin\frac{\theta+\varepsilon d}{2} &= \sin\frac{\theta}{2}+\varepsilon\frac{d}{2}\cos\frac{\theta}{2} \end{aligned}

という関係になるので，整理すると

\sigma = \cos\frac{\theta+\varepsilon d}{2} + \sin\frac{\theta + \varepsilon d}{2}\left(\vec{n}+\varepsilon \vec{m}\right).

となります．これがScLerpのデュアルクォータニオンです．パラメータ \alpha を使って補間をする場合は

\sigma(\alpha) = \cos\frac{(\theta+\varepsilon d)\alpha}{2} + \sin\frac{(\theta + \varepsilon d)\alpha}{2}\left(\vec{n}+\varepsilon \vec{m}\right).

となります．

📌 最後に

今回もなるべく新しい用語を出さずに，わかりやすくなるようにしてみたつもりです．また，数式も途中の式をなるべく表記するようにしました．今回の内容は最低限のものですので，興味があれば他の記事などを見てみてください．
少しでも参考になれば幸いです．

📌 参考文献

• 久保裕一郎, 「Enter the 3D Programming　第７回　回れクォータニオン」C Magazine 2001年4月号, ソフトバンク パブリッシング，2001
• Eric Lengyel著，狩野智英訳「ゲームプログラミングのための３Dグラフィックス数学」ボーンデジタル社，2006
• 広江克彦，「内積、外積の公式 」
• Kavan, L., Collins, S., Zara, J., O'Sullivan, C., Skinning with Dual Quaternions, 2007
• Yan-Bin Jia, Dual Quaternion