《Mathematics for 3D Game Programming and Computer Graphics》读书笔记 (上)

2019-07-24 todaylg读书流水

学习的过程中深感自己的地基不够牢固，决定补习巩固一波基础知识。

这本书涵盖了图形学相关的各方面数学知识，希望通过阅读这本书籍能补回丢掉的数学知识，也可以提高自己英语书籍的~~疯狂查词典~~阅读能力。

本着实践为主的目的，尽量能实践的就搂个Example，慢慢整吧～

公式编辑使用LateX，在线预览可在这里，养成好习惯从现在开始，强烈推荐。

Chapter 1：渲染管线

首先复习一波最重要的渲染管线（From：RTR4）：

定义：给定虚拟摄像机、三维物体、光源等的条件下，生成或者渲染出一幅二维图像的过程。

阶段划分：

应用程序阶段(Application) ：开发者能够完全控制的阶段，通常在此阶段实现视椎体裁剪、碰撞检测、输入检测、加速算法、动画等。
几何处理阶段(Geometry Processing)
- 顶点处理（不同坐标系的转换）
- 图元装配（顶点组成图元）
- 裁剪（对图元进行视体裁剪，舍弃无法成像的图元）
光栅化阶段(Rasterization)：由图元输出片元，三维顶点转二维像素
- 三角形设定阶段(Triangle Setup)：计算三角形表面到差异和三角形表面到其他相关数据
- 三角形遍历阶段(Triangle Traversal)：进行逐像素检查操作，找到哪些采样点或者像素位于三角形中到过程
像素处理阶段(Pixel Processing)
- 像素着色(Pixel Shading)：所有逐像素的着色计算都在像素着色阶段进行
- 融合阶段(Merging)：合成当前缓冲区中由之前的像素着色阶段产生的信息（颜色/深度/Alpha/模板/帧缓存）

图形处理

OpenGL支持绘制的图元有10种：

Points
Lines、Line Strip、Line Loop
Triangles、Triangle Strip、Triangle Fan
Quads、Quad Strip、Polygon

WebGL只支持其中的7种图元（除开Quads、Quad Strip、Polygon）

CPU与GPU间的通信：

CPU => App（Rendering Commands） => OpenGL/DX => Graphics Driver （Command Buffer）=> GPU

CPU内存与GPU内存的通信

CPU内存 => Vertex/Texture Data、Shader Parameters => GPU内存 VRAM(Video Random Access Memory)，存储例如Image Buffer、Depth/Stencil Buffer、Texture Maps、Vertex Buffers等数据，并支持front/back buffer swap

顶点变换

坐标空间：

局部空间(Local Space/Object Space)：局部坐标是对象相对于局部原点的坐标
世界空间(World Space)：世界空间坐标相对于世界的全局原点
观察空间(View Space/Eye Space)：使每个坐标都是从摄像机角度进行观察得到的
裁剪空间(Clip Space)：裁剪坐标会被处理至-1.0到1.0的范围内，并判断哪些顶点将会出现在屏幕上
屏幕空间(Window Space/Screen Space)：视口变换将位于-1.0到1.0范围的坐标变换到由glViewport函数所定义的坐标范围内

坐标变换：

模型(Model)矩阵：局部空间=>世界空间
观察(View)矩阵：世界空间=>观察空间
投影(Projection)矩阵：观察空间=>裁剪空间

MVP矩阵将Object Space转为齐次的Clip Space之后，最后经过视口变换（Viewport transform）最终转到 Window Space：

齐次的Clip Space中顶点需转换为NDC(Normalized Device Coordinates) 即x、y、z => [-1,1]
最终的viewPort Transform将规范化坐标映射到视图所覆盖的像素坐标实际范围，z坐标映射到浮点范围 => [0,1]

片元计算

片元集合了经过GPU计算的深度信息、插值后的顶点信息、插值后的UV坐标信息、像素本身的位置信息。

Face Culling => Rasterization => Fragment Shading => Fragment Operations

Fragment Operations：

