四维计算机图形学

四维计算机图形学：旋转篇

这次来看看如何实现最基础的坐标变换——旋转，它是渲染建模动画等一切技术的基石。这篇文章相当于一篇关于旋转的算法综述，虽然不需要读者熟悉计算机图形学，但会对数学水平要求稍微高一些，本文主要是罗列陈述算法，因此内容略显枯燥，若读者仅想学习使用这些工具构建四维交互场景，请关注后续的Tesserxel4D场景开发教程。

一般使用旋量表示N维空间中的旋转。特别地，3维与4维的旋量代数可以用四元数表示

特色内容

  • 三维四元数旋转相关算法
  • “几何代数”与“四元数”版四维旋量算法对比
  • 各种“lookAt”旋转生成算法
  • 四维相机的常见交互控制方式

前置知识

我本来打算先讲二维旋转、三维旋转的表示方法再讲到四维，逐渐引入旋转平面、2-向量、几何代数、旋量的，但由于低维旋转方面网上资料很多，且前面我也已经写过关于旋转平面、2-向量、几何代数、旋量的内容了，因此就不想啰嗦了（其实是懒~~），于是我整理了一下需要读者先了解的前置知识，大致知道这些内容后就可再继续阅读：
  1. 《四维空间（七）：N维的向量》 中开头至奇异2-向量小节。
  2. 《四维空间（十一）：几何代数、四元数与空间旋转》 中开头至四元数与四维空间旋转小节。

三维旋转的处理

首先空间中的点与向量都以三个浮点数来储存定义为三维向量类（class Vec3），我采用四元数表示三维旋转，因此需要定义个四元数类（class Quaternion），最后，矩阵（class Mat3、class Mat4）作为最经典的图形学工具也是必不可少的。

建议在计算机上采用四元数储存旋转，首先欧拉角有万向锁问题，其次矩阵有16个元素，远多于三维旋转的3个自由度，一是浪费内存，二是不便于纠正浮点数误差累积。若想利用显卡对矩阵乘法的并行优化，可在发送给显卡时再转换为矩阵。

基础数学类的运算

我们会定义上面类的成员变量和一些常用方法，下面以三维向量Vec3类为例，它的伪代码大概是这个样子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
class Vec3{

    // 成员变量：

    x: number; y: number; z: number;

    // 构造函数：

    constructor(x: number, y: number, z: number){
        this.x = x; this.y = y; this.z = z;
    }

    // 常用方法：

    // 反向
    neg():Vec3{
        return new Vec3(-this.x, -this.y, -this.z);
    }
    // 加法
    add(v3:Vec3):Vec3{
        return new Vec3(this.x + v3.x, this.y + v3.y, this.z + v3.z);
    }
    ...
    // 内积
    dot(v3:Vec3):number{
        return this.x * v3.x + this.y * v3.y + this.z * v3.z;
    }
    ...
}

如果你无法理解上述代码，提前学习typescript/javascript、java、C#或C++中某一个语言中关于类(class)的相关概念即可，它们并不复杂。本文采用类似typescript的伪代码来描述算法，即函数参数名冒号后面代表参数类型，函数返回值类型则在参数括号后以冒号标注。

四元数类的定义为：

1
2
3
4
5
class Quaternion{
    // r代表实数分量，i/j/k分别代表三个虚数前的系数
    r: number; i:number; j: number; k: number;
    ...
}

旋转坐标与四元数转矩阵

给定一个旋转$R$和一个向量$p$，如何计算旋转后的新向量$R(p)$是最基本的需求，我们将这个函数命名为apply。如果旋转用的是矩阵表示，则使用矩阵乘法完成计算：（矩阵与向量的乘法公式是明确的，本文假设读者有能力写出写出该函数的实现细节，故省略之，下同）

1
2
3
function apply(R:Mat3, p:Vec3):Vec3{
    return R.mul(p); // mul函数是矩阵与向量的乘法。
}

如果旋转用的是四元数表示，则我们采用这里的旋转坐标算法计算：

1
2
3
4
5
function apply(R:Quaternion, p:Vec3):Vec3{
    let p0 = new Quaternion(0, p.x, p.y, p.z); // 向量转四元数，注意实部为0。
    let p1 = R.mul(p0).mul(R.conj()); // mul、conj函数分别是四元数乘法与共轭。
    return new Vec3(p1.i, pi.j, p1.k); // 四元数转回向量。
}

