クォータニオンと回転行列の関係

はじめに

いま話題の3D Gaussian Splattingを最初に提案した論文「3D Gaussian Splatting for Real-Time Radiance Field Rendering」の中に、クォータニオン（Quaternion）を構成する4つの成分から回転行列を導出する式が記載されていた。この表式はクォータニオンの入門書には必ず示されているものである。今回はこの式の導出を行う。論文自体の紹介は機会を改める。

クォータニオンとは

クォータニオンとは、主に3Dコンピュータグラフィックス（3DCG）の世界で使われる手法である。クォータニオンを使うと回転を直感的に表現することができ、回転の関わる計算が容易になるため、3DCGの世界では頻繁に使われる。まず最初にクォータニオン\tilde{q}を定義しよう。

    \begin{align*} \tilde{q} = q_0+iq_1+jq_2+kq_3 \end{align*}

ここで、q_n\;(n=0,1,2,3)は実数であり、i,j,kに対しては次式が成り立つ。

    \begin{align*} & i^2=j^2=k^2=ijk=-1 \\ & ij=-ji=k,\;jk=-kj=i,\;ki=-ik=j \end{align*}

クォータニオンは複素数を拡張したような数である。さて、計算時の表記を簡略化するため以下の表記を導入する。

    \begin{align*} \tilde{q} &= q_0+q \\ q&=iq_1+jq_2+kq_3 \\ &=\boldsymbol{e}^T\cdot\boldsymbol{q} \\ \boldsymbol{e} &= (i,j,k)^T \\ \boldsymbol{q} &= (q_1,q_2,q_3)^T \end{align*}

q_0、qをそれぞれクォータニオンの実部、虚部と呼ぶ。また、Tは転置を表し、\boldsymbol{e}^T\cdot\boldsymbol{q}は2つのベクトル間の内積である。次に、2つのクォータニオンの積を計算してみよう。

    \begin{align*} \tilde{q}\tilde{p} &=\left(q_0+q\right)\left(p_0+p\right)\\ &=q_0p_0+q_0p+p_0q+qp \end{align*}

ここで

(1)   \begin{align*} qp=-\boldsymbol{q}^T\cdot\boldsymbol{p}+\boldsymbol{e}^T\cdot(\boldsymbol{q}\times\boldsymbol{p}) \end{align*}

が成り立つ。ただし、p=ip_1+jp_2+kp_3,\;\boldsymbol{p}=(p_1,p_2,p_3)^Tとおいた。qpはクォータニオンの虚部同士の積であり、その積が\boldsymbol{q}と\boldsymbol{p}の内積と外積から構成されることに注意する。また、\tilde{q}\tilde{p}は実部と虚部に分けることができるから、クォータニオン同士の積もクォータニオンになることが分かる。次に、\tilde{q}に共役なクォータニオン\tilde{q}^{*}を次式で定義する。

    \begin{align*} \tilde{q}^*=q_0-q \end{align*}

これを使うと

    \begin{align*} \tilde{q}\tilde{q}^*=q_0^2+q_1^2+q_2^2+q_3^2 \end{align*}

を得る。これより\tilde{q}の長さを次式で定義できる。

    \begin{align*} \|\tilde{q}\|=\sqrt{q_0^2+q_1^2+q_2^2+q_3^2} \end{align*}

回転行列

ここではクォータニオンから離れて回転行列について解説する。話を簡単にするため、z軸周りを\thetaだけ反時計回りに回転する行列を取り上げる。

右図において次式が成り立つ。

    \begin{align*} (x,y)&=(l\cos{\theta_0},l\sin{\theta_0})\\ (x^{'},y^{'})&=(l\cos{\theta_1},l\sin{\theta_1})\\ \end{align*}

\theta_1=\theta_0+\thetaを用いてx^'を変形すると

    \begin{align*} x^' &=l\cos{(\theta_0+\theta)}\\ &=l(\cos{\theta_0}\cos{\theta}-\sin{\theta_0}\sin{\theta})\\ &=\cos{\theta}\;x-\sin{\theta}\;y \end{align*}

を得る。同様にy^'を変形すると

    \begin{align*} y^' &=l\sin{(\theta_0+\theta)}\\ &=l(\sin{\theta_0}\cos{\theta}+\cos{\theta_0}\sin{\theta})\\ &=\cos{\theta}\;y+\sin{\theta}\;x\\ &=\sin{\theta}\;x+\cos{\theta}\;y \end{align*}

となる。この2式から回転行列を構成することができる。

    \begin{align*} \begin{pmatrix} x^' \\ y^' \end{pmatrix} = \begin{pmatrix} \cos{\theta} & -\sin{\theta}\\ \sin{\theta} & \cos{\theta} \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} \end{align*}

これを3次元に拡張すれば次式を得る。

    \begin{align*} \begin{pmatrix} x^' \\ y^' \\ z^' \end{pmatrix} = \begin{pmatrix} \cos{\theta} & -\sin{\theta} &0\\ \sin{\theta} & \cos{\theta} &0\\ 0 & 0 &1 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} \end{align*}

