1 前言

OpenGL渲染3D模型离不开空间几何的数学理论知识，而本篇文章的目的就是对空间几何进行简单的介绍，并对3D模型渲染需要使用到的几何知识，和其中各个坐标系，以及模型从3D空间投影到2D平面的过程进行详细说明。只有充分理解这些数学知识，我们才能熟练的渲染出如下有趣的场景。

下面的Demo会在后面的章节中提供完整的源码，在本篇文章中我们并不需要关注其具体实现。另外本文的知识不仅仅是OpenGL渲染3D图形所依赖的理论知识基础，同时也是其他图形渲染语言的基础知识。

2 向量和矩阵

向量(Vectors)和矩阵(Matrixs)是空间几何中的基本知识，也是3D图形渲染中最常用到的两个元素。

2.1 向量

顶点(Vertex)是OpenGL的主要输入对象，每个顶点可以看做是三维空间中的一个向量。OpenGL从顶点数组中建立三角形图元，通过在三维空间中连接各个顶点形成了一个有向矢量三角形。在着色器语言(GLSL)中二维、三维和四维向量分别被表示为vec2、vec3和vec4。对于向量OpenGL关心其方向(Direction)和长度(Magnitude)两个属性，另外在OpenGL中经常会用到标准化值以及标准化向量，标准化(Normalization)向量的定义是长度为1的向量，如(1, 0, 0)。

OpenGL用于表示三维空间中的方向向量通常只需要使用(x,y,z)三个坐标分量即可，如用于计算光照时使用的平面法向量。但是当处理三维物体变形时需要和4✖️4的矩阵做乘法运算，因此在做矩阵变化时必须使用顶点需要使用4维向量表示，这在本文稍后章节介绍齐次坐标系时会讲到。

2.1.1 基本运算

向量间的加减乘除要求进行运算的两个向量元素数量相等，计算时各个元素执行相应的运算即可。向量和标量的加减乘除基本元素直接将标量和每个向量元素执行相应运算即可。

2.1.2 点乘

向量的点乘(Dot Product)也称为向量的内积(Inner Product)，其运算结果等于各个元素相乘之和(v1 * v2 = v1.x * v2.x + v1.y * v2.y + v1.z * v2.z)，在几何意义上可以理解为两个向量的长度和其夹角的余弦值的积。在漫射光照(Diffuse Lighting)环境下，常需要计算平面的单位法向量和朝向光源单位向量点积，从而确定片段的颜色值。在着色器语言中，调用内建函数dot()计算两个向量的点乘。

2.1.3 叉乘

向量的叉乘(Cross Product)也成为矢量积，线性代数上运算时将两个向量和一个单位向量组合形成一个矩阵，借助部分行列式的运算规则得到新的向量。向量的叉乘也是一个向量。几何意义上表示为垂直于被计算两个向量组成平面的向量，另外对于向量v1和v2叉乘得到的向量v3，v3的长度等于v1的长度和v2的长度及它们之间的正弦乘积，即两个向量构成的平行四边形的面积。另外不同于点乘，向量的叉乘顺序会影响其最终的计算结果，即影响最终得到的向量方向。

向量a和向量b的叉乘可以通过下列行列式计算，通常I, j, k取单位向量(1,1,1)。a×b = ((ay×bz - az×by) × i, (az×bx - ax×bz) × j, (ax×by - ay × bx) × k)。在着色器语言中，调用函数cross()计算两个向量的叉乘。

向量v2×v1叉乘的几何意义可以表示为下图，其中v3垂直于平面v1-v2，其遵循右手螺旋法则。

2.1.4 长度

向量的长度(Length)也称向量的模，一个三维向量的长度计算公式如下。在着色器语言中，可以调用函数length()计算某个向量的长度。

2.1.5 折射和反射

在光照着色公式中常需要计算某个向量的折射(Refraction)和反射(Reflection)向量来确定片段最终的颜色，下图中Rin为入射标准向量，N为界面的标准法向量，Rreflect为反射标准向量。Rrefract为不同折射率的标准折射向量。

OpenGL计算反射向量时，入射向量和反射面的法向量都必须使用标准向量，这样，其计算结果只包含方向信息。反射向量的计算方法如下。在着色器语言中可以直接调用OpenGL函数reflect()计算反射向量。

OpenGL计算折射向量时需要指定折射率，其计算方法如下。在着色器语言中可以直接调用OpenGL函数refract()计算反射向量。

2.2 矩阵

在OpenGL中图形的变化几乎都是通过矩阵(Matrices)运算的方式实现。顶点v绕任意点以任意方向旋转得到的新顶点vn，vn的x坐标分量不仅与原顶点v的x坐标分量及旋转参数相关，还与v的y、z轴坐标分量相关。数学上的矩阵用行列式的形式表示，而OpenGL中矩阵以二维数组的形式表示。矩阵之间可以做加法和乘法运算，一个右向量表示的顶点乘以仿射矩阵后可以得到转换后新的点。

在4✖️4的矩阵中，即OpenGL中的数组中，前4个元素为矩阵第一列，紧接着4个元素为矩阵的第二列，以此类推。在将a点从当前坐标系转换到目标坐标系的场景中，4✖️4的矩阵还可以理解为是一个三维坐标系的3个轴单位向量及其原点，在目标坐标系中的坐标。例如下图的仿射矩阵，暂时不看最后一行，前三列分别代表了x，y，z轴的单位向量在目标坐标系中的坐标，最后1列代表当前坐标系原点在目标坐标系中的位置。左上角3*3矩阵可以实现物体的旋转以及缩放，最后一列可实现物体的平移。

在坐标系CoordinateSystemO，使用一个4✖️4矩阵Matrix表示一个它和坐标系CoordinateSystemN的关系，一个4维向量VecO表示CoordinateSystemO中的任意顶点，Matrix和VecO相乘得到新的4维向量VecN，则该向量为顶点在坐标系CoordinateSystemN中对应位置的顶点。