像素所有权测试（pixel ownership test）：像素所有权测试只是确定片元是否位于当前可见的viewport区域内，比如其他窗口挡住了当前viewPort的一部分的情况，此测试不可diable
裁剪测试（scissor test）：应用程序可以在视图端口中指定一个矩形，称为剪切矩形，任何落在剪切矩形外的碎片都将被丢弃。
透明度测试（alpha test）：alpha测试将片段的最终alpha值与应用程序预先设置的常量值进行比较
模板测试（stencil test）：根据模板缓冲的值决定是否通过
深度测试（depth test）：比对当前片元的深度值与深度缓冲中的深度值
混合（Blending）：当所有的测试通过后，片元的最终颜色经过混合（Blending）最终输出到 Image Buffer。Blending操作会通过结合当前片元的最终颜色和片段所在位置的图像缓冲区中已经存储的颜色来计算新的颜色（可以实现透明度效果之类的效果），Alpha值也是如此。

Chapter 2：向量

点乘

点乘公式：

$P \cdot Q = \sum_{i=1}^{n}P_{i}Q_{i}$

$P \cdot Q = |P||Q|cos(\alpha)$

证明过程：

由三角函数公式：

$c^2 = a^2+b^2-2abcos(\alpha)$

可得：

$|P - Q |^{2}= \left | P \right |^{2}+\left | Q \right |^{2}-2\left | P \right |\left | Q \right |cos\left [ \alpha \right ]$

化简即可得。

常用作用：

1.点乘结果正负表示两个向量的相对方向

PQ = 0 : 二者互相垂直
PQ > 0 : 二者同向
PQ < 0 : 二者反向

2.点乘结果表明两个向量角度的接近程度（即cos(a)的大小）

向量投影

P在Q上的投影(projection)向量：

${proj_{Q}}^{P}= \frac{P \cdot Q}{|Q|}\times \frac {Q}{|Q|} =\frac{P \cdot Q}{|Q|^2}Q$

其中Q/|Q|代表Q向量的单位长度。

向量相减可得垂线(perpendicular)向量：

${perp_{Q}}^{P}= P - {proj_{Q}}^{P} =P-\frac{P \cdot Q}{|Q|^2}Q$

叉乘

叉乘公式：

$P\times Q = \langle P_{y}Q_{z}-P_{z}Q_{y}, P_{z}Q_{x}-P_{x}Q_{z}, P_{x}Q_{y}-P_{y}Q_{x} \rangle$

$P \times Q = |P||Q|sin(\alpha)$

也可写作：

$P\times Q = \begin{bmatrix} 0 &-P_{z} &P_{y} \\ P_{z} &0 &-P_{x} \\ -P_{y} &P_{x} & 0 \end{bmatrix}\begin{bmatrix} Q_{x}\\ Q_{y}\\ Q_{z} \end{bmatrix}$

或伪行列式（好记）：

$P\times Q = \begin{vmatrix} i &j &k \\ P_{x} &P_{y} &P_{z} \\ Q_{x} &Q_{y} & Q_{z} \end{vmatrix}$

其中i=(1,0,0)、j=(0,1,0)、k=(0,0,1)

常用作用：

叉乘返回垂直于两个向量的新向量，常用与求解表面法线。

向量空间

施密特正交化（实践中暂未发现要用到施密特正交化的地方。。）

Chapter 3：矩阵

矩阵相乘：

$(FG)_{ij} = \sum_{k=1}^{m}F_{ik}G_{kj}$

列行：左列=右行为相乘前提条件
行列：左i行 * 右j列 = Cij

最简形矩阵条件：

对于每一个非零行，最左边的非零项（前导项）为 1
无零行在有零行之前
前导项所在的列上不再有其他的非零项
任意一对非零行i1，i2（i2>i1）对应的包含其前导项的列j1，j2需要满足j2>j1（即列位置在右方）

初等行变换：

交换矩阵的某两行
矩阵的某一行乘以非零数
将矩阵的某一行的倍数加到另一行

初等变换化简矩阵步骤

A.将第一列中最大数值（绝对值）的行交换至第一行
B.将第一行的先导项化为1（乘系数）
C.通过相乘相加第一列的方法将其他行的首个元素化为0
D.继续将第二行的先导项化为1
F.通过相乘相加第二列的方法将其他行的第二个元素化为0（包括第一行）
G.如此循环

