13. 回転とクォータニオン

第8章からrotateX、rotateY、rotateZで物を回して来ましたし、第11章では2つの角度からカメラを作りました。
其れで十分なのは、任意の向きを保存したり、沢山の回転を組み合わせたり、或る向きから別の向きへ滑らかに補間したりする必要が出て来る迄です。
そして其れこそ、アニメーションシステムが一日中やっている事です。
最後の此の章では、回転を表す3つの方法、オイラー角、軸角度、クォータニオンを見て行きます。
1つ目が何故問題を起こすのか、2つ目がどの様に一般的な回転行列に繋がるのか、そして1843年に発明された不思議な4次元の数である3つ目が、何故実際のゲームエンジンで使われているのかを確かめます。

13.1 オイラー角と其の問題

向きを表す最も人間らしい方法は、軸毎に1つずつ、3つの角度を使う事です。
此れをオイラー角と呼び、乗り物やカメラでは3つに名前が付いています：

此のシリーズでは$𝐑={𝐑}_y(\text{yaw}){𝐑}_x(\text{pitch}){𝐑}_z(\text{roll})$と組み合わせます。
第11章のカメラと同じ、ピッチの後にヨーという順番です。
3つの数だけで、エディターに打ち込むのも理解するのも簡単です。
では何が問題なのでしょうか？

問題1：順番で結果が変わる

行列の積は交換法則を満たさない（第7章）ので、同じ3つの角度でも順番が違えば違う向きになります。
第8章の行列で、$\cos 90\mathrm{°}=0$、$\sin 90\mathrm{°}=1$として、$𝐯=(1,0,0)$をX周りに90°、Y周りに90°回してみましょう：

Xが先、次にY。 ${𝐑}_x(1,0,0)=(1,0,0)$（ベクトルがX軸上にあるので動きません）、次に${𝐑}_y(1,0,0)=(0,0,-1)$です。

Yが先、次にX。 ${𝐑}_y(1,0,0)=(0,0,-1)$、次に${𝐑}_x(0,0,-1)=(0,0\cdot 0+(-1)(-1),0)=(0,1,0)$です。

同じ2つの角度なのに、片方ではベクトルが画面の奥（−Z）を、もう片方では上を向いています。
軸の順番は12通りあり、ソフトによって選ぶ順番が違うので、「ヨー30、ピッチ20、ロール10」は順番も言わない限り何の意味もありません。

問題2：ジンバルロック

3つの回転を、ジャイロスコープを支えるジンバルの様な、入れ子になった3つのリングだと考えて下さい：

では、丁度90°ピッチを上げます。
内側のリングの軸、つまり物体の前後の軸が真上を向き、外側のリングのY軸と重なります。
ヨーとロールが同じ直線の周りを回る様になり、3つ目の方向に物体を回せるリングがもう残っていません。
3つの自由度の内1つが、単純に消えてしまったのです。
此れがジンバルロックです。

数字でも確かめられます。
ピッチが90°の時、計算すると${𝐑}_y(\text{yaw}){𝐑}_x(90\mathrm{°}){𝐑}_z(\text{roll})={𝐑}_y(\text{yaw}-\text{roll}){𝐑}_x(90\mathrm{°})$となり、ヨー−ロールの差しか意味を持ちません。
例えば（ヨー30°、ピッチ90°、ロール0°）と（ヨー50°、ピッチ90°、ロール20°）は全く同じ行列になります：

90°の近くでも殆ど同じ位悪く、向きを少し変えるだけで角度が大きく跳ね、アニメーションが激しく暴れます。
此れはゲームだけの問題ではありません。
アポロ宇宙船の誘導装置は本物の3つのジンバルを使っていて、宇宙飛行士はジンバルロックを避けなければなりませんでした。
アポロ11号のマイケル・コリンズが「クリスマスに4つ目のジンバルを送ってくれないか」と冗談を言った話は有名です。

問題3：補間

敵がヨー170°からヨー−170°へ向きを変えるとします。
180°を通れば、たった20°しか離れていません。
でも角度の数字をlerp（第2章）で補間すると、角度は170、85、0、−85、−170と進み、逆向きに340°も回ってしまいます。
そして3つの角度が同時に変わると、補間の経路は最短の回転ではなく変な曲線になります。

従ってオイラー角は入力（エディターの回転欄や、ロールせずピッチを制限する一人称カメラ）には最適ですが、一般的な向きの保存、組み合わせ、補間には向きません。

13.2 軸角度とロドリゲスの回転公式

1775年、レオンハルト・オイラーは驚くべき事を証明しました。
3Dのどんな回転も、又どんな回転の連鎖も、或る軸周りの1回の回転と同じだという事です。
従って回転は、単位ベクトル$𝐤$（軸）と角度$\theta$で表せます。
正の向きは右手の法則に従います。
親指を$𝐤$の方向に向けると、他の指が曲がる向きに回ります。

公式

ロドリゲスの回転公式は、ベクトル$𝐯$を単位軸$𝐤$周りに$\theta$回転させます：

仕組み

第4章の射影を使って、$𝐯$を2つの部分に分けます：