当一个顶点经过多次矩阵变化到新的坐标空间时，也可以先将每次进行坐标变换的矩阵相乘，直接得到原始坐标系到最终坐标系的仿射矩阵。即A*(B*v) = (A*B)*v，但是这里需注意乘法顺序，因为矩阵的乘法只服从乘法结合律，而不服从交换律。

如下面的上3图中，正方体先围绕自己的z轴逆时针旋转一定角度，再沿着其x轴平移一段距离。而在下3图中，应用相反的变换过程，可以明显看到我们得到了完全不同的图像。

3 坐标空间

OpenGL等图形渲染语言都需要将物体从三维空间中经过一系列仿射变换和投影，最终得到二维的可显示图形，而在理解OpenGL中这些工作是如何进行之前，我们需要先认识OpenGL中的各个坐标系。在OpenGL中会使用到基于欧式几何建立模型空间、世界空间、相机空间和标准设备空间，以及基于投影几何建立的投影空间。

我们可以想象一个场景在我们眼睛中呈现的图像取决于场景中各个模型的相对位置，以及我们观察的姿势。同样的在OpenGL中图形渲染的第一步就是确定这些模型的相对位置，OpenGL每个模型自己的坐标系(模型坐标系)以确定每个模型的形状，定义了一个世界坐标系确定模型的相对位置，还定义了一个相机坐标系(视点坐标系)确定了我们的观察位置，这样将每个顶点的坐标转换到相机坐标系就能够确定它在我们视野中的位置。然后在裁剪空间中通过下文将会讲到的投影变换裁剪掉不可见部分，并将场景渲染到屏幕上。

3.1 模型空间

模型空间(Model Space)用于描述物体空间信息和决定模型自己的形态，通常也称为物体坐标空间(Object Space)。通常选择物体的重心或者某一特定的点建立三维空间坐标系，并用一系列顶点来描述一个三维模型，坐标原点通常不会在模型外部。在一个场景中，每个模型都存在自己的模型空间，如下图。

3.2 世界空间

世界空间(World Space)以一个全局的坐标原点为中心建立，在OpenGL中全局坐标原点是固定的，前文的模型空间都是在世界空间中以一个新的原点建立的3D坐标系，所有的模型都将被放在世界空间中。一些光照以及物理计算执行时都以世界空间为参照，或者相机空间。世界空间示意图如上。

3.3 视图空间

视图空间(View Space)是和观察者相关的空间，有时也称为相机或者视点空间(Camera or Eye Space)。同模型空间类似，视图空间也是基于世界坐标系建立的。其x，y正方向方向表示相机的右侧和上方，其z轴负方向表示相机的镜头朝向。

3.4 剪切空间

剪切空间(Clip Space)为一个平截椎体空间，由一个负责投影的矩形近平面和一个负责控制的矩形远平面组成，位于两个平面之间的区域即为剪切空间。 每个位于剪切空间中的顶点都包含4个分量，称为4维齐次坐标，其第四个分量和顶点到摄像机的距离成正比。在执行裁剪的时候，顶点的4个分量都会除以第四个分量w，当w不为1时，其余三个分量都会被放大或者缩小，这种近大远小的效果表现为透视和投影效果。

这里需要注意的是和前面的模型空间、世界空间、相机空间都是基于笛卡尔右手坐标系建立不同。剪切空间不是基于欧式几何建立的坐标系，它是基于投影几何建立的坐标系，和欧式坐标系中平行直线不会相交不同，在投影几何中，平行直线的投影在无穷远处相交。

3.5 标准设备空间

为了适配不同的设备，OpenGL采用标准设备空间(Normalized Device Space)，x轴正方向水平向右，y轴正方向垂直向上，z轴正方向垂直向内，它们的取值范围都是[-1, 1]。只有在真正显示在屏幕上的时候才会计算出每个顶点的真实像素坐标。经过剪切步骤的顶点都会被转化到标准设备空间内。标准设备坐标空间示意图如上。

另外这里需要注意的是，和前面的模型空间、世界空间、相机空间都是基于笛卡尔右手坐标系建立不同。标准设备空间是基于笛卡尔左手直角坐标系建立的。

4 仿射变换

数学中不能使用绝对值来描述物体的位置，只能描述其相对于参考系的位置，同样模型在空间中的变换也是相对位置的改变。描述模型变换主要有矩阵形式、欧拉角和四元数三种描述方式。在OpenGL中最常用的是矩阵，另外还可以使用四元素，但是欧拉角的表示方法存在死锁的问题，因此很少使用。

4.1 矩阵形式

OpenGL中模型的平移、选择和缩放操作都可以通过一个仿射矩阵和该模型的所有顶点向量相乘获得。iOS和MacOS中系统库GLKit提供了一系列矩阵的初始化和计算接口，在其他平台上也有开源的C++数学库可以使用，但是深入理解其原理对学习3D图形变换很有帮助。

4.1.1 单位矩阵

单位矩阵(Identity Matrix)的行列数相等，除右对角线元素为1外，其余元素为0。单位矩阵和任意模型顶点相乘都不会改变其位置。在GLKit框架中可以通过常量GLKMatrix4Identity获得一个4✖️4的单位矩阵。单位矩阵和向量的乘法计算如下图。

4.1.2 平移矩阵

前文说过4✖️4仿射矩阵最后一列表示当前坐标系原点在目标坐标系中的坐标，OpenGL中平移一个模型可以看做是将模型自己的坐标系平移，即将模型从旧坐标系转换到新坐标系的过程。一个模型应用平移矩阵(Translation Matrices)的结果如下图。

