四元数与旋转

本文章的很多图被掘金的logo覆盖了，大家直接看这篇原文吧： www.yxyy.name/blog/md.htm…

前言

源码

github.com/buglas/thre…

学习目标

• 掌握四元数运算规则
• 掌握绕轴旋转的方法
• 掌握四元数和矩阵相互转换的方法

知识点

• 四元数
• 矩阵
• 向量

1-四元数基础

四元数(quaternion)就是由四个元构成的数字，这个数字是复数，其写法为a+bi+cj+dk。

四元数在计算机图形学、计算机视觉、磁共振成像、晶体纹理分析等领域都有广泛的应用。

因为四元数是对复数的一种扩展，所以它有着一套不同于常规数字的运算规则，这种规则一度被人称为邪恶的黑魔法，它神奇、抽象、复杂，却又能高效的解决一些实际问题，比如三维旋转。

四元数在用于三维旋转的时候，比矩阵的计算速度更快，相比于欧拉角而言，还可以避免万向锁。

1-1-四元数的定义

四元数的定义如下：

• a, b, c, d 都是实数

• a 叫做标量(scalar part) 或实部(real part)

• b, c, d 叫做矢量(vector part) 或虚部(imaginary part)

• 1,i, j, k 是四元数的基向量(basis vector)，其乘法表如下：

  

四元数通常用粗体的H表示。所以大家看见某个变量属于H时，要知道它是属于四元数的意思：

如果四元数中，只有一个标量a，i,j,k都是0，则此四元数叫做标量四元数(scalar quaternion)。

如果四元数中，标量a 是0，i,j,k中至少有一个非零，则此四元数叫做矢量四元数(vector quaternion)，或纯四元数。

四元数也可以分成scalar part 和vector part 两部分，其写法也可以用这两部分表示：

• r 是scalar part，对应实数a
• v 是vector part，对应bi+cj+dk，可以理解为三维向量

在下面的课程里，大家需要对这两种四元数的书写方式自由切换。

比如当我说四元数(r,v)=1 的时候，大家不要觉得四元数可以等于一个实数，它等于的是一个scalar quaternion，即(r,v)=1+0i+0j+0k

一般而言，(r,v)的书写形式，我用得会多一些，因为它书写简洁、一目了然。

1-2-四元数的加法

四元数的加法规则如下：

1-3-四元数的数乘

四元数的数乘规则如下：

λ 是实数，它与虚数相乘的时候，符合结合律，即λa=aλ。

1-4-四元数的乘法

四元数的乘法符合结合律，但不符合分配率，其规则如下：

根据四元数的定义，简化上式：

四元数的乘法也可以用scalar part 和vector part 表示：

我在上面用不同的色块表示了两种四元数乘法的书写方式的对应关系。

通过上式可知，四元数的叉乘是不符合的交换律的，这是由上式的v1×v2 导致的。

只有两个四元数的vector part共线时，这2个四元数的乘法才满足交换律。因为两个共线的向量的叉乘结果是零向量。

纯四元数的乘积为：

四元数的平方为：

1-5-四元数的共轭

四元数的共轭(conjugate)就是将四元数的vector part 取反的四元数。

如四元数q：

其共轭四元数是：

四元数的共轭通常用以下*号表示。

2个四元数相乘后满足以下规则：

2个四元数相乘后的共轭满足以下规则：

四元数的共轭可以用四元数相乘和相加来表示：

共轭可用于提取四元数的标量部分(scalar part)和矢量部分( vector part).

scalar part：

vector part：

1-6-四元数的范数

四元数与其共轭的乘积的平方根称为其范数(norm)：

通过其定义可知，四元数q乘以其norm，等于其norm乘以四元数，其结果为四元数norm的平方：

实数与四元数数乘时，其norm为：

两个四元数乘积的norm等于这两个四元数norm的乘积：

1-7-单位四元数

单位四元数(Unit quaternion)是norm为1的非零四元数。

以四元数q为例：

非零四元数除以其norm 便可得到一个单位四元数Uq：

单位四元数的norm 为1：

纯单位四元数的平方为-1：

这里的-1 并不只是实数，它首先是scalar quaternion，然后才是一个实数。