• 軸に沿った部分：${𝐯}_{\parallel }=(𝐤\cdot 𝐯)𝐤$
• 軸に垂直な部分：${𝐯}_{\perp }=𝐯-{𝐯}_{\parallel }$

軸に沿った部分は、車輪の車軸の様に全く動きません。
垂直な部分は、$𝐤$に垂直な平面の中でくるくる回ります。
其の平面の中で回すには、第8章の2D回転の$𝐢$と$𝐣$の様に、同じ長さで互いに垂直な2つの方向が必要です。
1つ目は${𝐯}_{\perp }$自身です。
2つ目は$𝐤\times 𝐯$です。
此れは$𝐤$とも${𝐯}_{\perp }$とも垂直で、長さは${𝐯}_{\perp }$と同じです（第5章：長さは$\Vert 𝐤\Vert \Vert {𝐯}_{\perp }\Vert \sin 90\mathrm{°}$で、又$𝐤\times {𝐯}_{\parallel }=\mathrm{𝟎}$です）。

従って2D回転と全く同じ様に：

動かなかった部分を足して$𝐯\prime ={𝐯}_{\parallel }+{𝐯}_{\perp }\prime$とし、${𝐯}_{\perp }$を$𝐯-{𝐯}_{\parallel }$で置き換えます：

${𝐯}_{\parallel }=𝐤(𝐤\cdot 𝐯)$を入れれば、ロドリゲスの公式です。

計算例

$𝐯=(1,0,0)$を、対角線$𝐤=\frac{1}{\sqrt{3}}(1,1,1)$周りに120°回します。
必要なのは$\cos 120\mathrm{°}=-\frac{1}{2}$、$\sin 120\mathrm{°}=\frac{\sqrt{3}}{2}$、$1-\cos 120\mathrm{°}=\frac{3}{2}$です。

では3つの項です：

足し合わせます：

X軸がY軸に移りました。
対角線周りの120°回転は、軸を順番に入れ替えます。
xがyへ、yがzへ、zがxへ行きます。

行列にする

3つの項は、どれも行列×$𝐯$の形で書けます：

• $𝐯\cos \theta$は$\cos \theta 𝐈$×$𝐯$です。
• $𝐤\times 𝐯$は、下の「外積行列」$𝐊$を使って$𝐊𝐯$です（掛け算を展開すると第5章の外積の公式になります）。
• $𝐤(𝐤\cdot 𝐯)$は$(𝐤{𝐤}^{\mathrm{\top }})𝐯$です。$𝐤{𝐤}^{\mathrm{\top }}$は3×1の行列と1×3の行列の積で、3×3の行列になります（第7章）。

従って$𝐑=\cos \theta 𝐈+\sin \theta 𝐊+(1-\cos \theta )𝐤{𝐤}^{\mathrm{\top }}$です。
$c=\cos \theta$、$s=\sin \theta$、$t=1-c$と置き、要素毎に書くと：

此れがlinalg.hhのrotate(axis, angle)で、要素も1つ1つ同じです。
簡単な確認をしましょう。
Z軸$𝐤=(0,0,1)$なら$k_x$や$k_y$を含む物は全て消え、右下は$t+c=1$になり、第8章のrotateZ其の物が残ります。
上の例では$t{k}_i{k}_j$は全て$\frac{3}{2}\cdot \frac{1}{3}=\frac{1}{2}$、$s{k}_i$は全て$\frac{\sqrt{3}}{2}\cdot \frac{1}{\sqrt{3}}=\frac{1}{2}$なので：

$𝐢$の行き先である1列目は$(0,1,0)$で、さっきと同じです。

軸角度は回転を作るには素晴らしい方法で（「此の蝶番の周りに30°回す」）、ジンバルロックもありません。
でも2つの軸角度の回転を組み合わせる簡単な公式は無く、補間も軸が同じ時しか上手く行きません。
其の為にクォータニオンが必要です。

13.3 クォータニオン：公式

少し歴史を

1830年代、アイルランドの数学者ウィリアム・ローワン・ハミルトンは、複素数が2Dの点を掛け算する様に、3Dの点を掛け算しようとしていました。
何年も上手く行きませんでした。
1843年10月16日、妻とダブリンのロイヤル運河沿いを歩いていた時、彼は3つではなく4つの数が必要だと突然気付きました。
余りに興奮して、ナイフでブルーム橋（Brougham Bridge）の石に其の規則を刻み込んだ程です：

刻んだ文字はとっくに消えてしまいましたが、今は橋に記念の銘板があります。

クォータニオンとは

普通の数で2乗して−1になる物は無いので、$i$、$j$、$k$は上の規則に従うだけの新しい記号です。
クォータニオン（四元数）は次の物です：

$w,x,y,z$は4つの普通の数です。
$w$をスカラー部、$(x,y,z)$をベクトル部と呼び、良く$q=(w,𝐯)$と書きます。
（此の$i,j,k$は第2章の基底ベクトルではありませんが、其の祖先です。
3Dの単位ベクトルをi、j、kと書くのは、ハミルトンのクォータニオンから来ています。）