OpenGL中通常使用齐次坐标w分量(第4个分量)为1的四维向量表示空间中的一个点，当其w分量不为1.0时，其平移效果会在转换矩阵的基础上产生缩放效果。另外使用三维向量或者w分量为0的四维向量表示方向向量。平移仿射矩阵处理顶点的计算公式如下。在iOS和MacOS中可以通过函数GLKMatrix4MakeTranslation(float tx, float ty, float tz)获得一个平移矩阵。

4.1.3 旋转矩阵

绕坐标轴旋转的旋转矩阵
对于旋转矩阵(Rotation Matrices)，首先考虑模型绕x轴简单旋转，旋转后模型的x轴坐标不会变化。前文已经说明，4阶仿射矩阵的前三列代表当前坐标系轴的三个基向量在目标坐标系中的表示，如下图虚线是当前坐标系，而实线是目标坐标系。将当前坐标系绕X轴正方向旋转90度，当前坐标系的基向量在目标坐标系中的坐标如下。

根据上图可以直接写出这个仿射变换的坐标系如下，在库GLKit中可以使用函数直接获取旋转矩阵GLKMatrix4MakeXRotation(float radians)。

同理可以推出绕y轴和绕z轴旋转的矩阵如下，GLKit中可以使用函数GLKMatrix4MakeYRotation(float radians)、GLKMatrix4MakeZRotation(float radians)直接获取这些旋转矩阵。

当一个模型以物体坐标系为参考，分别实现绕x，y，z轴旋转后，也就实现了绕经过原点的任意向量旋转，GLKit内部代码会根据在三个轴上的旋转角度构建出类似下图的仿射矩阵。在GLKit中可以通过函数GLKMatrix4MakeRotation(float radians, float x, float y, float z)获得对应矩阵。

注意以世界坐标系作为参照进行仿射变换时，最后一次做的变换需要放在整个矩阵乘法公式的最右侧，这样才能将模型从最初的空间多次转换到最终的空间中。另外如果以物体坐标系为参照，则这种仿射变换等同于逆序的参照世界坐标系的变换，即需要将最第一次的仿射变换矩阵放在等式右侧。

绕任意轴旋转的旋转矩阵
如下图在空间中将向量v绕单位向量n旋转-θ得到向量v‘。

首先得出如下公式。

根据前4个公式可以推到除如下式子。

可以计算出对于基向量p = [1, 0, 0]，其旋转后的向量如下

可以计算出对于基向量q = [0, 1, 0]，其旋转后的向量如下

可以计算出对于基向量r = [0, 0, 1]，其旋转后的向量如下

前文已经说过仿射矩阵的前三列为三个单位基向量的方向，因此可以写出绕任意轴旋转的仿射矩阵如下。在GLKit中使用函数GLKMatrix4MakeRotation(float radians, float x, float y, float z)获得相应矩阵，该函数内部会对x-y-z进行标准化处理。

4.1.4 缩放矩阵和镜像矩阵

沿坐标轴缩放矩阵和镜像矩阵
沿坐标轴缩放矩阵(Scale Matrices)比较简单，只需将缩放后新的基向量直接组合成为仿射矩阵即可，但需注意其x-y-z轴上不等比缩放时会使模型发生变形。可使用函数GLKMatrix4MakeScale( float sx, float sy, float sz)获得。当缩放因子为-1则得到镜像变化的仿射矩阵。一个缩放矩阵表示如下，需要注意这里的基向量在其自己的坐标系内仍是单位向量。

沿任意方向缩放矩阵
对于空间中的向量v沿着单位向量n的方向以k为缩放因子进行缩放，如下图所示。注意使用n的反向量得到的是相同结果。

首先我们可以推断出如下4个公式。

根据上述公式我们可以计算出。

则对于基向量p = [1, 0, 0]，其缩放后可以表示为

则对于基向量q = [0, 1, 0]，其缩放后可以表示为

则对于基向量r = [0, 0, 1]，其缩放后可以表示为

上一节已经讲过物体变换使用到的仿射矩阵可以由当前坐标系的基向量在目标坐标系的坐标构建，因此可以得出其仿射矩阵如下。另外此处令缩放因子k=-1，则得到绕任意轴取镜像的仿射矩阵。

4.1.5 连接仿射矩阵

通常一个模型到目标多坐标需要多次仿射变换，此时我们只需要将每次仿射变换的矩阵相乘即可，只需要记住以世界坐标系作为参照进行仿射变换时，最后一次做的变换需要放在整个矩阵乘法公式的最右侧，以物体坐标系为参照，需要将最第一次的仿射变换矩阵放在等式右侧。

4.1.6 矩阵蠕变

使用矩阵时需注意矩阵蠕变现象，即当输入错误的数据，或者浮点运算产生的累积误差使得矩阵不再正交。此时需要矩阵正交化来消除误差，重新得到一个最接近原矩阵的正交矩阵。

正交矩阵
矩阵M和它的转置矩阵MT的乘积等于单位矩阵时，矩阵M称为正交矩阵。即M的转置等于M的逆。

在空间几何中，假定M是一个3✖️3矩阵，使用r1，r2，r3分别表示矩阵的1至3列，其中r1，r2，r3分别代表空间中的3个基向量，只要满足每个向量都是单位向量，同时任意两个向量垂直，此时M就可以称为正交矩阵。对于OpenGL中所有的仿射矩阵的左上角3✖️3子矩阵都满足上述两个条件，因此都为正交矩阵。

矩阵正交化
斯密特正交化的思想是逐列处理矩阵，第一列不变换，后面的每列表示的向量都减去前面所有列表示的向量方向上的分量。对于3✖️3阶矩阵M，r1，r2，r3分别表示其1至3列，正交化计算过程如下。

此时得到的新的基向量已经相互垂直，但是其并不是单位向量，此时还需对它们分别进行标准化操作。但是施密特正交化是有偏差的，因为其r1总不能得到校准。一种改进的方法是旋转一个小的因子k，每次只减去投影的k倍，通过多次迭代最终得到一个正交矩阵，其计算过程如下。

4.1.7 矩阵的优点