我们可以将单位四元数延伸到对于三角函数的思考。

单位向量可以写做以下形式：

• cosθ 对应四元数q里的scalar part，即a
• sinθ*u 对应四元数q里的vector part，即bi+cj+dk
• sinθ 是vector part 的长度
• u 是单位向量

根据三角函数的定义可知，这个四元数的norm 就是1。

以此类推，单位四元数的平方也可以用三角函数表示：

根据三角函数的倍角公式，上式还可以再做进一步简化：

这就是四元数可以应用于三维旋转的一个切入点，我们在后面会详解。

1-8-四元数的倒数

四元数的倒数(reciprocal)也可以理解为四元数逆，其计算方式如下：

可以简化上式：

用scalar part 和vector part 表示为：

四元数的倒数类似于矩阵里的逆矩阵。

由上面的定义可知，单位四元数的倒数等于其conjugate。

根据四元数的乘法定理可知，四元数与其倒数的乘积为1。

1-9-四元数的点积

四元数的叉乘跟vector part 无关，其叉乘结果是scalar quaternion。

以四元数p,q为例：

p和q的点积是：

由此可见，四元数的点积是符合交换律的。

点积也可以结合四元数的Conjugation来计算：

1-10-四元数的叉乘

四元数的叉乘结果跟scalar part 无关，四元数的叉乘结果是vector quaternion。

以四元数p,q为例：

p叉乘q是：

由此可见，四元数的叉乘是不符合交换律的。

叉乘也可以结合四元数来计算：

1-11-用矩阵表示四元数

以四元数q为例：

其对应的矩阵为：

我会在后面详解这个矩阵的推导过程。

2-四元数与旋转

一般当我们说起四元数的时候，首先就会想到绕某一轴旋转一定的度数，这是四元数在图形学中用得比较多的地方。

四元数定义了一种运算规则，这种规则为三维旋转提供了新的计算方式。这种计算方式快捷、高效，但是也很抽象。

所有在说四元数旋转之前，我们先说一些比较具象的绕轴旋转的解法，让大家对此有一个具象的认知。

2-1-几何解

已知：

• 单位向量 u(x,y,z)
• 点A绕u轴旋转θ°，落在点B
• 点O为坐标系原点

求：点B

解：

过点A，点B做它们在单位向量u所在直线上的共同的垂足P。

因为点O为坐标系原点，所以向量OB即点B；向量OA同理。

根据向量加法的三角形定理可知：

通过向量OA和单位向量u 的点积可以算出向量OP：

接下来再算出向量PB，便可通过它与向量OP的加法得向量OB。

做向量OA和单位向量u的垂线PC：

因为u是单位向量，所以PA和PC的长度相等。

将点B投影到向量PA、PC上，得点D和点E。

根据向量加法的三角形定理可知：

根据三角函数可知：

所以：

先在点B就已经算出来了，我们可以整体的回顾一下其计算流程。

2-2-矩阵解

在上面的公式里，我们知道了让任意点绕u轴旋转θ°的方法，所以我们也就可以知道让整个坐标系绕u轴旋转的矩阵。

我们可以分别让坐标系里的x,y,z 基向量绕u轴旋转，从而得到一个旋转矩阵。

让OA等于x轴基向量：

将向量u 分解为(x,y,z)，计算OB：

同理，当OA等于y轴基向量时：

OB如下：

当OA等于z轴基向量时：

OB如下：

最后，将三个基向量旋转后的结果合在一起，便是绕u轴旋转θ°的三维旋转矩阵：

2-3-基变换解

基变换的矩阵的结构就是矩阵M的两边各有一个互为逆矩阵的矩阵，如下所示：

基变换可以修改变换的基，比如把u 轴作为基向量进行变换。

接下来咱们具体看一下如何用基变换矩阵实现绕轴旋转。

已知：

• 单位向量 u(x,y,z)
• 点A绕u轴旋转θ°，落在点B
• 点O为坐标系原点

求：点B

思路：

让点A绕向量u旋转很难一步搞定。

但是如果我让点A绕当前坐标系里的基向量旋转，比如绕z轴，那就直接可以用三角函数实现。

所以接下来我们就把u轴，连带着点A，转到一个基向量上去。