求逆矩阵：

前提：含有全零行的矩阵不可逆

A.构造分块矩阵（A|E），E为单位阵
B.对矩阵（A|E）实施初等行变换，将其化为行最简形矩阵（即A=>E）
C.如果A不能行等价于E，则A不可逆。若A能行等价于E，则A可逆，且E就行等价于A^-1

行列式

计算公式：

$detM = \sum_{i=1}^{n}M_{ik}C_{ik}(M)$

其中：

$C_{ij}(M) = (-1)^{i+j}detM^{\{i,j\}}$

M^{i,j}表示从矩阵M中去掉第i行和第j列后的矩阵。

上面两个式子递归即可求得结果，比如2x2矩阵：

$\begin{vmatrix} a &b \\ c &d \end{vmatrix} = ad-bc$

3x3矩阵：

$\begin{vmatrix} a_{11} &a_{12} &a_{13} \\ a_{21} &a_{22} &a_{23} \\ a_{31} &a_{32} &a_{33} \\ \end{vmatrix} = a_{11}(a_{22}a_{33}-a_{23}a_{32})-a_{12}(a_{21}a_{33}-a_{23}a_{31})+a_{13}(a_{21}a_{32}-a_{22}a_{31})$

矩阵可逆的充分必要条件为：

$detM \neq 0$

方程组有非零解的充分必要条件是

使用行列式计算逆矩阵的方法：

$G_{ij} = \frac{C_{ij}(F)}{detF}$

初等行变换对行列式结果的影响：

交换两行 => 行列式结果取反
行乘系数a => 行列式结果也乘系数a
行与行之间相加对行列式结果没有影响

特征值与特征向量：

n阶矩阵M，若数λ和n维非零列向量α使下列的关系式成立：

$M\alpha = \lambda \alpha$

那么数λ称为矩阵M的特征值，非零向量α称为M对应于特征值λ的特征向量，比如：

$M\alpha = \begin{bmatrix} -1 &2 &0 \\ 0 &3 &0 \\ 2 &1 &-1 \\ \end{bmatrix}\begin{bmatrix} 1 \\ 2 \\ 1 \\ \end{bmatrix} = \begin{bmatrix} 3 \\ 6 \\ 3 \\ \end{bmatrix} = 3\begin{bmatrix} 1 \\ 2 \\ 1 \\ \end{bmatrix}$

那3就是矩阵M的特征值，α是M对应于特征值3的特征向量。

将式子移到左边可得：

$(M-\lambda I) \alpha= 0$

α为非零向量的话，(M-λI)必须不可逆，即：

$det(M-\lambda I) = 0$

求得非零解λ后重新代入方程：

$(M-\lambda I) \alpha= 0$

即可求得特征向量α。

举例说明，例如：

$M = \begin{bmatrix} 1 &1 \\ 3 &-1 \\ \end{bmatrix}$

矩阵M-λI即为：

$M-\lambda I = \begin{bmatrix} 1-\lambda &1 \\ 3 &-1-\lambda \\ \end{bmatrix}$

可得：

$det(M-\lambda I) = (1-\lambda)(-1-\lambda)-3 = \lambda ^2 -4 = 0$

所以矩阵M的特征值为λ1=2、λ2=-2。将其代入方程求特征向量可得：

$\begin{bmatrix} 3 &1 \\ 3 &1 \\ \end{bmatrix} \alpha_{1} = \begin{bmatrix} 0 \\ 0 \\ \end{bmatrix}$

$\begin{bmatrix} -1 &1 \\ 3 &-3 \\ \end{bmatrix} \alpha_{2} = \begin{bmatrix} 0 \\ 0 \\ \end{bmatrix}$

可求得特征向量α1与α2：

$\alpha_{1} = \begin{bmatrix} 1 \\ -3 \\ \end{bmatrix}$

$\alpha_{2} = \begin{bmatrix} 1 \\ 1 \\ \end{bmatrix}$

对角化