• 可以立即进行向量旋转：矩阵形式和四元数可以直接对向量进行运算得到旋转后的向量坐标，欧拉角不能。
• 矩阵形式被图形API使用：程序中保存方位可以使用欧拉角和四元数，但是OpenGL内部执行旋转时必须为其提供矩阵形式。
• 多个角位移连接：当对个坐标系转换时可以先合成一个转换矩阵，最后再对坐标进行转换。
• 矩阵的逆：使用逆矩阵可以撤销原来的操作，由于旋转矩阵是正交的，因此其逆矩阵等于其转置矩阵。

4.1.8 矩阵的缺点

• 矩阵占用了更多的内存：保存一个方位，矩阵需要使用9个数，欧拉角只是用3个数，四元数只使用4个数。
• 难于理解：矩阵不能直观表示物体的方位。
• 矩阵可能是病态的：矩阵会发生矩阵蠕变现象。
• 矩阵无法插值：给定两个矩阵，无法在其中进行曲面线性插值。

4.2 欧拉角形式

欧拉角的基本思想是3D物体的旋转分解为绕3个相互垂直轴的三个旋转组成的序列，这组描述三个角度及旋转信息的数称为欧拉角。

欧拉角并未指定旋转顺序和旋转轴，只要求两个连续的旋转轴必须垂直，因此旋转轴有多种组合。另外欧拉角还分静态定义和动态定义，静态定义指其旋转轴为世界坐标系，动态定义指其旋转轴为模型坐标系，静态定义不会发生万向节死锁（万向节死锁稍后解释），动态定义会。这里通过经典力学中使用的动态Z-X-Z规范做图形说明。

上图中蓝色部分xyz为世界坐标系，红色表示的XYZ为模型坐标系，N为XYO平面和xyO平面的交线。则该模型变换过程可以描述为，以模型坐标系为参考，首先绕Z轴旋转α度，再绕X轴旋转β度，最后绕Z轴旋转γ度。

但是对于3D图形旋转而言，通常使用的是Y-X-Z，即首先进行Yaw(偏航)，再进行倾覆(Pitch)，最后翻滚(Roll)的顺序进行。如下图所示。

4.2.1 万向节死锁

欧拉角有一个致命且无法避免的缺点，就是万向节死锁。即当第一次旋转的轴和第三次旋转的轴相对于世界坐标系重合时，欧拉角会丢失一个自由度，第一次和第三次相对于世界坐标系都绕同一个轴旋转。在上图中，假如第二次Pitch的角度为90°，此时第一次Yaw的角度和第三次Roll的角度产生的都是滚筒旋转。在处理这个问题时，通常在第二次旋转为90°时，第一次强制旋转为0°，让所有的Roll效果由第三次旋转实现。

称为万向节死锁的的原因是在早期航海时代需要使用类似下图陀螺仪来确定当前的方位，当船体倾覆的角度为90度时中间点两个圆环都在同一个竖直面上，此时再发生滚动运动时，陀螺仪已经失去调节作用了，因此讲将种情况称为万向节死锁。

4.2.2 欧拉角的优点

• 易于使用：欧拉角形象的表示出模型的朝向。
• 使用内存最少：欧拉角保存一个位置只需要使用三个数。
• 任意三个数都是合法的：任意三个角度都能找出一个对应的方位。

4.2.3 欧拉角的缺点

• 别名问题：同一个方位可以有无穷多个欧拉角组合，别名问题有三个原因。

第一，将一个角度加上360°后，方位不会改变，但是数值改变了。

第二，三个角度之间强耦合性，这使得如（Yaw 0, Pitch 135°, Roll 0）和（Yaw 180°, Pitch 45°, Roll 180°）表示同一个方位。

为了保证任意欧拉角都有唯一表示，通常对偏航、倾覆和翻滚的角度做出限制，令Yaw和Roll在180°到-180°之间，令Pitch在90°到-90°之间。

第三，由于万向节死锁使得（Yaw α, Pitch 90°, Roll β）和（Yaw η, Pitch 90°, Roll θ）欧拉角组合，只要满足α-β等于η-θ，这两个欧拉角表示同一个方位。为了解决这个问题，规定只要当pitch为±90°时，roll恒等于0。

• 插值问题：引发插值困难的因素有三个。

第一，例如在10°和300°之间插值会绕优弧插值，并不是最短路径，此问题可以通过限制欧拉角解决。

第二，在使用限制欧拉角后，对于170°和-170°之间插值，仍会绕优弧行走，解决办法是将其角度差α映射到-180°到180°之间。公式为wrap(α) = α - 360((α+180) / 360)。其中括号内除法得到的是整数。

第三，当发生万向节死锁时，无法在两个欧拉角之间进行插值，并且该问题无法解决。只能通过四元数的方式表示方位角实现球面插值。

4.3 四元数