但一般显卡中处理旋转还是用矩阵比较多，如何将四元数表示的旋转转为矩阵表示呢？下面提供两种方法：

  1. 硬推公式法：我们将上面算法中的四元数乘法按坐标展开，整理出线性变换的公式，将它写成矩阵即可，建议使用Mathematica等工具推导公式，推导结果可见我的tesserxel库中Quaternion类的toMat3函数。
  2. 变换基向量法：如果嫌公式推导麻烦，可直接将矩阵按列分块，会发现它是由所有坐标基向量旋转后得到的新向量拼成的。因此我们无需理会四元数如何旋转空间中的点，就能构造出矩阵（注意这些都是是列向量）：

$$M=\begin{pmatrix}R(e_x)&R(e_y)&R(e_z)\end{pmatrix}$$

翻译成伪代码则是：

1
2
3
4
5
6
7
8
9
10
function quaternion2mat3(R:Quaternion):Mat3{
    let ex = apply(R, new Vec3(1, 0, 0));
    let ey = apply(R, new Vec3(0, 1, 0));
    let ez = apply(R, new Vec3(0, 0, 1));
    return new Mat3(
        ex.x, ey.x, ez.x, 
        ex.y, ey.y, ez.y, 
        ex.z, ey.z, ez.z, 
    );
}

注意虽然第二种方法看似写出的代码优雅，但效率不如第一种快：基向量中要么是1要么是0，很多乘法与加法都是不必要的，将式子化简后效率更高，但其实就变成第一种方法了。

相反，把矩阵转换为四元数不太容易，我们将在后面解决这个问题。现在有了初步处理旋转的能力，下面先来看怎样生成旋转。

旋转的复合与逆

采用欧拉角或其它多步骤复合表示旋转：则使用该算法生成多次旋转的四元数，注意四元数$p$、$q$先后作用到点$x$上为$q(pxp^*)q^*=(qp)x(qp)^*$，相当于四元数$qp$单次作用，因此四元数跟矩阵一样表示复合旋转是直接按从右到左的顺序乘起来。求旋转的逆就更简单了，单位四元数乘以自己的共轭就是1，因此四元数共轭就是逆。

轴与角度生成旋转

  1. 直接给出旋转的轴与角度：其算法已在刚才的链接中与旋转坐标算法一起给出了，这里放一下伪代码：
    1
    2
    3
    4
    5
    6
    7
    // 给定单位向量axis表示的旋转轴，与旋转角度angle，生成四元数表示的旋转
    function axisAngle2quaternion(axis: Vec3, angle: number): Quaternion{
        angle *= 0.5; // 将角度减半，并计算其正余弦值s与c
        let s = Math.sin(angle);
        let c = Math.cos(angle);
        return new Quaternion(c, s * axis.x, s * axis.y, s * axis.z);
    }
  2. 直接给出旋转生成元2-向量：由于三维的2-向量可以霍奇对偶到普通的1-向量，因此其实这相当于给出了平行于旋转轴的向量，且向量的长度为旋转角度，相当于把旋转的轴与角度两个参数合二为一了，类似于角速度矢量这种东西，这种把大小方向合二为一的表示法又是比第一种更方便。它的算法其实是一样的，只是增加一些步骤：首先计算向量的长度的一半得到半角$\theta = ||l||/2$，然后将向量单位化后就跟第一种方法一样了。我们可进行化简：由于单位化即除以角度$2\theta$，最后四元数的虚部是单位向量乘以半角的正弦值$\sin(\theta)$，把这两步合并即向量乘以$\sin(\theta)/(2\theta)$。显然这个系数在角度接近零时分子分母都接近0，但极限为1/2，这并不是个数值稳定的方法，可以在角度小于$\epsilon$时将表达式用前两项泰勒展开$(1/2)(1-\theta^2/3!)$来避免。这相当于生成元到旋转的指数映射，因此该函数命名为exp。下面给出伪代码：
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    function exp(v: Vec3): Quaternion {
        let theta = v.length() * 0.5; // length函数求向量模长，得到半角
        let epsilon = 0.005;
        // 若角度的绝对值大于epsilon，则直接计算系数s，否则用泰勒展开
        let s:number;
        if (Math.abs(theta) > epsilon) {
            s = Math.sin(theta) / theta * 0.5;
        } else {
            s = 0.5 - theta * theta / 12;
        } 
        return new Quaternion(Math.cos(theta), s * v.x, s * v.y, s * v.z);
    }

lookAt生成旋转