ハミルトンの規則から、全ての積が求められます。
例えば$ijk=-1$の両辺に右から$k$を掛けると$ijkk=-k$で、$kk=-1$なので$-ij=-k$、つまり$ij=k$です。
同じ様な手順で表全体が出来ます（行×列）：

見て下さい。
$ij=k$なのに$ji=-k$です。
クォータニオンの掛け算は、回転と同じく交換法則を満たしません。
外積にも良く似ています：$𝐢\times 𝐣=𝐤$、$𝐣\times 𝐢=-𝐤$です。

ハミルトン積

2つのクォータニオンを掛けるとは、括弧を展開して表を使う事です。
一度一般の場合で其れをやって整理すると、内積と外積を使った簡潔な形になります：

成分毎に書き出すと、ヘッダーのoperator*(Quat, Quat)其の物です：

順番が大事になるのは外積の項の所為です。
${𝐯}_a\times {𝐯}_b=-{𝐯}_b\times {𝐯}_a$だからです。

計算例。 $(1+i+j)(2+k)$です。
表を使って括弧を展開すると：

ベクトルの形では、$q_a=(1,(1,1,0))$、$q_b=(2,(0,0,1))$として：

どちらの方法でも$2+3i+j+k$です。

長さ、共役、逆元

クォータニオンの長さは第3章の4D版で、長さ1の物を単位クォータニオンと呼びます：

共役はベクトル部の符号を反転します：

ベクトルの形でクォータニオンと其の共役を掛けると、スカラー部は$w^2+𝐯\cdot 𝐯$、ベクトル部は$-w𝐯+w𝐯+𝐯\times (-𝐯)=\mathrm{𝟎}$です。
従って$q{q}^{\ast }=\Vert q{\Vert }^2$で、単位クォータニオンなら1です。
つまり単位クォータニオンでは共役が逆元です。
回転行列の逆行列が転置だった（第10章）のと同じで、しかも更に安く、符号を3つ反転するだけです。

13.4 クォータニオンでベクトルを回す仕組み

公式

単位クォータニオンは回転です。
単位軸$𝐤$周りの$\theta$回転を表すクォータニオンは：

ベクトル$𝐯$を回すには、其れを「純粋な」クォータニオン$p=(0,𝐯)$として書き、挟み込みます：

$p\prime$のベクトル部が回転後のベクトルで、スカラー部は必ず0になります。

何故半分の角度なのか

$q$が左から1回、右から1回、合計2回掛かるからです。
実際に見てみましょう。
$𝐯=(1,0,0)$をZ周りに90°回します。
$h=\frac{\sqrt{2}}{2}\approx 0.7071$（$\cos 45\mathrm{°}$と$\sin 45\mathrm{°}$）とすると：

まず前の節の成分の公式で$qp$です：

従って$qp=(0,h,h,0)$です。
ベクトル$(h,h,0)$は$(1,0,0)$を45°だけ、つまり半分だけ回した物です。
次に其れに右から$q^{\ast }$を掛けます：

$p\prime =(0,0,1,0)$なので$𝐯\prime =(0,1,0)$です。
残りの45°で、合計90°です。
左右が半分ずつ仕事をしました。
（軸に垂直でないベクトルの場合、1回目の掛け算でスカラー部も0でなくなり、2回目の掛け算が其れを打ち消します。
其の為、1回の掛け算ではなく挟み込みが必要なのです。）

一般に、$q=(\cos \frac{\theta }{2},𝐤\sin \frac{\theta }{2})$で$qp{q}^{\ast }$を展開すると、13.2節のロドリゲスの公式になります。
但し$\cos \theta$の所が${\cos }^2\frac{\theta }{2}-{\sin }^2\frac{\theta }{2}$、$\sin \theta$の所が$2\sin \frac{\theta }{2}\cos \frac{\theta }{2}$、$1-\cos \theta$の所が$2{\sin }^2\frac{\theta }{2}$になっています。
此れらは三角関数の2倍角の公式なので、同じ数です。
$\theta =90\mathrm{°}$で確かめましょう：${\cos }^{2}45\mathrm{°}-{\sin }^{2}45\mathrm{°}=0.5-0.5=0=\cos 90\mathrm{°}$、$2\sin 45\mathrm{°}\cos 45\mathrm{°}=2\cdot 0.5=1=\sin 90\mathrm{°}$です。
挟み込みが角度を2倍にするので、半分の角度を入れておくのです。

13.5 組み合わせ、逆回転、変換

回転の組み合わせ

まず$q_1$で回し、次に$q_2$で回すには、2回挟み込みます：

（再び「靴下と靴」の規則$(q_2{q}_1{)}^{\ast }=q_1^{\ast }{q}_2^{\ast }$を使っています。）
従って組み合わせた回転は1つのクォータニオン$q_2{q}_1$で、行列と全く同じく右から左へ読みます。

計算例。 まずX周りに90°、次にY周りに90°、13.1節の「Xが先」の場合です。
$h^2=0.5$を使います：