让点A绕这个基向量转完之后，再把u轴，连带着点A，转回去。

解：

1.让向量u绕x轴逆时针旋转，转到Ozx平面中，设其旋转角度为α。

设相应的旋转矩阵为Mx

设向量u绕x轴旋转到Ozx平面上的向量为u'

连带着点A一起旋转，可得A'

2.让向量u' 绕y轴顺时针旋转，转到z轴上，设其旋转角度为β。

上面的旋转角度取负值是由顺时针旋转导致的。

设相应的旋转矩阵为My

设向量u' 绕y 轴旋转到z 轴上的向量为u''。

连带着点A'一起旋转，可得A''

3.此时，我们就可以让A'' 绕z 轴旋转θ°了，因为此时的u'' 与z轴共线。

设旋转矩阵为Mu''

设A'' 旋转后的位置为B''

4.将B'' 按照之前的My 和Mx 矩阵逆向旋转回去。

设My 和Mx 的逆矩阵分别是Myi 和Mxi。

因为正交旋转矩阵的逆矩阵就是其转置矩阵，所以：

通过Myi 和Mxi 将B'' 逆向旋转回去，即可得到OB

5.系统总结归纳一下这个过程。

通过上式可知，MxiMyi 和MyMx 是互逆的关系。

简化一下上式，设：

由此可得到一种经典的基变换模式：

解释一下其中的概念。

当前存在了两种坐标系的，我们可以暂且用世界坐标系和本地坐标系区分其关系。

详细解释一下其中的概念：

• 单位向量u：世界坐标系里的单位向量。当被变换到本地坐标系里时，与本地坐标系里的一条基向量共线，在我们当前的示例里，它与z 轴共线，当然它也可以与其它的基向量共线。
• 点A：世界坐标系里的点位。
• 矩阵N：将世界坐标系转本地坐标系的矩阵。
• 逆矩阵N：将本地坐标系转世界坐标系的矩阵。
• 矩阵M ：要对本地坐标系里的点A施加变换的矩阵。

由上可知，我们要将矩阵M作用于世界坐标系里的点A，首先要通过NA 方法将点A转入本地坐标系，然后再MNA。

MN*A 得到的是本地坐标系里的点，即点B''，而我们要求的是世界坐标系里的点B，所以需要再用矩阵N的逆矩阵将本地坐标系里的点B''转入世界坐标系。

2-4-四元数解

四元数旋转公式：

