ゲームやCGのための数学 Q -物体の回転について-

「オイラー回転は神の回転」「なにっ！？複素数？そんなもので回転が表現できるものか！貴様ーッ！！！オイラー先生を愚弄するかーッ！今すぐ取り消せーッ！」「し・・・四元数だと・・・？お、お前、変なクスリでもやってるのか？」

はい、すんません、このタイトルやりたかっただけです。
いつも通り、ちゃんとやりますんで許して

はじめに

この記事では物体の回転。それも３次元の物体の回転を取り扱います。これが2次元の話ならクッソ簡単なんですけどねえ…。
で、「物体の回転」と言いましたが、細かく言うとベクトルの回転…もしくは「原点を中心とした座標の回転」を取り扱います。

そして最終到達地点では「クオータニオンによる回転」の解説をするわけですが、CGプログラミング初心者はクオータニオンを道具として使いつつも「何であんなけったいな概念使わなあかんねん！」「そもそも計算が異次元なんじゃい！」「何であんな計算になるねん！」「マジで大丈夫なんか？実は騙されとるんちゃうんか？」と疑念を持ちつつ使っているかもしれません。

この記事ではそういう疑念を取っ払う一助になるかなと思って書いています。

多少、代数系の話が出てきますが、一応虚数とかが分かってれば大丈夫だと思います…たぶん。あと行列がしれっと出てくるので、分かんねー人はここでブラウザバックして行列はどっかで概念だけでも身につけてからこの記事を読み直してください。いいねとかストックはその後で構いません(笑)

それでは解説始めます。

オイラー回転

ゲームエンジン等で3D物体の回転を行うにはオイラー回転(ロールヨーピッチ)と、クオータニオンによる回転があります。
オイラー回転は分かりやすいですよね？物体の姿勢をX軸中心回転＊Y軸中心回転＊Z軸中心回転で表現するものですね。これさえあれば十分な気がします。

ご覧のようにそれぞれの軸で回転を組み合わせればどんな回転も表現できそうではありますし、最初に配置しておく静止物体ならそれで十分だと思います。あと地上の車や戦車なども大丈夫かなとは思います。でももうちょっと自由度の高い飛行機やら人体の関節やらが絡んでくるとちょ～～っと厳しくなってきます。それについてはもう少し後でお話しします。

オイラー回転の行列

仮にX軸中心にθ回転、Y軸中心にφ回転、Z軸中心にψ回転するとします。そうすると

R_x=\begin{pmatrix} 
1 & 0 & 0 \\ 
0 & cosθ & -sinθ \\
0 & sinθ & cosθ 
\end{pmatrix}
 , 
R_y=\begin{pmatrix} 
cosφ & 0 & sinφ \\ 
0 & 1 & 0 \\
-sinφ & 0 & cosφ 
\end{pmatrix}
 , 
 R_z=\begin{pmatrix} 
cosψ & -sinψ & 0 \\ 
sinψ & cosψ & 0 \\
0 & 0 & 1 
\end{pmatrix}

オイラー回転行列の乗算による操作の結合

そして、この回転行列というのは「回転操作」を意味しており、乗算することで回転操作を連続して行ったのと同じ結果が得られます。
というわけで３つ乗算してみますが、見た目がゴチャゴチャするので
$$A=cosθ , B=sinθ , C=cosφ , D=sinφ , E=cosψ , F=sinψ$$
としておきます。そうすると

R_xR_yR_z=\begin{pmatrix} 
CE & -CF & D \\  
BDE+AF & -BDF+AE & -BC \\
-ADE+BF & ADF+BE & AC
\end{pmatrix}

う～ん、ややこしい。一応展開してみます。

R_xR_yR_z=\begin{pmatrix} 
cosφcosψ & -cosφsinψ & sinφ \\  
sinθsinφcosψ+cosθsinψ & -sinθsinφsinψ+cosθcosψ & -sinθcosφ \\
-cosθsinφcosψ+sinθsinψ & cosθsinφsinψ+sinθcosψ & cosθcosφ
\end{pmatrix}

…やってられませんね。ひとまず「オイラー角でもまとめると超ヤヤコシイ形」というのがビジュアルで分かってればいいです。

ところで回転行列は転置行列を乗算すると「単位行列」になります

いきなり別の話をしているように感じると思いますが、これも回転においては大事な性質なので見ていきましょう
まず2Dから見てみましょうかね。２Dの回転行列は

R=\begin{pmatrix}
cosθ&-sinθ\\
sinθ&cosθ
\end{pmatrix}

で、これを転置(行列の縦と横を入れ替える)すると

R^T=\begin{pmatrix}
cosθ&sinθ\\
-sinθ&cosθ
\end{pmatrix}

この二つを乗算すると

RR^T=\begin{pmatrix}
cos^2θ+sin^2θ && cosθsinθ-sinθcosθ\\
sinθcosθ-cosθsinθ && sin^2θ+cos^2θ
\end{pmatrix}=
\begin{pmatrix}
1&&0\\
0&&1
\end{pmatrix}

になり、単位行列であることが分かりますね？

では3Dでも成り立つことを証明するために先ほどの3Dオイラー回転行列と、その転置行列を乗算してみます。見た目が大変な事になるためABCDEFで記述します。
ただし三角関数の性質より$$A^2+B^2=1 , C^2+D^2=1 , E^2+F^2=1$$とします。
ここで先ほどの

R_xR_yR_z=\begin{pmatrix} 
CE & -CF & D \\  
BDE+AF & -BDF+AE & -BC \\
-ADE+BF & ADF+BE & AC
\end{pmatrix}

と

(R_xR_yR_z)^T=\begin{pmatrix} 
CE & BDE+AF & -ADE+BF \\  
-CF & -BDF+AE & ADF+BE \\
D & -BC & AC
\end{pmatrix}

を乗算します(ところで転置しても左上から右下にかけては変化がなくて面白いですね)
一つ一つの計算が結構大変なので
$$R'=(R_xR_yR_z)(R_xR_yR_z)^T$$
として、それぞれの行,列をR₁₁ , R₁₂ , R₁₃ , R₂₁ …のように一つ一つ計算していきますね。
では記念すべきR₁₁から

\begin{align}
R_{11}&=(CE)^2+(CF)^2+D^2 \\
 &= C^2E^2 + C^2F^2+D^2 \\
 &= C^2(E^2 + F^2)+D^2 \\
 &= C^2 +D^2 = 1
\end{align}

おっ、幸先良いですね。この調子で残り8個やっつけちゃいましょう

\begin{align}
R_{12}&=(CE)(BDE+AF)-CF(-BDF+AE)-DBC \\
&=BCDE^2+ACEF+BCDF^2-ACEF-DBC \\
&=BCD(E^2+F^2)-DBC\\
&=BCD-DBC=0
\end{align}