従って$q_2{q}_1=(0.5,0.5,0.5,-0.5)$です。
此れはどんな回転でしょうか？
$\cos \frac{\theta }{2}=0.5$なので$\frac{\theta }{2}=60\mathrm{°}$、つまり$\theta =120\mathrm{°}$です。
軸はベクトル部を$\sin 60\mathrm{°}$で割った物で、$\frac{1}{\sqrt{3}}(1,1,-1)$です。
2回の90°回転が、対角線周りの1回の120°回転に纏まりました。
正にオイラーの定理の通りです。
此れで$(1,0,0)$を回すと$(0,0,-1)$になり、13.1節の行列と同じ答えです。

逆回転と二重被覆

単位クォータニオンの逆元は共役です。
$q^{\ast }$は同じ角度だけ逆向きに回します。

不思議な事があります。
$q$と$-q$は同じ回転なのです。
$(-q)p(-q{)}^{\ast }$では2つのマイナスが打ち消し合います。
半分の角度からも分かります。
$\theta +360\mathrm{°}$回すのは$\theta$回すのと同じ回転ですが、$\cos \frac{\theta +360\mathrm{°}}{2}=-\cos \frac{\theta }{2}$（sinも同様）なので、4つの数全ての符号が反転します。
従って全ての回転には丁度2つのクォータニオンがあり、此の性質を二重被覆と呼びます。
此れは補間で大事になります。

行列へ

GPUは行列を欲しがります。
回転行列の列は$𝐢$、$𝐣$、$𝐤$の行き先（第6章）なので、其々を$qp{q}^{\ast }$で回し、$w^2+x^2+y^2+z^2=1$を使って整理します。
結果がヘッダーのtoMat4です：

単位クォータニオンを前提にしています。
計算例。 13.2節の$(1,1,1)$周りの120°回転は、クォータニオンでは$(\cos 60\mathrm{°},\frac{1}{\sqrt{3}}(1,1,1)\sin 60\mathrm{°})=(0.5,0.5,0.5,0.5)$です。
2つの成分の積はどれも$0.25$なので：

ロドリゲスの公式で出したのと同じ行列です。
どの要素も2つの成分を掛けているので、$-q$でも同じ行列になります。
此れも二重被覆です。

再正規化

掛け算をする度に少しずつ丸め誤差が出ます（floatは約7桁しかありません）。
毎フレームq = step * qで向きを更新すると、其の誤差が溜まり、長さが少しずつ1からずれて行きます。
長さが1でないクォータニオンは、もう回転だけをする訳ではありません。
toMat4が物体を潰したり伸ばしたりもし始めます。
対策は安く、更新後にnormalize(q)を呼ぶだけです（割り算4回と平方根1回）。
3×3の回転行列では同じずれで列が垂直でなくなり、其れを直すには外積が何回も必要です。
此の安い修理も、エンジンが向きをクォータニオンで保存する理由の1つです。

13.6 補間：lerp、nlerp、slerp

アニメーションでは、$t$が0から1に進むにつれて、向き$a$から向き$b$へ滑らかに進む必要があります。
単位クォータニオンを4次元の球面上の点と考えて下さい。
$a$から$b$へ、球面に沿って歩きたいのです。

lerp。 4つの数を補間します（第2章）：$a+t(b-a)$。
此れは球の内側を真っ直ぐ突っ切るので、真ん中では長さが1より短くなります。
「回転無し」$(1,0,0,0)$と「Y周りに180°」$(0,0,1,0)$の真ん中は$(0.5,0,0.5,0)$で、長さは0.7071です。
もう回転ではありません。

nlerp。 lerpを正規化します：$\mathrm{normalize}(a+t(b-a))$。
此れで球面上に乗り、最短経路を通るので、此れで十分な事も多いです。
但し速さが一定ではありません。
両端では遅く、真ん中では速くなります。

slerp（球面線型補間）は球面に沿って一定の速さで歩きます。
4Dでの$a$と$b$の間の角度を$\mathrm{\Omega }$とし、4Dの内積（第4章）で求めます：$\cos \mathrm{\Omega }=a\cdot b=w_a{w}_b+x_a{x}_b+y_a{y}_b+z_a{z}_b$。
すると：

ヘッダーのslerpには、大事な工夫が2つ加わっています：

1. 短い方の道。 なら$b$の符号を反転します。二重被覆により$-b$は同じ回転ですが、球の近い側にあります。此れが無いと、0°から350°への補間は短い方の10°ではなく、遠回りで350°回ってしまいます。
2. ほぼ平行な時の代わり。 $a\cdot b>0.9995$なら$\mathrm{\Omega }$は殆ど0で$\sin \mathrm{\Omega }$も殆ど0なので、其れで割ると値が爆発します。其処迄近ければnlerpで殆ど同じ結果になるので、其れを使います。

計算例。 回転無し$a=(1,0,0,0)$から、Y周りに90°の$b=(0.7071,0,0.7071,0)$への真ん中です。
$a\cdot b=0.7071$なので$\mathrm{\Omega }=45\mathrm{°}$です。
（4Dでの角度は回転角度の半分で、此処でも半分の角度が出て来ます。）
$t=0.5$では2つの重みが等しくなります：