如果V1,V2,…Vn是一个n×n矩阵M的线性无关特征向量，那么相似变换矩阵A可以表示为：

$A = \begin{vmatrix} V_{1} &V_{2} &...&V_{n} \end{vmatrix}$

相似对角化矩阵M为：

$A^{-1}MA = \begin{vmatrix} \lambda_{1} &0 &...&0\\ 0 &\lambda_{2} &...&0\\ ... &... &...&...\\ 0 &0 &...&\lambda_{n}\\ \end{vmatrix}$

其中λ1、λ2…λn为矩阵M的特征值。

Chapter 4：变换

正交矩阵

正交矩阵的条件：

$M^{-1} = M^{T}$

正交矩阵保持其长度和角度，所以只受旋转和反射变换影响。？

平移(Translate)矩阵：

$\begin{bmatrix} 1 &0 &0, &a \\ 0 &1 &0, &b \\ 0 &0 &1, &c \\ 0 &0 &0, &1 \\ \end{bmatrix}$

缩放(Scale)矩阵：

$\begin{bmatrix} a &0 &0, &0 \\ 0 &b &0, &0 \\ 0 &0 &c, &0 \\ 0 &0 &0, &1 \\ \end{bmatrix}$

旋转(Rotation)矩阵：

首先求二维旋转矩阵，由球坐标系公式：

$x=rcos(\varphi ) \space y = rsin(\varphi )$

假设旋转的角度为θ，由三角函数正弦余弦公式可得：

${x}' = rcos(\theta+\varphi) = rcos(\theta)cos(\varphi)-rsin(\theta)sin(\varphi)$

${y}' = rsin(\theta+\varphi) = rsin(\theta)cos(\varphi)+rcos(\theta)sin(\varphi)$

代入x、y可得：

${x}' = xcos(\theta)-ysin(\theta)$

${y}' = xsin(\theta)+ycos(\theta)$

使用矩阵表示即为：

$\begin{bmatrix} {x}' \\ {y}' \\ \end{bmatrix} = \begin{bmatrix} cos(\theta) &-sin(\theta)\\ sin(\theta) &cos(\theta)\\ \end{bmatrix}\begin{bmatrix} {x} \\ {y} \\ \end{bmatrix}$

绕任意点的二维旋转可以通过级联两次位移矩阵计算（旋转点位移至原点 => 绕原点旋转 => 旋转点移回原位）

通过齐次坐标引入绕三个轴的三维旋转矩阵：

$R_x(\theta) = \begin{bmatrix} 1&0 &0\\ 0 &cos(\theta)&-sin(\theta)\\ 0&sin(\theta) &cos(\theta)\\ \end{bmatrix}$

$R_y(\theta) = \begin{bmatrix} cos(\theta) &0&-sin(\theta)\\ 0&1 &0\\ -sin(\theta)&0 &cos(\theta)\\ \end{bmatrix}$

$R_z(\theta) = \begin{bmatrix} cos(\theta) &-sin(\theta)&0\\ sin(\theta)&cos(\theta)&0\\ 0&0 &1\\ \end{bmatrix}$

围绕任意轴的旋转：

有两种思路：

通过将旋转轴旋转到与X/Y/Z某一轴重合，再执行基本旋转，之后再执行反向旋转恢复旋转轴方向。
计算旋转点在旋转轴上的平行和垂直分量，再分别计算分量绕A轴旋转角度后的所增加的分量

求点P绕任意轴A（单位向量，|A|=1）旋转的旋转矩阵：

首先求得在轴A上的投影及垂线分量：

${proj_{A}}^{P}= \frac{P \cdot A}{|A|}\times \frac {A}{|A|} = (A\cdot P)A$

${perp_{A}}^{P}= P - {proj_{A}}^{P} =P-(A \cdot P)A$

垂线分量（与A垂直）绕A旋转θ后的结果为：

$[P-(A \cdot P)A]cos(\theta) + (A \times P)sin(\theta)$

再加上A轴上的投影分量：

${P}' = A(A \cdot P) + [P-(A \cdot P)A]cos(\theta) + (A \times P)sin(\theta) = Pcos(\theta) + (A \times P)sin(\theta) + A(A \cdot P)(1-cos(\theta))$