\begin{align}
R_{13}&=CE(-ADE+BF)-CF(ADF+BE)+DAC \\
&=-ACDE^2+BCEF-ACDF^2-BCEF+DAC \\
&=-ACD(E^2+F^2)+BCEF-BCEF+DAC\\
&=-ACD+DAC=0
\end{align}

\begin{align}
R_{21}&=(BDE+AF)CE-(-BDF+AE)CF-BCD \\
 &= BCDE^2+ACEF+BCDF^2-ACEF-BCD \\
 &= BCD(E^2+F^2)+ACEF-ACEF-BCD \\
 &= BCD-BCD = 0
\end{align}

\begin{align}
R_{22}&=(BDE+AF)^2+(-BDF+AE)^2+B^2C^2 \\
&=B^2D^2E^2+2ABDEF+A^2F^2+B^2D^2F^2-2ABDEF+A^2E^2+B^2C^2 \\
&=B^2D^2(E^2+F^2)+2ABDEF-2ABDEF+A^2(F^2+E^2)+B^2C^2\\
&=B^2D^2+A^2F^2+A^2E^2+B^2C^2\\
&=B^2(D^2+C^2)+A^2(F^2+E^2)\\
&=B^2+A^2=1
\end{align}

\begin{align}
R_{23}&=(BFE+AF)(-ADE+BF)+(-BDF+AE)(ADF+BE)-ABC^2 \\
&=-ABDE^2F+B^2DEF-A^2DEF+ABF^2-ABD^2F^2-B^2DEF+ABE^2+A^2DEF-ABC^2 \\
&=-ABD^2E^2-ABD^2F^2+AB(F^2+E^2)-ABC^2\\
&=-ABD^2(E^2+F^2)+AB-ABC^2\\
&=-AB(D^2+C^2)+AB=0
\end{align}

\begin{align}
R_{31}&=-ACDE^2+BCEF-ACDF^2-BCEF+ACD \\
 &= -ACD(E^2+F^2)+ACD=0 
\end{align}

\begin{align}
R_{32}&=-ABD^2E^2-A^2DEF+B^2DEF+ABF^2-ABD^2F^2+A^2Def-B^2DEF+ABE^2-ABC^2 \\
&=-ABD^2(E^2+F^2)+AB(F^2+E^2)-ABC^2 \\
&=-ABD^2-ABC^2+AB\\
&=-AB(D^2+C^2)+AB=0\\
\end{align}

\begin{align}
R_{33}&=A^2D^2E^2-2ABDEF+B^2F^2+A^2D^2F^2+2ABDEF+B^2E^2+A^2C^2 \\
&=A^2D^2(E^2+F^2)+B^2(F^2+E^2)+A^2C^2 \\
&=A^2D^2+A^2C^2+B^2\\
&=A^2(D^2+C^2)+B^2\\
&=A^2+B^2=1
\end{align}

ハァ…ハァ…どうよ！？
「計算が疑わしいなぁ」と思った方は、ご自分で紙の上で計算してみてください。紙と鉛筆でやるのをお勧めします。相殺して消えるやつとか１にまとまるやつが多いので、鉛筆でシュッシュッと消していくのは快感です。

ともかく、複雑なオイラー回転でも要素が回転だけならば転置行列と乗算すれば単位行列になることが分かったと思います。
これは逆に言うと転置を乗算した結果単位行列になるのなら「オイラー回転のみの組み合わせ」で表せるということです。ちょっと覚えておいてください。後々の説明に使用します。

ところでオイラー回転による姿勢制御を図にするとこんな感じです。

X軸回転とY軸回転とZ軸回転が組み合わさって自由な姿勢を制御するわけです。ここでオイラー回転では致命的な弱点が出てきます。

図のような状態になった時、順番にもよりますがX-YもしくはY-ZもしくはZ-Xの軸が一致してしまい、この図の場合だと正面中心のロール、つまりZ軸回転をしようとしたときに「ギンッ！」となってしまい回転できなくなってしまいます。

これを「ジンバルロック」と言います。ジンバルロックで検索すると動画がたくさん出てくるので、分かりづらい人はそちらを参照してください。

さぁ、回転に対して万能であるかに見えたオイラー回転！どうしましょうか？
クォータニオンの出番か？

なぜ座標(ベクトル)の回転にクォータニオンが使用されてるん？クオータニオン以外に解決策はないん？

実際、UnityやUnrealEngineにおいて静的なオブジェクトの配置にはオイラー回転が使用されていますが、動的な姿勢制御にはクォータニオンが使用されています。なんでそんなけったいなものを使っているのか？

これ「オイラー回転だとジンバルロックを起こすから、ジンバルロックを起こさないクオータニオンが使われているんDA！」という説明をする人がいますが、う～ん、７０点。間違ってはないんですが、それだけが目的なら別にクオータニオンである必要はありません。要は「任意軸周りの回転」ができればいいわけです。

なので、一足飛びにクオータニオンやっちゃうと、理解ができないブラックボックスのまま使用することになり、ちょっと不安ですよね。なのでまずはクオータニオンより先に「ロドリゲスの方法」で、任意軸周りの回転をやってみましょう。

ロドリゲスの方法

任意軸回転を考える時に、ロドリゲスの公式というのがありますが、ひとまずそいつを知らないふりをして、自分で導出してみようと思います。
任意軸の回転というのは具体的には、特定の点Pを特定のベクトルVを中心に回転するという事です。ちなみにこの点Pも原点からベクトルを伸ばせばそれすなわちベクトルPとも言えます

ひとまず図のように、点PとそれをベクトルVを中心にθ回転したものをP'として考えます。XYZ軸のことは考えません(現段階では邪魔なだけです)
求めたいのはこの点P’なんですが、どうやって求めましょうかねえ？はい、もうちょっと補助線や情報を書いてみましょう

そうするとV'の長さはdですが、このdはPの射影と言えますので