そして$(0.9239,0,0.3827,0)=(\cos 22.5\mathrm{°},0,\sin 22.5\mathrm{°},0)$はY周りの45°回転で、丁度90°の半分です。

slerpとnlerpの比較。 回転無しからY周りに180°への、回転角度です：

slerpは均等に進みます。
nlerpは前半では遅れ、後半では先に進みますが、0、0.5、1では同じ所に来ます。
小さな幅（アニメーションの2つのキーフレームの間等）なら、差は僅かです。

数字で見る短い方の道。 Y周りに350°の回転は$b=(\cos 175\mathrm{°},0,\sin 175\mathrm{°},0)=(-0.9962,0,0.0872,0)$です。
するとなので、slerpは代わりに$-b=(0.9962,0,-0.0872,0)$、つまり−10°の回転を使います。
真ん中では−5°で、期待通り短い方を回ります。

13.7 グラフィックスプログラミングでの使い方

• スケルタルアニメーション。 各ボーンのキーフレームは回転をクォータニオンで保存します（平行移動と拡大縮小と一緒に）。再生する時、エンジンは現在の時刻を挟む2つのキーを探し、其の間をslerp又はnlerpします。歩きから走りへのブレンドも、ブレンドの重みで同じ事をするだけです。キー同士は近いので、安いnlerpを使うエンジンも多いです。
• 滑らかな旋回。 カメラや砲台はcurrent = slerp(current, target, amount)で毎フレーム少しずつ目標の向きへ回せます。ジンバルロックも無く、常に最短の弧を通ります。
• 物理エンジンは各剛体の向きをクォータニオンで保存します。角速度$𝝎$（ワールド空間、ラジアン毎秒）を使うと、1ステップは$q\leftarrow q+\frac{\mathrm{\Delta }t}{2}(0,𝝎)q$で、其の後にnormalizeします。
• ネットワーク。 単位クォータニオンは4つの数ですが、必要なのは3つだけです。$w^2+x^2+y^2+z^2=1$なので4つ目は復元出来ます。人気の方式「smallest three」は、一番大きい成分の番号を2ビットで、残りの3つと一緒に送ります。二重被覆のお陰で、送る側は符号を反転して、省いた成分を常に正に出来ます。

何故シェーダーには普通行列を渡すのか。 頂点シェーダーは既に平行移動、拡大縮小、投影も含んだ4×4行列を掛けていますし、回転だけでも行列の方がクォータニオンより演算が少なく済みます。
従って普通の方針は、CPU側では向きをクォータニオンで持ち、物体毎に1回toMat4で変換してモデル行列を作り、行列を送る事です。
クォータニオンがシェーダーに出て来るのは、メモリが大事な時です。
例えば巨大なインスタンスバッファのインスタンス毎の向きや、コンパクトなスキニングデータです。
標準的なシェーダー関数は次の通りです。
$qp{q}^{\ast }$を展開して、2回の外積に整理した物です：$𝐯\prime =𝐯+2𝐮\times (𝐮\times 𝐯+w𝐯)$（$𝐮$はベクトル部）。

// q = float4(x, y, z, w): xyz is the vector part, w is the scalar part
float3 rotateByQuat(float4 q, float3 v) {
  return v + 2.0 * cross(q.xyz, cross(q.xyz, v) + q.w * v);
}

// q = vec4(x, y, z, w): xyz is the vector part, w is the scalar part
vec3 rotateByQuat(vec4 q, vec3 v) {
  return v + 2.0 * cross(q.xyz, cross(q.xyz, v) + q.w * v);
}

// q = float4(x, y, z, w): xyz is the vector part, w is the scalar part
float3 rotateByQuat(float4 q, float3 v) {
  return v + 2.0 * cross(q.xyz, cross(q.xyz, v) + q.w * v);
}

保存順に注意。 シェーダーではクォータニオンはfloat4又はvec4として(x, y, z, w)の順に保存し、q.xyzがベクトル部、q.wがスカラー部になる様にします。
あたし達のC++のQuatは(w, x, y, z)の順で宣言されています。
従ってアップロードする時は並べ替えます：

float gpuQuat[4] = {q.x, q.y, q.z, q.w}; // shader order: x, y, z, w

忘れると、単位元$(w,x,y,z)=(1,0,0,0)$が$x=1$、$w=0$と読まれ、X周りに180°の回転になり、モデル全体が逆さまになります。

13.8 其の他の使い道

• 宇宙機。 人工衛星や宇宙ステーションの姿勢制御ソフトウェアは向きをクォータニオンで表し、スタートラッカーは宇宙機の姿勢をクォータニオンで報告します。ジンバルロックが無く、再正規化も安いので、何年も再起動せずに動き続けるコンピューターには大事です。
• ドローン及びスマートフォン。 IMUのセンサーフュージョンのアルゴリズム（MadgwickフィルターやMahonyフィルター等）は、ジャイロ、加速度センサー、磁気センサーを組み合わせて向きのクォータニオンを作ります。Androidの回転ベクトルセンサーは単位クォータニオンの成分を報告します。
• VR。 OpenXRは毎フレーム、ヘッドセットと各コントローラーの姿勢を位置と向きのクォータニオンで報告します。
• ロボット工学。 ROSのメッセージは向きをクォータニオンで、あたし達のシェーダーと同じ(x, y, z, w)の順で保存します。
• 分子モデリング。 ドッキングのシミュレーションは分子を絶えず回しますし、2つの分子を最も良く重ねる回転を求める問題には、クォータニオンを使った美しい解法があります。
• 3Dソフト。 Blenderでは、オブジェクトやボーン毎にオイラー回転モードとクォータニオン回転モードを切り替えられます。正に此の章で見た理由の為です。