• L(p')：四元数的左乘函数。
• p：要旋转的点。
• p' ：要旋转的点的四元数形态，p'=(0,P)。
• q：四元数，其中包含了旋转轴和旋转度数相关的数据。
• r：点p绕u 轴旋转θ°后的向量。

2-4-1-四元数旋转公式的推导

在几何解里，我们说过求OB的方式：

上式可以转化为四元数的形式：

因为：

所以：

对(0,PB)做分解，可得：

我可以通过逆推的方法验证一下这个分解过程：

因为单位向量u 和PA垂直，所以：

所以：

接下来我们继续推导四元数的旋转公式。

上式中的(cosθ,sinθ*u) 是一个四元数，为了方便运算，我们可以设其为w：

因此，(0,OB)可以简化为：

观察四元数w，我们可以发现它是一个单位四元数：

由此我们可以根据单位四元数的平方公式得出：

设：

根据单位四元数的平方公式可知：

将其代入(0,OB)：

在(0,OP)前面再写一个1

因为四元数与其倒数的乘积是1，所以：

因为四元数q 的倒数的vector part 的方向是u，向量OP 的方向也是u，所以两者共线，所以两者的相乘符合交换律，即：

根据四元数norm 的乘法规则可知：

因为四元数q 为单位四元数，所以：

根据上面的等式整理OB：

向量(0,OP) 和(0,OA) 可以合在一起：

现在，我们就看到了最纯粹的四元数变换公式：

• (0,OA) 就是我们要旋转的点的四元数形态。
• 四元数q及其倒数实现四元数旋转的方法。

设(0,OA)为P'，OA为点P，套上我们之前说过的四元数旋转公式，求出r值来，就是点A绕u轴旋转θ°后的点位。

2-4-2-四元数旋转与基变换的对比

对比四元数旋转与基变换矩阵的公式，你会发现它们长得都跟肉夹馍似的，都是中间夹一层肉肉(P',M)，两边是互为倒数的面饼。

它们的这种相似并不是偶然，因为其中都蕴含着一种反应其本质的基变换。

四元数的基变换是四维的，基变换矩阵是三维的。

与此同时，对比公式的长度，你还会发现，基变换矩阵的公式比四元数旋转公式多了一个A，这是四维和三维的差异。

四元数可以直接将变换作用于化身为四元数的三维向量。

而矩阵需要先算出三维变换矩阵，再将其应用于三维向量。

由此可见，矩阵旋转会造成运算的浪费。而且在3乘3的矩阵里有9个数字，它们的相乘，运算量会更大，远远不如四元数里4个数字的运算快捷。

3-四元数对象的封装

我们可以结合three.js 里相关方法来学习四元数的封装。

3-1-Quaternion 对象

在three.js 的Quaternion 对象是norm 为1 的单位四元数对象，它主要用于三维旋转。

构造函数

constructor( x = 0, y = 0, z = 0, w = 1 ) {
    this.isQuaternion = true;
    this._x = x;
    this._y = y;
    this._z = z;
    this._w = w;
}

x,y,z,w 对应下面的四元数：

setFromAxisAngle( axis, angle )

绕某一轴旋转特定弧度，返回四元数。

• axis 围绕旋转的轴，单位向量
• angle 弧度

代码如下：

setFromAxisAngle( axis, angle ) {
    // http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToQuaternion/index.htm
    // assumes axis is normalized

    const halfAngle = angle / 2, s = Math.sin( halfAngle );
    this._x = axis.x * s;
    this._y = axis.y * s;
    this._z = axis.z * s;
    this._w = Math.cos( halfAngle );
    this._onChangeCallback();
    return this;
}

通过代码，可以看出，它返回的就是我们之前算出的四元数q：

multiplyQuaternions( a, b )

四元数乘法，a乘以b。

multiplyQuaternions( a, b ) {

    // from http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/code/index.htm

    const qax = a._x, qay = a._y, qaz = a._z, qaw = a._w;
    const qbx = b._x, qby = b._y, qbz = b._z, qbw = b._w;

    this._x = qax * qbw + qaw * qbx + qay * qbz - qaz * qby;
    this._y = qay * qbw + qaw * qby + qaz * qbx - qax * qbz;
    this._z = qaz * qbw + qaw * qbz + qax * qby - qay * qbx;
    this._w = qaw * qbw - qax * qbx - qay * qby - qaz * qbz;

    this._onChangeCallback();

    return this;

}

公式如下：

conjugate()

计算四元数的conjugate。

conjugate() {
    this._x *= - 1;
    this._y *= - 1;
    this._z *= - 1;
    this._onChangeCallback();
    return this;
}

公式如下：

invert()

计算四元数的倒数，因为当前四元数默认就是单位四元数，所以其倒数就是其conjugate。

invert() {
    // quaternion is assumed to have unit length
    return this.conjugate();
}

公式如下：

dot( v )

四元数的点积。

dot( v ) {
    return this._x * v._x + this._y * v._y + this._z * v._z + this._w * v._w;
}

公式如下：

length()

四元数的norm。

length() {
    return Math.sqrt( this._x * this._x + this._y * this._y + this._z * this._z + this._w * this._w );
}

公式如下：

normalize()

四元数归一化，即单位四元数。

normalize() {
    let l = this.length();
    if ( l === 0 ) {
        this._x = 0;
        this._y = 0;
        this._z = 0;
        this._w = 1;
    } else {
        l = 1 / l;
        this._x = this._x * l;
        this._y = this._y * l;
        this._z = this._z * l;
        this._w = this._w * l;
    }
    this._onChangeCallback();
    return this;
}

公式如下：

setFromRotationMatrix( m )

从旋转矩阵中获取四元数，我们稍后会细说这个方法。

其它的方法我就不说了，一般大家只要具备了四元数基础，就能看得懂。

其中还有跟欧拉相互转换的方法，我暂且不说，等以后说欧拉的时候一起说。

3-2-Vector3 对象里的applyQuaternion(q) 方法

Vector3 对象的applyQuaternion(q) 方法就是用四元数旋转三维点的。

applyQuaternion( q ) {
    const x = this.x, y = this.y, z = this.z;
    const qx = q.x, qy = q.y, qz = q.z, qw = q.w;
  
    // calculate quat * vector
    const ix = qw * x + qy * z - qz * y;
    const iy = qw * y + qz * x - qx * z;
    const iz = qw * z + qx * y - qy * x;
    const iw = - qx * x - qy * y - qz * z;

    // calculate result * inverse quat
    this.x = ix * qw + iw * - qx + iy * - qz - iz * - qy;
    this.y = iy * qw + iw * - qy + iz * - qx - ix * - qz;
    this.z = iz * qw + iw * - qz + ix * - qy - iy * - qx;

    return this;
}

这个方法看着挺复杂的，但对于理解了四元数旋转公式的同学而言，便可以一眼看穿其本质。

我结合四元数旋转公式来解释上面的方法。

1.计算qp'

const ix = qw * x + qy * z - qz * y;
const iy = qw * y + qz * x - qx * z;
const iz = qw * z + qx * y - qy * x;
const iw = - qx * x - qy * y - qz * z;

其算法就是根据四元数的乘法公式写的：

2.用qp'的结果乘以q的倒数。

this.x = ix * qw + iw * - qx + iy * - qz - iz * - qy;
this.y = iy * qw + iw * - qy + iz * - qx - ix * - qz;
this.z = iz * qw + iw * - qz + ix * - qy - iy * - qx;

整体逻辑就这么简单。

当然，你也可以用四元数旋转公式里计算r 的公式一气呵成。

这个公式的代码就留给大家自己写了，算是给大家练练手，有啥问题再随时跟我说(wx:1051904257)。

4-四元数转矩阵

通常，为了快捷高效的变换物体，物体的位移、旋转和缩放数据会一起存储的到一个模型矩阵里，如three.js里的Object3D对象的matrix 属性。

为了方便操作，我们通常还会把物体的位移、旋转和缩放再进行单独封装，这样就可以更简单的控制物体的某一种变换状态。

物体的位移、和缩放数据会通过向量存储，而旋转数据会通过四元数存储。

因此，当我们把四元数所代表的旋转数据和其它的变换数据一起存储到模型矩阵中时，就需要将其转为矩阵。

接下来我们就说一下将四元数转矩阵的方法。

之前我们说了绕u轴旋转θ°的矩阵：

其中的x,y,z 是单位向量u的x,y,z，为了区分后面要说的四元数，我将x,y,z 改成ux,uy,uz

假设代表了旋转的四元数如下所示：

那我接下来想办法把绕u轴旋转θ°的矩阵里的cosθ、sinθ、ux、uy、uz 替换成四元数里的w、x、y、z 即可。

首先，我们要知道cosθ、sinθ、ux、uy、uz 和w、x、y、z 的关系：

接下来以矩阵里的第一行第一列为例说一下替换方法：

以此类推，其它的矩阵因子也做类似的转换，你就会得到这样的一个由四元数里的4个数字构成的矩阵：

在three.js的Matrix4对象的makeRotationFromQuaternion(q) 方法里，就是按照这个原理将四元数转矩阵的。

大家可以先默认position为(0,0,0)，scale为(1,1,1)，只关注quaternion 部分。

compose( position, quaternion, scale ) {
    const te = this.elements;
    const x = quaternion._x, y = quaternion._y, z = quaternion._z, w = quaternion._w;
    const x2 = x + x,    y2 = y + y, z2 = z + z;
    const xx = x * x2, xy = x * y2, xz = x * z2;
    const yy = y * y2, yz = y * z2, zz = z * z2;
    const wx = w * x2, wy = w * y2, wz = w * z2;
    const sx = scale.x, sy = scale.y, sz = scale.z;

    // 1-2y²+2z²
    te[ 0 ] = ( 1 - ( yy + zz ) ) * sx;
    // 2xy+2wz
    te[ 1 ] = ( xy + wz ) * sx;
  // 2xz-2wy
    te[ 2 ] = ( xz - wy ) * sx;
    te[ 3 ] = 0;

    te[ 4 ] = ( xy - wz ) * sy;
    te[ 5 ] = ( 1 - ( xx + zz ) ) * sy;
    te[ 6 ] = ( yz + wx ) * sy;
    te[ 7 ] = 0;

    te[ 8 ] = ( xz + wy ) * sz;
    te[ 9 ] = ( yz - wx ) * sz;
    te[ 10 ] = ( 1 - ( xx + yy ) ) * sz;
    te[ 11 ] = 0;

    te[ 12 ] = position.x;
    te[ 13 ] = position.y;
    te[ 14 ] = position.z;
    te[ 15 ] = 1;
    return this;
}

5-矩阵转四元数

我们可以根据之前所说的四元数旋转矩阵推导出其中的四个元数。

将上面的矩阵简写为以下结构：

观察四元数旋转矩阵，我们可以发现一些计算四元数的规律。

举一个计算四元数中的w的例子。

观察m00+m11+m22的和：

你会发现w和m00+m11+m22的规律：

这样，我们通过四元数旋转矩阵计算w的方法就有了。

同理，我们可以再找出求x,y,z的规律：

然而，这个思路是对的，但这里面还埋了个坑。

根号下的1+m00+m11+m22不能为负，x,y,z中的除数4w 也不能等于0的，所以1+m00+m11+m22 要大于0。

在有的旋转矩阵里，1+m00+m11+m22 是不满足上述条件的，如下面的旋转矩阵：

这个问题的解决方法是：

在w,x,y,z 中寻找满足计算条件的数，在其中肯定会有满足计算条件的数。

在找到此数后，就像上面的w一样，以此数的4倍为除数，算出其它的数。

c++ 代码如下：

float tr = m00 + m11 + m22

if (tr > 0) { 
  float S = sqrt(tr+1.0) * 2; // S=4*qw 
  qw = 0.25 * S;
  qx = (m21 - m12) / S;
  qy = (m02 - m20) / S; 
  qz = (m10 - m01) / S; 
} else if ((m00 > m11)&(m00 > m22)) { 
  float S = sqrt(1.0 + m00 - m11 - m22) * 2; // S=4*qx 
  qw = (m21 - m12) / S;
  qx = 0.25 * S;
  qy = (m01 + m10) / S; 
  qz = (m02 + m20) / S; 
} else if (m11 > m22) { 
  float S = sqrt(1.0 + m11 - m00 - m22) * 2; // S=4*qy
  qw = (m02 - m20) / S;
  qx = (m01 + m10) / S; 
  qy = 0.25 * S;
  qz = (m12 + m21) / S; 
} else { 
  float S = sqrt(1.0 + m22 - m00 - m11) * 2; // S=4*qz
  qw = (m10 - m01) / S;
  qx = (m02 + m20) / S;
  qy = (m12 + m21) / S;
  qz = 0.25 * S;
}

第一个if 条件之所以是tr > 0，而不是tr+1 > 0，这是为了避免计算机浮点误差，其余的if 条件也是同理。

three.js中Quaternion 对象的setFromRotationMatrix(m)方法也是按照这个逻辑走的：

setFromRotationMatrix( m ) {
    // http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
    // assumes the upper 3x3 of m is a pure rotation matrix (i.e, unscaled)
    
  const te = m.elements,
        m11 = te[ 0 ], m12 = te[ 4 ], m13 = te[ 8 ],
        m21 = te[ 1 ], m22 = te[ 5 ], m23 = te[ 9 ],
        m31 = te[ 2 ], m32 = te[ 6 ], m33 = te[ 10 ],
        trace = m11 + m22 + m33;

    if ( trace > 0 ) {
        const s = 0.5 / Math.sqrt( trace + 1.0 );
        this._w = 0.25 / s;
        this._x = ( m32 - m23 ) * s;
        this._y = ( m13 - m31 ) * s;
        this._z = ( m21 - m12 ) * s;
    } else if ( m11 > m22 && m11 > m33 ) {
        const s = 2.0 * Math.sqrt( 1.0 + m11 - m22 - m33 );
        this._w = ( m32 - m23 ) / s;
        this._x = 0.25 * s;
        this._y = ( m12 + m21 ) / s;
        this._z = ( m13 + m31 ) / s;
    } else if ( m22 > m33 ) {
        const s = 2.0 * Math.sqrt( 1.0 + m22 - m11 - m33 );
        this._w = ( m13 - m31 ) / s;
        this._x = ( m12 + m21 ) / s;
        this._y = 0.25 * s;
        this._z = ( m23 + m32 ) / s;
    } else {
        const s = 2.0 * Math.sqrt( 1.0 + m33 - m11 - m22 );
        this._w = ( m21 - m12 ) / s;
        this._x = ( m13 + m31 ) / s;
        this._y = ( m23 + m32 ) / s;
        this._z = 0.25 * s;
    }
    this._onChangeCallback();
    return this;
}

6-代码测试

这里的代码测试，主要是避免有的同学学了之后心里没底，不确定我这是不是瞎扯，所以咱们画出来看看。

当然，若你将前面的内容都理解透彻了，这里可以略过。

当前测试使用的是three.js。

6-1-几何解测试

根据几何解的公式写代码：

整体代码如下：

import React, { useRef, useEffect } from 'react'
import {
  AxesHelper,
    BufferGeometry,
    Color,
    Float32BufferAttribute,
    Line,
    LineBasicMaterial,
    Vector3,
} from 'three'
import Stage from '../component/Stage'
import './fullScreen.css'

/* 快速初始化项目 */
const stage = new Stage(1, 2, 4)
const { scene, renderer } = stage
renderer.shadowMap.enabled = true

/* 已知条件 */
const u=new Vector3(1,1,1).normalize()
const A=new Vector3(0.75,0.4,0.4)
const OP=u.clone().multiplyScalar(A.dot(u))
const PA=A.clone().sub(OP)
const uxA=u.clone().cross(A)

/* u轴 */
const lineGeo = new BufferGeometry()
lineGeo.setAttribute('position', new Float32BufferAttribute([0,0,0,u.x,u.y,u.z], 3))
lineGeo.setAttribute('color', new Float32BufferAttribute([0,1,0,1,1,0], 3))
const lineMat = new LineBasicMaterial({ vertexColors:true})
const line = new Line(lineGeo, lineMat)
scene.add(line)

/* 坐标系 */
const axes=new AxesHelper(1)
const axesColor= new Color(0x666666)
axes.setColors(axesColor,axesColor,axesColor)
scene.add(axes)

/* 绕u旋转的轨迹 */
const points=[A.x,A.y,A.z]
const trackPointNum=20
const trackGeo = new BufferGeometry()
const trackColors=[]
for(let i=0;i<trackPointNum;i++){
  const n=i/trackPointNum
  trackColors.push(1,n,0)
}
trackGeo.setAttribute('color', new Float32BufferAttribute(trackColors, 3))
const trackMat = new LineBasicMaterial({ vertexColors:true})
const track = new Line(trackGeo, trackMat)
scene.add(track)

/* 时间 */
const time=new Date().getTime()
/* 连续旋转 */
stage.beforeRender=function(){
  // 旋转弧度
  const theta=(new Date().getTime()-time)*0.006
  const c=Math.cos(theta)
  const s=Math.sin(theta)
  // 旋转后的位置
  const B=OP.clone()
  .add(PA.clone().multiplyScalar(c))
  .add(uxA.clone().multiplyScalar(s))
  // 过滤极小距离点
  if(new Vector3(points[0],points[1],points[2]).sub(B).lengthSq()<0.001){
    return
  }
  // 在数组前面添加点
  points.unshift(B.x,B.y,B.z)
  // 删除最后一个点
  const pointsLen=points.length
  if(pointsLen/3>trackPointNum){
    points.splice(pointsLen-4,pointsLen-1)
  }
  // 更新position
  trackGeo.setAttribute('position', new Float32BufferAttribute(points, 3))
}

const Rotate01: React.FC = (): JSX.Element => {
    const divRef = useRef<HTMLDivElement>(null)
    useEffect(() => {
        const { current } = divRef
        if (!current) {
            return
        }
        current.append(renderer.domElement)
        stage.animate()
        return () => {
            renderer.domElement.remove()
        }
    }, [])
    return <div ref={divRef} className="canvasWrapper"></div>
}

export default Rotate01

效果如下：

在上面的代码里，我通过连续绕u轴旋转点A，得到了一条旋转路径。

在这条路径里，当顶点超过一定数量时，我会将路径最后面的顶点删除，以确保它不会无限长。

其旋转的核心代码如下：

  // 旋转弧度
  const theta=(new Date().getTime()-time)*0.006
  const c=Math.cos(theta)
  const s=Math.sin(theta)
  // 旋转后的位置
  const B=OP.clone()
  .add(PA.clone().multiplyScalar(c))
  .add(uxA.clone().multiplyScalar(s))

6-2-矩阵解测试

根据矩阵解的矩阵写代码：

在上一个示例的基础上写代码：

/* 连续旋转 */
stage.beforeRender=function(){
  // 旋转弧度
  const theta=(new Date().getTime()-time)*0.006
  const c=Math.cos(theta)
  const s=Math.sin(theta)
  
  // 旋转后的位置
  // const B=OP.clone()
  // .add(PA.clone().multiplyScalar(c))
  // .add(uxA.clone().multiplyScalar(s))

  // 旋转矩阵
  const {x,y,z}=u
  const n=1-c
  const rotateMatrix=new Matrix3().set(
    x*x*n+c,x*y*n-z*s,x*z*n+y*s,
    x*y*n+z*s,y*y*n+c,y*z*n-x*s,
    x*z*n-y*s,y*z*n+x*s,z*z*n+c,
  )
  // 旋转后的位置
  const B=A.clone().applyMatrix3(rotateMatrix)
  ……
}

效果和之前都是一样的。

6-3-基变换解测试

根据矩阵解的矩阵写代码：

代码如下：

/* 已知条件 */
const u=new Vector3(1,1,1).normalize()
const A=new Vector3(0.75,0.4,0.4)

const alpha=Math.atan2(u.z,u.y)
const ca=Math.cos(alpha)
const sa=Math.sin(alpha)
const Mx=new Matrix3().set(
  1,0,0,
  0,ca,-sa,
  0,sa,ca
)
const u1=u.clone().applyMatrix3(Mx)
const beta=-Math.atan2(u1.x,u1.z)
const cb=Math.cos(beta)
const sb=Math.sin(beta)
const My=new Matrix3().set(
  cb,0,sb,
  0,1,0,
  -sb,0,cb
)
const N=new Matrix3().multiplyMatrices(My,Mx)
const Ni=N.clone().transpose()

//……

/* 连续旋转 */
stage.beforeRender=function(){
  // 旋转弧度
  const theta=(new Date().getTime()-time)*0.006
  const c=Math.cos(theta)
  const s=Math.sin(theta)
  
  // 旋转矩阵
  const M=new Matrix3().set(
    c,-s,0,
    s,c,0,
    0,0,1
  )
  const rotateMatrix=Ni.clone().multiply(M).multiply(N)

  // 旋转后的位置
  const B=A.clone().applyMatrix3(rotateMatrix)

  // ……
}

在计算N的逆矩阵的时候，我对它进行的是转置操作，这是因为正交矩阵的逆矩阵就是其转置矩阵，转置矩阵的运算效率要比逆矩阵高。

6-4-四元数解测试

四元数旋转公式是这样的：

我们可以通过计算r 值得到旋转后的点位，但我就不算了。

接下来我直接通过three.js 里的Quaternion 对象转一下看看。

/* 连续旋转 */
stage.beforeRender=function(){
  // 旋转弧度
  const theta=(new Date().getTime()-time)*0.006
  // 四元数
  const quaternion=new Quaternion().setFromAxisAngle(u,theta)
  // 旋转后的位置
  const B=A.clone().applyQuaternion(quaternion)
  // ……
}

其效果和之前都是一样的。

总结

这一篇我们主要说了四元数的定义、四元数的基本运算规则、多种绕轴旋转的实现方法、四元数与矩阵的相互转化。

四元数是用于存储三维物体旋转量的首选，因为其性能优于矩阵和欧拉。

参考链接：

en.wikipedia.org/wiki/Quater…

www.euclideanspace.com/maths/geome…

zhuanlan.zhihu.com/p/56587491

blog.csdn.net/loongkingwh…