最后代入前面的点乘叉乘的矩阵表示可得：

${P}' = \begin{bmatrix} 1&0 &0\\ 0 &1&0\\ 0&0 &1\\ \end{bmatrix}Pcos(\theta) + \begin{bmatrix} 0&-A_{z} &A_{y}\\ A_{z}&0 &-A_{x}\\ -A_{y}&A_{x} &0\\ \end{bmatrix}Psin(\theta) + \begin{bmatrix} A_{x}^2&A_{x}A_{y} &A_{x}A_{z}\\ A_{x}A_{y}&A_{y}^2 &A_{y}A_{z}\\ A_{x}A_{z}&A_{y}A_{z} &A_{z}^2\\ \end{bmatrix}P(1-cos(\theta))$

最终围绕任意旋转轴A的旋转矩阵：

$R_A(\theta) = \begin{bmatrix} c+(1-c)A_{x}^2&(1-c)A_{x}A_{y}-sA_{z} &(1-c)A_{x}A_{z}+sA_{y}\\ (1-c)A_{x}A_{y}+sA_{z}&c+(1-c)A_{y}^2 &(1-c)A_{y}A_{z}-sA_{x}\\ (1-c)A_{x}A_{z}-sA_{y}&(1-c)A_{y}A_{z}+sA_{x} &c+(1-c)A_{z}^2\\ \end{bmatrix}$

转换为Shader代码在顶点着色器中使用：

mat4 rotationMatrix(vec3 axis, float angle) {
    axis = normalize(axis);
    float s = sin(angle);
    float c = cos(angle);
    float oc = 1.0 - c;

    return mat4(oc * axis.x * axis.x + c,           oc * axis.x * axis.y - axis.z * s,  oc * axis.z * axis.x + axis.y * s,  0.0,
                oc * axis.x * axis.y + axis.z * s,  oc * axis.y * axis.y + c,           oc * axis.y * axis.z - axis.x * s,  0.0,
                oc * axis.z * axis.x - axis.y * s,  oc * axis.y * axis.z + axis.x * s,  oc * axis.z * axis.z + c,           0.0,
                0.0,                                0.0,                                0.0,                                1.0);
}

vec3 rotate(vec3 v, vec3 axis, float angle) {
    mat4 m = rotationMatrix(axis, angle);
    return (m * vec4(v, 1.0)).xyz;
}

齐次坐标：

三维表示的点和向量可能会被搞混，齐次坐标可以解决这个问题，其对三维空间中的点和向量的表示都是四维的：

任意点(w分量为1)：

$\begin{bmatrix} a_{1}\\ a_{2}\\ a_{3}\\ 1\\ \end{bmatrix}$

任意向量(w分量为0)：

$\begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\\ 0\\ \end{bmatrix}$

两点相减即表示向量（w分量变为0，不受平移矩阵影响）。

所有仿射变换（具有保持直线的特性）都可以借助齐次坐标表示成矩阵相乘的形式。

比如三维表示的平移变换为矩阵加法，采样齐次坐标的话就可以统一表示为矩阵乘法。

3X3的变换矩阵M和3维平移矩阵T可以结合为4X4的变换矩阵F：

$F=\begin{bmatrix} M_{11}&M_{12}&M_{13}&T_{x}\\ M_{21}&M_{22}&M_{23}&T_{y}\\ M_{31}&M_{32}&M_{33}&T_{z}\\ 0&0&0&1\\ \end{bmatrix}$

法向量矩阵：

在经过非正交矩阵变换之后，法向量通常会不再垂直于表面。

法向量(N)与切向量(T)互相垂直，而切向量经过变换后仍是切向量，则可得：

${N}'\cdot{T}' = (GN)\cdot(MT) = (GN)^{T}(MT) = N^{T}G^{T}MT = 0$

由N·T=0（互相垂直）=> G·MT = T => G·M = I(I为单位阵) =>

$G=(M^{-1})^{T}$

因此法线变换矩阵应该使用变换矩阵的逆转置矩阵。

四元数：

