上篇文章 介绍了使用欧拉角表示定向，但是欧拉角有万向节死锁问题，我们可以使用的四元数来解决这个问题。四元数比较复杂，在介绍四元数之前让我们来看看和它相关的一些内容。
# 根据欧拉旋转定理，任何三维旋转都可以找到一个轴 n ^ \hat{n} n ^ θ \theta θ n ^ \hat{n} n ^ θ \theta θ
# 上面轴-角的旋转轴使用的是单位矢量，我们可以让旋转角度乘上旋转轴而不丢失任何信息，来产生一个矢量 e = θ n ^ e=\theta \hat{n} e = θ n ^ θ = ∥ e ∥ \theta = \|e \| θ = ∥ e ∥ e e e
复数（Complex Number）表示为 z = a + b i z=a+bi z = a + b i a a a b b b i 2 i^2 i 2 − 1 -1 − 1
# 复数加减法是将的对应分量相加减。
z 1 ± z 2 = ( a 1 ± a 2 ) + ( b 1 ± b 2 ) i z1 \pm z2=(a1 \pm a2) + (b1 \pm b2)i z 1 ± z 2 = ( a 1 ± a 2 ) + ( b 1 ± b 2 ) i 用分配律来计算两个复数乘积。
z 1 z 2 = ( a 1 + b 1 i ) ( a 2 + b 2 i ) = a 1 a 2 + a 1 b 2 i + a 2 b 1 i + b 1 b 2 i 2 = a 1 a 2 − b 1 b 2 + ( a 2 b 1 + a 1 b 2 ) i \begin{aligned} z1z2&=(a1+b1i)(a2+b2i) \\ &=a1a2 + a1b2i + a2b1i + b1b2i^2 \\ &=a1a2-b1b2+(a2b1+a1b2)i \end{aligned} z 1 z 2 = ( a 1 + b 1 i ) ( a 2 + b 2 i ) = a 1 a 2 + a 1 b 2 i + a 2 b 1 i + b 1 b 2 i 2 = a 1 a 2 − b 1 b 2 + ( a 2 b 1 + a 1 b 2 ) i # 复数的模长等于 ∥ z ∥ = a 2 + b 2 \|z\| =\sqrt{a^2+b^2} ∥ z ∥ = a 2 + b 2 z = a − b i z=a-bi z = a − b i
∥ z ∥ 2 = ( a + b i ) ( a − b i ) = a 2 − a b i + a b i + b 2 = a 2 + b 2 \begin{aligned} \| z \|^2 &= (a+bi)(a-bi) \\ &=a^2-abi+abi+b^2 \\ &=a^2+b^2 \end{aligned} ∥ z ∥ 2 = ( a + b i ) ( a − b i ) = a 2 − a b i + a b i + b 2 = a 2 + b 2 # 我们可以将复数画在二维坐标轴上。
从上图 z = a + b i z=a+bi z = a + b i a = c o s ( θ ) ∗ ∥ z ∥ a=cos(\theta)*\|z\| a = c o s ( θ ) ∗ ∥ z ∥ b = s i n ( θ ) ∗ ∥ z ∥ b=sin(\theta)*\|z\| b = s i n ( θ ) ∗ ∥ z ∥ z = ∥ z ∥ ( c o s ( θ ) + i s i n ( θ ) ) z=\|z\|(cos(\theta)+isin(\theta)) z = ∥ z ∥ ( c o s ( θ ) + i s i n ( θ ) )
根据欧拉公式 。
e i θ = c o s ( θ ) + i s i n ( θ ) e^{i\theta}=cos(\theta)+isin(\theta) e i θ = c o s ( θ ) + i s i n ( θ )
我们令 r = ∥ z ∥ r=\|z\| r = ∥ z ∥ z = r e i θ z=re^{i\theta} z = r e i θ r r r 1 1 1 v = x + y i v=x+yi v = x + y i z = e i θ z=e^{i\theta} z = e i θ θ \theta θ
v ′ = e i θ v v'=e^{i\theta}v v ′ = e i θ v # 四元数（Quaternion）由 4 个数字组成，这也是它名字的由来。和复数一样它也是由实部和虚部组成，只不过它有 3 个虚部。
q = a + b i + c j + d k q = a + bi + cj + dk q = a + b i + c j + d k i 2 = j 2 = k 2 = i j k = − 1 i^2 = j^2 = k^2 = ijk = -1 i 2 = j 2 = k 2 = i j k = − 1 四元数是由哈密顿在 1843 年爱尔兰发现的。当时他正研究扩展复数到更高的维次（复数可视为平面上的点）。他不能做到三维空间的例子，但四维则造出四元数。他于 10 月 16 日跟他的妻子在都柏林的皇家运河（Royal Canal）上散步时突然想到 的方程解。之后哈密顿立刻将此方程刻在附近布鲁穆桥（Brougham Bridge，现称为金雀花桥 Broom Bridge）。
我们还可以将它看成是 1 个标量和 1 个矢量组成。标量用 w 表示，v 或 [x, y, z] 表示矢量。
[ w v ] [w \quad \bold{v}] [ w v ] [ w ( x , y , z ) ] [w \quad (x, y, z)] [ w ( x , y , z ) ] 它和上面介绍的轴-角的关系如下。
[ w v ] = [ c o s ( θ 2 ) s i n ( θ 2 ) ∗ n ^ ] [w \quad \bold{v}] = [cos(\frac{\theta}{2}) \quad sin(\frac{\theta}{2}) * \hat{n}] [ w v ] = [ c o s ( 2 θ ) s i n ( 2 θ ) ∗ n ^ ] [ w ( x , y , z ) ] = [ c o s ( θ 2 ) ( s i n ( θ 2 ) ∗ n x , s i n ( θ 2 ) ∗ n y , s i n ( θ 2 ) ∗ n z ) ] [w \quad (x, y, z)] = [cos(\frac{\theta}{2}) \quad (sin(\frac{\theta}{2}) * n_x, sin(\frac{\theta}{2}) * n_y, sin(\frac{\theta}{2}) * n_z)] [ w ( x , y , z ) ] = [ c o s ( 2 θ ) ( s i n ( 2 θ ) ∗ n x , s i n ( 2 θ ) ∗ n y , s i n ( 2 θ ) ∗ n z ) ] 四元数也可以变负，只要将的 4 个数字都变为负数就行 [-w (-x, -y, -z)]，比较有意思的是的正 q 和负 -q四元数表示的是相同的角位移，因为变负意味着沿轴的另一端反方向旋转。
这篇文章只是简单介绍下四元数，如果想深入了解四元数请阅读这篇文章 。通过这个视频 可视化四元数。这个网站 可以在线交互四元数。
# 纯四元数是只有虚部的四元数。因为纯四元数仅由虚部，我们可以将任意的 3D 矢量转换为纯四元数。
q = [ 0 , v ] \bold{q} = [0, \bold{v}] q = [ 0 , v ] # 和矢量一样四元数也有大小并且计算方式也一样。
∥ q ∥ = ∥ [ w ( x , y , z ) ] ∥ = w 2 + x 2 + y 2 + z 2 \| \bold{q} \| = \| [w \quad (x, y, z)] \| = \sqrt{w^2 + x^2 + y^2 + z^2} ∥ q ∥ = ∥ [ w ( x , y , z ) ] ∥ = w 2 + x 2 + y 2 + z 2 我们只关心旋转量，也就是矢量是单位矢量。
∥ q ∥ = ∥ [ w , v ] ∥ = c o s ( θ 2 ) 2 + ( s i n ( θ 2 ) n ) 2 = 1 = 1 \begin{aligned} \| q \| &= \| [w, v] \| \\ &=\sqrt{cos(\frac{\theta}{2})^2 + (sin(\frac{\theta}{2})n)^2} \\ &=\sqrt{1} \\ &=1 \end{aligned} ∥ q ∥ = ∥ [ w , v ] ∥ = c o s ( 2 θ ) 2 + ( s i n ( 2 θ ) n ) 2 = 1 = 1 我们使用四元数的目的是为了表示定向，所以我们这里使用的四元数都是长度为 1 的四元数。
# 如果一个四元数的大小是 1，那么这个四元数就被称为单位四元数。我们可以将一个四元数除以它的大小得到一个单位四元数。
U q = q ∥ q ∥ U_q = \frac{q}{\| q \| } U q = ∥ q ∥ q # 和复数一样四元数的共轭表示为 q ∗ \bold{q^*} q ∗ [ w ( − x , − y , − z ) ] [w \quad (-x, -y, -z)] [ w ( − x , − y , − z ) ]
四元数的逆表示为 q − 1 \bold{q^{-1}} q − 1
q − 1 = q ∗ ∥ q ∥ \bold{q^{-1}} = \frac{q^*}{\| q\| } q − 1 = ∥ q ∥ q ∗ 因为我们只关心旋转，旋转轴是单位矢量四元数大小也是 1，所以对于我们来说四元数的逆和共轭是相等的。
# 和复数一样，四元数的加减法只需要将相应分量相加或相减就可以了。
q 1 + q 2 = a + b i + c j + d k + e + f i + g j + h k = ( a + e ) + ( b + f ) i + ( c + g ) j + ( d + h ) k \begin{aligned} q1+q2&=a+bi+cj+dk+e+fi+gj+hk \\ &=(a+e)+(b+f)i+(c+g)j+(d+h)k \end{aligned} q 1 + q 2 = a + b i + c j + d k + e + f i + g j + h k = ( a + e ) + ( b + f ) i + ( c + g ) j + ( d + h ) k q 1 − q 2 = ( a − e ) + ( b − f ) i + ( c − g ) j + ( d − h ) k q1-q2=(a-e)+(b-f)i+(c-g)j+(d-h)k q 1 − q 2 = ( a − e ) + ( b − f ) i + ( c − g ) j + ( d − h ) k 四元数和标量相乘，结果是四元数的每个分量都乘以标量。
s q = s ( a + b i + c j + d k ) = s a + s b i + s c j + s d k \begin{aligned} sq&=s(a+bi+cj+dk) \\ &=sa+sbi+scj+sdk \end{aligned} s q = s ( a + b i + c j + d k ) = s a + s b i + s c j + s d k 四元数与四元数相乘，其结果还是四元数。
q 1 q 2 = ( a + b i + c j + d k ) ( e + f i + g j + h k ) = a e + a f i + a g j + a h k + b e i + b f i 2 + b g i j + b h i k + c e j + c f j i + c g j 2 + c h j k + d e k + d f k i + d g k j + d h k 2 \begin{aligned} q1q2&=(a+bi+cj+dk)(e+fi+gj+hk) \\ &=ae+afi+agj+ahk+ \\ &\quad\: bei+bfi^2+bgij+bhik+ \\ &\quad\: cej+cfji+cgj^2+chjk+ \\ &\quad\: dek+dfki+dgkj+dhk^2 \end{aligned} q 1 q 2 = ( a + b i + c j + d k ) ( e + f i + g j + h k ) = a e + a f i + a g j + a h k + b e i + b f i 2 + b g i j + b h i k + c e j + c f j i + c g j 2 + c h j k + d e k + d f k i + d g k j + d h k 2 我们可以利用 i j k ijk i j k
q 1 q 2 = ( a e − b f − c g − d h ) + ( b e + a f − d g + c h ) i + ( c e + d f + a g − b h ) j + ( d e − c f + b g + a h ) k \begin{aligned} q1q2&=(ae-bf-cg-dh)+ \\ &\quad\: (be+af-dg+ch)i+ \\ &\quad\: (ce+df+ag-bh)j+ \\ &\quad\: (de-cf+bg+ah)k \end{aligned} q 1 q 2 = ( a e − b f − c g − d h ) + ( b e + a f − d g + c h ) i + ( c e + d f + a g − b h ) j + ( d e − c f + b g + a h ) k 观察上面简化后的等式，我们还可以将乘法写成矩阵形式。
q 1 q 2 = [ a − b − c − d b a − d c c d a − b d − c b a ] [ e f g h ] q1q2=\begin{bmatrix} a & -b & -c & -d \\ b & a & -d & c \\ c & d & a & -b \\ d & -c & b & a \end{bmatrix} \begin{bmatrix} e \\ f \\ g \\ h \end{bmatrix} q 1 q 2 = ⎣ ⎢ ⎢ ⎢ ⎡ a b c d − b a d − c − c − d a b − d c − b a ⎦ ⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ ⎡ e f g h ⎦ ⎥ ⎥ ⎥ ⎤ 四元数乘法是不满足交换律的，如果是右乘 q 1 q1 q 1
q 2 q 1 = [ a − b − c − d b a d − c c − d a b d c − b a ] [ e f g h ] q2q1=\begin{bmatrix} a & -b & -c & -d \\ b & a & d & -c \\ c & -d & a & b \\ d & c & -b & a \end{bmatrix} \begin{bmatrix} e \\ f \\ g \\ h \end{bmatrix} q 2 q 1 = ⎣ ⎢ ⎢ ⎢ ⎡ a b c d − b a − d c − c d a − b − d − c b a ⎦ ⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ ⎡ e f g h ⎦ ⎥ ⎥ ⎥ ⎤ 如果我们将四元数看成矢量形式 v = [ b , c , d ] \bold{v}=[b,c,d] v = [ b , c , d ] u = [ f , g , h ] \bold{u}=[f,g,h] u = [ f , g , h ]
q 1 q 2 = [ a e − v ⋅ u , a u + e v + v × u ] = [ a e − b f − c g − d h , ( a f + e b + c h − d g , a g + e c + d f − b h , a h + e d + b g − c f ) ] \begin{aligned} q1q2&=[ae - v \cdot u, au+ev+v\times u] \\ &=[ae-bf-cg-dh, \\ &\quad\: (af+eb+ch-dg, \\ &\quad\: ag+ec+df-bh, \\ &\quad\: ah+ed+bg-cf)] \end{aligned} q 1 q 2 = [ a e − v ⋅ u , a u + e v + v × u ] = [ a e − b f − c g − d h , ( a f + e b + c h − d g , a g + e c + d f − b h , a h + e d + b g − c f ) ] 四元数乘法可以结合但是不可以交换，四元数乘积的大小等于，它们大小的乘积，四元数乘积的倒数等于相反顺序倒数的乘积。
( a b ) c = a ( b c ) (ab)c=a(bc) ( a b ) c = a ( b c ) ∥ a b ∥ = ∥ a ∥ ∗ ∥ b ∥ \| ab \| =\| a \| * \| b \| ∥ a b ∥ = ∥ a ∥ ∗ ∥ b ∥ ( a b ) − 1 = b − 1 a − 1 (ab)^{-1}=b^{-1}a^{-1} ( a b ) − 1 = b − 1 a − 1 # 假设有一个矢量或点，我们将它定义为纯四元数 p = [ 0 , v ] p = [0, v] p = [ 0 , v ] θ \theta θ q = [ c o s ( θ 2 ) s i n ( θ 2 ) ∗ u ^ ] q = [cos(\frac{\theta}{2}) \quad sin(\frac{\theta}{2}) * \hat{u}] q = [ c o s ( 2 θ ) s i n ( 2 θ ) ∗ u ^ ] p p p θ \theta θ p ′ p' p ′
p ′ = q p q ∗ = q p q − 1 = [ 0 , c o s ( θ ) v + ( 1 − c o s ( θ ) ) ( u ⋅ v ) u + s i n ( θ ) ( u × v ) ] p'=qpq^*=qpq^{-1}=[0, cos(\theta)v+(1-cos(\theta))(u \cdot v)u+sin(\theta)(u \times v)] p ′ = q p q ∗ = q p q − 1 = [ 0 , c o s ( θ ) v + ( 1 − c o s ( θ ) ) ( u ⋅ v ) u + s i n ( θ ) ( u × v ) ] 多个旋转也可以组合就和矩阵一样，假设现在有两个旋转四元数，先旋转 a a a b b b b a ba b a
p ′ = b ( a p a − 1 ) b − 1 = ( b a ) p ( b a ) − 1 \begin{aligned} p'&=b(apa^{-1})b^{-1} \\ &=(ba)p(ba)^{-1} \end{aligned} p ′ = b ( a p a − 1 ) b − 1 = ( b a ) p ( b a ) − 1 # 我们还可以将 v ′ = q v q ∗ v'=qvq^* v ′ = q v q ∗ v v v q ∗ q^* q ∗ q q q v ′ v' v ′ q q q q ∗ q* q ∗ q q q q = a + b i + c j + d k q=a+bi+cj+dk q = a + b i + c j + d k a 2 + b 2 + c 2 + d 2 = 1 a^2+b^2+c^2+d^2=1 a 2 + b 2 + c 2 + d 2 = 1
q v q ∗ = L ( q ) R ( q ∗ ) v = [ a − b − c − d b a − d c c d a − b d − c b a ] [ a b c d − b a − d c − c d a − b − d − c b a ] v = [ 1 0 0 0 0 1 − 2 c 2 − 2 d 2 2 b c − 2 a d 2 a c + 2 b d 0 2 b c + 2 a d 1 − 2 b 2 − 2 d 2 2 c d − 2 a b 0 2 b d − 2 a c 2 a b + 2 c d 1 − 2 b 2 − 2 c 2 ] v \begin{aligned} qvq^*&=L(q)R(q^*)v \\ &=\begin{bmatrix} a & -b & -c & -d \\ b & a & -d & c \\ c & d & a & -b \\ d & -c & b & a \end{bmatrix} \begin{bmatrix} a & b & c & d \\ -b & a & -d & c \\ -c & d & a & -b \\ -d & -c & b & a \end{bmatrix} v \\ &=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1-2c^2-2d^2 & 2bc-2ad & 2ac+2bd \\ 0 & 2bc+2ad & 1-2b^2-2d^2 & 2cd-2ab \\ 0 & 2bd-2ac & 2ab+2cd & 1-2b^2-2c^2 \end{bmatrix} v \end{aligned} q v q ∗ = L ( q ) R ( q ∗ ) v = ⎣ ⎢ ⎢ ⎢ ⎡ a b c d − b a d − c − c − d a b − d c − b a ⎦ ⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ ⎡ a − b − c − d b a d − c c − d a b d c − b a ⎦ ⎥ ⎥ ⎥ ⎤ v = ⎣ ⎢ ⎢ ⎢ ⎡ 1 0 0 0 0 1 − 2 c 2 − 2 d 2 2 b c + 2 a d 2 b d − 2 a c 0 2 b c − 2 a d 1 − 2 b 2 − 2 d 2 2 a b + 2 c d 0 2 a c + 2 b d 2 c d − 2 a b 1 − 2 b 2 − 2 c 2 ⎦ ⎥ ⎥ ⎥ ⎤ v 矩阵的最外圈不会对 v v v 3x3 旋转的部分，就是四元数 3D 旋转矩阵形式了。
[ 1 − 2 c 2 − 2 d 2 2 b c − 2 a d 2 a c + 2 b d 2 b c + 2 a d 1 − 2 b 2 − 2 d 2 2 c d − 2 a b 2 b d − 2 a c 2 a b + 2 c d 1 − 2 b 2 − 2 c 2 ] \begin{bmatrix} 1-2c^2-2d^2 & 2bc-2ad & 2ac+2bd \\ 2bc+2ad & 1-2b^2-2d^2 & 2cd-2ab \\ 2bd-2ac & 2ab+2cd & 1-2b^2-2c^2 \end{bmatrix} ⎣ ⎢ ⎡ 1 − 2 c 2 − 2 d 2 2 b c + 2 a d 2 b d − 2 a c 2 b c − 2 a d 1 − 2 b 2 − 2 d 2 2 a b + 2 c d 2 a c + 2 b d 2 c d − 2 a b 1 − 2 b 2 − 2 c 2 ⎦ ⎥ ⎤ # 我们还可以利用四元数乘法计算两个四元数之间的差值，也就是一个方向旋转到另一个方向的角位移。 d a = b da=b d a = b d d d a a a b b b
( d a ) a − 1 = b a − 1 d ( a a − 1 ) = b a − 1 d [ 1 , ( 0 , 0 , 0 ) ] = b a − 1 d = b a − 1 \begin{aligned} (da)a^{-1}&= ba^{-1} \\ d(aa^{-1})&=ba^{-1} \\ d[1,(0,0,0)]&=ba^{-1} \\ d&=ba^{-1} \end{aligned} ( d a ) a − 1 d ( a a − 1 ) d [ 1 , ( 0 , 0 , 0 ) ] d = b a − 1 = b a − 1 = b a − 1 = b a − 1 和矢量一样，四元数点积的结果也是一个标量。
q 1 ⋅ q 2 = [ w 1 v 1 ] ⋅ [ w 2 v 2 ] = w 1 w 2 + v 1 ⋅ v 2 = w 1 w 2 + x 1 x 2 + y 1 y 2 + z 1 z 2 \begin{aligned} q1 \cdot q2&=[w1 v1] \cdot [w2 v2] \\ &=w1w2 + v1 \cdot v2 \\ &=w1w2+x1x2+y1y2+z1z2 \end{aligned} q 1 ⋅ q 2 = [ w 1 v 1 ] ⋅ [ w 2 v 2 ] = w 1 w 2 + v 1 ⋅ v 2 = w 1 w 2 + x 1 x 2 + y 1 y 2 + z 1 z 2 对于单位四元数，它们的大小 − 1 < = a ⋅ b < = 1 -1 <= a \cdot b <= 1 − 1 < = a ⋅ b < = 1 q 1 ⋅ q 2 q1 \cdot q2 q 1 ⋅ q 2 q 1 q1 q 1 q 2 q2 q 2
# 和复数一样四元数，单位四元数 u = [ 0 , u ] u=[0,u] u = [ 0 , u ] e u θ e^{u\theta} e u θ q = [ c o s ( θ ) , s i n ( θ ) u ] q=[cos(\theta), sin(\theta)u] q = [ c o s ( θ ) , s i n ( θ ) u ] l o g ( q ) = l o g ( e u θ ) = [ 0 , u θ ] log(q)=log(e^{u\theta})=[0,u\theta] l o g ( q ) = l o g ( e u θ ) = [ 0 , u θ ] t t t
q t = ( e u θ ) t = e u ( t θ ) = [ c o s ( t θ ) , s i n ( t θ ) u ] q^t=(e^{u\theta})^t=e^{u(t\theta)}=[cos(t\theta), sin(t\theta)u] q t = ( e u θ ) t = e u ( t θ ) = [ c o s ( t θ ) , s i n ( t θ ) u ] 对于标量 a a a a 0 = 1 ; a 1 = a a^0=1; a^1=a a 0 = 1 ; a 1 = a t t t a t a^t a t a a a q t q^t q t [ 1 , 0 ] [1,\bold{0}] [ 1 , 0 ] q q q t t t q 8 q^8 q 8 q q q q 8 q^8 q 8 q t q^t q t
function interpolate ( [ w , x , y , z ] , t ) { if ( Math . abs ( w ) < 0.999 ) { const alpha = Math . acos ( w ) const newAlpha = alpha * t w = Math . cos ( newAlpha ) const mult = Math . sin ( newAlpha ) / Math . sin ( alpha ) x *= mult y *= mult z *= mult } return [ w , x , y , z ] } Copy # Lerp 表示线性插值它会沿着一条直线进行插值。
L e r p ( q 0 , q 1 , t ) = q 0 + t ( q 1 − q 0 ) = ( 1 − t ) q 0 + t q 1 Lerp(q_0, q_1, t)=q_0+t(q_1-q_0)=(1-t)q_0+tq_1 L e r p ( q 0 , q 1 , t ) = q 0 + t ( q 1 − q 0 ) = ( 1 − t ) q 0 + t q 1
可以看到它是沿着一条直线进行插值的，但是这样插值出来的四元数并不是单位四元数。
# 线性插值的结果不是单位四元数，但是我们可以在插值完毕后再将它变成单位四元数，这个被称为正规化线性插值 Nlerp。
class Quat { static lerp ( a , b , t , out = [ ] ) { const ax = a [ 0 ] , ay = a [ 1 ] , az = a [ 2 ] , aw = a [ 3 ] ; const bx = b [ 0 ] , by = b [ 1 ] , bz = b [ 2 ] , bw = b [ 3 ] ; const cosom = ax * bx + ay * by + az * bz + aw * bw ; const k0 = 1 - t ; const k1 = cosom >= 0 ? t : - t ;
out [ 0 ] = k0 * ax + k1 * bx out [ 1 ] = k0 * ay + k1 * by out [ 2 ] = k0 * az + k1 * bz out [ 3 ] = k0 * aw + k1 * bw
const s = 1 / Math . hypot ( out [ 0 ] , out [ 1 ] , out [ 2 ] , out [ 3 ] ) out [ 0 ] *= s out [ 1 ] *= s out [ 2 ] *= s out [ 3 ] *= s
return out } } Copy 上面代码可以对两个单位四元数进行线性插值，返回的也是单位四元数。上面代码中当两个四元数的点积小于 0 时会将其中一个四元数变负，原因请看下面的 Slerp。
Nlerp 插值仍然有一定问题，当需要插值的弧比较大时，角速度会有显著的变化。下图中线性插值是等长的，但是它们的弧长完全不相等。
# 四元数的一个优势就是球面线性插值 Slerp，它可以提供两个定向之间平滑的插值。它是对两个四元数的夹角进行插值，这样就可以不改变四元数的大小。
和线性插值一样我们首先需要求出两个四元数之间的夹角差，上面我们已经求过了，两个四元数的差如下。
Δ q = q 1 q 0 − 1 \Delta q=q_1q_0^{-1} Δ q = q 1 q 0 − 1 然后需要计算这个差的倍数，也就是旋转量。
然后通过四元数乘法组合来旋转 q 0 q_0 q 0
S l e r p ( q 0 , q 1 , t ) = ( q 1 q 0 − 1 ) t q 0 Slerp(q_0,q_1,t)=(q_1q_0^{-1})^tq_0 S l e r p ( q 0 , q 1 , t ) = ( q 1 q 0 − 1 ) t q 0 我们还可以想象两个二维矢量，将 v t v_t v t v 0 v_0 v 0 v 1 v_1 v 1
也就是存在 a a a b b b v t = a v 0 + b v 1 v_t=av_0+bv_1 v t = a v 0 + b v 1 v 0 v_0 v 0 v 1 v_1 v 1
通过上图我们可以看出来 s i n ( θ ) = s i n ( t θ ) b sin(\theta)=\frac{sin(t\theta)}{b} s i n ( θ ) = b s i n ( t θ ) b = s i n ( t θ ) s i n ( θ ) b=\frac{sin(t\theta)}{sin(\theta)} b = s i n ( θ ) s i n ( t θ ) a = s i n ( ( 1 − t ) θ ) s i n ( θ ) a=\frac{sin((1-t)\theta)}{sin(\theta)} a = s i n ( θ ) s i n ( ( 1 − t ) θ )
v t = a v 0 + b v 1 = s i n ( ( 1 − t ) θ ) s i n ( θ ) v 0 + s i n ( t θ ) s i n ( θ ) v 1 v_t=av_0+bv_1=\frac{sin((1-t)\theta)}{sin(\theta)}v_0+\frac{sin(t\theta)}{sin(\theta)}v_1 v t = a v 0 + b v 1 = s i n ( θ ) s i n ( ( 1 − t ) θ ) v 0 + s i n ( θ ) s i n ( t θ ) v 1 同样的思路我们扩展到四元数空间，可以重新得到 Slerp 函数。
S l e r p ( q 0 , q 1 , t ) = s i n ( ( 1 − t ) θ ) s i n ( θ ) q 0 + s i n ( t θ ) s i n ( θ ) q 1 Slerp(q_0,q_1,t)=\frac{sin((1-t)\theta)}{sin(\theta)}q_0+\frac{sin(t\theta)}{sin(\theta)}q_1 S l e r p ( q 0 , q 1 , t ) = s i n ( θ ) s i n ( ( 1 − t ) θ ) q 0 + s i n ( θ ) s i n ( t θ ) q 1 其中 θ \theta θ a r c c o s ( q 0 ⋅ q 1 ) arccos(q_0 \cdot q_1) a r c c o s ( q 0 ⋅ q 1 )
之前说过 q q q − q -q − q q 0 q_0 q 0 q 1 q_1 q 1 q 0 ⋅ q 1 q_0 \cdot q_1 q 0 ⋅ q 1 q 0 q_0 q 0 q 1 q_1 q 1
class Quat { static slerp ( a , b , t , out = [ ] ) { const ax = a [ 0 ] , ay = a [ 1 ] , az = a [ 2 ] , aw = a [ 3 ] ; let bx = b [ 0 ] , by = b [ 1 ] , bz = b [ 2 ] , bw = b [ 3 ] ; let theta , cosom , absCosom , sinom , k0 , k1 ;
cosom = ax * bx + ay * by + az * bz + aw * bw ; absCosom = Math . abs ( cosom ) if ( 1 - absCosom > 0.000001 ) { theta = Math . acos ( absCosom ) sinom = 1 / Math . sin ( theta ) k0 = Math . sin ( ( 1 - t ) * theta ) * sinom k1 = Math . sin ( t * theta ) * sinom } else { k0 = 1 - t k1 = t }
k1 = ( cosom >= 0 ) ? k1 : - k1 ; out [ 0 ] = k0 * ax + k1 * bx out [ 1 ] = k0 * ay + k1 * by out [ 2 ] = k0 * az + k1 * bz out [ 3 ] = k0 * aw + k1 * bw
return out ; } } Copy 我们将式子变成了上面函数，需要注意的是如果两个四元数夹角非常接近我们使用线性插值，不然 Slerp 会除 0。
# 除普通的四元数外，几何代数中还衍生出来了一个二重四元数（Dual Quaternion）。它不仅能够表示 3D 旋转，还能够表示 3D 空间中的任何的刚体运动（Rigid Motion），即旋转、平移、反射和均匀缩放。和普通四元数一样，二重四元数同样可以表示为 q = a + b i + c j + d k q=a+bi+cj+dk q = a + b i + c j + d k a , b , c , d a, b, c, d a , b , c , d 二元数 。
# 上面我们已经求出四元数 3D 旋转矩阵形式。
[ 1 − 2 y 2 − 2 z 2 2 x y − 2 w z 2 w y + 2 x z 2 x y + 2 w z 1 − 2 x 2 − 2 z 2 2 y z − 2 w x 2 x z − 2 w y 2 w x + 2 y z 1 − 2 x 2 − 2 y 2 ] \begin{bmatrix} 1-2y^2-2z^2 & 2xy-2wz & 2wy+2xz \\ 2xy+2wz & 1-2x^2-2z^2 & 2yz-2wx \\ 2xz-2wy & 2wx+2yz & 1-2x^2-2y^2 \end{bmatrix} ⎣ ⎢ ⎡ 1 − 2 y 2 − 2 z 2 2 x y + 2 w z 2 x z − 2 w y 2 x y − 2 w z 1 − 2 x 2 − 2 z 2 2 w x + 2 y z 2 w y + 2 x z 2 y z − 2 w x 1 − 2 x 2 − 2 y 2 ⎦ ⎥ ⎤ 所以直接将这个矩阵变成代码就行了。
class Mat4 { static fromQuat ( [ x , y , z , w ] , out = [ ] ) { const x2 = x + x , y2 = y + y , z2 = z + z ; const xx = x * x2 , yy = y * y2 , zz = z * z2 ; const yx = y * x2 , zx = z * x2 , zy = z * y2 ; const wx = w * x2 , wy = w * y2 , wz = w * z2 ;
out [ 0 ] = 1 - yy - zz out [ 1 ] = yx + wz out [ 2 ] = zx - wy out [ 3 ] = 0 out [ 4 ] = yx - wz out [ 5 ] = 1 - xx - zz out [ 6 ] = zy + wx out [ 7 ] = 0 out [ 8 ] = zx + wy out [ 9 ] = zy - wx out [ 10 ] = 1 - xx - yy out [ 11 ] = 0 out [ 12 ] = 0 out [ 13 ] = 0 out [ 14 ] = 0 out [ 15 ] = 1 return out } } Copy # 矩阵转四元数同样还是利用上方求出的四元数矩阵，有多种方式通过矩阵的项求解四元数的分量，我们来看一种常用的方式，因为是单位四元数，所以 x 2 + y 2 + z 2 + w 2 = 1 x^2+y^2+z^2+w^2=1 x 2 + y 2 + z 2 + w 2 = 1
m 11 + m 22 + m 33 = 3 − 4 x 2 − 4 y 2 − 4 z 2 = 4 w 2 − 1 m_{11}+m_{22}+m_{33}=3-4x^2-4y^2-4z^2=4w^2-1 m 1 1 + m 2 2 + m 3 3 = 3 − 4 x 2 − 4 y 2 − 4 z 2 = 4 w 2 − 1 通过将矩阵对角线的 3 项相加可以得到的 w w w w w w
m 32 − m 23 = 2 w x + 2 y z − 2 y z + 2 w x = 4 w x m 13 − m 31 = 2 w y + 2 x z − 2 x z + 2 w y = 4 w y m 21 − m 12 = 2 x y + 2 w z − 2 x y + 2 w z = 4 w z \begin{aligned} m_{32}-m_{23} &=2wx+2yz-2yz+2wx=4wx \\ m_{13}-m_{31} &=2wy+2xz-2xz+2wy=4wy \\ m_{21}-m_{12} &=2xy+2wz-2xy+2wz=4wz \end{aligned} m 3 2 − m 2 3 m 1 3 − m 3 1 m 2 1 − m 1 2 = 2 w x + 2 y z − 2 y z + 2 w x = 4 w x = 2 w y + 2 x z − 2 x z + 2 w y = 4 w y = 2 x y + 2 w z − 2 x y + 2 w z = 4 w z 我们可以发现将上面式子的结果除去 4 w 4w 4 w w w w 0 0 0 w w w
m 11 − m 22 − m 33 = 4 x 2 − 1 − m 11 + m 22 − m 33 = 4 y 2 − 1 − m 11 − m 22 + m 33 = 4 z 2 − 1 \begin{aligned} m_{11}-m_{22}-m_{33} &=4x^2-1 \\ -m_{11}+m_{22}-m_{33} &=4y^2-1 \\ -m_{11}-m_{22}+m_{33} &=4z^2-1 \end{aligned} m 1 1 − m 2 2 − m 3 3 − m 1 1 + m 2 2 − m 3 3 − m 1 1 − m 2 2 + m 3 3 = 4 x 2 − 1 = 4 y 2 − 1 = 4 z 2 − 1 为了获取稳定的四元数，我们应该选择 x , y , z , w x,y,z,w x , y , z , w
class Quat { static fromMat3 ( m , out = [ ] ) { const fTrace = m [ 0 ] + m [ 4 ] + m [ 8 ] ; let fRoot ;
if ( fTrace > 0 ) { fRoot = Math . sqrt ( fTrace + 1.0 ) ; out [ 3 ] = 0.5 * fRoot ; fRoot = 0.5 / fRoot ; out [ 0 ] = ( m [ 5 ] - m [ 7 ] ) * fRoot out [ 1 ] = ( m [ 6 ] - m [ 2 ] ) * fRoot out [ 2 ] = ( m [ 1 ] - m [ 3 ] ) * fRoot } else { let i = 0 ; if ( m [ 4 ] > m [ 0 ] ) i = 1 ; if ( m [ 8 ] > m [ i * 3 + i ] ) i = 2 ; let j = ( i + 1 ) % 3 ; let k = ( i + 2 ) % 3 ;
fRoot = Math . sqrt ( m [ i * 3 + i ] - m [ j * 3 + j ] - m [ k * 3 + k ] + 1.0 ) ; out [ i ] = 0.5 * fRoot ; fRoot = 0.5 / fRoot ; out [ 3 ] = ( m [ j * 3 + k ] - m [ k * 3 + j ] ) * fRoot ; out [ j ] = ( m [ j * 3 + i ] + m [ i * 3 + j ] ) * fRoot ; out [ k ] = ( m [ k * 3 + i ] + m [ i * 3 + k ] ) * fRoot ; }
return out } } Copy 上面代码将 3x3 的旋转矩阵转换成四元数。
# 将欧拉角转换成四元数和欧拉角转换成矩阵类似，我们找到绕 3 个基本轴的旋转四元数，然后将它们相乘就行了。四元数可以写成这种形式 [ c o s ( θ 2 ) ( s i n ( θ 2 ) x , s i n ( θ 2 ) y , s i n ( θ 2 ) z ) ] [cos(\frac{\theta}{2}) \quad (sin(\frac{\theta}{2})x, sin(\frac{\theta}{2})y, sin(\frac{\theta}{2})z)] [ c o s ( 2 θ ) ( s i n ( 2 θ ) x , s i n ( 2 θ ) y , s i n ( 2 θ ) z ) ]
X = [ c o s ( x / 2 ) , ( s i n ( x / 2 ) , 0 , 0 ) ] Y = [ c o s ( y / 2 ) , ( 0 , s i n ( y / 2 ) , 0 ) ] Z = [ c o s ( z / 2 ) , ( 0 , 0 , s i n ( z / 2 ) ) ] \begin{aligned} X&=[cos(x/2),(sin(x/2),0,0)] \\ Y&=[cos(y/2),(0,sin(y/2),0)] \\ Z&=[cos(z/2),(0,0,sin(z/2))] \end{aligned} X Y Z = [ c o s ( x / 2 ) , ( s i n ( x / 2 ) , 0 , 0 ) ] = [ c o s ( y / 2 ) , ( 0 , s i n ( y / 2 ) , 0 ) ] = [ c o s ( z / 2 ) , ( 0 , 0 , s i n ( z / 2 ) ) ] 然后将这三个四元数相乘就可以将欧拉角转换成四元数了，不过不同顺序的欧拉角四元数相乘的顺序也不同。
class Quat { static fromEuler ( x , y , z , out = [ ] , order = 'yxz' ) { x *= 0.5 ; z *= 0.5 ; y *= 0.5 ; const sx = Math . sin ( x ) , cx = Math . cos ( x ) ; const sy = Math . sin ( y ) , cy = Math . cos ( y ) ; const sz = Math . sin ( z ) , cz = Math . cos ( z ) ;
switch ( order ) { case "xyz" : out [ 0 ] = sx * cy * cz + cx * sy * sz ; out [ 1 ] = cx * sy * cz - sx * cy * sz ; out [ 2 ] = cx * cy * sz + sx * sy * cz ; out [ 3 ] = cx * cy * cz - sx * sy * sz ; break ;
case "xzy" : out [ 0 ] = sx * cy * cz - cx * sy * sz ; out [ 1 ] = cx * sy * cz - sx * cy * sz ; out [ 2 ] = cx * cy * sz + sx * sy * cz ; out [ 3 ] = cx * cy * cz + sx * sy * sz ; break ;
case "yxz" : out [ 0 ] = sx * cy * cz + cx * sy * sz ; out [ 1 ] = cx * sy * cz - sx * cy * sz ; out [ 2 ] = cx * cy * sz - sx * sy * cz ; out [ 3 ] = cx * cy * cz + sx * sy * sz ; break ;
case "yzx" : out [ 0 ] = sx * cy * cz + cx * sy * sz ; out [ 1 ] = cx * sy * cz + sx * cy * sz ; out [ 2 ] = cx * cy * sz - sx * sy * cz ; out [ 3 ] = cx * cy * cz - sx * sy * sz ; break ;
case "zxy" : out [ 0 ] = sx * cy * cz - cx * sy * sz ; out [ 1 ] = cx * sy * cz + sx * cy * sz ; out [ 2 ] = cx * cy * sz + sx * sy * cz ; out [ 3 ] = cx * cy * cz - sx * sy * sz ; break ;
case "zyx" : out [ 0 ] = sx * cy * cz - cx * sy * sz ; out [ 1 ] = cx * sy * cz + sx * cy * sz ; out [ 2 ] = cx * cy * sz - sx * sy * cz ; out [ 3 ] = cx * cy * cz + sx * sy * sz ; break ;
default : throw new Error ( 'Unknown angle order ' + order ) ; }
return out ; } } Copy # 我们可以结合上篇文章 求出旋转矩阵和这篇文章上面求出的四元数旋转矩阵来讲四元数转换成欧拉角。
[ c o s ( H ) c o s ( B ) + s i n ( H ) s i n ( P ) s i n ( B ) s i n ( H ) s i n ( P ) c o s ( B ) − s i n ( B ) c o s ( H ) s i n ( H ) c o s ( P ) s i n ( B ) c o s ( P ) c o s ( P ) c o s ( B ) − s i n ( P ) c o s ( H ) s i n ( P ) s i n ( B ) − s i n ( H ) c o s ( B ) s i n ( B ) s i n ( H ) + c o s ( H ) s i n ( P ) c o s ( B ) c o s ( H ) c o s ( P ) ] \begin{bmatrix} cos(H)cos(B)+sin(H)sin(P)sin(B) & sin(H)sin(P)cos(B)-sin(B)cos(H) & sin(H)cos(P) \\ sin(B)cos(P) & cos(P)cos(B) & -sin(P) \\ cos(H)sin(P)sin(B)-sin(H)cos(B) & sin(B)sin(H)+cos(H)sin(P)cos(B) & cos(H)cos(P) \end{bmatrix} ⎣ ⎢ ⎡ c o s ( H ) c o s ( B ) + s i n ( H ) s i n ( P ) s i n ( B ) s i n ( B ) c o s ( P ) c o s ( H ) s i n ( P ) s i n ( B ) − s i n ( H ) c o s ( B ) s i n ( H ) s i n ( P ) c o s ( B ) − s i n ( B ) c o s ( H ) c o s ( P ) c o s ( B ) s i n ( B ) s i n ( H ) + c o s ( H ) s i n ( P ) c o s ( B ) s i n ( H ) c o s ( P ) − s i n ( P ) c o s ( H ) c o s ( P ) ⎦ ⎥ ⎤ [ 1 − 2 y 2 − 2 z 2 2 x y − 2 w z 2 w y + 2 x z 2 x y + 2 w z 1 − 2 x 2 − 2 z 2 2 y z − 2 w x 2 x z − 2 w y 2 w x + 2 y z 1 − 2 x 2 − 2 y 2 ] \begin{bmatrix} 1-2y^2-2z^2 & 2xy-2wz & 2wy+2xz \\ 2xy+2wz & 1-2x^2-2z^2 & 2yz-2wx \\ 2xz-2wy & 2wx+2yz & 1-2x^2-2y^2 \end{bmatrix} ⎣ ⎢ ⎡ 1 − 2 y 2 − 2 z 2 2 x y + 2 w z 2 x z − 2 w y 2 x y − 2 w z 1 − 2 x 2 − 2 z 2 2 w x + 2 y z 2 w y + 2 x z 2 y z − 2 w x 1 − 2 x 2 − 2 y 2 ⎦ ⎥ ⎤ 可以发现 − s i n ( P ) = 2 y z − 2 w x -sin(P) = 2yz-2wx − s i n ( P ) = 2 y z − 2 w x p = a r c s i n ( − 2 ( y z − w x ) ) p=arcsin(-2(yz-wx)) p = a r c s i n ( − 2 ( y z − w x ) )
class Quat { static toEuler ( [ x , y , z , w ] , out = [ ] ) { const sp = - 2 * ( y * z - w * x )
if ( Math . abs ( sp ) < 0.99999 ) { out [ 0 ] = Math . asin ( sp ) out [ 1 ] = Math . atan2 ( x * z + w * y , 0.5 - x * x - y * y ) out [ 2 ] = Math . atan2 ( x * y + w * z , 0.5 - x * x - z * z ) } else { out [ 0 ] = Math . PI * ( sp < 0 ? - 0.5 : 0.5 ) out [ 1 ] = Math . atan2 ( w * y - x * z , 0.5 - y * y - z * z ) out [ 2 ] = 0 }
return out } } Copy 