13.9 コード

linalg.hhの回転の部分です。
まず13.2節のrotate(axis, angle)、其の後にクォータニオンの全てが続きます。
Quatのデフォルト値に注目して下さい。
w = 1で残りが0なので、Quat{}は「回転無し」です。
fromAxisAngleは軸を正規化してくれますし、半分の角度を使っています。

// 任意軸周りの回転（ロドリゲスの回転公式） / rotation around any unit axis (Rodrigues) (chapter 13)
inline Mat4 rotate(Vec3 axis, float angle) {
  Vec3 a = normalize(axis);
  float c = std::cos(angle), s = std::sin(angle), t = 1.0f - c;
  Mat4 m = Mat4::identity();
  m.at(0, 0) = t * a.x * a.x + c;
  m.at(0, 1) = t * a.x * a.y - s * a.z;
  m.at(0, 2) = t * a.x * a.z + s * a.y;
  m.at(1, 0) = t * a.x * a.y + s * a.z;
  m.at(1, 1) = t * a.y * a.y + c;
  m.at(1, 2) = t * a.y * a.z - s * a.x;
  m.at(2, 0) = t * a.x * a.z - s * a.y;
  m.at(2, 1) = t * a.y * a.z + s * a.x;
  m.at(2, 2) = t * a.z * a.z + c;
  return m;
}

// q = w + xi + yj + zk
struct Quat {
  float w = 1.0f, x = 0.0f, y = 0.0f, z = 0.0f;

  // 単位軸axis周りにangleラジアン回転 / rotation of angle radians around a unit axis
  static Quat fromAxisAngle(Vec3 axis, float angle) {
    Vec3 a = normalize(axis);
    float s = std::sin(angle * 0.5f);
    return {std::cos(angle * 0.5f), a.x * s, a.y * s, a.z * s};
  }
};

// ハミルトン積 / Hamilton product
inline Quat operator*(Quat a, Quat b) {
  return {
    a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
    a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
    a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
    a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
  };
}

inline Quat conjugate(Quat q) { return {q.w, -q.x, -q.y, -q.z}; }

inline Quat normalize(Quat q) {
  float len = std::sqrt(q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z);
  return {q.w / len, q.x / len, q.y / len, q.z / len};
}

// ベクトルを回転：v' = q v q* / rotate a vector: v' = q v q*
inline Vec3 rotate(Quat q, Vec3 v) {
  Quat p{0.0f, v.x, v.y, v.z};
  Quat r = q * p * conjugate(q);
  return {r.x, r.y, r.z};
}

// 単位クォータニオンから回転行列へ / unit quaternion to rotation matrix
inline Mat4 toMat4(Quat q) {
  float w = q.w, x = q.x, y = q.y, z = q.z;
  Mat4 m = Mat4::identity();
  m.at(0, 0) = 1 - 2 * (y * y + z * z);
  m.at(0, 1) = 2 * (x * y - w * z);
  m.at(0, 2) = 2 * (x * z + w * y);
  m.at(1, 0) = 2 * (x * y + w * z);
  m.at(1, 1) = 1 - 2 * (x * x + z * z);
  m.at(1, 2) = 2 * (y * z - w * x);
  m.at(2, 0) = 2 * (x * z - w * y);
  m.at(2, 1) = 2 * (y * z + w * x);
  m.at(2, 2) = 1 - 2 * (x * x + y * y);
  return m;
}

// 球面線型補間 / spherical linear interpolation
inline Quat slerp(Quat a, Quat b, float t) {
  float d = a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z;
  if (d < 0.0f) { // 短い方の道を選ぶ / take the shorter path
    b = {-b.w, -b.x, -b.y, -b.z};
    d = -d;
  }
  if (d > 0.9995f) { // ほぼ同じ向き：線型補間で十分 / nearly equal: plain lerp is enough
    Quat r{a.w + t * (b.w - a.w), a.x + t * (b.x - a.x), a.y + t * (b.y - a.y), a.z + t * (b.z - a.z)};
    return normalize(r);
  }
  float theta = std::acos(d);
  float sa = std::sin((1.0f - t) * theta) / std::sin(theta);
  float sb = std::sin(t * theta) / std::sin(theta);
  return {sa * a.w + sb * b.w, sa * a.x + sb * b.x, sa * a.y + sb * b.y, sa * a.z + sb * b.z};
}

第3章のnormalize(Vec3)と同じく、normalize(Quat)もゼロのクォータニオンは直せません（ゼロ除算になります）。
でも回転のクォータニオンがゼロになる事は無いので、問題ありません。