右辺に現れる行列が、z軸周りを反時計回りに\thetaだけ回転する行列R_z(\theta)である。

クォータニオンによる回転行列の表記

上で求めた回転行列R_z(\theta)を、クォータニオンを用いて書くことができる。最初に次のクォータニオンを用意する。これは単位ベクトル\boldsymbol{n}の周りを反時計回りに\thetaだけ回転するクォータニオンを表す（右図参照）。

    \begin{align*} \tilde{q}&=\alpha+\beta n\\ \alpha&=\cos{\frac{\theta}{2}}\\ \beta&=\sin{\frac{\theta}{2}}\\ n&=in_1+jn_2+kn_3\\ \boldsymbol{n}&=(n_1,n_2,n_3)^T\\ \|\boldsymbol{n}\|&=1 \end{align*}

クォータニオンで回転を表現する際は、\|\tilde{q}\|=1が要求される。上式がこれを満たしていることは以下のように確認できる。

    \begin{align*} \tidle{q}\tilde{q}^{*} &=(\alpha+\beta n)(\alpha-\beta n)\\ &=\alpha^2-\beta^2nn \\ &=\alpha^2-\beta^2\left(-\boldsymbol{n}^T\cdot\boldsymbol{n}+\boldsymbol{e}^T\cdot\left(\boldsymbol{n}\times\boldsymbol{n}\right)\right)\\ &=\alpha^2+\beta^2\\ &=\cos^2{\frac{\theta}{2}}+\sin^2{\frac{\theta}{2}}\\ &=1 \end{align*}

さて

    \begin{align*} r&=ix+jy+kz\\ r^{'}&=ix^{'}+jy^{'}+kz^{'}\\ \boldsymbol{r}&=(x,y,z)^T \\ \boldsymbol{r}^{'}&=(x^{'},y^{'},z^{'})^T \end{align*}

