はじめに
この記事は検証用です。直交行列を回転行列で表現します。回転行列とは、単位ベクトル$(x,y,z)$と角度$\theta$に対して次の形式で表現されるものです。
\begin{pmatrix}
\cos\theta + (1-\cos\theta)x^2 & (1-\cos\theta)xy - (\sin\theta)z & (1-\cos\theta)xz + (\sin\theta)y \\
(1-\cos\theta)xy + (\sin\theta)z & \cos\theta + (1-\cos\theta)y^2 & (1-\cos\theta)yz - (\sin\theta)x \\
(1-\cos\theta)xz - (\sin\theta)y & (1-\cos\theta)yz + (\sin\theta)x & \cos\theta + (1-\cos\theta)z^2
\end{pmatrix}
単位行列の場合は$(x,y,z)$に意味はなく、$\theta$は$0$で、そういう表現になります。たとえば軸が$(0,0,1)$の場合、これを列ベクトル$(a,b,c)$に左から掛けると、$(\cos\theta, \sin\theta, 0)$になるので、ちゃんと回転の表現になっています。
この記事の目的は、任意の直交行列をこの形式で表現することです。
直交行列の定義
実数成分の$3\times 3$行列は、自身の転置との積が単位行列の場合に直交行列と呼ばれます。
V = \begin{pmatrix}
a & b & c \\
d & e & f \\
g & h & i
\end{pmatrix}, ~~~~~VV^{t} = I
ただし
I = \begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{pmatrix}
なお、その定義の中に行列式が1であることは含まれていません。この記事では行列式が1のものを直交行列と呼んで議論することとします。個別の名前が欲しい所です...
\mathrm{det}(V)=1.
定義により、
a^2+b^2+c^2=1,~~~ d^2 + e^2 + f^2 = 1, ~~~g^2 + h^2+i^2=1,
ad+be+cf=0,~~~ag+bh+ci=0,~~~dg+eh+fi=0
が成立します。しかし、逆行列の性質により$(V^{t})V=I$もまた成立するゆえ、
a^2+d^2+g^2=1,~~~b^2+e^2+h^2=1,~~~c^2+f^2+i^2=1,
ab+de+gh=0,~~~bc+ef+hi=0,~~~ac+df+gi=0
も成立します。さらに、余因子行列の理論により、行列式が1であることから、
a=ei-fh,~~~b=fg-di,~~~c=dh-eg,
d=ch-bi,~~~e=ai-cg,~~~f=bg-ah,
g=bf-ce,~~~h=cd-af,~~~i=ae-bd
が成立します。行列式が1でないとこの一連の式は成立しないので注意してください。（この最後の式はベクトルの外積でも出すことができます。要は正規直交基底の話なので、2つのベクトルの外積でもう一つが出る、それをあらわしています。が、厳密性を重んじる立場からこのように述べています。）
以上の公式を用いて、回転行列の形式で表現されることを証明します。
トレースが3のとき単位行列になること
次の事実が成立します。
a+e+i=3 ~~\Leftrightarrow ~~ V=I.
すなわち直交行列において、それが単位行列であることと、トレースが3であることは同値になります。証明は非常にエレガントで、要するに単位行列の定義をダイレクトに使うだけですね。
\begin{align}
& ~~~~~(a-1)^2 + b^2+c^2+d^2+(e-1)^2+f^2+g^2+h^2+(i-1)^2 \\
&= 3+1+1+1-2a-2e-2i \\
&= 2(3-a-e-i).
\end{align}
平方の和に関する先ほどの式を使って変形するとこうなります。これは$0$以上ですが、$0$になることと$a+e+i$が$3$になることが同値です。そして$V$が単位行列になることとも同値です。ゆえに同値であり、さらに$V$が単位行列でない場合は$a+e+i<3$であることも分かりました。
各種変数の定義
以下、$V$は単位行列でないとし、議論を進めます。$\theta,x,y,z$は次のように定義します。まず、
p = \frac{a+e+i-1}{2} ~~(<1)
とおきます。不等式は仮定から従います。これは今後$\cos\theta$となるものです。このとき、
\begin{align}
& a+e+i = 2p+1, \\
& (a-p)+(e-p)+(i-p) = 1-p > 0
\end{align}
ですから、$a-p,~~e-p,~~i-p$のうち少なくともひとつは正でなければなりません。$a-p>0$としても一般性を失わないので、そのようにします。そこで、
x = \sqrt{\frac{a-p}{1-p}}~~ (>0)
と定義します。そして、
q = \frac{h-f}{2x}
とします。このとき、
p=\cos\theta,~~~q=\sin\theta
を満たすように$\theta$を定義します。また、
y=\frac{b+d}{2x(1-p)},~~~z=\frac{c+g}{2x(1-p)}
で定めます。これで$\theta,x,y,z$が決定されました。あとは式を確かめるだけです。
うそです。$p^2+q^2=1$を証明してないですね。これを示さないと、$\theta$が取れることを確かめることができないですね。
\begin{align}
1-p^2-q^2 &= 1-p^2-\frac{(h-f)^2}{4x^2} = 1-p^2-\frac{(h-f)^2(1-p)}{4(a-p)} \\
&= (1-p)((1+p) - \frac{(h-f)^2}{4(a-p)}) \\
&= \frac{1-p}{4(a-p)}\cdot (4(a-p)(1+p) - (h-f)^2).
\end{align}
ここで、
4(a-p)(1+p) = (a-e-i+1)(a+e+i+1) = (a+1)^2 - (e+i)^2
なので、$d^2+e^2+f^2=1,~~g^2+h^2+i^2=1,~~a^2+d^2+g^2=1$に注意して、
\begin{align}
(a+1)^2 - (e+i)^2 - (h-f)^2 &= a^2+2a+1 - e^2-2ei-i^2 - h^2+2hf-f^2 \\
&= a^2+d^2-1 + g^2-1+2a+1-2ei+2hf \\
&= 2(a-ei+hf) \\
&= 0.
\end{align}
最後は余因子の式を使いました。$p^2+q^2=1$なので、確かに$\theta$を取ることが出来るようです。良かったです。
よくないですね。$x^2+y^2+z^2=1$が証明できてないですね。これを示さないと、$(x,y,z)$が単位ベクトルになりませんね。
\begin{align}
1-x^2-y^2-z^2 &= 1-\frac{a-p}{1-p} - \frac{(b+d)^2}{4x^2(1-p)^2} - \frac{(c+g)^2}{4x^2(1-p)^2} \\
&= \frac{1-a}{1-p} - \frac{(b+d)^2+(c+g)^2}{4(a-p)(1-p)} \\
&= \frac{4(1-a)(a-p) - (b+d)^2 - (c+g)^2}{4(a-p)(1-p)}.
\end{align}
ここで分子について、$2(a-p) = a-e-i+1$なので、
\begin{align}
&~~~~~~ 4(1-a)(a-p) - (b+d)^2 - (c+g)^2 \\
&= 2(1-a)(a-e-i+1) - b^2-2bd-d^2-c^2-2cg-g^2 \\
&= 2(1-a^2) - 2(1-a)(e+i) - (b^2+c^2) - (d^2+g^2) - 2bd-2cg \\
&= -2e-2i+2ae+2ai-2bd-2cg \\
&= (ae-bd-i) + (ai-cg-e) \\
&= 0.
\end{align}
ふぅ。ちゃんと単位ベクトルになっていますね。よかったです。
証明
- $a=p+(1-p)x^2.$
定義から明らか。
- $e=p+(1-p)y^2.$
\begin{align}
(e-p)-(1-p)y^2 &= (e-p) - \frac{(b+d)^2}{4x^2(1-p)^2} = (e-p) - \frac{(b+d)^2}{4(a-p)} \\
&= \frac{4(a-p)(e-p) - (b+d)^2}{4(a-p)}.
\end{align}
ここで、
\begin{align}
4(a-p)(e-p) &= (a-e-i+1)(-a+e-i+1) = (1-i)^2 - (a-e)^2 \\
&= 1-2i+i^2-a^2+2ae-e^2.
\end{align}
ゆえに、$i^2+c^2+f^2=1$に注意して、
\begin{align}
4(a-p)(e-p) - (b+d)^2 &= 1-2i+i^2-a^2+2ae-e^2-b^2-2bd-d^2 \\
&= 2-2i-c^2-f^2-a^2+2ae-e^2-b^2-2bd-d^2 \\
&= -2i+2ae-2bd.
\end{align}
縦の行のノルム$1$と横の行のノルム$1$を両方使うわけです。で、最後のこれは余因子のところで紹介した
i=ae-bd
によりゼロになります。これで証明できました。
- $i=p+(1-p)z^2.$
$e$の証明とほぼ同じなので省略。
- $b = (1-p)xy- qz.$
\begin{align}
(b+qz)-(1-p)xy &= b + \frac{(h-f)(c+g)}{4x^2(1-p)} - \frac{b+d}{2} \\
&= \frac{b-d}{2} + \frac{(h-f)(c+g)}{4(a-p)} \\
&= \frac{2(a-p)(b-d) + (h-f)(c+g)}{4(a-p)}.
\end{align}
ここで$2(a-p) = a-e-i+1$なので、
\begin{align}
&~~~~~2(a-p)(b-d) - (h-f)(c+g) \\
&= (a-e-i+1)(b-d) + (h-f)(c+g) \\
&= ab-be-bi+b-ad+de+di-d +ch-cf+gh-fg \\
&= (ab+de+gh) - (ad+be+cf) + b-d-bi+di +ch-fg \\
&= (b-fg+di) - (d-ch+bi) \\
&= 0.
\end{align}
内積がゼロになったり余因子がゼロになったりしています。大変だ。こんな感じですが、ともあれ証明はできています。
- $d=(1-p)xy+qz.$
これは先ほどの結果と$y$の定義から明らかですね。$b+d=2xy(1-p)$なので。
あと$c,f,g,h$が残っていますが、$c,g$は$z$の定義から片方だけ調べればOKですね。$f,h$に関しても、$q$の定義から片方だけでOKです。$c,g$の証明はさっきと大体一緒なので、あとは$f$だけ説明することとします。
- $f=(1-p)yz-qx.$
\begin{align}
f+qx - (1-p)yz &= f+ \frac{h-f}{2} - \frac{(b+d)(c+g)}{4x^2(1-p)} \\
&= \frac{h+f}{2} - \frac{(b+d)(c+g)}{4(a-p)} \\
&= \frac{2(a-p)(h+f) - (b+d)(c+g)}{4(a-p)}.
\end{align}
ここで$2(a-p) = a-e-i+1$なので、
\begin{align}
&~~~~~2(a-p)(h+f) - (b+d)(c+g) \\
&= (a-e-i+1)(h+f) - (b+d)(c+g) \\
&= ah-eh-ih+h+af-ef-if+f - bc-cd-bg-dg \\
&= (ah-bg+f) + (af-cd+h) - (eh+if+dg) - (bc+ef+ih) \\
&= 0.
\end{align}
わぁすごい。手品みたい。そういうわけで成立します。これで証明終わりとします。
回転行列の導出について
いわゆるロドリゲスの回転公式ですが、全く言及しないのもあれなので、一応述べておきます。まず単位ベクトル$v=(x,y,z)$を用意します。ベクトル$p$に対し、これを$v$の周りに$\theta$だけ回転させることを考えます。具体的には$v$を上から見るとして、これが$z$軸正方向だとした場合に、$x$軸を最短で$y$軸に重ねる、その向きに回転させます。一言で言うと「正の向き」です。時計を持ち出すと余計な突っ込みをされるので敢えて厳密に説明しました。
p = (p\cdot v)v + (p-(p\cdot v)v),~~~p' = p - (p\cdot v)v.
第一項はいわゆる$v$成分です。これは軸に平行なので不変です。第二項を$p'$とおきます。これが回転します。これは$v$に垂直です。$v,p'$の順で外積を取ると、$p'$と同じ大きさのベクトルを用意でき、これを使えば線形結合で目的のベクトルを出すことができます。
p'' = (\cos\theta)p' + (\sin\theta)(v\times p').
これとさっきの$(p\cdot v)v$を足すと回転の結果が得られます（$q$とおく）。
\begin{align}
q &= (p\cdot v)v + (\cos\theta)p' + (\sin\theta)(v\times p') \\
&= (p\cdot v)v + (\cos\theta)(p-(p\cdot v)v) + (\sin\theta)(v\times (p-(p\cdot v)v)) \\
&= (\cos\theta)p + (1-\cos\theta)(p\cdot v)v + (\sin\theta)(v\times p).
\end{align}
実際に$p=(a,b,c)$として適用すると、これを列ベクトルとして、一番最初に用意した行列を左から掛けた場合の結果と一致します。
一般性について
$a-p>0$としても一般性を失わない、としましたが、要は$a,e,i$のうち$a$が最大としているようなものです。他の成分が最大の場合であっても同じような議論になります。実装の際には、まさに最大の対角成分を使って回転行列の形式を割り出しました。それも大小比較なので、負荷は無いに等しいです。
回転行列が直交行列になることについて
直交行列が回転行列形式になることは証明できましたが、回転行列そのものの形に対して、それが直交行列であることを証明しないと片手落ちなので、一応補足の意味で計算します。証明することは、角度$\theta$と単位ベクトル$(x,y,z)$からできる行列：
V=\begin{pmatrix}
\cos\theta + (1-\cos\theta)x^2 & (1-\cos\theta)xy - (\sin\theta)z & (1-\cos\theta)xz + (\sin\theta)y \\
(1-\cos\theta)xy + (\sin\theta)z & \cos\theta + (1-\cos\theta)y^2 & (1-\cos\theta)yz - (\sin\theta)x \\
(1-\cos\theta)xz - (\sin\theta)y & (1-\cos\theta)yz + (\sin\theta)x & \cos\theta + (1-\cos\theta)z^2
\end{pmatrix}
について、その行列式が1であることと、$VV^t=I$（単位行列）が成立することです。
- 行列式が1になること
各成分を次のようにおきます。
V=\begin{pmatrix}
a & b & c \\ d & e & f \\ g & h & i \end{pmatrix}
察してください。で、たとえば次のようになります。
\begin{align}
& ~~~~~ae-bd \\ &= \cos^2\theta + (1-\cos\theta)^2x^2y^2 + \cos\theta(1-\cos\theta)(x^2+y^2) - (1-\cos\theta)^2x^2y^2 + (\sin^2\theta)z^2 \\
&= \cos^2\theta + \cos\theta(1-\cos\theta)(1-z^2) + (\sin^2\theta)z^2 \\
&= \cos\theta + (1-\cos\theta)z^2 \\
&= i.
\end{align}
こんな感じで、$dh-ge = c, ~~bg-ah = f$を示したうえで、行列式の定義（3列目による展開）を適用すると、
\begin{align}
&~~~~~c^2+f^2+i^2 \\
&= \sin^2\theta(x^2+y^2) + \cos^2\theta + (1-\cos\theta)^2(x^2+y^2+z^2)z^2 +2\cos\theta(1-\cos\theta)z^2 \\
&= \sin^2\theta(1-z^2) + \cos^2\theta + (1-2\cos\theta+\cos^2\theta)z^2 + 2\cos\theta(1-\cos\theta)z^2 \\
&= (\sin^2\theta+\cos^2\theta) + z^2(\cos^2\theta-1+1-2\cos\theta+\cos^2\theta+2\cos\theta-2\cos^2\theta) \\
&= 1.
\end{align}
ゆえに、行列式は1になります。
- $VV^t=I$となること
$V$を次のように分解します。
V = (\cos\theta)I + (1-\cos\theta)ww^t + (\sin\theta)J.
ただし、$w=(x,y,z)$であるとし（列ベクトル）、
J = \begin{pmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{pmatrix}
です。交代行列ですね。$ww^t$はいわゆるダイアド積です。ベクトルの積で出来る三次正方行列です。
これを切り崩すには$JJ^t$を計算すればいいです。
JJ^t = \begin{pmatrix}
y^2+z^2 & -xy & -xz \\ -xy & x^2+z^2 & -yz \\ -xz & -yz & x^2+y^2
\end{pmatrix} = I - ww^t.
単位ベクトルの性質を使いました。$J^t=-J$（交代行列）なので、
ww^t = I-JJ^t = I+J^2
であり、これを代入すれば、
V = I + (\sin\theta)J + (1-\cos\theta)J^2
と書けます。ゆえに
V^t = I - (\sin\theta)J + (1-\cos\theta)J^2
であり、両者の積は、
\begin{align}
&~~~~VV^t \\
&= (I + (\sin\theta)J + (1-\cos\theta)J^2)(I - (\sin\theta)J + (1-\cos\theta)J^2) \\
&= I + (1-\cos\theta)^2J^4 +2(1-\cos\theta)J^2 - (\sin\theta)^2J^2 \\
&= I + (1-\cos\theta)^2J^4 + (1-2\cos\theta+1-\sin^2\theta)J^2 \\
&= I + (1-\cos\theta)^2(J^4 + J^2)
\end{align}
となります。行列代数というやつで、単位行列とそうでない一種類の正方行列のべき和は多項式のように扱うことができ、すべての積が可換となります。とても便利ですね。最後は
J^4+J^2=O
（ゼロ行列）を示すだけですが、計算で、
-J^2 = JJ^t = \begin{pmatrix}
1-x^2 & -xy & -xz \\ -xy & 1-y^2 & -yz \\ -xz & -yz & 1-z^2
\end{pmatrix},~~~~~(-J^2)J = \begin{pmatrix}
0 & -z & y \\ z & 0 & -x \\ -y & x & 0
\end{pmatrix} = J
が分かるので、従います。証明終わり。
（テイラー展開みたいだなって思うでしょう。実はその通りで、$\theta$が十分小さいとき、これは$\theta$について二次の項までのテイラー展開のように見えます。さしずめ$J$は導関数といったところでしょうか。すみません、わかるのはここまでなのであとは自分で勉強してみてください）
（追記：調べてたらありました。リー代数の話で、これですね...
CV・CG・ロボティクスのためのリー群・リー代数入門: (5) 回転ベクトル
テイラー展開の変形で出るんですね...すごい。いずれきちんと読んでみたいと思います。なるほど...なるほど...）
おわりに
見ての通り、軸にしても角度にしても割と一瞬で出るので、計算速度も割と速いです。実際にカメラのビューの補間をこのメソッドで実装したところ爆速でした。ビューの補間まででいいならクォータニオンは要らないわけです。あれが必要とされるのはもっと別の文脈なのでしょうが、知識とか色々足りないのでわかりません。機会があったらもうちょっと勉強してみようと思います。
そもそも運用の段階ではこんな面倒な計算も、検証も不要なので、速いのはある意味当然ですね...
ここまでお読みいただいてありがとうございました。（不備がありましたらコメントでよろしくお願いします）