そしてテストプログラムです。
此の章の計算例を再現し、扱いにくい軸と角度でクォータニオンとロドリゲスの行列を比べ、13.5節のずれも見せます。

#include <cstdio>
#include "linalg.hh"

static void printQuat(const char *name, Quat q) {
  printf("%-8s (w %.4f, x %.4f, y %.4f, z %.4f)\n", name, q.w, q.x, q.y, q.z);
}

static void printVec(const char *name, Vec3 v) {
  printf("%-8s (%.4f, %.4f, %.4f)\n", name, v.x, v.y, v.z);
}

int main() {
  // 120 degrees around (1, 1, 1): x -> y -> z -> x
  Quat q = Quat::fromAxisAngle({1, 1, 1}, radians(120.0f));
  printQuat("q", q);
  printVec("q x", rotate(q, Vec3{1, 0, 0}));

  // compare with Rodrigues' matrix for a less friendly axis and angle
  Vec3 axis{1, 2, 3};
  float angle = 0.7f;
  Vec3 v{0.3f, -1.0f, 2.0f};
  Mat4 r = rotate(axis, angle);
  Vec4 a = r * Vec4{v.x, v.y, v.z, 0.0f};
  Vec3 b = rotate(Quat::fromAxisAngle(axis, angle), v);
  printVec("matrix", {a.x, a.y, a.z});
  printVec("quat", b);

  Mat4 m = toMat4(Quat::fromAxisAngle(axis, angle));
  float maxDiff = 0.0f;
  for (int row = 0; row < 4; ++row)
    for (int col = 0; col < 4; ++col)
      maxDiff = std::fmax(maxDiff, std::fabs(m.at(row, col) - r.at(row, col)));
  printf("toMat4 == rotate: %s\n", maxDiff < 1e-5f ? "yes" : "no");

  // combining: first 90 degrees around X, then 90 degrees around Y
  Quat q1 = Quat::fromAxisAngle({1, 0, 0}, radians(90.0f));
  Quat q2 = Quat::fromAxisAngle({0, 1, 0}, radians(90.0f));
  Quat both = q2 * q1;
  printQuat("q2*q1", both);
  printVec("x", rotate(both, Vec3{1, 0, 0}));

  // slerp halfway from "no rotation" to 90 degrees around Y
  Quat half = slerp(Quat{}, q2, 0.5f);
  printQuat("half", half);
  printf("angle    %.1f degrees\n", degrees(2.0f * std::acos(half.w)));

  // drift: apply a tiny rotation a million times
  Quat step = Quat::fromAxisAngle({0, 1, 0}, 0.001f);
  Quat drift{}, fixed{};
  for (int i = 0; i < 1000000; ++i) {
    drift = step * drift;
    fixed = normalize(step * fixed);
  }
  auto len = [](Quat x) { return std::sqrt(x.w * x.w + x.x * x.x + x.y * x.y + x.z * x.z); };
  printf("length without normalize %.6f, with %.6f\n", len(drift), len(fixed));
  return 0;
}

期待される出力：

q        (w 0.5000, x 0.5000, y 0.5000, z 0.5000)
q x      (-0.0000, 1.0000, 0.0000)
matrix   (1.5069, -0.8098, 1.4709)
quat     (1.5069, -0.8098, 1.4709)
toMat4 == rotate: yes
q2*q1    (w 0.5000, x 0.5000, y 0.5000, z -0.5000)
x        (0.0000, 0.0000, -1.0000)
half     (w 0.9239, x 0.0000, y 0.3827, z 0.0000)
angle    45.0 degrees
length without normalize 1.002127, with 1.000000

$(1,1,1)$周りの120°のクォータニオンは13.5節の通り$(0.5,0.5,0.5,0.5)$、組み合わせた回転は$(0.5,0.5,0.5,-0.5)$、slerpの真ん中は13.6節の45°回転です。
ずれの行の細かい桁はコンパイラーやCPUによって違いますが、normalizeしないと長さは必ず1から離れて行きます。
（-0.0000はまた負のゼロです。第11章を参照して下さい。）

13.10 振り返りと此れから

此れでシリーズはおしまいです。
どれだけ進んで来たか振り返ってみましょう。
矢印と数の並び（第2章）から始めて、其れを測り（第3章）、「どれだけ同じ方向を向いているか」と「どの方向が垂直か」に答える2つの積を学びました（第4章と第5章）。
線型結合を行列に詰め込み（第6章と第7章）、其れで拡大縮小、回転、せん断、鏡映をし（第8章）、4つ目の座標を加えて移動も出来る様にし（第9章）、其の全てを元に戻す方法を学びました（第10章）。
そしてワールドにカメラを置き（第11章）、ワールドを画面に平らにし（第12章）、最後に回転をもっと上手く扱う方法を見付けました（此の章）。

では、全ての頂点シェーダーの心臓部を見て下さい：

float4 toClip(float4x4 proj, float4x4 view, float4x4 model, float3 position) {
  return mul(proj, mul(view, mul(model, float4(position, 1.0))));
}