介绍可看之前LGL总结的Math部分，这里不再重复，了解一下就好..

绕A轴旋转角度θ的四元数表示（没有对比就没有伤害…）：

$q = cos\frac{\theta}{2}+Asin\frac{\theta}{2}$

应用到点P：

${P}' = qPq^{-1}$

相当于3x3矩阵：

$R_{q} = \begin{bmatrix} 1-2y^2-2z^2&2xy-2wz&2xz+2wy\\ 2xy+2wz&1-2x^2-2z^2&2yz-2wx\\ 2xz-2wy&2yz+2wx&1-2x^2-2y^2\\ \end{bmatrix}$

球内插值（0<=t<=1）：

$q(t) = \frac{sin(1-\theta)}{sin(\theta)}q_{1}+\frac{sin(\theta t)}{sin(\theta)}q_{2}$

Chapter 5：几何体

点线距离(直接计算其垂直分量即可)：

$d^2 = (Q-S)^2 - [proj_{v}(Q-S)]^2 = (Q-S)^2 - [\frac{(Q-S)\cdot V}{V^2}V]^2$

面：

给定法线向量N和点P，则由二者确定的平面需满足（Q为面内点）：

$N\cdot (Q-P) = 0$

再由平面方程：

$Ax+By+Cz+D=0$

可知：A、B、C为法向量N的xyz分量，D=-N·P：

$d=N\cdot Q+D$

当：

d=0时：点Q在平面中
d>0时，点Q在正面上
d<0时，点Q在背面下

线面交点：

设直线过点S，平行于向量V，则直线方程可以表示为：P(t) = S + tV，可得与面交点方程：

$N\cdot P(t) + D = 0$

可得：

$t=\frac{-(N \cdot S+D)}{N\cdot V}$

视椎体:

假设水平视野的角度为α，相机到视椎体近平面的距离e为：

$e= \frac{1}{tan(\frac{a}{2})}$

垂直视野角度β可以表示为：

$\beta = 2tan^{-1}(\frac{\alpha}{e})$

透视插值校正：

深度插值公式：

$\frac{1}{z_3} = \frac{1}{z_1}(1-y) + \frac{1}{z_2}t$

结果表明在透视投影中，深度值z1和z2需要通过对它们的倒数进行线性插值来正确地插值:。

顶点插值校正：

$\frac{b_3}{z_3} = [\frac{b_1}{z_1}(1-t) + \frac{b_2}{z_2}t]$

透视矩阵：

假设点P(Px,Py,Pz,1)是在视椎体中的一个点（视椎体中n为负)，可得在视椎体近平面的坐标x,y：

$\frac{x}{P_x} = -\frac{n}{P_z} \Rightarrow x = -\frac{n}{P_z}P_x$

同理y：

$\frac{y}{P_y} = -\frac{n}{P_z} \Rightarrow y = -\frac{n}{P_z}P_y$

透视矩阵需要将视椎体映射成齐次裁剪空间的立方体，即需要将X轴[l,r] => [-1,1]、Y轴[b,t]=>[-1,1]，可得关系式：

${x}' = 2\frac{(x-l)}{r-l}-1$

${y}' = 2\frac{(y-b)}{t-b}-1$

代入近平面上点xy的关系式并化简后可得：

${x}' = \frac{2n}{r-l}(-\frac{P_{x}}{P_{z}})-\frac{r+l}{r-l}$

${y}' = \frac{2n}{t-b}(-\frac{P_{y}}{P_{z}})-\frac{t+b}{t-b}$

Z轴坐标的变换复杂一些，为了保证线性，我们需要构造z的倒数关系式：

${z}' = A\frac{1}{z}+B$

由-n => -1、-f => 1可得：

$-1 = \frac{A}{-n}+B \space and\space 1 = \frac{A}{-f}+B$

可计算出A、B：

$A = \frac{2nf}{f-n} \space and \space B = \frac{f+n}{f-n}$

最终z坐标的变换公式：

${z}' = -\frac{2nf}{f-n}(-\frac{1}{P_z})+\frac{f+n}{f-n}$

对xyz化简后可得：

