3次直交行列の補間を行う話

p5.jsというライブラリの関数を作る一環で、3次直交行列の補間が必要になったので、その際の計算について触れたいと思います。Mathlogを触るのが久しぶりなので、うまくできるかわかりませんが...

扱う行列の定義

　ここで扱うのは行列式が1である3次の実直交行列です。いわゆる$SO(3)$と呼ばれるクラスになります。このクラスの行列の補間を行います。以下、単に$SO(3)$と書いてこのような行列の全体とします。いちいち実であることに言及するのが面倒なので、そのようにします。

性質

　$SO(3)$の行列には特徴がいくつかあります。固有値のひとつは必ず1であること、残りの2つは互いに共役な絶対値が1の複素数であること、です。
　この行列は線形変換としてどのような処理を実行するのかというと、実は球面の回転として表現することができます。たとえば$SO(2)$の行列はすべて原点中心の回転（もしくは恒等写像）として表現できます。$SO(3)$の行列は、単位行列でない限り、なんらかの軸となるベクトルを持っています。加えて角度を持っており、軸の周りにその角度だけ回転させる作用となっています。軸をどちらから見るかにより、回転の方向に2通りの任意性があるので、具体的に確定するわけではなく、いずれかを選ばなければなりません。

軸と角度による表現

　軸となる単位ベクトルを$(a,b,c)$とし、角度を$\theta$とします。このときの$SO(3)$の行列の表現は次のようになります。もちろん、単位行列ではないことを前提とします。軸が定まらないからです。

$$V = \threeMat{\cos\theta + (1-\cos\theta)a^2}{(1-\cos\theta)ab-(\sin\theta) c}{(1-\cos\theta)ca+(\sin\theta) b}{(1-\cos\theta)ab+(\sin\theta) c}{\cos\theta + (1-\cos\theta)b^2}{(1-\cos\theta)bc-(\sin\theta) a}{(1-\cos\theta)ca-(\sin\theta) b}{(1-\cos\theta)bc+(\sin\theta) a}{\cos\theta + (1-\cos\theta)c^2}$$
なお、$\theta$が0の場合、$a,b,c$の値によらずこれは単位行列となります。このように軸によらないので、一意的な表現がありません。
　作用を確かめると、まず$(a,b,c)$が$(0,0,1)$の場合、$(1,0,0)$を列ベクトルとみなして行列を左から掛けると、$(\cos\theta, \sin\theta, 0)$になります。確かに、z軸を上から見た場合の$\theta$の回転となっていますね。

軸と角度の算出

　軸と角度をいかにして求めるかについて説明します。$SO(3)$の行列$V$を取ってきます。単位行列ではないとし、成分を次のように書きます。
$$V = \threeMat{a_{00}}{a_{10}}{a_{20}}{a_{01}}{a_{11}}{a_{21}}{a_{02}}{a_{12}}{a_{22}}$$
対角成分の和を考えます。いわゆるトレースですが、これは$SO(3)$の行列の性質により、$1+2\cos\theta$になります。$\theta$は上記の表現における角度ですが、$\pm$の任意性があります。単位行列ではないため、常に$\cos\theta$は1より小さいです。固有値の内訳が1と、互いに共役な絶対値が1の複素数2つなので、実部を考えるとそのようになります。
先にこれにより$\cos\theta$を出しておきます。
$$\cos\theta = \frac{a_{00}+a_{11}+a_{22}-1}{2}.$$
すると、$\cos\theta<1$なので、次の式が成立します。
$$a_{00} + a_{11} + a_{22} =1 + 2\cos\theta ,$$
$$(a_{00} - \cos\theta) + (a_{11}-\cos\theta) + (a_{22}-\cos\theta) = 1 - \cos\theta > 0.$$
これより、いずれかの$a_{ii}-\cos\theta$は正であることがわかります。$a_{00}-\cos\theta$が正であるとしても一般性を失いません。このとき、目標の行列の左上の成分から、
$$a^2 = \frac{a_{00}-\cos\theta}{1-\cos\theta},~~~~a = \sqrt{\frac{a_{00}-\cos\theta}{1-\cos\theta}}$$
として$a$が定まります。これは二つありますが、正の方を取ります。これが軸ベクトルの任意性で、常に2通りの取り方が存在します。$b,c$が成分の和を比較することで次のように求まります。
$$a_{10}+a_{01} = 2ab(1-\cos\theta),~~~b=\frac{a_{10}+a_{01}}{2a(1-\cos\theta)}.$$
$$a_{20}+a_{02} = 2ca(1-\cos\theta),~~~c=\frac{a_{20}+a_{02}}{2a(1-\cos\theta)}.$$
また、$\sin\theta$も定まります。
$$a_{12} - a_{21} = 2a\sin\theta,~~~\sin\theta = \frac{a_{12}-a_{21}}{2a}.$$
$\cos\theta$と$\sin\theta$から具体的な角度が定まります。これで上記の計算をすれば、元の行列が復元されます。
　$a_{11}-\cos\theta$が正の場合、$a_{22}-\cos\theta$が正の場合も同様です。

補間するには

　上記の方法で角度を算出したうえで、角度に適当な0～1の係数を掛ければ、0のときは単位行列、1のときは元の行列ですから、0から1まで動かすと元の行列に滑らかに変化していきます。これで補間することができました。なお$\theta$を取る際に$-\pi$から$\pi$までの範囲の値を取れば、常に短い方で補間できます。その辺りは自由に選ぶことができます。

一般の場合

　一般に$SO(3)$の行列$V_0,V_1$を取ってこれらを補間するにはどうするのかというと、片方の転置（というか逆行列）を取って掛け算します。
$$V=V_1 {V_0}^t$$
とおきます。この$V$はもちろん$SO(3)$の行列です。これを先ほどのように補間して単位行列から元の行列まで動かします。このとき、$VV_0$を考えれば、これは$V$が単位行列から$V$まで動く間、$V_0$から$V_1$まで動きます。これで滑らかに補間することができました。

補足

　ここで計算しているのはクォータニオンの成分です。が、自分はこれについては詳しくないです。クォータニオンのまま扱うと独特の積の法則などがあり割とややこしいので、直交行列のまま扱いたいと思ってこのような方法を採用しました。自分が開発のお手伝いをしているp5.jsというライブラリがカメラの取り扱いにおいてクォータニオンを採用していないので、このような方法を考案することになりました。Three.jsのこれが似たような計算をしているかもしれません。
　 setFromRotationMatrix
ただこのやり方だと$\sin\theta$が0のとき使えないので、成分のまま計算するには不向きなように感じました。それでこうしました。

実装の実際

　実際にカメラの補間に用いた例がこちらになります：
　 p5.Camera.set()_DEMO_2
この中でカメラの行列を補間するのに実際にこれを用いています。$a_{ii}-\cos\theta$が正になるものを取るくだりでは、確実に正になるよう、最も大きい対角成分を取ることで対処しています。角度の算出にはMath.atan2()を用いています。なので$-1$に対しては常に$\pi$が返ります。
　以上です。ここまでお読みいただいてありがとうございました。