vec4 toClip(mat4 proj, mat4 view, mat4 model, vec3 position) {
  return proj * view * model * vec4(position, 1.0);
}

float4 toClip(float4x4 proj, float4x4 view, float4x4 model, float3 position) {
  return mul(proj, mul(view, mul(model, float4(position, 1.0))));
}

全ての章が此の中に入っています。
ベクトル、同次座標のw = 1、右から左へ掛ける3つの行列の積、拡大縮小とクォータニオンから作った回転と平行移動で出来たモデル行列、lookAtのビュー行列、そして透視除算の準備をする投影行列です。
此のシリーズの前は、此の1行は魔法でした。
今は手で計算出来ます。

次はSlang、GLSL、Vulkan、OpenGLのシリーズです。
其処では此れら全てをGPUの上で使います。
C++から行列をアップロードし、内積と法線でライティングをし、カメラ、シャドウマップ、スキニングされたキャラクターを作ります。
linalg.hhは取っておいて下さい。
其方でも使います。
最後迄読んでくれて、ありがとうございました！

13.11 まとめ

• オイラー角（ヨー、ピッチ、ロール）は人間には分かり易いですが、順番に依存し、ジンバルロックが起き（ピッチ90°でヨーとロールが同じ軸周りに回る）、補間が苦手です。
• どんな回転も、単位軸$𝐤$周りの角度$\theta$の1回の回転です。ロドリゲスの公式$𝐯\cos \theta +(𝐤\times 𝐯)\sin \theta +𝐤(𝐤\cdot 𝐯)(1-\cos \theta )$でベクトルを回し、其の行列がrotate(axis, angle)です。
• クォータニオンは$i^2=j^2=k^2=ijk=-1$を満たす$w+xi+yj+zk$です。ハミルトン積は$(w_a{w}_b-{𝐯}_a\cdot {𝐯}_b,w_a{𝐯}_b+w_b{𝐯}_a+{𝐯}_a\times {𝐯}_b)$で、交換法則を満たしません。
• $𝐤$周りの$\theta$回転は$q=(\cos \frac{\theta }{2},𝐤\sin \frac{\theta }{2})$で、$qp{q}^{\ast }$として使います。挟み込みで$q$が2回掛かるので、半分の角度になります。
• $q_2{q}_1$は「まず$q_1$、次に$q_2$」です。単位クォータニオンでは共役が逆元です。$q$と$-q$は同じ回転です。
• toMat4でGPU用の行列に変換します。丸め誤差のずれと戦う為、定期的に正規化します。
• lerpは単位球面から外れ、nlerpは安いが不均等、slerpは均等です。slerpはなら$b$を反転し、$a\cdot b>0.9995$ならnlerpに切り替えます。
• シェーダーはクォータニオンを(x, y, z, w)で保存しますが、C++のQuatは(w, x, y, z)です。

13.12 課題

1. ★ $(1+2i)(3+j)$と$(3+j)(1+2i)$を、掛け算の表とハミルトン積のベクトルの形の両方で手計算して下さい。順番を入れ替えると答えのどの部分が変わりますか？其れは何故ですか？
2. ★ X軸周りに180°の単位クォータニオンを書き、$𝐯=(0,1,0)$を$qp{q}^{\ast }$で手計算で回して下さい。第8章のrotateXの結果と一致しますか？
3. ★ ロドリゲスの公式で、$𝐯=(1,2,3)$を$𝐤=(0,0,1)$周りに90°回して下さい。rotateZの行列で答えを確かめて下さい。
4. ★ $q=(0.5,0.5,0.5,0.5)$の回転軸と角度を求めて下さい。次に$-q$のtoMat4を計算し、$q$と同じ行列になる事を確かめて下さい。
5. ★★ slerpと同じ短い方の道の反転を入れたQuat nlerp(Quat a, Quat b, float t)を書いて下さい。Quat{}とY周りの180°の間で、t = 0、0.25、0.5、0.75、1の回転角度$2\arccos (w)$を度で表示して下さい。0、36.87、90、143.13、180になり、slerpでは0、45、90、135、180になる筈です。
6. ★★ 3つのfromAxisAngleのクォータニオンをY、X、Zの順で掛けたQuat fromEuler(float yaw, float pitch, float roll)を書き、其のtoMat4がrotateY(yaw)、rotateX(pitch)、rotateZ(roll)を此の順で掛けた積と等しい事を確かめて下さい。次にジンバルロックを数字で見せて下さい。（30°、90°、0°）と（50°、90°、20°）は、どちらもクォータニオン$(0.6830,0.6830,0.1830,-0.1830)$を表示する筈です。
7. ★★★ ハミルトン積のベクトルの形を使って、純粋なクォータニオン$p=(0,𝐯)$と単位クォータニオン$q$について、$qp{q}^{\ast }$のスカラー部が常に0になる事を示して下さい。次に、$(q_2{q}_1{)}^{\ast }=q_1^{\ast }{q}_2^{\ast }$こそが「$q_2{q}_1$は$q_1$が先」を成り立たせている理由を説明して下さい。