仅给出物体的初始朝向与最终朝向，找到符合要求的旋转：比如在一个射击游戏里，大炮需要始终自动对焦到移动的敌机，设大炮的局部坐标系中炮管的方向为$x$轴，则要求一个旋转变换，它能把$x$轴旋转到从炮到敌机连线的方向$v$。几乎所有游戏引擎中都有现成的函数供程序员使用，一般叫“lookAt”函数，下面我们来看如何实现它。
给定两个单位向量$u$与$v$，我们希望构造旋转$R$，使得$R(u)=v$。由于两个向量确定一个平面，最简便的方式就是在$u$-$v$平面中进行旋转，或绕轴$u\times v$旋转，旋转的角度则可以通过计算向量的夹角得到，此时已知转轴和夹角，使用前面的算法就能生成旋转了：

1
2
3
4
5
6
function lookAt(u:Vec3, v:Vec3):Quaternion{
    let axis = u.cross(v); // 通过向量叉乘函数cross计算旋转轴
    let e_axis = axis.normalize(); // normalize函数将转轴缩放为单位向量
    let angle = Math.acos(u.dot(v)); // 内积dot函数求夹角余弦，反余弦函数Math.acos求角度
    return axisAngle2quaternion(e_axis, angle); // 已知转轴和夹角，使用前面的算法生成旋转
}

其实上面的代码漏掉了一个特殊情况：若$u$与$v$方向相同或相反，叉乘得到的轴为$0$，这些特殊情况的处理留作习题。
注意，满足$R(u)=v$的旋转$R$是不唯一的：$R$上再任意复合一个以$v$为轴的旋转也符合要求。上面算法求得的旋转在某种意义下是“距离”最近的旋转，但有时候会像下图那样出现不好的“倾斜”感，比如摄像机运镜时可能就希望画面不要左右倾斜，我们马上来解决这个问题。将蓝色方向对齐的“最短”旋转路径会让绿色水平轴不再水平

给出物体的初始朝向与最终朝向，与物体的上方，找到保持物体不侧倾的旋转。方法有两种，一是将旋转拆成水平与垂直的两步可以避免物体侧倾，二是可以先用lookAt函数旋转到位后，再叠加一个以目标为轴的旋转纠正倾斜。这两种方法其实都相当于欧拉角，只是旋转的次序存在差异，具体代码实现留作习题。避免绿轴“侧倾”的旋转生成方法之一

旋转矩阵转四元数

前面说过旋转矩阵的自由度超过了旋转的自由度，数值误差导致不一定是精确的旋转矩阵，因此很可能问题没有解，但有了“lookAt”函数，我们就可以通过类似“Schmidt正交化”的方式来逐步逼近得到四元数：旋转矩阵按列分块后能够告诉我们原来的各坐标轴旋转到的新位置，因此可先使用lookAt函数生成对齐x轴前后位置的旋转r1，此时y轴与z轴还没对齐，然后再来一次lookAt对齐y轴，剩下的z轴就自动对齐了，将两次lookAt得到的四元数乘起来就得到了最终的转换结果。

1
2
3
4
5
6
7
8
function mat32quaternion(m:mat3): Quaternion{
    // 第一次旋转：对齐x轴
    let r1 = lookAt(new Vec3(1, 0, 0), m.x_()); // x_函数返回矩阵m的第一个列向量
    // 第二次旋转：对齐y轴
    let newY = apply(r1, new Vec3(0, 1, 0)); // 注意y轴已经被r1旋转到了newY
    let r2 = lookAt(newY, m.y_()); // y_函数返回矩阵m的第二个列向量
    return r2.mul(r1); // 将两步旋转复合
}

已知旋转找转轴与角度

有些时候会遇到旋转的生成的逆问题：给定旋转找出它的生成元，由于生成元生成旋转是指数映射exp，因此这个函数叫它对数映射log。求逆映射并不难，很容易通过四元数实部值求出半角的余弦，下面是伪代码：

1
2
3
4
5
log(R: quaternion): Vec3 {
    let theta = Math.acos(R.r); // 通过四元数实部求半角
    let s = 2 * theta / Math.sin(theta); // 恢复旋转轴，再乘上角度值
    return new Vec3(R.i * s, R.j * s, R.k * s);
}

旋转的插值

在动画制作中常会遇到已知物体始末初始状态，需补出中间帧以实现动画效果。若已知$0$时刻的旋转对应四元数$R_a$，$1$时刻对应四元数$R_b$，求时刻$t$的旋转对应的四元数$R_t$。旋转是由单位四元数表示的，相当于在超球面$S^3$上插值，算法可参考这里，代码见这里的slerp函数。

数值稳定性

