はじめに
最近四元数について勉強したのでまとめました
本記事の対象者
- 四元数とは何か知りたい方
- 四元数で回転を表す方法を知りたい方
- 四元数でなぜ回転を表せるのか知りたい方
- LookAt, Slerp, Quaternion.idnetity の意味を知りたい方
1. 複素数
この項は複素数が２次元の回転を表せることを確認するためのものです。
複素数平面がわかっている方は飛ばしてください。
1-1. 虚数
実数の世界では、次の方程式は解なしとなります。
$$x^2 = -1$$
ここで、2乗すると $-1$ になる数である虚数単位 $i$ を導入することで、解を得ることができるようになります
$$x = \pm i$$
1-2. 複素数
実数と虚数を足し合わせて実数を拡張したもの
$x =a + bi$
1-3. 複素数平面
実部を横軸、虚部を縦軸として、平面を作ります。
複素数 $x = a + bi$ は、複素数平面上で $(a, b)$ に対応します。
大きさ $r$ , 実軸からの偏角 $\theta$ を用いる極形式で表すと以下のようになります。
$x = r (\cos\theta + i \sin\theta )$
ノルム(大きさ)
大きさを三平方の定理を用いて求めます。
$|a + bi| = \sqrt{a^2 + b^2}$
$|r (\cos\theta + i \sin\theta )| = \sqrt{r^2(\cos^2\theta + \sin^2\theta)} = r$
共役な複素数
複素数 $x = a+bi$ の虚部を反転した複素数 $x^* = a-bi$ を $x$ の共役な複素数と呼びます。
$x^*$ を用いると、ノルムや実部・虚部を $a, b$ を用いずに表すことができます。
ノルム
$$|x| = \sqrt{xx^*}$$
実部
$$Re(x) = \frac{x + x^*}{2}$$
虚部
$$Im(x) = \frac{x - x^*}{2i}$$
1-4. オイラーの公式
オイラーの公式
$$e^{i\theta} = \cos\theta + i\sin\theta$$
を用いると、複素数を指数で表せます。
$x = r e^{i \theta}$
1-5. 掛け算を用いた平面上の回転
極形式の複素数 $x, y$ を掛け算してみます
$x = r_1 e^{i \theta_1}, y = r_2 e^{i \theta_2}$
$x y = r_1 r_2 e^{i (\theta_1 + \theta_2)}$
計算結果は、
- 大きさ : 元の複素数の大きさの積
- 偏角 : 元の複素数の偏角の和
の複素数となります。
この性質によって、複素数を使用して平面上の回転を表せすことができます。
特にノルムが1の複素数 $e^{i\theta}$ を掛ければ大きさを変えずに回転できます。
2. 四元数(Quaternion)とは？
四元数について、次のように認識している方は多いのではないでしょうか？
- 複素数を拡張したもの
- 三次元の回転に使用できる
この項では、四元数の数学的な定義と性質を確認していきます。
2-1. 四元数の定義
四元数とは、$i, j, k$ の３つの虚数単位を用いる数の体系です。
$$q = w + i x + j y + k z$$
$i , j , k$ の3要素をベクトルとしてみて、スカラーとベクトルの和として扱います。
$$q = w + (x, y, z) = w + \boldsymbol{v}$$
虚数単位 $i, j, k$ については、以下の関係が成り立ちます。
\begin{align}
i^2 = j^2 = k^2 = -1 \\
ij = -ji = k \\
jk = -kj = i \\
ki = -ik = j
\end{align}
2-2. 四元数の性質
加法の結合法則
$(p + q) + r = p + (q + r)$
加法の交換法則
$q + p = p + q$
乗法の結合法則
$(p \times q) \times r = p \times (q \times r)$
乗法の交換法則は成立しない
$ij = -ji$
これは３次元の回転の計算には順番が重要であることを示しています。
分配法則
$p \times (q + r) = p \times q + p \times r$
2-3. そのほかのルール
共役
$q = w + i x + j y + k z = w + \boldsymbol{v}$
に対して、共役な四元数 $q^*$ を次のように定義します。
$q^* = w - i x - j y - k z = w - \boldsymbol{v}$
回転を考える上でうれしい性質が成り立ちます（ → 6-1. 連続回転 ）
$$q_1^* q_2^* = (q_2q_1)^*$$
ノルム
複素数をそのまま拡張する形で定義します。
$|q| = \sqrt{q q^{*}} = \sqrt{w^2 + x^2 + y^2 + z^2}$
逆元
$$\frac{a + ib + jc + kd}{w + ix + jy + kz}$$
上のような四元数同士の割り算はどう考えるでしょうか？
結論としては、
$$q q^{-1} = 1$$
となるような $q^{-1}$ を掛けることで実現します。
$$q q^{*} = |q|^2$$
$$q\frac{q^{*}}{|q|^2} = 1$$
$$\boldsymbol{q^{-1}=\frac{q^{*}}{|q|^2}}$$
よって最初の割り算は次のように計算できます。
$$\frac{a + ib + jc + kd}{w + ix + jy + kz} = \frac{(a + ib + jc + kd)(w -ix-jy-kz)}{w^2 + x^2 + y^2 + z^2}$$
3. なぜ三元数ではないのか
複素数では、実部と虚部の２つの変数を平面に対応させる形で表現できました。
なぜ３次元の回転に四元数を用いるのでしょうか？
3-1. i×j 問題
三元数(Ternion) $\hat{t} = a + ib + jc$ を考えてみます。
実部を$x$軸 , $i$を$y$軸 , $j$を$z$軸 にそれぞれ対応させて、$\cos \theta + i \sin \theta$ を掛けてみます。
これが $z$ 軸周りの回転を表せていると嬉しそうです。
$e^{i \theta} \hat{t} = (a\cos\theta - b\sin\theta) +\boldsymbol{i}(b\cos\theta + a \sin\theta) +\boldsymbol{j}(c\cos\theta) + \underline{\boldsymbol{ij}} (c\sin\theta)$
実部と $i$ は$z$軸まわりに回転できましたが、動かしたくない $j$ も回転してしまった上、 $ij$ という新たな数が現れてしまいました。この$ij$をどうにか解釈する必要があります。
3-2. i j を無理やり解釈してみる
i j = 0
$e^{i \theta} \hat{t} = (a\cos\theta - b\sin\theta) +\boldsymbol{i}(b\cos\theta + a \sin\theta) +\boldsymbol{j}(c\cos\theta)$
一応計算結果に意味を持たせることはできましたが、ノルムを計算してみると...
$|\hat{t}| = a^2 + b^2 + c^2$
$|e^{i \theta} \hat{t}|^2 = a^2 + b^2 + c^2\cos^2\theta$
大きさ1の要素 $e^{i\theta}$ を掛けたのに大きさが変わってしまっています。
これでは回転を考えるには扱いづらいです。
i j = ±1
$e^{i \theta} \hat{t} = (a\cos\theta - b\sin\theta \pm c\sin\theta) +\boldsymbol{i}(b\cos\theta + a \sin\theta) +\boldsymbol{j}(c\cos\theta)$
$|e^{i \theta} \hat{t}|^2 = a^2 + b^2 + c^2 + \mp bc\sin^2\theta + \pm ca\sin\theta\cos\theta$
この場合も同様の問題が発生します。
i j = k
$ij$を新たな要素としてそのまま解釈してノルムを計算してみます。
\begin{split}
e^{i \theta} \hat{t} &= (a\cos\theta - b\sin\theta) +\boldsymbol{i}(b\cos\theta + a \sin\theta) +\boldsymbol{j}(c\cos\theta) + \boldsymbol{ij} (c\sin\theta) \\
|e^{i \theta} \hat{t}| &= (a\cos\theta - b\sin\theta)^2 +(b\cos\theta + a \sin\theta)^2 + (c\cos\theta)^2 + (c\sin\theta)^2 \\
&= a^2 + b^2 + c^2
\end{split}
これは、大きさ1の要素をかけても大きさを変えないことを実現できています。
以上の理由から、$\boldsymbol{ij}$をあらたな軸 $\boldsymbol{k}$ として考えた方が性質が良く、三元数ではなく四元数となっています。
$$q = w + ix + jy + kz$$
ハミルトンはこれが成り立ってほしかったらしい
|p q| = |p||q|
4. どう回転に用いるのか
四元数が良い性質をもつことは確認できたましたが、
ただ掛けるだけでは、回転させたくない要素も変わってしまうことも確認できました。
この項では、四元数をどのように用いると回転を表せるのかを示します。
4-1. 単位四元数
回転に用いるということで、複素数を参考に四元数を指数表記してみます。
参考 : 複素数の指数表記
$$e^{i\theta} = \cos\theta + i\sin\theta, |e^{i\theta}| = 1$$
単位四元数を$q_u$として
$$q_u = \cos\theta + \boldsymbol{v_u}\sin\theta, |q_u| = 1$$
と書いたとき、これを満たすような $\boldsymbol{v_u} = i v_x + j v_ y + k v_z$ はどんなものか考えます。
$$|q_u| = \sqrt{\cos^2\theta + |\boldsymbol{v_u}|^2\sin^2\theta} = 1$$
$$|\boldsymbol{v_u}| = 1$$
ノルムをとると、$\boldsymbol{v_u}$は単位ベクトルである必要があることがわかります。
$\boldsymbol{v_u}$ を $\boldsymbol{I}$ と書き直して次のように四元数を指数表記します。
$$q_u = \cos\theta + \boldsymbol{I}\sin\theta = e^{\boldsymbol{I}\theta}$$
$q_u$ のようなノルムが$1$の回転に用いる四元数を単位四元数と呼びます。
4-2. 回転の公式
いきなりですが、ゴールである回転の公式を見てみます。
点の座標 : $\boldsymbol{p} = (x, y, z)$
回転軸 : $\boldsymbol{I}$
回転する角度(右手系) : $\theta$
回転前座標四元数 $P = 0 + \boldsymbol{p} = ix + jy + kz$
回転後座標四元数 $P^{'} = 0 + \boldsymbol{p^{'}} = ix^{'} + jy^{'} + kz^{'}$
とするとき、原点から $\boldsymbol{p}$ 離れた点を軸 $\boldsymbol{I}$ まわりに $\theta$ 回転した後の座標は次のように表せます。
$$P^{'} = e^{\boldsymbol{I}\frac{\theta}{2}} P e^{\boldsymbol{I}(-\frac{\theta}{2})}$$
次の項から、なぜこの式が回転を表せるのか見ていきます。
5. なぜ回転が表せるのか？
この項では、前項で登場した $P^{'} = e^{\boldsymbol{I}\frac{\theta}{2}} P e^{\boldsymbol{I}(-\frac{\theta}{2})}$ がなぜ回転を表せるのかを考えていきます。
5-1. 視覚化してみる
いきなり式を見てもイメージしにくいので、Unityを使って視覚化しました。
左が $\frac{\pi}{2}$ 回転、右が $\pi$ 回転の計算を２段階に分けて行う様子をGIFにしたものです。
見てみると、1段階目の不要な回転を2段階目で直しながら回してそうです。
5-2. z軸まわりの回転を理解する
視覚化から立てた予想が正しいことを確認していきます。
最初に簡単な例として、$z$軸まわりで $\theta$ 回転させることを考えます。
点の座標 : $\boldsymbol{p} = (p_x, p_y, p_z) = (r\cos\phi, r\sin\phi, p_z)$
回転軸 : $\boldsymbol{I}=(0, 0, 1)$
ここで座標四元数 $P$ を後の計算結果が見やすくなるように工夫して書いておきます
\begin{split}
P &= 0 + ip_x + jp_y + kp_z \\
&= p_z \sin0 + i r\cos\phi + j r\sin\phi + k p_z \cos0
\end{split}
計算を進めると
e^{\boldsymbol{I}\frac{\theta}{2}} P
= p_z \sin(0 - \frac{\theta}{2}) +
\begin{pmatrix}
r\cos(\phi + \frac{\theta}{2}) \\
r\sin(\phi + \frac{\theta}{2}) \\
p_z\cos(0 - \frac{\theta}{2})
\end{pmatrix}
１段階目では $xy$ 平面では目的の方向に回転し、$wz$ 平面では不要な回転が発生しています。
(e^{\boldsymbol{I}\frac{\theta}{2}} P) e^{\boldsymbol{I}(-\frac{\theta}{2})}
= p_z \sin((0 - \frac{\theta}{2}) +\frac{\theta}{2}) +
\begin{pmatrix}
r\cos(\phi + \frac{\theta}{2} +\frac{\theta}{2}) \\
r\sin(\phi + \frac{\theta}{2} +\frac{\theta}{2}) \\
p_z \cos((0 - \frac{\theta}{2}) +\frac{\theta}{2})
\end{pmatrix}
２段階目では $xy$ 平面では目的の方向に回転し、$wz$ 平面では不要な回転が戻されています。
$z$ 軸まわりの回転では予想通りの結果になっています。
5-3. ロドリゲスの回転公式を理解する
任意の軸の回転を考える前に参考になる考え方を見ておきます。
ロドリゲスの回転公式とは、
点 $\boldsymbol{p}$ を $\boldsymbol{I}$ 軸まわりに $\theta$ 回転させた点 $\boldsymbol{p^{'}}$ が次の式で表せるというものです。
$$\boldsymbol{p^{'}} = (\boldsymbol{p}\cdot\boldsymbol{I})\boldsymbol{I}+(\boldsymbol{p} - (\boldsymbol{p}\cdot\boldsymbol{I})\boldsymbol{I})\cos\theta + (\boldsymbol{I}\times \boldsymbol{p})\sin\theta$$
つまり四元数の回転の式 $e^{\boldsymbol{I\frac{\theta}{2}}}Pe^{-\boldsymbol{I\frac{\theta}{2}}}$ のベクトル版なのですが、こちらは以外と簡単に導出することができます。
$\boldsymbol{I}$ を$\hat{z}$軸とみなして、$\boldsymbol{p}$ から$\boldsymbol{I}$ に下した垂線を$\hat{x}$軸とし、さらにこの２つと直交する軸を$\hat{y}$軸とします。
この座標変換によって $\boldsymbol{p}$ の $\hat{z}$軸成分は変化せず、$\hat{x}\hat{y}$ 平面上で回転させることだけを考えればよくなります。後は三角関数を使ってxy平面上で回転させ、ベクトルを足し合わせることで回転後の点 $\boldsymbol{p^{'}}$ を求めることができます。
ロドリゲスの回転公式に関する詳しい内容や、導出に関しては、
が大変参考になります。
5-4. ロドリゲス回転公式の座標変換方式を利用して
任意の軸まわりの回転を理解する
ロドリゲスの回転公式の座標の取り方を利用して、任意の軸$\boldsymbol{I}$まわりの回転を$xy$平面上の回転に変換していきましょう。
$\boldsymbol{I}$の向きを$\hat{z}$軸として、直交座標系をとります。
$\hat{x}$ 軸と $\hat{y}$ 軸は $\hat{z}$ 軸と直交するようにとればなんでもいいですが、ロドリゲスの回転公式のときと合わせて $\boldsymbol{p}$ から、$\boldsymbol{I}$ におろした垂線を $\hat{x}$ 軸、この２つと直交する軸を $\hat{y}$ 軸とします。
この座標系の基本ベクトル $(\hat{\boldsymbol{e_1}}, \hat{\boldsymbol{e_2}}, \hat{\boldsymbol{e_3}})$ を次のように表します
$\hat{\boldsymbol{e_1}} = (e_{11}, e_{12}, e_{13})$
$\hat{\boldsymbol{e_2}} = (e_{21}, e_{22}, e_{23})$
$\hat{\boldsymbol{e_3}} = (e_{31}, e_{32}, e_{33})$
各成分については、$\boldsymbol{I}$ や $\boldsymbol{p}$ を用いて求めることができますが、今回欲しい結論には必要ありません。
この $\hat{x}\hat{y}\hat{z}$ 座標系と元の $xyz$ 座標系の変換は次のような行列 $T$ を用いて表せます。
\begin{split}
\boldsymbol{v} &= T \hat{\boldsymbol{v}} , \\
T &=
\begin{pmatrix}
e_{11} & e_{12} & e_{13} \\
e_{21} & e_{22} & e_{23} \\
e_{31} & e_{32} & e_{33} \\
\end{pmatrix}
\end{split}
上の式に、元の座標系の虚数単位 $i, j, k$ と、座標変換した座標系の虚数単位 $\hat{i}, \hat{j}, \hat{k}$ を代入すると、
\begin{split}\begin{pmatrix}
\hat{i} \\
\hat{j} \\
\hat{k} \\
\end{pmatrix}
&= T
\begin{pmatrix}
i \\
j \\
k \\
\end{pmatrix} = \begin{pmatrix}
ie_{11} + je_{12} + ke_{13} \\
ie_{21} + je_{22} + ke_{23} \\
ie_{31} + je_{32} + ke_{33} \\
\end{pmatrix}\end{split}
次のような関係がわかります。
\begin{align}
\hat{i} = ie_{11} + je_{12} + ke_{13} \\
\hat{j} = ie_{21} + je_{22} + ke_{23} \\
\hat{k} = ie_{31} + je_{32} + ke_{33}
\end{align}
ここで、単位四元数と座標四元数を $\hat{i}, \hat{j}, \hat{k}$ を使って書き直すと、
\begin{split}
e^{\boldsymbol{I}\frac{\theta}{2}} &= \cos\frac{\theta}{2} + \boldsymbol{I}\sin\frac{\theta}{2}
&= \cos\frac{\theta}{2} + \hat{k}\sin\frac{\theta}{2} (\because \boldsymbol{I} = \hat{k})
\end{split}
$$P = 0 + \hat{i}\hat{x} + \hat{j}\hat{y} + \hat{k}\hat{z}$$
のようになり、元の座標系の $z$ 軸まわりの回転と非常に似ています。
$$
e^{k\frac{\theta}{2}} = \cos\frac{\theta}{2} + k\sin\frac{\theta}{2}
$$
$$P = 0 + ix + jy + kz$$
そのため、座標変換した後の $\hat{i}, \hat{j}, \hat{k}$ も法則
\begin{align}
\hat{i}^2 = \hat{j}^2 = \hat{k}^2 = -1 \\
\hat{i}\hat{j} = -\hat{j}\hat{i} = \hat{k} \\
\hat{j}\hat{k} = -\hat{k}\hat{j} = \hat{i} \\
\hat{k}\hat{i} = -\hat{i}\hat{k} = \hat{j}
\end{align}
を満たすとすると、
$\hat{z}$軸まわりの $\hat{x}\hat{y}$ 平面上の回転として考えることができます。
実際に計算して法則が成り立つことを示します。
$$
\hat{i} ^2 = (ie_{11} + je_{12} + ke_{13} )(ie_{11} + je_{12} + ke_{13} ) = -( e_{11}^2 + e_{12}^2 + e_{13}^2 ) = -1
$$
$\hat{j}^2 = -1, \hat{k}^2 = -1$ も同様に成り立ちます。
$\hat{i} \hat{j}$ の値を求めるために、$\hat{\boldsymbol{e_1}}, \hat{\boldsymbol{e_2}}, \hat{\boldsymbol{e_3}}$ を別の形で表しておきます。
\begin{split}
\hat{\boldsymbol{e_3}} =
\hat{\boldsymbol{e_1}} \times \hat{\boldsymbol{e_2}} =
\begin{pmatrix}
e_{12} e_{23} - e_{13} e_{22} \\
e_{13} e_{21} - e_{11} e_{23} \\
e_{11} e_{22} - e_{12} e_{21} \\
\end{pmatrix}
\end{split}
\begin{split}
\hat{\boldsymbol{e_1}} =
\hat{\boldsymbol{e_2}} \times \hat{\boldsymbol{e_3}} =
\begin{pmatrix}
e_{22} e_{33} - e_{23} e_{32} \\
e_{23} e_{31} - e_{21} e_{33} \\
e_{21} e_{32} - e_{22} e_{31} \\
\end{pmatrix}
\end{split}
\begin{split}
\hat{\boldsymbol{e_2}} =
\hat{\boldsymbol{e_3}} \times \hat{\boldsymbol{e_1}} =
\begin{pmatrix}
e_{32} e_{13} - e_{33} e_{12} \\
e_{33} e_{11} - e_{31} e_{13} \\
e_{31} e_{12} - e_{32} e_{11} \\
\end{pmatrix}
\end{split}
準備ができたので、$\hat{i} \hat{j}$ を計算します。
\begin{split}
\hat{i} \hat{j} &=
(ie_{11} + je_{12} + ke_{13} )(ie_{21} + je_{22} + ke_{23} ) \\
&=
\begin{pmatrix}
-( e_{11} e_{21} + e_{12} e_{22} + e_{13} e_{23} ) \\
e_{12} e_{23} - e_{13} e_{22} \\
e_{13} e_{21} - e_{11} e_{23} \\
e_{11} e_{22} - e_{12} e_{21} \\
\end{pmatrix} \\
&=\begin{pmatrix} - \boldsymbol{\hat{e_x}} \cdot \boldsymbol{\hat{e_y}} = 0 \\
e_{31} \\
e_{32} \\
e_{33} \\
\end{pmatrix} \\
&= \hat{k}
\end{split}
$\hat{j}\hat{k} = \hat{i} , \hat{k}\hat{i} = \hat{j}$ についても同様に成り立ちます。
さらに、$\hat{i}, \hat{j}, \hat{k}$ は実部を含まないので、 $q_1 q_2 = - q_2 q_1$ が成り立つため、
\begin{align}
\hat{i}\hat{j} = -\hat{j}\hat{i}\\
\hat{j}\hat{k} = -\hat{k}\hat{j}\\
\hat{k}\hat{i} = -\hat{i}\hat{k}
\end{align}
も成立します。
以上より虚数単位の性質が変換後も成り立つため、任意の軸周りの回転を表せていることが確認できました。
この項は以下の記事を参考にさせていただきました。
クォータニオンによる回転の原理をスッキリ理解したい
6. Tips
ここからは四元数や回転にまつわる細かい話を見ていきます。
6-1. 連続回転
四元数の共役について、
q_{1}^{*}q_{2}^{*} = (q_{2}q_{1})^{*}
が成り立ちます。
これを用いると、$q_1$ 回転させた後に $q_2$ 回転させる動作は次のように変形できます。
q_2 (q_1 p q_1^*) q_2^* = (q_2q_1) p (q_1^*q_2^*) = (q_2q_1) p (q_2q_1)^*
これによって、各軸ごとに回転させるオイラー角から四元数に変換したり、役割の異なる四元数を順番に適用することを１つの四元数で表すことができます。
6-2. Quaternion.identity
Unity等を触っていると、よく見る無回転を表すQuaternion.identityですが、実際はどのような値で、どのように使用されているのでしょうか？
まず値について、.NET System.Numerics の資料を参考にすると、
値
(0, 0, 0, 1)が である四元数。
とあります。
$$
q = 1 + (0, 0, 0) = 1
$$
そもそもUnityのtransform.rotationは、「自分がそのように回転していたら子要素や頂点はどこに移動するか」という計算に使用される値です。自分は点のようなものです
例えば、(0, 0, 0)にある長方形のオブジェクトがz軸まわり$\frac{\pi}{2}$に回転させる表現は、頂点の位置をオブジェクトの原点を中心に回転させた位置に移動させることで実現します。
そのため、回転とは位置の移動または座標変換のことになります。
自分のrotation $q$
子要素や頂点の回転前位置 $p$
子要素や頂点の回転後位置 $p'$
$$p' = q p q^*$$
ここで、Quaternion.identity = 1 なので、
$$p' = p$$
となるため、子要素や頂点の位置を移動しないことになり、無回転を表します。
実際は頂点の回転(座標変換)はシェーダーで行うため、回転行列を用います。
6-3. LookAt
UnityにはLookAtという、指定したtargetの方を向くようにrotationを変更するメソッドがあります。
これを自作するにはどうすればよいでしょうか？
transform.LookAt(target, Vector3.up);
targetの方向へ回転させる q を求める
targetの方向を向くということは、ある四元数 $q$ を使用して回転したときに、オブジェクトの前方向がtargetのある方向に変換されているということです。
オブジェクトの前方向 : $\boldsymbol{f}=(1, 0, 0)$
回転後の前方向 : $\boldsymbol{F} = \frac{\boldsymbol{v_{target}}}{|\boldsymbol{v_{target}}|}$
回転に用いる四元数 : $q$
$$\boldsymbol{F} = q \cdot \boldsymbol{f} \cdot q^{*}$$
ただ、下の動画の通り、targetの方向を向いている状態は無数にあり、上の式だけでは $q$ を一意に求めることはできません。
解を絞る
UnityのLookAt関数の動作をよく見てみると、回転後のオブジェクトの左右方向が平面内に収まっていることがわかります。
これを条件として加えることで、$q$について解くことができるようになります。
オブジェクトの上方向 : $\boldsymbol{u}=(0, 0, 1)$
オブジェクトの左方向 : $\boldsymbol{l}=(0, 1, 0)$
回転後の左方向 : $\boldsymbol{L} = \boldsymbol{F} \times \boldsymbol{u}$
$$\boldsymbol{L} = q \cdot \boldsymbol{l} \cdot q^{*}$$
もっと簡単に q を求める
求めることができるといっても方程式を解くのは大変です。
そこで、オイラー角の考え方を用いることで目的の $q$ を求めます。
オイラー角とは、$x, y, z$軸周りの回転を繰り返すことで目的の回転を表す方式です。この考えを用いると、LookAtは下の画像のように回転させることで実現できます。
①②の操作を順番に行うことは、 6-1. 連続回転 より、
①を行う四元数 $q_1=e^{(0, 1, 0) \phi}$ と、
②を行う四元数 $q_2=e^{(0, 0, 1) \theta}$ を掛けた
$$q = q_2 q_1$$
を使用することで実現できます。
6-4. Slerp
Slerpとは、2つのrotation($q_1, q_2$)間を等速で回転するように補完する関数です。
まず、Lerp(Liner interpolation, 線形補完)という２点間を直線で補完する方式があります。
$$p_t = (1-t)p_1 + tp_2 , (0 \leq t \leq 1)$$
Lerpに対して、Slerp(Spherical linear interpolation, 球面線形補完)とは、２点間の補完を球面上で行う方式です。
$$p_t = \frac{\sin((1-t)\Omega)}{\sin\Omega} p_1 + \frac{\sin(t\Omega)}{\sin\Omega} p_2 , (0 \leq t \leq 1, \cos\Omega = \frac{p_1 \cdot p_2}{|p_1||p_2|})$$
2点間移動ではLerpを使用しますが、回転にLerpを使用すると、一定速度の回転にはなりません。
下の動画は、赤がLerpを用いた補完、青がSlerpを用いた回転の補完です。
Lerpの方では、加減速を繰り返していることがわかります。
そのため、2つのrotationを等速で移動したい場合は、基本的にはSlerpを使用していくことになりますが、２つの四元数間の角度が小さい場合には、LerpとSlerpはほぼ一致します。($\Omega \to 0$ で一致)
下の動画は $\Omega = \frac{\pi}{3}$ の場合です。上の場合($\Omega = \pi$) と比べるとLerpとSlerpの差はかなり小さいですね。
.NET System.NumericsのSlerpの実装では、これを利用して $\Omega$ が小さいときにはLerpが使われています。
これは三角関数の計算の有無によってLerpの方が計算量が小さくなるためです。
src/libraries/System.Private.CoreLib/src/System/Numerics/Quaternion.cs 616行目
private const float SlerpEpsilon = 1e-6f;
// 省略
if (cosOmega > (1.0f - SlerpEpsilon))
{
// Too close, do straight linear interpolation.
s1 = 1.0f - t;
s2 = (flip) ? -t : t;
}
Slerpの詳しい内容についてはこちらの記事が大変参考になります。
ねこじゃらシティ 【Unity】Vector3.Lerp／Slerpの使い方と内部挙動
視覚化に使用したUnityプロジェクト
本記事のGIF作成に使用したプロジェクトです。
GitHub URL
参考文献
第８章 四元数
四元数を発見法的に導く
高校数学の美しい物語 ロドリゲスの回転公式
クォータニオンによる回転の原理をスッキリ理解したい
回転行列、クォータニオン(四元数)、オイラー角の相互変換
かみのメモ 回転ベクトル・回転行列・クォータニオン・オイラー角についてまとめてみた
Quaternion.identity
src/libraries/System.Private.CoreLib/src/System/Numerics/Quaternion.cs
Mathlog 四元数による回転の記述