はじめに
クォータニオンの挟み込みで回転が表現されるのを自分なりに説明してみる試み。
複素数の場合
クォータニオンのベースになってるのは複素数で、これが理解できていないとクォータニオンも理解できないと思う。$i,j,k$はいずれも$1$と一緒に複素数体と同型な部分体を作るので。そこでまず複素数の方で考える。ベクトル$(a,b)$を複素数平面上で考える。
v=a+bi.
これの原点を中心とする角度$\theta$の回転は、
\zeta(\theta) = \cos\theta + i\sin\theta,~~~~~v \mapsto \zeta(\theta) v = (a\cos\theta - b\sin\theta) + i(a\sin\theta + b\cos\theta)
で表現される。実際、これにより$(1,0)$は$(\cos\theta,\sin\theta)$に移る。高校で習うはず...今やってるかどうかは知らない。そもそも複素数についても既知としたいところ。ざっくり言うと、
i^2=-1
を満たす$i$について$a+bi$（$a,b$は実数）の全体を考えて、文字列計算で、上記の式を適用して、常に$a+bi$という形になるよう変形する、そういう体系である。体になる。体と言ってもわかんないだろうからざっくり言うと足し算、引き算、掛け算、$0$でない数での割り算ができるということ。
なお、ゼータとしたのはこれが$z$軸周りの回転をイメージしてるからである。このあと同じことを$3$次元で考える。
クォータニオン
クォータニオン、というか四元数とは、$1,i,j,k$の実数係数の線形結合の全体に、
i^2=j^2=k^2=-1,~~~~ij=k=-ji,~~~~jk=i=-kj,~~~~ki=j=-ik
によって積を入れたものである（すごくざっくり）。例えばこんな感じで計算する。$a,b,c,a',b'$は実数とする。
\begin{align}
& ~~~~~(a+bi+ck)(a'+b'j) \\
&=aa' + ba'i + ca'k + ab'j + bb'ij + cb'kj \\
&= aa' + ba'i + ca'k + ab'j + bb'k -cb'i \\
&= aa' + (ba'-cb')i + ab'j + (ca'+bb')k.
\end{align}
定義を見てもわかる通り、非可換である。つまり、一般には$qq' = q'q$とは限らない（$q,q'$はクォータニオン）。
重要なことは、このうち$1$と$i$だけからなる四元数の全体が、すなわち複素数体と一緒（体論的には「同型」）であること。もちろん$1$と$j$でもよく、$1$と$k$でもよい。
ベクトルをクォータニオンで表す
3次元のベクトルをクォータニオンであらわすことができる。それには、$(a,b,c)$だとして、
v=ai+bj+ck
とする。純虚数のようなものである。純四元数？（そんな感じ）
共役も一応触れておく。
q = w+xi+yj+zk,~~~~~~\bar{q} = w-xi-yj-zk.
q\bar{q} = \bar{q}q = w^2+x^2+y^2+z^2
はクォータニオン$q$のノルムと言われる。単純計算で示せるので、練習だと思って実行してみるといいと思う。これが複素数のノルムに相当するようなものである。先ほどのベクトルも、ノルムを取れる。で、
v\bar{v}=\bar{v}v = a^2+b^2+c^2
となるわけ。これはベクトルのノルムである。
z軸周りの回転
クォータニオンが複素数っぽいなら、こんな感じで各軸の周りの回転を定義できそうな気がする。
\begin{align}
\xi(\theta) = \cos\theta + i\sin\theta, \\
\eta(\theta) = \cos\theta + j\sin\theta,\\
\zeta(\theta) = \cos\theta + k\sin\theta.
\end{align}
そんで、これを掛けるのが回転だと考えるのは、とても自然。実際、$(a,b,0)$に相当する$ai+bj$を考えると、$j=ki$に注意すれば、
\begin{align}
\zeta(\theta)(ai+bj) &= \zeta(\theta)(a+bk)i = (\cos\theta+k\sin\theta)(a+bk)i \\
&= (a\cos\theta-b\sin\theta + k(a\sin\theta+b\cos\theta))i \\
&= (a\cos\theta - b\sin\theta)i + (a\sin\theta + b\cos\theta)j.
\end{align}
となる。上の方でやった計算と同じですね。なんだ回転これでいいじゃんと思うでしょう。
ダメです。
今取り扱っているのは3次元のベクトルなので、これじゃまずい。どういうことかというと、ベクトル$(a,b,c)$に対して$z$軸周りの回転を適用した結果というのは、
\begin{pmatrix}
\cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}a \\ b \\ c \end{pmatrix} =
\begin{pmatrix} a\cos\theta-b\sin\theta \\ a\sin\theta+b\cos\theta \\ c \end{pmatrix}
なんだから、$(a,b,c)$に相当する$ai+bj+ck$に対して$\zeta(\theta)$を適用した結果が
\zeta(\theta)(ai+bj+ck) = (a\cos\theta - b\sin\theta)i + (a\sin\theta + b\cos\theta)j + ck
になるんだったらいいでしょう。しかし実際は...
\zeta(\theta)(ai+bj+ck) = -c\sin\theta + (a\cos\theta - b\sin\theta)i + (a\sin\theta + b\cos\theta)j + kc\cos\theta
となってしまうのでだめですね。単純な掛け算では、実部が$0$になってくれない。
困ったね。
逆から掛ける
四元数は複素数体と違って非可換なので、逆から掛けるという選択肢が存在する。交代規則：
ij=k=-ji,~~~~jk=i=-kj,~~~~ki=j=-ik
によれば、
\zeta(\theta)i = i\cos\theta + ki\sin\theta = i\cos\theta-ik\sin\theta = i\zeta(-\theta).
同様に、
\zeta(\theta)j = j\zeta(-\theta),~~~\zeta(\theta)k = k\zeta(\theta).
つまり$\zeta(\theta)$は$i,j$をかざすと角度がマイナス、というか共役になり、$k$をかざす場合はそのままである。この性質を使って、左から$\zeta(\theta)$を掛ける代わりに、左から$\zeta(\theta/2)$を掛けて、右から$\zeta(-\theta/2)$を掛ける、という表現を考えることができる。そうすれば、最初に考えた$ai+bj$の場合、$i,j$をかざすと符号が逆になるので、
\zeta(\theta/2)(ai+bj)\zeta(-\theta/2) = \zeta(\theta/2)\zeta(\theta/2)(ai+bj) = \zeta(\theta)(ai+bj)
で結果は同じである。この方法の偉い所は、$ck$がひっついていてもびくともしないところである。
\zeta(\theta/2)(ck)\zeta(-\theta/2) = \zeta(\theta/2)\zeta(-\theta/2)(ck) = \zeta(0)(ck) = ck.
そういうわけで、片側から$\zeta(\theta)$を掛けるのではなく、左から$\zeta(\theta/2)$を掛けて、右から$\zeta(-\theta/2)$を掛けるのが正しいように思える。これでベクトルの回転を表現できた。$z$軸の場合、だけ。
もちろん同じことを$x$軸と$y$軸に対しても行うことができる。すなわち$x$軸周りの回転が欲しければ、
\xi(\theta/2)(ai+bj+ck)\xi(-\theta/2)
を計算すればいい。$\xi(\theta)$は$i$をかざしても変化しないが、$j,k$に対しては符号が変化する。そして$k=ij$だから$(b+ci)j$に左から掛ける形になり、これは回転である。最後に$\eta$の場合、これも同様である。そういうわけで、ベクトル$ai+bj+ck$の$x,y,z$軸周りの角度$\theta$の回転はそれぞれ、
\begin{align}
\xi(\theta/2)(ai+bj+ck)\xi(-\theta/2),\\
\eta(\theta/2)(ai+bj+ck)\eta(-\theta/2),\\
\zeta(\theta/2)(ai+bj+ck)\zeta(-\theta/2)
\end{align}
で表現できることが分かった。
一般の軸周りの回転
いよいよクライマックス、一般のベクトルの周りの回転に移る。記号法として、
\begin{align}
\xi_\theta = \cos(\theta/2)+i\sin(\theta/2), \\
\eta_\theta = \cos(\theta/2)+j\sin(\theta/2),\\
\zeta_\theta = \cos(\theta/2)+k\sin(\theta/2)
\end{align}
とおく。これの意味はもう分かると思う。すなわち、ベクトル$p=ai+bj+ck$（$a,b,c$は実数）に対し、
\xi_\theta p \xi_{-\theta},~~~~\eta_\theta p \eta_{-\theta},~~~~\zeta_\theta p \zeta_{-\theta}
がそれぞれ$p$の$x$軸、$y$軸、$z$軸周りの角度$\theta$の回転を施されたベクトルになるというわけ。
今、単位ベクトル$e$を考える。$e$の成分は敢えて明示しない。実は不要である。$e$は単位ベクトルであるから、同じく単位ベクトルである$(0,0,1)$、というか$k$に対する、上記の軸周りの回転の連続適用で作れるはず。今それが、まず$x$軸周りに$\phi$だけ回転させた後、$y$軸周りに$\psi$だけ回転させて得られるとする。球面上の点を想像すればわかる。これを式で書くと、
e = \eta_\psi \xi_\phi k \xi_{-\phi} \eta_{-\psi}
となる。$k$を$e$に重ねられるなら、どんな二本の軸でもいい。ベクトル$p$を$e$のまわりに角度$\theta$だけ回転させたいとする。そのために次のように考える。まず、$p$と$e$に同時に軸周りの回転を適用して、$e$が$k$と重なるようにする。つまり、$p$に$y$軸周りの$-\psi$の回転をしたあと、$x$軸周りの$-\phi$の回転をする。そうして$e$が$k$と重なったら、そのときの$p$を$z$軸の周りに$\theta$だけ回転させる。そのうえで、再び$k$が$e$に重なるように逆回転を施す。つまり、$x$軸の周りに$\phi$だけ回転させた後、$y$軸の周りに$\psi$だけ回転させる。このとき、$p$は正しく$e$の周りに角度$\theta$だけ回転しているはずである。式で書くと、こう。変換後を$p'$として、
p' = (\eta_\psi \xi_\phi \zeta_\theta \xi_{-\phi} \eta_{-\psi})\cdot p \cdot (\eta_\psi \xi_\phi \zeta_{-\theta} \xi_{-\phi} \eta_{-\psi})
となる。ところで、共役操作は積の順序を変える。たとえば$ij$の共役は$-k$なので$(-j)(-i)$である。すなわち、
q = \eta_\psi \xi_\phi \zeta_\theta \xi_{-\phi} \eta_{-\psi}
とおくと、
p' = q\cdot p \cdot \bar{q}
である。ところが、
\zeta_\theta = \cos(\theta/2) + k\sin(\theta/2)
であるゆえ、
\begin{align}
q &= (\eta_\psi \xi_\phi)(\cos(\theta/2) + k\sin(\theta/2))(\xi_{-\phi} \eta_{-\psi}) \\
&= \cos(\theta/2) + \sin(\theta/2)(\eta_\psi \xi_\phi k \xi_{-\phi} \eta_{-\psi}) \\
&= \cos(\theta/2) + \sin(\theta/2)e.
\end{align}
これが適用される。
成分表示
実際に適用してみる。それには真ん中の$p$として$i,j,k$をそれぞれ採用すればいい。そうすれば成分が自然に表れる。線型性とはそういうものである。実際に確かめるため、
e=xi+yj+zk
とおく。このとき、
q=\cos(\theta/2) + (\sin(\theta/2))(xi+yj+zk)
である。$qi\bar{q},~~qj\bar{q},~~qk\bar{q}$をそれぞれ見ればいい。簡単のため、
C=\cos(\theta/2),~~~~S=\sin(\theta/2)
とおくと、
\begin{align}
qi\bar{q} &= (C + S(xi+yj+zk))\cdot i \cdot (C - S(xi+yj+zk)) \\
&= (C + S(xi+yj+zk)) \cdot (C - S(xi-yj-zk)) \cdot i \\
&= (C+Syj+Szk + Sxi)(C+Syj+Szk -Sxi)i \\
&= (C^2-S^2y^2-S^2z^2 +S^2x^2)i + (2CSy -2S^2xz)ji + (2CSz+2S^2xy)ki \\
&= (C^2-S^2+2S^2x^2)i + (2CSz + 2S^2xy)j + (2S^2xz - 2CSy)k.
\end{align}
計算では$i,j,k$の規則を用いました。$i$の前の計算において$i$の項が消滅するところでは、たとえば$jk+kj=0$を用いているのでその辺りは注意が必要でしょう。ここで、
C^2-S^2=cos\theta,~~~~2S^2=1-\cos\theta,~~~~2CS = \sin\theta
なので、結局どうなるかというと、
任意の行列式が1の直交行列を回転行列で表現する話（計算編）
の冒頭の式の一列目が現れます。ようやくつながった！
\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}.
二列目、三列目も$j,k$で同様に計算すると得られるので演習問題とします。
なお勘違いしやすいのですが、四元数では可換規則が成り立たないので、
(a+b)(a-b)=a^2-b^2
という、おそらく中学高校以来お世話になっている等式が（必ずしも）成立しません。きちんと順序を考えて計算しましょう。
そういうわけで、この変換は回転行列と等価なのでした。
残りもやってみる。
\begin{align}
qj\bar{q} &= (C + S(xi+yj+zk))\cdot i \cdot (C - S(xi+yj+zk)) \\
&= (C + S(xi+yj+zk))\cdot (C - S(-xi+yj-zk)) \cdot j \\
&= (C+Sxi+Syj+Szk)(C+Sxi-Syj+Szk)j \\
&= (C^2 - S^2x^2 + S^2y^2 - S^2z^2)j + (2CSx + 2S^2yz)ij + (2CSz - 2S^2xy)kj \\
&= (2S^2xy - 2CSz)i + (C^2-S^2 + 2S^2y^2)j + (2CSx + 2S^2yz)k.
\end{align}
できてますね。じゃあ$qk\bar{q}$よろしくお願いします...
単位クォータニオン
ノルムが$1$のクォータニオンのことを単位クォータニオンと呼びます。$q=w+xi+yj+zk$として、
q\bar{q} = w^2+x^2+y^2+z^2 = 1.
このうち実部$w$が$0$以上のものを考えると、$w=\cos(\theta/2)$（ただし$-\pi\leq\theta\leq\pi$）と書けて、
x^2+y^2+z^2 = 1-w^2 = \sin^2(\theta/2)
となるので、
x=x'\sin(\theta/2),~~~~y=y'\sin(\theta/2),~~~~z=z'\sin(\theta/2),~~~~{x'}^2+{y'}^2+{z'}^2=1
と書けます。これは先ほどの形式なので、実部が$0$以上の単位クォータニオンで軸周りの回転を表現できる、という「短絡」が成立します（軸ベクトル$(x',y',z')$のまわりの角度$\theta$の回転）。短絡は使いどころを間違えると混乱の原因になりやすいので、きちんと理解する必要があります。
（※もちろん$w$が$1$の場合は普通に$\theta$は$0$ですし、$x',y',z'$を取るような操作は存在しません（そもそも軸が無い）。念のため...クォータニオンの場合そういうことを考えなくても普通にすべて$0$なので分かりやすくていいですね。後述しますが、クォータニオンは$0$の周りを扱うのに相性がいいです。逆に$\pm\pi$のまわりについてきちんと扱うことに長けているのが回転行列の形式です。どちらか一方がソリッドに決まるからです。）
正規直交基底との関係
軸周りの回転をあらわすクォータニオン$q$に対し、上で述べたように、
qi\bar{q},~~~qj\bar{q},~~~qk\bar{q}
は正規直交基底になります。逆に正規直交基底、というか行列式が1の直交行列は（列ベクトルとして正規直交基底を扱っている）、以前の記事により、回転行列の形式で表現できます。そしてそこから出る角度と軸ベクトルでクォータニオンを作れば、それによる挟み込み積でベクトルへの作用を同じ回転行列として表現できます。
今、正規直交基底とそれをあらわすクォータニオン$q$があるとき、正規直交基底のすべてのベクトルを軸ベクトル$u$の周りに角度$\phi$だけ回転させることを考えると、クォータニオン$s=\cos(\phi/2)+\sin(\phi/2)u$を作れば（$u$は純四元数として扱っています）、変換後の各ベクトルは、
sqi\bar{q}\bar{s},~~~~sqj\bar{q}\bar{s},~~~~sqk\bar{q}\bar{s}
となります。これに相当するクォータニオンが$sq$なのは明らかです。すなわち、$s$を左から掛けることで、正規直交基底の回転を表現できることが分かります。掛け算は結合律が成立するので、複数回掛ける場合は複数回掛けたものを扱えばいいです。
ところで、いちいち角度や軸ベクトルを出さなくても、直接成分からベクトルを出せた方が便利なので、それを計算しておきます。やり方は簡単で、さっきの計算結果で$C,Sx,Sy,Sz$をそれぞれ$w,x,y,z$で置き換えればダイレクトに得られますね。最後の
qk\bar{q} = (2S^2xz + 2CSy)i + (2S^2yz - 2CSx)j + (C^2-S^2+2S^2z^2)k
も追加します。また、$C^2-S^2=2C^2-1$としておきます。これでクォータニオン$q$に相当する正規直交基底が、
\begin{align}
(2w^2-1+2x^2, ~~2(xy+zw), ~~2(xz-yw)) \\ (2(xy-zw),~~2w^2-1+2y^2,~~2(yz+xw)) \\ (2(xz+yw),~~2(yz-xw),~~2w^2-1+2z^2))
\end{align}
となることが分かりました。計算量はともかく、これで「クォータニオンの乗算」と「そこからのベクトル抽出」という形で軸周りの正規直交基底の回転を整理できることが分かりました。
正規直交基底からクォータニオンを出す
これも可能ですね。やるには成分を使えばいいです。まず、行列式が1の直交行列を
V=\begin{pmatrix}a & b & c \\ d & e & f \\ g & h & i\end{pmatrix}
とおいて、
\begin{pmatrix}2w^2-1+2x^2 & 2(xy-zw) & 2(xz+yw) \\ 2(xy+zw) & 2w^2-1+2y^2 & 2(yz-xw) \\ 2(xz-yw) & 2(yz+xw) & 2w^2-1+2z^2\end{pmatrix}
と比較します。さっきのベクトルを列ベクトルで並べただけです。対角成分の和が$x^2+y^2+z^2+w^2=1$により$4w^2-1$となるので、
w = \sqrt{(a+e+i+1)/4}
で定めればいいです。あとは成分から、
x=\frac{h-f}{4w},~~~~~y=\frac{c-g}{4w},~~~~~z=\frac{d-b}{4w}
とすればいいです。これらについて$x^2+y^2+z^2+w^2=1$が成り立つことを確かめる必要があります。
\begin{align}
&~~~~16w^2(x^2+y^2+z^2+w^2-1)\\
&=(h-f)^2 + (c-g)^2 + (d-b)^2 - 16w^2(1-w^2) \\
&= h^2+f^2+c^2+g^2+b^2+d^2-2(hf+cg+bd) - (1+a+e+i)(3-a-e-i) \\
&= h^2+f^2+c^2+g^2+b^2+d^2-2(hf+cg+bd) - (3-a^2-e^2-i^2+2(a+e+i)-2(ae+ei+ia)) \\
&= -2(hf+cg+bd) -2(a+e+i) + 2(ae+ei+ia) \\
&= 0.
\end{align}
余因子の式などを使いました。これで単位クォータニオンになりますね！めでたし、めでたし、...
うそです。
はじめに、根号内が負にならないことですが、これは前回の記事で確かめたような内容が根拠になっています。$p^2+q^2=1$を示したと思いますが、あれにより結局$a+e+i$は$-1$から$3$の間の値を取ることが分かるからです。具体的には$(a+e+i-1)/2$をニ乗すればわかるかと思いますが冗長なのでやりません（※固有値を使えば見通しよく証明できます...）。
次に、成分比較をやっていません。$x,y,z$の定義から検証すべき箇所は6か所ですが、どれも前回の記事でやったような余因子の法則に基づいた単純計算です。そこまでやってしまうとさすがに冗長の極みなので、読者にお任せします（演習問題！）。
本題ですが、$w$が$0$になる場合、コーディング的には$0$に非常に近い場合なども含みますが、この議論ができないです。それは角度が$\pi$に近い場合です。割り算ができません。以前の記事で直交行列からダイレクトに角度と軸を出しましたが、あの方法は逆に単位行列の近くでうまくいかないです（軸が定まらないので）。クォータニオン的には、積で扱うので、これは問題になりません。なので、一番最適な方法は両者を組み合わせることです。つまり、$w$に相当する値が$0$に近くなりそうなら角度と軸を出すやり方で扱い、$1$に近くなりそうならこっちの方法でダイレクトアタックすれば完璧ですね。
コメントがありましたが、Threeはまさにその方針でクォータニオンを求めているようです。
https://github.com/mrdoob/three.js/blob/r172/src/math/Quaternion.js#L294
まずトレースが正の場合、としていますが、これはこの記事でいうところの$w$が$0$にならず安心して計算できる場合を示しており、そのように計算されています。見ての通り、全く同じ計算です。次にそうでない場合ですが、その場合は以前書いた記事のように対角成分のうちどれが一番大きいかということを考えています。あの記事ではクォータニオンまでは出しませんでしたが、データを使えばこのソースになるようなやり方で出せます。
$p^2+q^2=1$を使います。$\epsilon = \mathrm{sgn}(h-f)$（符号）とおきます。$q=\epsilon\sqrt{1-p^2}$なので、
\begin{align}
\cos\frac{\theta}{2} &= \sqrt{\frac{p+1}{2}} = \sqrt{\frac{1-p^2}{2(1-p)}} = \frac{\epsilon q}{\sqrt{2}\sqrt{1-p}} \\
&= \frac{\epsilon(h-f)}{\sqrt{2}\sqrt{1-p}\cdot 2x} = \frac{\epsilon(h-f)}{2\sqrt{2}\sqrt{a-p}}.
\end{align}
ここで最後に出てきた式の分母が先ほどのソースにある$s$です。これで$w$を求めたことになっており、また他の成分については、以前書いた記事における$x,y,z$に$\sin(\theta/2)$を掛ければ得られるので、たとえば
\sin\frac{\theta}{2} = \epsilon\sqrt{\frac{1-p}{2}},~~~~x\sin\frac{\theta}{2} = \epsilon\sqrt{\frac{a-p}{1-p}}\cdot \sqrt{\frac{1-p}{2}} = \epsilon\sqrt{\frac{a-p}{2}} = \epsilon\frac{s}{4}
となります。$y,z$も同様です。これをすべての場合に計算しているわけですね。繰り返しになるかもしれませんが、単位行列になってしまうケースはトレースが正でないという前提から該当しないので、安心して計算できます。
全体に符号$(\epsilon)$が付いていますが、これが負の場合であっても全体がマイナスになるだけで、クォータニオンとしての挙動は同じなので、これでよしとします（Threeの方でも$w$が負になるかどうかについては特に気にしていないようです。こちらでは角度範囲の都合上必然的に$\cos(\theta/2)$は$0$以上になりますし、そのように全体の符号が決まっています。）。
しかし、クォータニオンを算出するのが目的であれば、$x^2,y^2,z^2$のうち最大のものを取ればそれは正であることが保証される、という方針でやった方が見通しよく求められます。ここでは前回の結果を再利用したかったので、そのようにしました。式については成分をよく見れば結果的にソースと同じ式が出てきます。その場合の検証については、読者にお任せします。
長くなりましたが、単に場合分けするだけなら前回の記事のように角度とベクトルを各場合に出してから適切に半角だの掛け算だのやって終わりです。ここで追記した意図はThreeの計算式の説明のためで、ただのおまけだと思っていただけると幸いです。それでもいい加減な事はできないので、厳密にやらせていただきました。
右から掛ける場合
単位クォータニオンで正規直交基底を表す場合、これに左から掛け算すると、すべての軸がそのクォータニオンのあらわす角度、軸に従って回転します。これはグローバルな回転と呼ぶことができます。実は右から掛ける場合にも意味があります。言ってしまうと、この場合はローカルな回転となります。たとえば$\zeta_\theta$を使う場合、左から掛ける場合はグローバルの$z$軸の周りに角度$\theta$だけ回転しますが、右から掛ける場合、ローカルの$z$軸（要するに$qk\bar{q}$）の周りの角度$\theta$の回転となります。これを見るには、単位ベクトル$e$について、
q(\cos\frac{\theta}{2} + e\sin\frac{\theta}{2}) = (\cos\frac{\theta}{2} + \sin\frac{\theta}{2}\cdot (qe\bar{q}))q
という計算をすればいいです。$qe\bar{q}$はローカル空間におけるベクトル$e$なので、その周りに回転するというわけですね。
以上より、正規直交基底に対して、それを表現するクォータニオンを用意すると、それのグローバル、ないしはローカルな軸周りの回転を、然るべきクォータニオンを用いて実行できることが分かりました。さらにその結果から正規直交基底は容易に算出できます。クォータニオン、とても便利ですね。
実部が負のクォータニオン
単位クォータニオンにおいて実部が負の場合、実部を$\cos(\theta/2)~~(-\pi\leq \theta \leq\pi)$の形で書くことができません。
この辺はちょっと理解が足りていません。事実を述べると、まず単位クォータニオン$q$について、$q$と$-q$は相似変換における作用が全く同一です（相殺するので）。次に、$q$でも、$-q$でも、正規直交基底の計算の結果は変わりません（二次の項しかないので）。そういう事情から、おそらく適宜$q$を$-q$ですりかえたとしても特に問題ないように思われます。クォータニオンの重要な仕事である相似変換と正規直交基底の算出において同じ挙動をするのであれば、符号については特に問題ないということでしょうか。複素数でもそうですが、実部が$0$以上の元の積もまた実部が$0$以上であることは保証されていません。なので、掛け算をするたびに、そこら辺の調整が必要かもしれないです。あるいは、
\cos((\theta+2\pi)/2)=-\cos(\theta/2),~~~~~\sin((\theta+2\pi)/2)=-\sin(\theta/2)
ですから、マイナスを掛けることは角度を$2\pi$ずらすことに相当します。挙動が同じなのはそういうわけです。どうやら適宜$-q$で置き換えても問題なさそうです。
（追記：先ほどThreeのQuarternionのslerpの実装：
Quaternion.slerp()
を斜め読みしていたのですが、ある場合にマイナスのクォータニオンで取り換える処理をしています。きちんと理解していないのですが、ともかく「$q$と$-q$の挙動は一緒である」という事実を使っていることは間違いなさそうです。）
もしかすると単位クォータニオンのなす群を$1,-1$の作る群で割ったものを考えるといいのかもしれません。積で作用するので、群として考えても差し支えなく、$-1$の作用が恒等なので、割っても問題ないですね。そこら辺どうなんだろう...（分からない人は読み飛ばしてください）
slerpの計算
次いでと言っては何ですが、slerpについても扱っておきます。自分は計算馬鹿なので、計算ベースの説明にはなってしまいますが、お付き合いください。厳密ではないです。あくまで概念的な説明です。
まず$q_0,q_1$を単位クォータニオンとします。$q=q_1{q_0}^{-1}$とおくとき、これの実部は$0$以上とは限りません。そこで必要なら$q_0$を$-q_0$で置き換えて、実部が$0$以上となるようにしておきます。この単位クォータニオンについて、
q=\cos\frac{\theta}{2}+e\sin\frac{\theta}{2}
と書きます。なお$e$は単位ベクトルです。このとき補間は、
q(t) = \cos\frac{t\theta}{2}+e\sin\frac{t\theta}{2}~~~(0\leq t \leq 1)
となります。$q(0)$を$q_0$に掛ければ$q_0$になり、$q(1)$を$q_0$に掛ければ$q_1$になります（もちろん左から）。最終的に欲しいのは、$q(t)q_0$というわけです。ここで、
q(t)q(1-t)=q(1-t)q(t)=q(1),~~~q(1)q_0=q_1
ですから、逆元が共役であることより、
q(t)q_0 = (\cos\frac{(1-t)\theta}{2} - e\sin\frac{(1-t)\theta}{2})q_1
が成立します。$q_t$で結果を表現すると、
\begin{align}
q_t &= (\cos\frac{t\theta}{2} + e\sin\frac{t\theta}{2})q_0, \\
q_t &= (\cos\frac{(1-t)\theta}{2} - e\sin\frac{(1-t)\theta}{2})q_1
\end{align}
というわけですね。ここでふたたび逆元が共役であることを用いると、
\begin{align}
q_0 &= (\cos\frac{t\theta}{2} - e\sin\frac{t\theta}{2})q_t, \\
q_1 &= (\cos\frac{(1-t)\theta}{2} + e\sin\frac{(1-t)\theta}{2})q_t
\end{align}
となります。$e$の項が邪魔ですね。消えないかと思うでしょう。消えないかと願うのではなく、消えるように計算するのです。計算にはより良い見通しを持とうとする意志が大切です。
\begin{align}
&~~~~~(\sin\frac{(1-t)\theta}{2})q_0 + (\sin\frac{t\theta}{2})q_1 \\
&= ((\cos\frac{t\theta}{2})(\sin\frac{(1-t)\theta}{2}) + (\cos\frac{(1-t)\theta}{2})(\sin\frac{t\theta}{2}))q_t \\
&= q_t\sin\frac{\theta}{2}.
\end{align}
これを除算で調整すると例の式ができます。もちろん$\sin\frac{\theta}{2}$が$0$の場合、というか$q$が1の場合は除算ができないので、場合分けで工夫することになります。この議論では最終的に$e$が消えます。$e$は何かというと、$q_0{q_1}^{-1}$のダイレクトな計算で現れるわけですね。その負荷を軽減するための、この工夫です。なおThreeでは角度の算出にはサインの方も計算してatan関数を使っています。サインは平方を成分から出して平方根を取っていますね。
（追記：あの後きちんと読んだのですが、サインについてはコサインから公式で出してますね...つまり外積の計算をしていません。そしてそれが非常に小さい場合には（判定にはNumber.EPSILONを使っていますね）線形補間で出しています。意地でも厳密な$q_0{q_1}^{-1}$のダイレクト計算はしないぞという意思を感じます。そのうえでatanにもっていってるので実質acosとあんま変わんないかもしれないですね。）
（すごく細かいことを言うと、cosとsinが両方得られる場合、両方の情報を使う方が角度の精度は高いです。片方だとどうしても偏ってしまいます。しかしそこまで問題にはならないのでまあ計算スピードの方を優先しているんだと思います。トレードオフというやつですね。）
これも商群で考えれば、符号が負にならないように代表を取って、それにより作られる元の同値類を得る、という考え方でもいいのかもしれないですね。作用は同じですからね。同値類同士の演算ではよくあることです。
おわりに
以上、お話として単位クォータニオンによるサンドイッチが回転行列と等価になることの簡単な説明でした。ここまでお読みいただいてありがとうございました。