在对场景的计算中，随着三维物体的不断旋转，四元数一直在做连乘，这是数值不稳定的，因为只要有点误差，四元数的模长不为1，就会导致一段时间后要么衰减至零或发散至无穷。解决该问题的方法很简单，隔一段时间（如刷新每一帧时）强行把四元数除以它的长度单位化即可。

四维旋转的处理

计算机表示四维空间中的点、向量跟三维向量是大同小异的，比如有四维向量类（class Vec4），和矩阵（class Mat4），甚至你可以定义仿射矩阵（class Mat5）。四维旋转同样不建议用矩阵表示，它采用旋量（Rotor）表示旋转，且旋量还分四元数与几何代数两种版本。无论用哪个版本，都会用到表示平面的2-向量（class Bivec），我们先来看它的定义。

2-向量类

四维空间的2-向量则用六个浮点数来储存，定义为2-向量类（class Bivec4）。由于我不打算开发通用的高维图形库，因此不会出现Bivec5、Bivec6……又由于霍奇对偶，三维中的Bivec3其实就是Vec3、因此也就只有Bivec4才有用，于是我把Bivec4简写为Bivec不会有歧义。它的伪代码定义如下：

1
2
3
4
class Bivec{
    xy: number; xz:number; xw: number;
    yz: number; yw:number; zw: number;
}

其实它很像一个6维的向量，比如可以定义加法(add)、减法(sub)、取反(neg)、数乘(mul)、内积(dot)、单位化(normalize)等方法。同时，我们添加一些2-向量独有的方法，比如霍奇对偶(dual)、混合积（也叫交换积）(cross)，并给四维向量类添加楔积运算(wedge)得以生成2-向量等。

四维旋量

四维旋量分两种表示法：一是四元数版的旋量（class Rotorq），二是几何代数版的旋量（class Rotorg），其中几何代数版最原始且易懂，而四元数版处理旋转算法更多更方便。由于几乎没人会直接输入旋量的各分量来描述旋转，它们仅为程序内部的中间数据，所以搭建渲染引擎时一般仅选其一，比如Marc ten Bosch的四维引擎就仅支持几何代数版本的旋量，而我的tesserxel引擎就只支持四元数版本的旋量而不支持几何代数版的，下面我把这两种版本旋量各自的相关旋转算法分为两节进行罗列，请按需展开阅读：


+旋量类（几何代数版）



+旋量类（四元数版）


随机方向/旋转的生成

随机单位向量的生成

随机单位向量的生成问题其实等价于在球面上生成均匀分布的点。这个问题在数学上已经研究得十分透彻了：首先直接均匀生成一段区间上的随机数，再根据坐标映射的“伸缩”系数通过一个非线性函数映射到球坐标上的经纬度值。下面直接列出伪代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
function randVec3(){
    let a = Math.random() * 2.0 * Math.PI;
    let c = Math.random() * 2.0 - 1.0;
    let b = Math.sqrt(1.0 - c * c);
    return new Vec3(b * Math.cos(a), b * Math.sin(a), c);
}
function randVec4(){
    let a = Math.random() * 2.0 * Math.PI;
    let b = Math.random() * 2.0 * Math.PI;
    let c = Math.random();
    let sc = Math.sqrt(c);
    let cc = Math.sqrt(1.0 - c);
    return new Vec4(sc * Math.cos(a), sc * Math.sin(a), cc * Math.cos(b), cc * Math.sin(b));
}

顺便提一句，这两个函数在上一篇文章提到的路径追踪算法生成照片级渲染中大有用处。

随机旋转的生成

由于三维旋转可以由同构于超球面的单位四元数表示，因此可以使用前面生成四维随机单位向量的算法生成四元数。若采用四元数版本的旋量，则四维旋转由左、右两个四元数表示，随机生成两个超球上的两个点即可；若采用几何代数版本的旋量，我目前还不知道有很方便的生成算法。

1
2
3
4
5
6
function randQuaternion():Quaternion{
    return new Quaternion(randVec4());
}
function randRotorq():Rotorq{
    return new Rotorq(randQuaternion(), randQuaternion());
}

随机简单旋转的生成

上面的随机旋转算法生成的朝向是任意的，如果要求必须生成简单旋转则必须先生成随机朝向的简单单位2-向量，再随机生成一个角度来实现。最简单生成随机朝向的简单单位2-向量的做法为先生成两个随机向量，通过楔积张成2-向量将它单位化即可。其实还有一种更快捷的生成方法：简单2-向量的自对偶与反自对偶部分大小相等，因此它们的长度都是$\sqrt 2/2$。而每个部分都相当于一个二维球面，因此可以用过随机生成二维球面上的点来生成这些单位自对偶/反自对偶2-向量，将它们乘以$\sqrt 2/2$加在一起就得到了单位化的简单2-向量。