とおくと、\boldsymbol{r}を軸\boldsymbol{n}の周りに回転した\boldsymbol{r}^{'}を、r^{'} =\tilde{q}r\tilde{q}^{*}で求めることができる。これを計算すると

    \begin{align*} r^{'} &=\tilde{q}r\tilde{q}^{*}\\ &=(\alpha+\beta n)r(\alpha-\beta n)\\ &=\alpha^2 r-\alpha\beta(rn-nr)-\beta^2nrn \end{align*}

となる。ここで、rnとnrはクォータニオンの虚部同士の積であり式(1)を用いて計算することができる。

    \begin{align*} rn&=-\boldsymbol{r}^T\cdot\boldsymbol{n}+\boldsymbol{e}^T\cdot(\boldsymbol{r}\times\boldsymbol{n}) \\ nr&=-\boldsymbol{n}^T\cdot\boldsymbol{r}+\boldsymbol{e}^T\cdot(\boldsymbol{n}\times\boldsymbol{r}) \end{align*}

さらに、上の2番目の式を使うと

    \begin{align*} nrn &=\left(-\boldsymbol{n}^T\cdot\boldsymbol{r}+\boldsymbol{e}^T\cdot(\boldsymbol{n}\times\boldsymbol{r})\right)n \\ &=-(\boldsymbol{n}^T\cdot\boldsymbol{r})n+\left[\boldsymbol{e}^T\cdot(\boldsymbol{n}\times\boldsymbol{r})\right]n\\ &=-(\boldsymbol{n}^T\cdot\boldsymbol{r})n-(\boldsymbol{n}\times\boldsymbol{r})^T\cdot\boldsymbol{n}+\boldsymbol{e}^T\cdot \left[(\boldsymbol{n}\times\boldsymbol{r})\times\boldsymbol{n}\right]\;\;\;(\because (1))\\ &=-(\boldsymbol{n}^T\cdot\boldsymbol{r})n-(\boldsymbol{n}\times\boldsymbol{r})^T\cdot\boldsymbol{n}+\boldsymbol{e}^T\cdot \left[-(\boldsymbol{n}^T\cdot\boldsymbol{r})\cdot\boldsymbol{n}+\boldsymbol{r} \right]\\ &=-(\boldsymbol{n}^T\cdot\boldsymbol{r})n+\boldsymbol{e}^T\cdot \left[-(\boldsymbol{n}^T\cdot\boldsymbol{r})\cdot\boldsymbol{n}+\boldsymbol{r} \right]\;\;\;(\because (\boldsymbol{n}\times\boldsymbol{r})^T\cdot\boldsymbol{n}=0)\\ &=\boldsymbol{e}^T\cdot\left\{ -(\boldsymbol{n}^T\cdot\boldsymbol{r})\boldsymbol{n}+ \left[-(\boldsymbol{n}^T\cdot\boldsymbol{r})\cdot\boldsymbol{n}+\boldsymbol{r} \right]\right\}\;\;\;(\because n=\boldsymbol{e}^T\cdot\boldsymbol{n}) \end{align*}

を得る。これらを代入すると

    \begin{align*} r^{'} =\boldsymbol{e}^T\cdot\left[ (\alpha^2-\beta^2)\boldsymbol{r}-2\alpha\beta(\boldsymbol{r}\times\boldsymbol{n})+2\beta^2(\boldsymbol{n}^T\cdot\boldsymbol{r})\boldsymbol{n} \right] \end{align*}

となる。r^{'}=\boldsymbol{e}^T\cdot\boldsymbol{r}^{'}であるから

    \begin{align*} \boldsymbol{r}^{'} =(\alpha^2-\beta^2)\boldsymbol{r}-2\alpha\beta(\boldsymbol{r}\times\boldsymbol{n})+2\beta^2(\boldsymbol{n}^T\cdot\boldsymbol{r})\boldsymbol{n} \end{align*}

を得る。\alpha=\cos{\theta/2},\beta=\sin{\theta/2}を代入して少し計算すると

    \begin{align*} \boldsymbol{r}^{'} =\left(\boldsymbol{r}^T\cdot\boldsymbol{n}\right)\boldsymbol{n} +\left[\boldsymbol{r}-\left(\boldsymbol{r}^T\cdot\boldsymbol{n}\right)\boldsymbol{n}\right]\cos{\theta}+\left(\boldsymbol{n}\times\boldsymbol{r}\right)\sin{\theta} \end{align*}

となる。これが、\boldsymbol{r}を単位ベクトル\boldsymbol{n}の周りに反時計回りに\thetaだけ回転させる式である。いま、\boldsymbol{n}=(0,0,1)^Tとすれば上で求めた回転行列R_z(\theta)が導出される。

クォータニオンによる回転行列の表記（一般化）

論文「3D Gaussian Splatting for Real-Time Radiance Field Rendering」に載っていた式はもっと一般的な回転行列であったので最後にこれを求めてお終いにしたい。回転を表すクォータニオン\tilde{q}=q_0+qをrに作用させる。ここで\|\tilde{q}\|=1が満たされていることに注意する。

    \begin{align*} r^{'} &=\tilde{q}r\tilde{q}^{*}\\ &=\left(q_0+q\right)r\left(q_0-q\right) \end{align*}

これをゴリゴリと計算すると

    \begin{align*} r^{'}= \boldsymbol{e}^T\cdot\left[ \left(q_0^2-\|\boldsymbol{q}\|^2\right)\boldsymbol{r}-2q_0\left(\boldsymbol{r}\times\boldsymbol{q}\right)+2\left(\boldsymbol{q}^T\cdot\boldsymbol{r}\right)\boldsymbol{q} \right] \end{align*}

となり

    \begin{align*} \boldsymbol{r}^{'}= \left(q_0^2-\|\boldsymbol{q}\|^2\right)\boldsymbol{r}-2q_0\left(\boldsymbol{r}\times\boldsymbol{q}\right)+2\left(\boldsymbol{q}^T\cdot\boldsymbol{r}\right)\boldsymbol{q} \end{align*}

を得る。これを行列形式で表記すると

    \begin{align*} \begin{pmatrix} x^' \\ y^' \\ z^' \end{pmatrix} = \begin{pmatrix} q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2-q_0q_3) &2(q_1q_3+q_0q_2)\\ 2(q_1q_2+q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 &2(q_2q_3-q_0q_1)\\ 2(q_1q_3-q_0q_2) & 2(q_2q_3+q_0q_1)&q_0^2-q_1^2-q_2^2+q_3^2 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} \end{align*}

となる。拘束条件q_0^2+q_1^2+q_2^2+q_3^2=1を用いて書き直すと

    \begin{align*} \begin{pmatrix} x^' \\ y^' \\ z^' \end{pmatrix} =2 \begin{pmatrix} \frac{1}{2}-(q_2^2+q_3^2) & q_1q_2-q_0q_3 &q_1q_3+q_0q_2\\ q_1q_2+q_0q_3 & \frac{1}{2}-(q_1^2+q_3^2) &q_2q_3-q_0q_1\\ q_1q_3-q_0q_2 & q_2q_3+q_0q_1&\frac{1}{2}-(q_1^2+q_2^2) \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} \end{align*}

となる。この行列が論文内の式(10)である。4つの未知数に1つの拘束条件なので上の行列の自由度は3となる。3次空間内の回転の自由度は3なので矛盾はない。

まとめ

今回は、クォータニオンを用いた回転行列の表記について解説した。一番最後に示した行列が最近読んだ論文に出てきた式である。論文中では、ある目的を満たすように回転行列を最適化しているのであるが、クォータニオンを使うとその最適化が容易になるそうである。クォータニオンに関しては、3DCGをかじっていた頃に勉強したことがあるが、それ以来出会うことはなかった。今回改めて勉強し直したのでブログにまとめた。

参考文献

Kumada Seiya

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

最近の記事

  • 関連記事
  • おすすめ記事
  • 特集記事

アーカイブ

カテゴリー

PAGE TOP