$-{x}'P_z = \frac{2n}{r-l}P_{x}-\frac{r+l}{r-l}P_{z}$

$-{y}'P_z = \frac{2n}{t-b}P_{y}-\frac{t+b}{t-b}P_{z}$

$-{z}'P_z = \frac{f+n}{f-n}P_z-\frac{2nf}{f-n}$

使用矩阵表示即为：

$P' = M_{frustum}P = \begin{bmatrix} \frac{2n}{r-l}&0&\frac{r+l}{r-l}&0\\ 0&\frac{2n}{t-b}&\frac{t+b}{t-b}&0\\ 0&0&-\frac{f+n}{f-n}&-\frac{2nf}{f-n}\\ 0&0&-1&0\\ \end{bmatrix}\begin{bmatrix} P_{x}\\ P_{y}\\ P_{z}\\ 1\\ \end{bmatrix}$

正交矩阵：

与投影矩阵同理，只不过x’/y/’z’与x/y/z之间的关系在正交投影下的关系式是直接相等的：

${x}' = \frac{2}{r-l}{x}-\frac{r+l}{r-l}$

${y}' = \frac{2}{t-b}{y}-\frac{t+b}{t-b}$

Z轴也是同理：

${z}' = \frac{-2}{f-n}{z}-\frac{f+n}{f-n}$

最终的变换矩阵即为：

$P' = M_{ortho}P = \begin{bmatrix} \frac{2}{r-l}&0&0&-\frac{r+l}{r-l}\\ 0&\frac{2}{t-b}&0&\frac{t+b}{t-b}\\ 0&0&-\frac{2}{f-n}&-\frac{f+n}{f-n}\\ 0&0&0&1\\ \end{bmatrix}\begin{bmatrix} P_{x}\\ P_{y}\\ P_{z}\\ 1\\ \end{bmatrix}$

斜视锥体深度投影和裁剪：

涉及反射场景的情况一般是通过一个反射平面，反射主相机生成一个虚拟相机对场景进行渲染。

需要注意的是当有物体横跨反射平面的时候，渲染反射的虚拟相机中的场景就会出现问题（会把平面以下的部分渲染到反射结果中）。

最简单的解决方法是让GPU暴露一个自定义选定剪切平面的方法，截断反射表面下的所有几何体。但是很不幸的是并不是所有GPU都支持这个自定义的剪切选项，并且即使支持改起来也会很麻烦（因为要改动代码）。

这时候斜视锥体裁剪的办法就派上用场了，对视椎体进行裁剪，用反射平面替代原来的近裁剪面（即新的斜视椎体），这样反射平面以下的物体都不会渲染（因为在视椎体以外了），其中变换后的新投影矩阵M’与原投影矩阵M的关系：

$M' = \begin{bmatrix} M_1\\ M_2\\ \frac{-2Q_{z}}{C\cdot Q}C+<0,0,1,0>\\ M_4\\ \end{bmatrix}$

其中:

$Q = M^{-1}Q'$

$Q' = <sgn(C_x',C_y',1,1)>$

$C' = (M^{-1})^{T}C$

代码实现：

//reflectorPlane.setFromNormalAndCoplanarPoint
reflectorPlane.normal.copy(normal);
reflectorPlane.constant = -reflectorWorldPosition.dot(reflectorPlane.normal);

//reflectorPlane.applyMatrix4
reflectorPlane.m1.getNormalMatrix(virtualCamera.viewMatrix);
reflectorPlane.v1.copy(reflectorPlane.normal).multiply(-reflectorPlane.constant);
reflectorPlane.v1.applyMatrix4(virtualCamera.viewMatrix);
reflectorPlane.normal.applyMatrix3(reflectorPlane.m1).normalize();
reflectorPlane.constant = -reflectorPlane.v1.dot(reflectorPlane.normal);

//C'
clipPlane.set(reflectorPlane.normal.x, reflectorPlane.normal.y, reflectorPlane.normal.z, reflectorPlane.constant);