\displaylines{
d=\hat{V} \cdot \vec{P} \\
\vec{V'}=d\hat{V}
}

そうすれば、Vの直線に対するPからの垂線ベクトルは$$\vec{V'P}$$と言えますね。
また、操作が回転であるという事は

\|\vec{V'P'}\|=\|\vec{V'P}\|

なので、P'からV'Pに垂線を下した座標をQとするとここに直角三角形△V'P'Qができるわけですから

\displaylines{
\|\vec{V'Q}\|=\|\vec{V'P}\|cosθ \\ 
\|\vec{QP'}\|=\|\vec{V'P}\|sinθ \\ 
}

となります。長さがわかったものの方向が分からないですね。ここで外積を使います。
ベクトルQ'Pは△OV'Pで作られる面と直交しているわけですから、つまりその面の法線ベクトルの定数倍という事になります。

\vec{QP'}=\frac{\vec{V'} \times \vec{P}}{\|\vec{V'} \times \vec{P}\|}\|\vec{QP'}\|

ここでQP'の大きさはV'Psinθなので

\vec{QP'}=\frac{\vec{V'} \times \vec{P}}{\|\vec{V'} \times \vec{P}\|}\|\vec{V'P}\|sinθ

なお、

\frac{\vec{V'} \times \vec{P}}{\|\vec{V'} \times \vec{P}\|}=\frac{\hat{V} \times \vec{P}}{\|\hat{V} \times \vec{P}\|}

です。外積の部分は約分できるため、Vを正規化Vにしておきます。納得できない人はとりあえず手計算して試してください。
ここで、|V'P|はベクトルＰに対する|P|sinφなので

\|\vec{V'P}\|=\|\vec{P}\|sinθ

また、外積とsinの性質である

\|\vec{V}\times\vec{P}\|=\|\vec{V}\|\|\vec{P}\|sinθ

より

\|\vec{V'P}\|=\|\vec{P}\|sinθ=\frac{\vec{V}\times\vec{P}}{\|\vec{V}\|}=\|\hat{V}\times\vec{P}\|

となりますので、前述のQPの式は

\vec{QP'}=\frac{\hat{V'} \times \vec{P}}{\|\hat{V'} \times \vec{P}\|}\|\hat{V}\times\vec{P}\|sinθ=(\hat{V'} \times \vec{P})sinθ=(\hat{V} \times \vec{P})sinθ

となります。最終的に知りたいP'の座標(ベクトル)は

\vec{P'}=\vec{V'}+\vec{V'Q}+\vec{QP'}

であるため、それぞれをVとPとθで表すと

\begin{align}
\vec{P'}&=\hat{V}(\hat{V}\cdot\vec{P})+\vec{V'P}cosθ+(\hat{V} \times \vec{P})sinθ \\
&=\hat{V}(\hat{V}\cdot\vec{P})+(\vec{P}-\hat{V}(\hat{V}\cdot\vec{P}))cosθ+(\hat{V} \times \vec{P})sinθ \\
&=\vec{P}cosθ+\hat{V}(\hat{V}\cdot\vec{P})(1-cosθ)+(\hat{V} \times \vec{P})sinθ
\end{align}

これがロドリゲスの式です。
次にこれを行列化しましょう

ロドリゲスの方法の行列化

なお、内積と外積の行列化について最初に述べておこうと思います。内積は簡単ですね。そもそも行列の演算自体が内積の集合体みたいなもんなんですから。例えば、$V'=\hat{V}(\hat{V}\cdot{P}$)を計算したいとすると

\begin{align}
(V'_x,V'_y,V'_z)&=(V_xP_x+V_yP_y+V_zP_z)(V_x,V_y,V_z)\\
&=(V_x(V_xP_x+V_yP_y+V_zP_z),V_y(V_xP_x+V_yP_y+V_zP_z),V_z(V_xP_x+V_yP_y+V_zP_z))
\end{align}

となり、行列は乗算結果の1行目がX、2行目がY、3行目がZに割り当たるため

(P_x,P_y,P_z)
\begin{pmatrix}
V_x^2 & V_xV_y & V_xV_z \\
V_xV_y & V_y^2 & V_yV_z \\
V_xV_z & V_yV_z & V_z^2
\end{pmatrix}
=(V'_x.V'_y,V'_z)

といった具合になります。また、外積$\hat{V}\times\vec{P}$を計算したいとすると

(V_yP_z-V_zP_y , V_zP_x-V_xP_z , V_xP_y-V_yP_x)

ですから、

(P_x,P_y,P_z)
\begin{pmatrix}
0 & V_z & -V_y \\
-V_z & 0 & V_x \\
V_y & -V_x & 0
\end{pmatrix}
=(V_yP_z-V_zP_y , V_zP_x-V_xP_z , V_xP_y-V_yP_x)

になりますね？では、これらの武器を使って、ロドリゲスの方法を行列化しましょう。
まず、Pをそのまま出力したい時は単位行列を使用します。ここで単位行列をEとして、T=$\hat{V}\times{P}$として、$D=\hat{V}(\hat{V}\cdot\vec{P})$とします(アルファベットに特に意味はありません)

\displaylines{
E=\begin{pmatrix}
1&0&0\\
0&1&0\\
0&0&1\end{pmatrix}\\

T=\begin{pmatrix}
0 & V_z & -V_y \\
-V_z & 0 & V_x \\
V_y & -V_x & 0
\end{pmatrix}\\

D=\begin{pmatrix}
V_x^2 & V_xV_y & V_xV_z \\
V_xV_y & V_y^2 & V_yV_z \\
V_xV_z & V_yV_z & V_z^2
\end{pmatrix}
}

はい、これを用いて行列化してみます。Pは乗算対象だからこの式からは消えます

\vec{P}cosθ+
\hat{V}(\hat{V}\cdot\vec{P})(1-cosθ)+
(\hat{V} \times \vec{P})sinθ=
Ecosθ+D(1-cosθ)+Tsinθ

\begin{align}
R&=Ecosθ+D(1-cosθ)+Tsinθ\\
&=
\begin{pmatrix}
1&0&0\\
0&1&0\\
0&0&1\end{pmatrix}cosθ+
\begin{pmatrix}
V_x^2 & V_xV_y & V_xV_z \\
V_xV_y & V_y^2 & V_yV_z \\
V_xV_z & V_yV_z & V_z^2
\end{pmatrix}(1-cosθ)+
\begin{pmatrix}
0 & V_z & -V_y \\
-V_z & 0 & V_x \\
V_y & -V_x & 0
\end{pmatrix}sinθ \\
&=
\begin{pmatrix}
cosθ&0&0\\
0&cosθ&0\\
0&0&cosθ\end{pmatrix}+
\begin{pmatrix}
V_x^2(1-cosθ) & V_xV_y(1-cosθ) & V_xV_z(1-cosθ) \\
V_xV_y(1-cosθ) & V_y^2(1-cosθ) & V_yV_z(1-cosθ) \\
V_xV_z(1-cosθ) & V_yV_z(1-cosθ) & V_z^2(1-cosθ)
\end{pmatrix}+
\begin{pmatrix}
0 & V_zsinθ & -V_ysinθ \\
-V_zsinθ & 0 & V_xsinθ \\
V_ysinθ & -V_xsinθ & 0
\end{pmatrix}\\
&=
\begin{pmatrix}
cosθ+V_x^2(1-cosθ)& V_xV_y(1-cosθ) +  V_zsinθ &V_xV_z(1-cosθ)-V_ysinθ\\
V_xV_y(1-cosθ)-V_zsinθ&cosθ+V_y^2(1-cosθ)&V_yV_z(1-cosθ)+V_xsinθ\\
V_xV_z(1-cosθ)+V_ysinθ&V_yV_z(1-cosθ)-V_xsinθ &cosθ+ V_z^2(1-cosθ)\end{pmatrix}
\end{align}

うわぁ～～～！気が狂いそうなくらいめんどくせえ～！！でもこの後もっと面倒な事をやります。
これがオイラー回転行列と同等のものであることを証明します

ロドリゲスの行列と逆行列を乗算すると…？

さぁ、クッソ面倒だけどやっていきましょう。

R=\begin{pmatrix}
cosθ+V_x^2(1-cosθ)& V_xV_y(1-cosθ) +  V_zsinθ &V_xV_z(1-cosθ)-V_ysinθ\\
V_xV_y(1-cosθ)-V_zsinθ&cosθ+V_y^2(1-cosθ)&V_yV_z(1-cosθ)+V_xsinθ\\
V_xV_z(1-cosθ)+V_ysinθ&V_yV_z(1-cosθ)-V_xsinθ &cosθ+ V_z^2(1-cosθ)\end{pmatrix}

として、$H=R\times R^T$を計算していきましょう。Hなのは特に理由はありません。なんとなくです。
また、この計算はややこしく長くなるため、簡略化のために$C=cosθ , S=sinθ , X=V_x , Y=V_y , Z=V_z$とします。
また$C^2+S^2=1 , X^2+Y^2+Z^2=1$から$S^2=1-C^2 , Y^2+Z^2=1-X^2$などを利用して計算をしていきますのでご了承ください。

H=\begin{pmatrix}
C+X^2(1-C)& XY(1-C)+ZS &XZ(1-C)-YS\\
XY(1-C)-ZS&C+Y^2(1-C)&YZ(1-C)+XS\\
XZ(1-C)+YS&YZ(1-C)-XS &C+ Z^2(1-C)\end{pmatrix}
\begin{pmatrix}
C+X^2(1-C)& XY(1-C)-ZS &XZ(1-C)+YS\\
XY(1-C)+ZS&C+Y^2(1-C)&YZ(1-C)-XS\\
XZ(1-C)-YS&YZ(1-C)+XS &C+ Z^2(1-C)\end{pmatrix}

はい、それではひとつずつ

\begin{align}
H_{11}&=(C+X^2(1-C))^2+(XY(1-C)+ZS)^2+(XZ(1-C)-YS)^2\\
&=C^2+2X^2C(1-C)+X^4(1-C)^2+X^2Y^2(1-C)^2-2XYZ(1-C)S+Z^2S^2+X^2Z^2(1-C)^2+2XYZ(1-C)S+Y^2S^2\\
&=C^2+2X^2C(1-C)+X^4(1-C)^2+X^2(Y^2+Z^2)(1-C)^2+(Y^2+Z^2)S^2\\
&=C^2+2X^2C(1-C)+X^4(1-C)^2+X^2(1-X^2)(1-C)^2+(1-X^2)S^2\\
&=C^2+2X^2C(1-C)+X^4(1-C)^2+X^2(1-C)^2-X^4(1-C)^2+S^2-X^2S^2\\
&=C^2+S^2+2X^2C(1-C)+X^2(1-C)^2-X^2S^2\\
&=1+2X^2C-2X^2C^2+X^2-2CX^2+C^2X^2-X^2S^2\\
&=1+2X^2C-2X^2C^2+X^2-2CX^2+C^2X^2-X^2(1-C^2)\\
&=1+2X^2C-2X^2C^2+X^2-2CX^2+C^2X^2-X^2+X^2C^2\\
&=1\\
\end{align}

はい、１つめだけは丁寧に計算過程を書きましたが、次以降はもうちょっと途中を簡略化して書いていきます。

\begin{align}
H_{12}&=(C+X^2(1-C))(XY(1-C)+ZS)+(XY(1-C)+ZS)(C+Y^2(1-C))+(XZ(1-C)+YS)(YZ(1-C)-XS)\\
&=XYC(1-C)+X^3Y(1-C)^2+ZCS+X^2ZS(1-C)+XYC(1-C)+XY^3(1-C)^2-ZCS-Y^2ZS(1-C)+XYZ^2(1-c)^2-X^2ZS(1-C)+Y^2ZS(1-C)-XYS^2\\
&=XYC(1-C)+XY(X^2+Y^2+Z^2)(1-C)^2+XYC(1-C)-XYS^2\\
&=XY(2C-2C^2+(1-C)^2-S^2)\\
&=XY(2C-2C^2+1-2C+C^2-S^2)\\
&=0\\
\end{align}

\begin{align}
H_{13}&=(C+X^2(1-C))(XY(1-C)-YS)+(XY(1-C)-ZS)(YZ(1-C)+XS)+(XZ(1-C)+YS)(C+Z^2(1-C))\\
&=XZC(1-C)-YCS+X^2Z(1-C)^2-X^2YS(1-C)+XY^2Z(1-C)^2+X^2YS(1-C)-YZ^2S(1-C)-XZS^2+XZC(1-C)+YCS+XZ^3(1-C)^2+YZ^2S(1-C)\\
&=XZ(2C-2C^2+X^2(1-C)^2+Y^2(1-C)^2+Z^2(1-C)^2-S^2)-YCS+YCS-X^2YS(1-C)+X^2YS(1-C)-YZ^2S(1-C)+YZ^2S(1-C)\\
&=XZ(2C-2C^2+1-2C+C^2-S^2)\\
&=0\\
\end{align}

なお、転置との乗算において$H_{nm}=H_{mn}$同じ内積になるため、$H_{21}=H_{12} , H_{31}=H_{13} , $$H_{32}=H_{23}$となりますので、計算を省略できます。

\begin{align}
H_{21}&=H_{12}\\
&=0
\end{align}

\begin{align}
H_{22}&=(XY(1-C)+XS)^2+(C+Y^2(1-C))^2+(YZ(1-C)-XS)^2\\
&=Y^2(1-C)^2(X^2+Y^2+Z^2)+S^2(1-Y^2)+C^2+2Y^2C(1-C)\\
&=Y^2-C^2Y^2-S^2Y^2+S^2+C^2\\
&=1
\end{align}

\begin{align}
H_{23}&=(XY(1-C)+ZS)(XZ(1-C)-YS)+(C+Y^2(1-C))(YZ(1-C)+XS)+(YZ(1-C)-XS)(C+Z^2(1-C))\\
&=X^2YZ(1-C)^2-XY^2S(1-C)+XZ^2S(1-C)-YZS^2+YZC(1-C)+XCS+Y^2Z(1-C)^2+XY^2S(1-C)+YZC(1-C)+YZ^3(1-C)^2-XCS-XZ^2S(1-C)\\
&=YZ(1-C)^2(X^2+Y^2+Z^2)-YZS^2+2YZC(1-C)\\
&=YZ(1-2C+C^2-S^2+2C-2C^2)\\
&=0
\end{align}

\begin{align}
H_{31}&=H_{13}\\
&=0
\end{align}

\begin{align}
H_{32}&=H_{23}\\
&=0
\end{align}

\begin{align}
H_{33}&=(XZ(1-C)-YS)^2+(YZ(1-C)+XS)^2+(C+Z^2(1-C))^2\\
&=X^2Z^2(1-C)^2+Y^2Z^2(1-C)^2+Z^4(1-C)^2-2XYZS(1-C)+2XYZS(1-C)+X^2S^2+Y^2S^2+C^2+2CZ^2-2C^2Z^2\\
&=Z^2(1-C)^2(X^2+Y^2+Z^2)+C^2+(X^2+Y^2)S^2-2C^2Z^2\\
&=Z^2-Z^2S^2+Z^2C^2-2Z^2C^2+C^2+S^2\\
&=1
\end{align}

うお～～～終わった～！とてもアホくさいんだけど、全部手計算したよ

つまり、

H=R\times R^T=
\begin{pmatrix}
1&0&0\\
0&1&0\\
0&0&1\end{pmatrix}

さて、これで何が言えるのかというと、ロドリゲス行列はオイラー回転と同等と言えるわけです。

さて、ここまでロドリゲス行列を示してきた理由は「クオータニオンを使うのはジンバルロックを起こさないため」というのが、半分正解で半分不正解という事をいうためだったわけです。
別にそれだけならクオータニオンは必須ではありません。

では、いよいよ真打ち、「クオータニオン」の解説をしていきます

クオータニオンとは…

さて、クオータニオンとは何でしょうか？
「ハミルトン」って数学者が定義した数の体系です。
一言で言うと、複素数(実数+虚数)の次元を拡張したものと言えます。
複素数同様、もともと別に回転のためものではないのです。

というわけで、まず虚数というか複素数を回転に利用する話をします

複素数による回転

まず、２つの複素数を考えます。

\displaylines{
Z_F=a+bi\\
Z_G=c+di
}

なお、この$Z_F , Z_G$は正規化というか、大きさが１とします。$a^2+b^2=1 , c^2+d^2=1$
そしてこの二つを乗算します

Z_F Z_G = (a+bi)(c+di)=ac-bd+(ad+bc)i

で？まぁ、とりあえずこの合成複素数の大きさを測ってみましょうか

\begin{align}
L&=(ac-bd)^2+(ad+bc)^2
&=a^2c^2-2abcd+b^2d^2+a^2d^2+2abcd+b^2c^2
&=a^2c^2+b^2d^2+a^2d^2+b^2c^2
&=a^2(c^2+d^2)+b^2(d^2+c^2)
&=1
\end{align}

はい１になりますね。乗算してもベクトルの大きさ変わらない…回転行列に似てますね？

で？
うん、実数と虚数は別世界なので一緒くたに演算はできないんだけど、ただ2乗するとマイナスの実数になる…ということで、これ2次元の行列として表せます

1=\begin{pmatrix}
1&0\\
0&1
\end{pmatrix}

i=\begin{pmatrix}
0&-1\\
1&0
\end{pmatrix}

i^2=\begin{pmatrix}
0&-1\\
1&0
\end{pmatrix}
\begin{pmatrix}
0&-1\\
1&0
\end{pmatrix}=\begin{pmatrix}
-1&0\\
0&-1
\end{pmatrix}=-1

実数部と虚数部の場所が異なるため加算は混ざりませんが
実数と実数の乗算は実数、虚数と虚数の乗算は実数のように、乗算の時は影響を及ぼすのが分かります。
ところで、回転行列は

\begin{pmatrix}
cosθ&-sinθ\\
sinθ&cosθ
\end{pmatrix}

ですが、これを分解して複素数で表すと

\begin{pmatrix}
cosθ&-sinθ\\
sinθ&cosθ
\end{pmatrix}=
cosθ
\begin{pmatrix}
1&0\\
0&1
\end{pmatrix}=sinθ\begin{pmatrix}
0&-1\\
1&0
\end{pmatrix}=cosθ+isinθ

う～ん、でもこれで回転…できる？
では具体的な数字で確認してみましょう。
30°回転して60°回転してそこからさらに120°回転させてみましょう
まず30°回転なら$Z_{30}=\frac{\sqrt{3}+i}{2}$、60°回転なら$Z_{60}=\frac{1+\sqrt{3}i}{2}$、120°回転なら$Z_{120}=\frac{-1+\sqrt{3}i}{2}$
まず、30°と60°を乗算します

\begin{align}
Z_{30}Z_{60}&=\frac{\sqrt{3}+i}{2}\frac{1+\sqrt{3}i}{2}\\
&=\frac{(\sqrt{3}+i)(1+\sqrt{3}i)}{4}\\
&=\frac{(\sqrt{3}+i+3i-\sqrt{3})}{4}\\
&=\frac{4i}{4}\\
&=i=Z_{90}
\end{align}

となり、乗算が角度の加算に対応しているのが分かります。さらに120°を乗算しましょう

\begin{align}
Z_{90}Z_{120}&=i\frac{-1+\sqrt{3}i}{2}\\
&=\frac{-\sqrt{3}-i}{2}\\
&=Z_{210}
\end{align}

となり、90°+120°=210°になるのがわかります。少なくとも行列の代わりにしてもいいのが分かると思います。

複素共役(複素共軛)と転置行列

ちなみに$Z=a+bi$に対して、虚数側の符号を逆転させたものを複素共役と言い、$\overline{Z}=a-bi$と表します。
これは乗算すると回転がなかったことになりZの大きさの二乗になります。

Z\overline{Z}=(a+bi)(a-bi)=a^2+b^2=\|Z\|^2

複素共役を行列で表すと転置行列となり

Z\overline{Z}=
\begin{pmatrix}
a&-b\\
b&a
\end{pmatrix}
\begin{pmatrix}
a&b\\
-b&a
\end{pmatrix}
=\begin{pmatrix}
a^2+b^2&ab-ba\\
ba-ab&b^2+a^2
\end{pmatrix}=
(a^2+b^2)\begin{pmatrix}
1&0\\
0&1
\end{pmatrix}=\|Z\|^2

となります。ともかく四元数やりましょ、四元数！

四元数(クオータニオン)の体系

ここからしばらくは一旦、回転のこととか忘れて代数的な所に注目して考えてください。

ハミルトンさんは、「あ～、この複素数を三次元に拡張してえなあ～俺もな～」と考えてこのようなルールを作りました。

虚数単位を３つ(i,j,k)とし、それと実数を組み合わせたものを四元数とする
虚数単位なので当然$i^2=-1 , j^2=-1 , k^2=-1$とする
加えて$ijk=-1$とする
i,j,kはそれぞれ別次元とし、それぞれ互いの加算の影響は受けないが互いの乗算によって次元を変える
- $ij=k , ji=-k$
- $jk=i , kj=-i$
- $ki=j , ik=-j$
これらをまとめたものを$Q=a+bi+cj+dk$とする
こういう体系を作りました。

さて、四元数も複素数と同様に行列で表せるのでしょうか？考えてみましょう。
いったん虚数単位のiをまた使います。これ自体は四元数のiではなく、ただ単に虚数単位と思ってください
加算で被らないように、2乗すると-単位行列になるように…
そうすると

\displaylines{
Q_R=\begin{pmatrix}
1&0\\
0&1
\end{pmatrix},Q_I=\begin{pmatrix}
0&-1\\
1&0
\end{pmatrix}\\
Q_J=\begin{pmatrix}
0&-i\\
-i&0
\end{pmatrix},Q_K=\begin{pmatrix}
-i&0\\
0&i
\end{pmatrix}
}

となります。
こう書くと、クオータニオンの虚数単位を入れ替えると符号が逆になるのもわかると思います。

なお、あまり意味はありませんが、一応このiをさらに行列化すると4x4行列となり

\displaylines{
Q_R=\begin{pmatrix}
1&0&0&0\\
0&1&0&0\\
0&0&1&0\\
0&0&0&1\\
\end{pmatrix},Q_I=\begin{pmatrix}
0&0&-1&0\\
0&0&0&-1\\
1&0&0&0\\
0&1&0&0\\
\end{pmatrix}\\
Q_J=\begin{pmatrix}
0&0&0&-1\\
0&0&1&0\\
0&-1&0&0\\
1&0&0&0\\
\end{pmatrix},Q_K=\begin{pmatrix}
0&1&0&0\\
-1&0&0&0\\
0&0&0&-1\\
0&0&1&0\\
\end{pmatrix}
}

となります。
見事に加算で影響しあわないようになっていますね。

ということはクオータニオン全体としては

Q=a\begin{pmatrix}
1&0\\
0&1
\end{pmatrix}
+b\begin{pmatrix}
0&-1\\
1&0
\end{pmatrix}
+c\begin{pmatrix}
0&-i\\
i&0
\end{pmatrix}
+d\begin{pmatrix}
-i&0\\
0&i
\end{pmatrix}
=\begin{pmatrix}
a-di&-b-ci\\
b+ci&a+di
\end{pmatrix}

と表せます

まぁ、それはともかく、回転のための材料集めをしていきましょう。

スカラー三重積、ベクトル三重積

今から言うのは、クオータニオンでも行列でもありませんが、後で使うのでスカラー三重積とベクトル三重積について
お話しします。

スカラー三重積

これは、$\vec{A}\cdot(\vec{B}\times \vec{C})$のことで、実はこれは平行6面体の体積になっています。
ベクトルの外積って、2つのベクトルに直交するベクトルであり、かつその大きさはBとCで作る平行四辺形の
体積になっています。
で、BとCで作った直交ベクトルとAのベクトルの内積をとる事で、この内積結果が射影になるため平行四辺形を底辺とした平行六面体の体積になるわけです。
このため

ベクトル三重積

これは$\vec{A}\times \vec{B}\times \vec{C}$のことです。
これはむしろ成分表示したほうが分かりやすいかもしれませんね。

\begin{align}
\vec{A}\times \vec{B}\times \vec{C}&=\vec{A} \times (\vec{B}\times \vec{C})\\
&=\vec{A}\times (B_yC_z-B_zC_y,B_zC_x-B_xC_z,B_xC_y-B_yC_x)\\
&=(A_x,A_y,A_z)\times (B_yC_z-B_zC_y,B_zC_x-B_xC_z,B_xC_y-B_yC_x)\\
&=(A_y(B_xC_y-B_yC_x)-A_z(B_zC_x-B_xC_z),
A_z(B_yC_z-B_zC_y)-A_y(B_xC_y-B_yC_x),
A_x(B_zC_x-B_xC_z)-A_z(B_yC_z-B_zC_y))\\
&=(A_y(B_xC_y-B_yC_x)-A_z(B_zC_x-B_xC_z),
A_z(B_yC_z-B_zC_y)-A_y(B_xC_y-B_yC_x),
A_x(B_zC_x-B_xC_z)-A_z(B_yC_z-B_zC_y))\\
&=(A_yB_xC_y-A_yB_yC_x-A_zB_zC_x+A_zB_xC_z,
A_zB_yC_z-A_zB_zC_y-A_yB_xC_y+A_yB_yC_x,
A_xB_zC_x-A_xB_xC_z-A_zB_yC_z+A_zB_zC_y)\\
&=(B_x(A_yC_y+A_zC_z)-C_x(A_yB_y+A_zB_z),
B_y(A_zC_z+A_yC_x)-C_y(A_zB_z+A_yB_x),
B_z(A_xC_x+A_zC_y)-C_z(A_xB_xC_z-A_zB_y))\\
\end{align}

ここで

B_x(A_yC_y+A_zC_z)-C_x(A_yB_y+A_zB_z)

の両項にそれぞれ $A_xB_xC_x$ を足してもプラスマイナスゼロになるわけだから

\begin{align}
B_x(A_yC_y+A_zC_z)-C_x(A_yB_y+A_zB_z)&=(B_x(A_xC_x+A_yC_y+A_zC_z)-C_x(A_xB_x+A_yB_y+A_zB_z)\\
&=B_x(\vec{A}\cdot{C})-C_x(\vec{A}\cdot \vec{B})
\end{align}

つまり

\vec{A}\times \vec{B}\times \vec{C}=\vec{B}(\vec{A}\cdot \vec{C})-\vec{C}(\vec{A}\cdot \vec{B})

はい、まともにやると面倒なので、ちょっと覚えておきましょう

クオータニオンの共役

複素数と同様に共役クオータニオンというものがあります。複素数の時と同様に虚数部分にマイナスをつければいいだけです。
$Q=a+bi+cj+dk$ならば共役クオータニオンは$\overline{Q}=a-bi-cj-dk$になります。

クオータニオンの積

はい、ではまずクオータニオン同士の乗算について考えてみましょう
とりあえず、２つのクオータニオンをP,Qとあらわしてみます

\displaylines{
P=P_a+P_bi+P_cj+P_dk \\ 
Q=Q_a+Q_bi+Q_cj+Q_dk \\ 
}

この二つを乗算します

\begin{align}
PQ&=(P_a+P_bi+P_cj+P_dk)(Q_a+Q_bi+Q_cj+Q_dk)\\
&=P_aQ_a-P_bQ_b-P_cQ_c-P_dQ_d\\
&+(P_aQ_b+P_bQ_a+P_cQ_d-P_dQ_c)i\\
&+(P_aQ_c-P_bQ_d+P_cQ_a+P_dQ_b)j\\
&+(P_aQ_d+P_bQ_c-P_cQ_b+P_dQ_a)k
\end{align}

3Dベクトルのクオータニオン表記

3Dベクトルをクオータニオンで表すには、xyzそれぞれをijkに割り当てればいいので

V=V_xi + V_yj + V_zk

になります。このベクトルクオータニオン同士の乗算は、例えばVとWがあるとすると

\begin{align}
VW&=(V_xi + V_yj + V_zk)(W_xi + W_yj + W_zk)\\
&=-(V_xW_x+V_yW_y+V_zW_z)\\
&+(V_yW_z-V_zW_y)i\\
&+(V_zW_x-V_xW_z)j\\
&+(V_xW_y-V_yW_x)k
\end{align}

になります。
ん？よく見たらこれって…

VW=-\vec{V}\cdot \vec{W}+(\vec{V}\times \vec{W})_xi
+(\vec{V}\times \vec{W})_yj
+(\vec{V}\times \vec{W})_zk

の形になっていますね。面白いですね。ひとまず記憶にとどめておきましょう

ロドリゲスの式とクオータニオン

ところで、ロドリゲスの式を覚えているでしょうか？ベクトルVまわりの点P回転のことでした。
式は

\vec{P}cosθ+\hat{V}(\hat{V}\cdot\vec{P})(1-cosθ)+(\hat{V} \times \vec{P})sinθ

でしたね。
うーん、外積と内積とベクトルPの組み合わせ…前述のベクトルクオータニオンどうしの乗算がなんか惜しい感じですね。

どうにか利用できないものでしょうか？

ちょっとまだ武器が足りない気がします。
ベクトルクオータニオンと、実数部付きのクオータニオンを乗算してみましょうか

\begin{align}
VQ&=(0+V_xi + V_yj + V_zk)(Q_r+Q_xi + Q_yj + Q_zk)\\
&=-V_xQ_x-V_yQ_y-V_zQ_z\\
&+(V_xQ_r+V_yQ_z-V_zQ_y)i\\
&+(-V_xQ_z+V_yQ_r+V_zQ_x)j\\
&+(V_xQ_y-V_yQ_x+V_zQ_r)k\\

&= Q_r\vec{V}-\vec{V}\cdot Q_{xyz}+(\vec{V}\times \vec{Q})_xi
+(\vec{V}\times \vec{Q})_yj
+(\vec{V}\times \vec{Q})_zk

\end{align}

ん－、近い。$Q_r=cosθ$にしたらなんかいい感じになりそうな気もしますね
あと一歩って感じ…
そういえば、共役クオータニオンってのがありましたね。乗算してみましょう(わざとらしいｗ)

\begin{align}
V\overline{Q}&=(0+V_xi + V_yj + V_zk)(Q_r-Q_xi - Q_yj - Q_zk)\\
&=V_xQ_xV_yQ_yV_zQ_z\\
&+(V_xQ_r-V_yQ_z+V_zQ_y)i\\
&+(V_xQ_z+V_yQ_r-V_zQ_x)j\\
&+(-V_xQ_y+V_yQ_x+V_zQ_r)k\\

&= Q_r\vec{V}+\vec{V}\cdot Q_{xyz}
-(\vec{V}\times \vec{Q})_xi
-(\vec{V}\times \vec{Q})_yj
-(\vec{V}\times \vec{Q})_zk

\end{align}

外積部分がマイナスになってるのが邪魔だな…あれ？そういえば、クオータニオンって乗算の順序入れ替えるとプラスマイナス逆になるよね？
じゃあ…

\begin{align}
\overline{Q}V&=(Q_r-Q_xi - Q_yj - Q_zk)(0+V_xi + V_yj + V_zk)\\
&=V_xQ_x+V_yQ_y+V_zQ_z\\
&+(V_xQ_r+V_yQ_z-V_zQ_y)i\\
&+(V_xQ_z-V_yQ_r+V_zQ_x)j\\
&+(V_xQ_y+V_yQ_x-V_zQ_r)k\\
&= Q_r\vec{V}+\vec{V}\cdot Q_{xyz}
+(\vec{V}\times \vec{Q})_xi
+(\vec{V}\times \vec{Q})_yj
+(\vec{V}\times \vec{Q})_zk
\end{align}

お！？キタよ？キタよ？キタキタキタよ？
ここで、$|Q_{xyz}|^2=1$の単位ベクトル状態だとします
で、さらにクオータニオン自体を単位クオータニオンだとすると、sin,cosを使って

\displaylines{
Q_r^2+\|Q_{xyz}\|^2=1
(cosθ)^2+(sinθ\|Q_{xyz}\|)^2=1
}

とも言えるわけです。さて、これを共役クオータニオンとベクトルの式に当てはめてみます


\overline{Q}V= cosθ\vec{V}+sinθ\vec{V}\cdot Q_{xyz}
+sinθ((\vec{V}\times \vec{Q})_xi
+(\vec{V}\times \vec{Q})_yj
+(\vec{V}\times \vec{Q})_zk)

んんんんん～～！！！惜しい！内積のsinが邪魔！！

どうにかしてコイツを(1-cos)とかにできませんかねえ！
そこで加法定理を思い出してみましょう

\begin{align}
sin(α+β)&=sinαcosβ+sinβcosα\\
cos(α+β)&=cosαcosβ-sinαsinβ\\
sin2α=sin(α+α)&=sinαcosβ+sinβcosα\\
&=sinαcosα+sinαcosα\\
&=2sinαcosα\\
cos2α=cos(α+α)&=cosαcosβ-sinαsinβ\\
&=(cosα)^2-(sinα)^2\\
&=1-2(sinα)^2=2(cosα)^2-1
\end{align}

どちらか使っていい感じになりませんかねぇ…？
ひとまず右側に(共役でない)クオータニオンをかけてみよう
チョット式がめんどくさくなるので

(\vec{V}\times \vec{Q})_xi+(\vec{V}\times \vec{Q})_yj+(\vec{V}\times \vec{Q})_zk

の部分を$(\vec{V}\times Q_{ijk})$って書き換えときますね？

また、ここで言うベクトルは基底を(i,j,k)としているので実数部は内積部分だけで、あとはベクトルになっている部分は虚数部になります。
ここら辺に注目しつつ計算していきましょう。仮に

S=\overline{Q}V=cosθ\vec{V}+sinθ\vec{V}\cdot Q_{xyz}+sinθ(V\times Q_{ijk})

とおくと
Sの実数部分は$S_r=sinθ\vec{V}\cdot Q_{xyz}$であり、虚数部分は$S_{xyz}=cosθ\vec{V}+sinθ(V\times Q_{ijk})$です
ここに右側から$Q=cosθQ_r+sinθ Q_{xyz}$を乗算します。さぁ、また地味な計算頑張りましょう。
$Q_r=1$として考えるため$Q=cosθ+sinθQ_{xyz}$であることから

\begin{align}
\overline{Q}V{Q}&=SQ\\
&=(cosθ\vec{V}+sinθ\vec{V}\cdot Q_{xyz}+sinθ(V\times Q_{ijk}))(cosθ+Q_{xyz}sinθ)\\
&=(cosθ\vec{V}+sinθ\vec{V}\cdot Q_{xyz}+sinθ(V\times Q_{ijk}))cosθ\\
&+(cosθ\vec{V}+sinθ\vec{V}\cdot Q_{xyz}+sinθ(V\times Q_{ijk}))Q_{xyz}sinθ\\
&=\vec{V}cos^2θ+\vec{V}\cdot Q_{xyz}sinθcosθ+(V\times Q_{ijk})sinθcosθ\\
&+\vec{V}Q_{xyz}sinθcosθ+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ+(V\times Q_{ijk})Q_{xyz}sin^2θ\\
&=\vec{V}cos^2θ+\vec{V}\cdot Q_{xyz}sinθcosθ+(V\times Q_{ijk})sinθcosθ\\
&+\vec{V}Q_{xyz}sinθcosθ+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ+(V\times Q_{ijk})Q_{xyz}sin^2θ\\
\end{align}

うーん、ぱっと見どうにもならなさそう！
でも、前の計算で$\vec{V}Q_{xyz}sinθcosθ$の部分って、$-(\vec{V}\cdot Q_{xyz}) sinθcosθ +(V\times Q_{ijk})sinθcosθ$になってたので、

\begin{align}
\overline{Q}V{Q}&=\vec{V}cos^2θ+\vec{V}\cdot Q_{xyz}sinθcosθ+(V\times Q_{ijk})sinθcosθ\\
&+(-\vec{V}\cdot Q_{xyz}+\vec{V}\times Q_{ijk})sinθcosθ\\
&+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ+(V\times Q_{ijk})Q_{xyz}sin^2θ\\
&=\vec{V}cos^2θ+\vec{V}\cdot Q_{xyz}sinθcosθ+(V\times Q_{ijk})sinθcosθ\\
&-\vec{V}\cdot Q_{xyz}sinθcosθ+\vec{V}\times Q_{ijk}sinθcosθ\\
&+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ+(V\times Q_{ijk})Q_{xyz}sin^2θ\\

&=\vec{V}cos^2θ+(V\times Q_{ijk})sinθcosθ+(\vec{V}\times Q_{ijk})sinθcosθ\\
&+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ+(V\times Q_{ijk})Q_{xyz}sin^2θ\\

&=\vec{V}cos^2θ+(\vec{V}\times Q_{ijk})2sinθcosθ\\
&+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ+(V\times Q_{ijk})Q_{xyz}sin^2θ
\end{align}

前述の$sin2θ=2sinθcosθ$より

\begin{align}
\overline{Q}V{Q}
&=\vec{V}cos^2θ+(V\times Q_{ijk})sin2θ\\
&+(\vec{V}\cdot Q_{xyz})Q_{xyz}sin^2θ + (\vec{V}\times Q_{ijk})Q_{xyz}sin^2θ
\end{align}

さぁ！もうひと頑張りだ！
頑張れ💛頑張れ💛
で、$(\vec{V}\times Q_{ijk})Q_{xyz}$この部分がまた「虚数のみのクオータニオンの乗算」になっているので
ちょっといったん$H=(\vec{V}\times Q_{ijk})$とおくと

HQ_{xyz}=- \vec{H_{xyz}}\cdot \vec{Q_{xyz}}+
(\vec{H_{xyz}}\times \vec{Q_{xyz}})_xi + 
(\vec{H_{xyz}}\times \vec{Q_{xyz}})_yj + 
(\vec{H_{xyz}}\times \vec{Q_{xyz}})_zk

段々ややこしくなってきましたがここでHを戻します

\begin{align}
(\vec{V}\times Q_{ijk})Q_{xyz}&=- (\vec{V}\times \vec{Q_{ijk}})\cdot \vec{Q_{xyz}}\\
&+((\vec{V}\times \vec{Q_{xyz})}\times \vec{Q_{ijk}})_xi \\ 
&+((\vec{V}\times \vec{Q_{xyz})}\times \vec{Q_{ijk}})_yj \\
&+((\vec{V}\times \vec{Q_{xyz})}\times \vec{Q_{ijk}})_zk 
\end{align}

はい、出ましたベクトル三重積。

\vec{A}\times \vec{B}\times \vec{C}=\vec{B}(\vec{A}\cdot \vec{C})-\vec{C}(\vec{A}\cdot \vec{B})

でしたね？だから

\vec{V}\times \vec{Q_{xyz}}\times \vec{Q_{ijk}}=\vec{Q_{xyz}}(\vec{V}\cdot \vec{Q_{ijk}})-\vec{Q_{ijk}}(\vec{V}\cdot \vec{Q_{xyz}})

なお、$Q_{ijk}$と$Q_{xyz}$は同じものです。同じベクトル同士の内積は1ですこのため

\begin{align}
\vec{Q_{xyz}}\times \vec{V}\times \vec{Q_{ijk}}&=\vec{V}(\vec{Q_{xyz}}\cdot \vec{Q_{ijk}})-\vec{Q_{ijk}}(\vec{Q_{xyz}}\cdot \vec{V})\\
&=\vec{V}-\vec{Q_{ijk}}(\vec{Q_{xyz}}\cdot \vec{V})
\end{align}

また、$-\vec{Q_{xyz}}\cdot (\vec{V}\times \vec{Q_{ijk}})$この部分ですが、中の外積結果は$Q_{ijk}$に直交しており、それと同じ$Q_{xyz}$を内積しているわけだから$\vec{A}\perp \vec{B} \Leftrightarrow \vec{A}\cdot\vec{B}=0$であるため

\begin{align}
Q_{xyz}(\vec{V}\times Q_{ijk})&=\vec{V}-\vec{Q_{ijk}}(\vec{Q_{xyz}}\cdot \vec{V})
\end{align}

すみません！計算が…あわない…ッ！

あ…あれ？
計算が合わない…

ここまで読んでくれた人、すみません！アドベントカレンダーに間に合わないため、結論だけ書きます。
計算が合うようになったらまた更新しますので、気になる人はストックしておいてください！！！

で、結論から言うと

クオータニオンが$\overline{Q}VQ$なんていうけったいな計算になっているのは、ロドリゲスの公式に合わせるため
で、その計算過程で2sinθcosθが出てきてsin2θになるため、回転のつじつまを合わせるためにθ/2を用いている

という事です。

マジですんません！！あとで必ず更新しますんで！！！！

で、残りは
「こんなけったいなクオータニオンを使うと何が嬉しいのか」
なので、そこまでやったら完成です。