はじめに

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

扱う行列の定義

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

性質

　$SO(3)$の行列には特徴がいくつかあります。固有値のひとつは必ず$1$であること、残りの2つは互いに共役な絶対値が1の複素数であること、です。

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

軸と角度による表現

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

V = \begin{pmatrix}
\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
\end{pmatrix}

なお、$\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 = \begin{pmatrix}
a_{00} & a_{10} & a_{20} \\
a_{01} & a_{11} & a_{21} \\
a_{02} & a_{12} & a_{22}
\end{pmatrix}

対角成分の和を考えます。いわゆるトレースですが、これは$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$までの範囲の値を取れば、常に短い方で補間できます。その辺りは自由に選ぶことができます。

計算による証明の手引き

　面倒なので具体的な成分比較による証明はしませんが、コツを述べておきます。余因子行列の概念を使うと、行列式が1なので、各行、各列のノルムが1だとか直交するだとかで得られる式とは別に、たとえば

  a_{00}a_{11} - a_{01}a_{10} = a_{22}

のような一連の式が得られます。そういうのを駆使すれば、頑張れば、できます。
　それと、証明の際は添え字ではなくアルファベットの方が書きやすいです。ではなぜこの記事で添え字をつかっているかというと、これがMathlogに投稿した記事のコピペを編集したものだからです。
　

一般の場合

　一般に$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$が返ります。

おわりに

　その後、p5.js/1.7.0は無事にリリースされ、カメラの補間の関数もめでたく実装されました。良かったです。
　ここまでお読みいただいてありがとうございました。