let projectionMatrix = virtualCamera.projectionMatrix;
q.x = (Math.sign(clipPlane.x) + projectionMatrix[8]) / projectionMatrix[0];
q.y = (Math.sign(clipPlane.y) + projectionMatrix[9]) / projectionMatrix[5];
q.z = - 1.0;
q.w = (1.0 + projectionMatrix[10]) / projectionMatrix[14];

// Calculate the scaled plane vector
let res = 2.0 / clipPlane.dot(q);
clipPlane.multiply(res);

// Replacing the third row of the projection matrix
projectionMatrix[2] = clipPlane.x;
projectionMatrix[6] = clipPlane.y;
projectionMatrix[10] = clipPlane.z + 1.0 - clipBias;
projectionMatrix[14] = clipPlane.w;

实践反射场景：

Examples: Reflector

Paper

reflect向量计算

Chapter 6：光线追踪

Chapter 7：光照与着色

Light Source

环境光(Ambient Light):

环境光从各个方向以相同的强度照射过来，并且均匀地照亮物体的每一部分。在场景中环境光的系数通常是一个常量（也可以读纹理）。

定向光(Directional Light ):

定向光源是一种从无限远的地方向单一方向照射的光源。定向光通常用于模拟太阳，它的光线可以认为是平行的，且其光强不会随着距离的增加而减弱。

点光源(Point Light)：

$C=\frac{1}{k_c+k_{l}d+k_qd^2}C_0$

其中光源强度为C，d为距离光源的距离，C0为灯的颜色，常量Kc、Kl、Kq控制衰减。

聚光灯(Spot Light)：

$C=\frac{max\{-R\cdot L,0\}^p}{k_c+k_{l}d+k_qd^2}C_0$

其中R为聚光灯方向，L为照射点指向光源的单位向量，指数p控制聚光灯的集中度，控制光照强度随R和- L夹角的增大而减小的速率。

环境光与漫反射光(Ambient and Diffuse Light)：

漫反射(Lambertian)指入射光在表面某一点上的光的一部分向随机方向散射，因为光在每个方向上都是均匀反射的，所以朗伯反射的计算不要考虑观察者的位置因素。

考虑环境光强度与n个光源的贡献，可以将漫反射分量表示为：

$K_{diffuce} = DA+D\sum_{i=1}^{n}C_imax\{N\cdot L_i,0\}$

其中D是表面的漫反射颜色，A代表环境光的颜色，N是反射表面的法向量，L是指向第i个光源的单位向量，Ci是第i个光源的光照强度。

镜面光照(Specular Lighting)：

$K_{specular} = S\sum_{i=1}^{n}C_imax\{N\cdot H_i,0\}^m\space(N\cdot L_i > 0)$

其中S表面的镜面反射颜色，Hi观察者与光源Li的半角向量，m控制反射的锐度，N与Li的点乘结果决定表面是正对光源还是背对光源（背对的话就不用计算了）。

实践可见：Exmaple

纹理映射(Texture Mapping)

Projective Texture Maps：

投射纹理映射的一个应用是模拟聚光灯将图像投射到环境中，随着与聚光灯距离的增加，投影图像变大。这种效果是通过使用一个4×4纹理矩阵将一个物体的顶点位置映射到纹理坐标（s, t, 0, q），这样除以q就得到了正确的二维纹理坐标（s，t），用于对投影图像进行采样。

首先需要乘以转换矩阵，将每个顶点的坐标先转换为以聚光灯位置为原点的坐标系下的坐标，再乘以投影矩阵（仅xy）

//Todo: 待实践

Cube Texture Maps：

立方体纹理贴图通常用于近似环境周围的反射，纹理的采样坐标为三维向量，表示从立方体中心发出的指向要采样的texel的方向向量。

着色模型(Shading Models)：

Blinn-Phong Shading：

实践可见：Example

凹凸映射(Bump Mapping)：

Bump Map Construction
Tangent Space
Calculating Tangent Vectors
Implementation

//Todo: 实践

物理反射模型

Bidirectional Reflectance Distribution Functions（BRDF）
Cook-Torrance Illumination
The Fresnel Factor
The Microfacet Distribution Function
The Geometrical Attenuation Factor
Todo：Implementation

PBR实践可见：Example