四维相机的常见控制方式

轨道球（TrackBall）控制模式

该模式类似一般3D软件中的视图旋转，对应到四维则类似Jenn3D的操作方式。该控制器存储一个旋转中心点和距离，再加一个相机朝向的旋量，该模式中相机始终面对中心点，设相机坐标系下前方为w轴，则可按以下方式更新它的位置：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class TrackBallController{
    // 旋转中心坐标
    center: Vec4;
    // 相机距中心的距离
    distance: number;
    // 相机朝向
    rotor: Rotorq;
    update(camera: Object4D){
        // 更新相机旋转
        camera.rotation = this.rotor;
        // 初始位置设为相机前方给定长度的向量，通过相机朝向求出相机相对中心的位置，最后叠加上中心
        camera.position = center.add(apply(this.rotor,new Vec4(0,0,0,distance)));
    }
}

用户使用键鼠交互更新距离与旋量，比如鼠标左键在屏幕上上下左右拖拽时在每帧的移动量被转换成一个在相机坐标系下仅有yw和xw分量的2-向量$B$，它生成的旋转$V’=\exp(B)$首先变换至世界坐标$V$，即$V=RV’R^{-1}$，然后再叠加到目前相机的朝向$R$上得到总旋转$VR = RV’R^{-1}R=RV’=R\exp(B)$，这可以总结为若后面叠加的旋转是在旋转物体的局部坐标系中表示的就用右乘，若都是同样的全局坐标系则使用通常的左乘。伪代码如下：

1
2
3
4
5
6
7
8
9
10
11
class TrackBallController{
    ...
    onLeftButtonDrag(dx:number, dy:number){
        // 根据鼠标移动距离生成2-向量 dx*exw + dy * eyw
        let localBivec = new Bivec(0, 0, dx, 0, dy, 0);
        // 生成旋转R
        let localRotation = exp(localBivec);
        // 将其转换至世界坐标系，叠加旋转后更新新的相机朝向，注意叠加的旋转在局部坐标，乘法顺序颠倒
        this.rotor = localRotation.mul(this.rotor);
    }
}

自由飞行（FreeFly）控制模式

该模式一般模拟无重力下的第一人称游走，它不用单独储存信息，直接读取相机的位置与朝向来根据用户的键鼠输入更新即可，需注意用户的旋转平移都应在相机的局部坐标系中定义再转换至世界坐标，这里不再给出伪代码。

保持竖直（KeepUp）模式

该模式一般模拟在重力环境下的第一人称游走，要求相机的侧前后方与左右方张成的平面始终保持在水平胞中不倾斜。有两种方法，第一种是使用类似欧拉角的技术，用两个旋量分别储存水平旋转与垂直旋转，第二种则是在每次旋转后纠偏。实际使用中我发现第二种方法稳定性差一些，因此tesserxel引擎中采用第一种方式。该旋转的自由度仅为4，它由水平胞内的三维自由旋转加一个垂直方向的俯仰角构成，因此我们用一个旋量与一个俯仰角来储存旋转，我采用的是鼠标拖拽控制水平胞内的旋转，使用滚轮控制俯仰角。最后每帧按下面方式更新相机朝向:

1
2
3
4
5
6
7
8
9
10
11
12
class KeepUpController{
    // 水平方向的朝向由一个旋量表示
    rotorH: Rotorq;
    // 垂直方向的俯仰角由一个标量数字表示
    rotorV: number;
    update(camera: Object4D){
        // 计算仅考虑水平胞旋转的相机局部坐标系中产生俯仰角的垂直旋转，它由yw平面上的2-向量生成
        let verticalRotation = planeAngle2rotorq(new Bivec(0,0,0,0,1,0), this.rotorV);
        // 本来该按水平旋转再垂直旋转的顺序合成最终旋转，但叠加的旋转在局部坐标，乘法顺序颠倒
        camera.rotation = rotorH.mul(this.verticalRotation);
    }
}

控制相机移动可以有两种模式，一是飞行模式，即在旋转时除了使用上述方法保持水平胞不倾斜外，其余的体验跟自由飞行模式相当，即若相机有俯仰角时，前进后退会让相机高度随之改变，二是地面漫游模式，这种模式下相机的前进后退不会考虑俯仰角，始终在地平面内运动（比如Minecraft在地面行走的玩家），并额外增加直接控制竖直升降相机的键盘交互。