四元数起源于复数，复数对(a，b)定义了复数P(a+bi)，其共轭复数定义为P’(a-bi)。复数的模定义为||P|| = square(P*P')，||P|| = square(a^2 + b^2)。将复数的实部作为二维笛卡尔坐标系的横轴，虚部作为其纵轴，构建一个复平面，任意复数都可以表示为该平面上的点。此时引入新的复数q，通过复数乘法能够表示平面上的旋转。

爱尔兰数学家William Hamilton一直致力于将复数从2D扩展到3D，最后他扩展了复数系统，定义复数P(w, xi, yj, zk)。其中i、j、k关系如下。该复数可以用四元数M[w, (x, y, z)]表示。类似于复数能够旋转2维向量，四元数能旋转3D向量。

一个四元数包含一个标量分量和一个3D向量分量，有如下两种表示方式。

4.3.1 四元数和轴-角对

3D中任意角位移都能表示为绕单一轴旋转，这种描述方式称为轴-角描述法，绕轴n旋转θ度的角位移可以用轴角对(n, θ)表示。四元数能够被解释为角位移的的轴角对方式。即上述轴角对可以用下面的四元数表示。

4.3.2 负四元数

四元数求负即将其各个分量分别取负，在空间几何中四元数M和-M表示相同的角位移。

4.3.3 单位四元数

数学上的单位四元数为M [1, 0]，任何四元数乘以单位四元数都等于自身。在空间几何中还有单位四元数N [-1, 0]，表示不发生任何角位移。

4.3.4 四元素运算

四元数的模
对于任意四元数，其求模公式如下。

四元数的共轭和逆
四元数的共轭表示为q*，将四元数的向量部分取负即可得到。四元数的逆表示为q^-1。q^-1 = q的共轭除以q的模。因此对于单位四元数，其共轭和逆四元数是相等的。

四元数的对数、指数
对于如下四元数

其对数公式定义如下

指数以相反的方式定义，假设存在如下四元数

则其指数公式表示如下

四元数的幂运算
当四元数q的指数从0变化为1时，其值从 [1, 0] 变化为q。这里需注意当q代表的角位移角度为30°时，q²代表60°，q⁴代表的是240°，(q⁴)^1/2表示的是120°，可以看出(q⁴)^1/2不等于q²，因为数学上关于指数运算的代数公式对于四元数的求幂运算都不适用。

四元数和标量乘法运算
四元数和标量的乘法运算较为简单，直接将标量和四元素的各个分量相乘即可。

四元数间的叉乘
四元数的叉乘公式如下，需注意的是四元数的乘法满足结合律，不满足交换律。另外由乘法公式可以推导四元数乘积的模等于模的乘积。另外如未写两个四元数乘法的符号，默认都是叉乘。

四元数的点乘
四元数的点乘运算方式类似于向量的点乘，其公式如下。其几何意义可以表示为四元数a和b的点乘绝对值越大，它们代表的角位移越相似。

四元数的差
四元数的差被定义为一个方位到另一个方位的角位移，给定方位a和b，从a经过角位移d变化到b可用公式表示为ad = b。此时可计算出d = a^-1b。

4.3.5 四元数旋转空间向量

对于空间任意向量P，在用四元素计算将其绕单位向量n旋转θ度得到向量p‘时，先将其扩充到四维向量p = [0, P]，再引入角轴对四元数q [cos(θ/2), nsin(θ/2)]。其计算公式为p’ = qpq^-1，可以通过乘法展开，并将p'和使用旋转矩阵方式得到的结果比较从而证明该等式成立。实际上，该方法的时间成本和将四元数转换为旋转矩阵后再旋转向量的时间成本相当。

四元数的乘法可以连接多个角位移，当p分别经过轴角对a和b旋转，其旋转公式表示如下。需要注意这里的四元数乘法是叉乘。

4.3.6 四元数的插值

Slerp插值
欧拉角和一般线性插值都不能用于空间球面插值，四元数的Slerp(Spherical Linear Interpolation)可以用于空间球面平滑插值。两个标量之间的一般线性插值可以如下定义。
d = a₁ - a₀； a_t = a₀ + td

同样的对于四元数的插值也可以定义如下。
lerp(a₀, a₁, t) = a₀ + td

首先计算两个值的差：上文讲过计算四元数q0到q1的差可以由角位移d = q0^-1q1得出。
然后计算差的一部分：四元数的幂运算可以将四元数等分，dt = d^t。
然后在初值上加上偏移量：通过四元数的乘法可以组合角位移。
最后，得出四元数的插值公式为。
slerp(q₀, q₁, t) = q0 (q₀^-1 q₁)^t

在实际应用中，可以通过空间几何推导出四元数插值一般公式。对于向量v₀，旋转w可以得到向量v₁，对于任意0到1之间的时间点t，其旋转后的向量表示为v_t。

首先可以得到等式：v_t = k₀v₀ + k₁v₁
根据上图由空间几何可以推导出k₁ = sin(tw) / sin(w)
再推导出k₀ = sin((1-t)w) / sin(w)
最后可以推导出v_t = (sin((1-t)w) / sin(w)) v₀ + (sin(tw) / sin(w)) v₁

将该公式推导至四元数，因此四元数的插值一般公式可以表示为
slerp(q₀, q₁, t) = (sin((1-t)w) / sin(w)) q₀ + (sin(tw) / sin(w)) q₁

这里需要注意，当q₀和q₁点乘为负时，其角位移方向为优弧表示的方向，因此需要对其中一个取负，另外当q₀和q₁非常接近时sinw的值会引起不必要的计算误差，此时需使用四元数插值的理论公式。

Squad插值
Slerp提供了在两个方位之间插值，Squad允许在一个方位序列中进行平滑插值，方位序列由一组四元数控制点构成。

此外，还需要引入辅助四元数序列si，将它作为临时控制点，为了为每个控制点都扩展一个对应的临时控制点，另外控制点前后分别增加虚拟控制点q₀ = q₁和q_n+1 = q_n。

引入插值变量h，当h从0变化到1时，squad描述了q_i到q_i+1之间的曲线，整条插值曲线分别用如下公式获得。此处并未深入讨论Squad的细节，如需了解，参考书籍Quaternions, Interpolation and Animation。

4.3.7 四元数的优点

• 平滑插值：Slerp和Squad提供了方位间的平滑插值。
• 快速连接和角位移求逆：四元数叉乘能够将角位移序列转换为单个角位移，比矩阵连接更快。单位四元数的逆可以通过共轭四元数得到，相较于正交矩阵通过转置求逆更快。
• 可以和矩阵快速互相转换
• 内存占用较小：四元数存储一个方位只需要使用四个数。

4.3.8 四元数的缺点

• 四元数可能不合法：坏数据的输入或者浮点数的累积误差都可能引起四元数不合法，可用四元数标准化解决这个问题，即将四元数除以它的模，确保其为单位大小。

四元数对于空间方位描述有着不可替代的作用，在GLKit中可以使用GLKQuaternion相关方法对四元数进行运算。

4.4 方位表达方式之间的互相转换

前文已经说过，OpenGL的底层是通过矩阵变化来旋转物体的，因此无论我们用的是哪一种形式表达刚体的方位，最终都需要转化为矩阵形式。

4.4.1 从欧拉角到仿射旋转矩阵

欧拉角转化到矩阵比较简单，上文旋转矩阵中已经分别介绍了绕各个轴单独的旋转矩阵，只需按规定欧拉角的旋转顺序将各个矩阵相乘即可得到，如x-y-z旋转顺序动态欧拉角对应的模型坐标系向世界坐标系转化的仿射旋转矩阵表示如下。由于是以物体坐标系为参考，因此第一次的绕x轴变换在等式最右侧。从世界坐标系向模型坐标系转换的仿射矩阵使用相反的乘法，即使用下述矩阵的转置矩阵。

4.4.2 从仿射旋转矩阵到欧拉角

根据上一小结中的旋转矩阵，可以直接使用反三角函数求出三个欧拉角，得出Ψ = arcTan(-m₁₂ / m₁₁)，θ = arcSin(m₁₃)，Φ = arcTan(-m₂₃ / m₃₃)。当θ = ±90°时，上述4个分量都为0，此时计算出θ = arcSin(m₁₃)，再令Ψ = 0，在通过矩阵计算出Φ的值。

4.4.3 从四元数到矩阵

前文已经得到绕任意轴旋转的矩阵表达式如下。

表示绕任意轴旋转的四元数q = [w, x, y, z]可以表示如下。

此时通过数学上的一些技巧可以将矩阵中各个元素用四元数分量表示出来。对于对角线上的元素如m₁₁，可以通过下面的推导公式求出：

进一步求解如下

对于非对角线的元素如m₁₂，可以通过以下推导获得。

因此最后可以写出用四元数转换为的旋转仿射矩阵表示如下

4.4.4 从矩阵到四元数

根据上述公式，各个元素之间代数运算可得到下述关系

首先解平方根确定一个分量w，再根据如下等式计算剩余分量，这样会有两个不同的解，但是它们代表了相同的角位移。

4.4.5 从欧拉角到四元数

将绕x-y-z轴旋转p、h、b度的四元数表示如下。

假定欧拉角为y-x-z的动态旋转顺序，则表示描述物体向世界坐标系转换的四元数q = BPH表示如下。

4.4.6 从四元数到欧拉角

以物体坐标系为参照，分别绕x-y-z轴旋转ϕ-θ-ψ度的四元数可以通过上面公式解方程求出p-h-b(即ϕ-θ-ψ)的值，但是通常利用矩阵转换到欧拉角到公式如下。

再使用如下四元数到矩阵的转换公式得出矩阵元素的值，最后计算出欧拉角。注意，不同旋转顺序其欧拉角有不同的解。

5 OpenGL中的空间变换

OpenGL在渲染3D模型时需要将其从模型坐标系转换到窗口坐标系，其转换过程如下。

5.1 模型坐标系向世界坐标系转换

一个向量经过多个连续的空间变换可以通过矩阵相乘合成一个变换，然后右乘列向量。但是这里需要注意的是，矩阵乘法的顺序非常重要，当基于模型坐标系做动态变化时，矩阵的乘法顺序和变化顺序相反，当基于世界坐标系做静态变化时，矩阵的乘法顺序和变化顺序相同。例如基于世界坐标系绕y轴旋转一个模型，然后再将其平移(4,10,-20)的代码应该如下表示。

GLKVector4 inputVector;
GLKMatrix4 rotationMatrix = GLKMatrix4MakeRotation(M_PI_2, 0, 1, 0);
GLKMatrix4 translationMatrix = GLKMatrix4MakeTranslation(4, 10, -20);
GLKMatrix4 compositeMatrix = GLKMatrix4Multiply(rotationMatrix, translationMatrix);
GLKVector4 outputVector = GLKMatrix4MultiplyVector4(compositeMatrix, inputVector);

5.2 模型坐标系向视图坐标系转换

在OpenGL中，模型首先被转换到世界坐标系中，再将其转换到视图坐标系中。将模型从模型空间转换到视图空间称为模型-视图转换(Model-View Transform)，对应的矩阵成为模型视图矩阵(Model-View Matrix)。由模型-世界转换(Model-World Transform)以及世界-视图转换(World-View Transform)合成，可以通过世界视图矩阵(World-View Matrix)右乘模型世界矩阵(Model-World Matrix)计算得到。

5.2.1 观察矩阵

世界视图矩阵也可以称为观察矩阵(Lookat Matrix)，通过定义相机位置，即观察者位置构建。在世界坐标系定义一个相机即观察点Eye，观察中心Center，相机向上参考向量up，它可以不严格和Y轴正方向重合，它们在有xyz构建的世界坐标系中表示如下。

则视图坐标系的基向量XYZ在世界坐标系中的表示可以由如下公式推导。

根据前面讲到的仿射矩阵的构建知识，从相机坐标系到世界坐标系的3✖️3转换矩阵R_{相机->世界}可以由相机坐标系3个轴的基向量在世界坐标系总的坐标构建如下。

而由世界坐标系向相机坐标系的3✖️3转换矩阵R_{世界->相机}等于R_{相机->世界}的转置矩阵，可以表示如下。

最后我们需要使用的4✖️仿射矩阵可以定义如下。

由世界坐标系中相机的位置在相机坐标系中是原点，因此通过以下计算能够计算出abc的值，从而得到最终的由世界坐标系向相机坐标系转换的仿射矩阵T_{世界->相机}。

在GLKit中可以调用如下函数构造该仿射矩阵。

GLKMatrix4MakeLookAt(float eyeX, float eyeY, float eyeZ,
                     float centerX, float centerY, float centerZ,
                     float upX, float upY, float upZ)

5.3 视图坐标系向标准设备坐标系转换

当模型顶点坐标被转换到视图坐标系后，OpenGL会进行投影变换操作，该操作后模型的坐标就会被映射到裁剪空间。投影变换(Projection Transformations)可以分为正交投影和透视投影两种，通常我们在项目中为了体现出模型近小远大的自然特征而采用后者。

5.3.1 齐次坐标系

在正式了解投影变换之前我们需要首先了解再投影变换中使用到的齐次坐标系。首先我们考虑一个描述“两条平行的直线没有交点”，这个描述在高中学习的基于欧式几何建立的笛卡尔坐标系下是正确的。但是几何分为欧式几何和非欧几何两个大类，其中非欧几何泛指非欧式几何。而欧式几何不同于其他几何的地方就是平行公理，即过直线外一点有且仅有一条平行线。

我们在投影变换中使用的齐次坐标系就是非欧式几何的一种，考虑实际生活中的一个场景，假如我们站在🛤️上并望向远方，我们会发现平行的两根轨道会在无穷远处相交，这也是投影几何中的特点，即平行线相交于其直线方向的无穷远点。

那么我们在笛卡尔坐标系中如何定义无穷点了，一种可能的表示方式为P(∞, ∞, ∞)。这样表示有两个问题，第一符号∞没有实际意义，第二该坐标不能表示无穷远点的方向。为了解决这两个问题，我们对笛卡尔坐标进行扩展，新增第四维度来表示三维空间中的点。

在笛卡尔坐标系中，直线的参数函数可以表示如下。

令t=1/w，则上述函数可以转换为如下形式

则直线上的任意点P的笛卡尔坐标为(m/w, n/w, q/w)，齐次坐标为(m, n, q, w)，可以发现当w趋于0时，其笛卡尔坐标趋于无穷大，因此当使用齐次坐标点P(m, n, q, 0)表示在向量V(m, n, q)方向上的无穷远点。这种表示方法使用了有意义的符号，并同时描述了无穷远点的方向，从而弥补了笛卡尔坐标的劣势。由各个无穷远点组成的面称为无穷远面。

刚刚讲到在笛卡尔坐标系中两条平行线无交点，即下面的函数无解。

使用齐次坐标系改写函数如下。

则上述函数存解P(x, y, z, 0)，其中xyz需要满足函数Ax+By+Cz=0。几何意义上就是在直线两端无穷远处它们相交。

5.3.2 正交投影

在正交投影(Orthographic Projection)中，模型将被垂直投影到屏幕上，即无论观察者距离物体的距离有多远，它投影后的大小都相同。正交投影通常用于二维图形和文字的处理。

正交矩阵
构建好的正交矩阵如下图所示

GLKit中使用如下函数构建正交矩阵

GLKMatrix4 GLKMatrix4MakeOrtho(float left, float right,
                               float bottom, float top,
                               float nearZ, float farZ)

5.3.3 透视投影

透视投影(Perspective Projection)的特点是对于同样尺寸的模型，其距离观察者的距离越远，它被投影到屏幕上的尺寸越小。

透视矩阵(Perspective Matrices)
透视矩阵又被称为平截椎体矩阵(Frustum Matrix)，在OpenGL中通过指定其近矩形平面的四个方向距离中心的距离，以及4个剪切平面（及平截椎体除去近、远剪切面）的坐标可以直接构造如下透视矩阵。这里需要注意left、right、top、bottom是坐标，即有正负，而near、far是观察点到近投影面和远投影面的距离，都为正值。

GLKit中使用如下函数构建透视矩阵

GLKMatrix4MakeFrustum(float left, float right,
                      float bottom, float top,
                      float nearZ, float farZ)

另外，还可以通过指定垂直方向上的可视角度(以pi为单位)，指定水平和垂直宽度比，以及远近剪切面的距离直接构建对称的透视矩阵

GLKMatrix4 GLKMatrix4MakePerspective(float fovyRadians, 
                                     float aspect, 
                                     float nearZ, float farZ)

5.4 从投影坐标系转换到标准设备坐标系

经过投影变换，顶点的坐标就呗转换到投影空间中，前面讲过将齐次坐标系转换为笛卡尔坐标系的方法是讲其xyz分量都除以w分量即，P点点齐次坐标(x, y, z, w)转换得到的笛卡尔坐标为(x/w, y/w, z/w)，这个操作称为透视除法。经过这个变换后模型的顶点坐标就被转换到了标准设备坐标系中。

标准设备坐标系以设备中心为原点，水平向右为x轴正方向，垂直向上为y轴正方向，垂直屏幕朝内为z轴正方向，它们的取值范围都为[-1, 1]。前面章节将图形渲染管道的时候还讲过，这一阶段OpenGL将会裁剪掉不在视野范围内的图元以提升图形渲染性能，即会裁剪掉坐标超出[-1, 1]的图元。

5.5 从标准设备坐标系转换到窗口坐标系

图像最终是需要显示在窗口上的，这意味着我们在渲染的最终阶段必须得到每个片段的窗口坐标，经过前面的几节我们已经讲模型的各个顶点从模型坐标系转换到标准设备坐标系之中了，接下来我们需要讲其转换到窗口坐标系中，如下图。

这个过程在OpenGL中称为视口变换(Viewport Transformation)，其计算过程如下。

其中(x_d, y_d, z_d)表示点P在标准设备坐标系中的坐标，(x_w, y_w, z_w)表示该点在窗口坐标系中的坐标。O_xO_yp_xp_y分别是在设置窗口位置时调用的函数glViewport(_ x: GLint, _ y: GLint, _ width: GLsizei, _ height: GLsizei)时设置的四个参数，n和f则是设置深度区间时调用函数glDepthRange(_ zNear: GLclampd, _ zFar: GLclampd)时设置的两个参数。

6 插值和曲线

6.1 插值

在OpenGL中常用到插值计算(Interpolation)来确定中间状态的值，如前文讲到的使用四元数的插值公式实现曲面插值。此外在绘制一个颜色过渡的三角形图元时，通常我们设置了三个顶点的颜色值，而在片段着色器中我们需要确定每个像素的颜色，此时就需要进行插值计算。插值计算分为线性插值计算和非线性插值计算。

线性插值：对于因变量Y和自变量t的函数，t的取值范围为[0, 1]，Y从A逐渐变换到B，可以表示为y=f(t)，如果这个函数是线性变换的则称将这个计算过程称为线性插值计算。OpenGL在颜色混合等场景中会用到，着色器中可以使用内置函数vec4 mix(vec4 A, vec4 B, float t);计算。

非线性插值：非线性插值和线性插值的区别在于插值函数y=f(t)是非线性的。OpenGL在片段着色器中得到的输入变量默认情况下都是图像渲染管线前端计算得到的顶点属性在投影平面非线性插值的结果，具体逻辑会在片段着色器章节详细讲解。

6.2 曲线

OpenGL中常用的曲线类型是贝塞尔曲线和样条曲线，通过这些光滑的曲线我们可以模拟出一些有趣的现象，如迎风招展的旗子。

6.2.1 贝塞尔曲线

贝塞尔曲线(Bezier Curve)由一个起点一个终点以及其中一系列控制点组成。由起点A，终点C，控制点B描绘出贝塞尔曲线如下图。

当t从0变化到1时，对于P点，其满足下列公式：
D = A + t(B - A)
E = B + t(C - B)
P = D + t(E -D)
可以得到：
P = A + 2t(B - A) + t²(C - 2B + A)

在着色器中可以使用多个线性插值函数绘制由1个控制点构成的2次贝塞尔曲线：

vec4 quadratic_bezier(vec4 A, vec4 B, vec4 C, float t) {
    vec4 D = mix(A, B, t); 
    vec4 E = mix(B, C, t);
    vec4 P = mix(D, E, t); 
    return P;
}

由起点A，终点D，控制点BC描绘的贝塞尔曲线如下图。

当t从0变化到1时，对于P点，其满足下列公式：
E = A + t(B - A)
F = B + t(C - B)
G = C + t(D - C)
H = E + t(F - E)
I = F + t(G - F )
P = H + t(I - H)

在着色器中可以使用多个线性插值函数绘制由两个控制点构成的3次贝塞尔曲线：

vec4 cubic_bezier(vec4 A, vec4 B, vec4 C, vec4 D, float t) {
    vec4 E = mix(A, B, t);      //E=A+t(B-A) 
    vec4 F = mix(B, C, t);      //F=B+t(C-B)
    vec4 G = mix(C, D, t);      //G=C+t(D-C)
    vec4 H = mix(E, F, t);      //H=E+t(F-E) 
    vec4 I = mix(F, G, t);      //I=F+t(G-F)
    vec4 P = mix(H, I, t);      //P=H+t(I-H)
    return P;
}

可以简写为：

vec4 cubic_bezier(vec4 A, vec4 B, vec4 C, vec4 D, float t) {
    vec4 E = mix(A, B, t); 
    vec4 F = mix(B, C, t); 
    vec4 G = mix(C, D, t);
    return quadratic_bezier(E, F, G, t); 
}

更高阶的贝塞尔曲线都可以按照这样的方式构建，但是通常高于包含4个点的时曲线的描绘就不再使用贝塞尔曲线了，而使用样条曲线。含有5个点的贝塞尔曲线构建方式如下：

vec4 quintic_bezier(vec4 A, vec4 B, vec4 C, vec4 D, vec4 E, float t) {
    vec4 F = mix(A, B, t); 
    vec4 G = mix(B, C, t); 
    vec4 H = mix(C, D, t); 
    vec4 I = mix(D, E, t);
    return cubic_bezier(F, G, H, I, t); 
}

6.2.2 样条曲线

样条曲线(Splines)是由多个小的曲线(如贝塞尔曲线)组成，每个曲线称为片段(Segment)，片段两头的控制点称为焊点(Welds)，片段内部的点称为节点(Knots)，焊点一定被相邻片段共享，有时节点也会被相邻片段共享。

一个由三个贝塞尔曲线组成的B类立方样条曲线(Cubic B-spline)(ABCD)-(DEFG)-(GHIJ)如下图所示。注意这里CDE共线，FGH共线，其中D、G都为各自线段的中点。

当t从0分别变化到1、2、3时，P分别变化到D、G和J。其样条曲线生成函数可以如下构建。使用改方程时，匀速变化的t会得到匀速变化的p，反之p将不会匀速插值。

vec4 cubic_bspline_10(vec4 CP[10], float t) {
    if (t <= 0.0) {
        return CP[0];
    }
    if (t >= 1.0) {
        return CP[9];
    }

    float f = t * 3.0;
    int i = int(floor(f));   //  取整数
    float s = fract(f);      //  取小数

    vec4 A = CP[i * 3]; 
    vec4 B = CP[i * 3 + 1]; 
    vec4 C = CP[i * 3 + 2]; 
    vec4 D = CP[i * 3 + 3];
    return cubic_bezier(A, B, C, D, s); 
}

如果一个函数具有连续的一阶导数(First Derivative)则被称为C¹连续，如果它具有连续的二阶导数(Second Derivative)称为C²连续。贝塞尔曲线同是C¹和C²连续，为了保证焊点位置是连续的，每个焊点出发到下一个节点的方向和t变化速率都必须同上一个节点到该焊点的对应方向和速率(焊点为两个相邻节点的中点)相同。满足上述条件样条函数也满足C¹和C²连续。这样的B类立方体样条曲线也称为三次埃米尔特样条插值(Cubic Hermite Spline)简称为(Cspline)。