Home

工欲善其事，必先利其器。

# 齐次坐标和仿射变换 (Homogeneous Coordinates and Affine Transformations)

在以前的三维空间中，我们常常用向量v=(x,y,z)T\textbf{v}=(x,y,z)^T 表示一个点，或者一个向量。
其中，为什么能表示一个点呢？其实我们省略了一个前提假设：坐标系原点是0\textbf{0}。
对于点（用大写字母表示）和向量（用小写加粗字母表示），其实它们之间可以进行运算。原先 向量集合 和 向量间的运算 构成了 向量空间（vector space） ，现在想办法把点也加进去。
例子很好理解，譬如：

P+vP−QP−Q+v\begin{aligned} & P+\textbf{v}\\ & P-Q\\ & P-Q+\textbf{v} \end{aligned}

第一个表示点加向量等于点，第二个表示两个点相减等于向量，第三个表示点减点加向量等于向量。很显然，只需要把点当作 1，向量当作 0，然后做加减法。如果结果是 1 那么就是点，如果结果是 0 就是向量，否则就是不合法的运算。
那么根据这个性质，我们就可以尝试拓展 向量空间 。

定义：我们定义一个三维欧氏空间的坐标系形如：

b=(e1e2e3O)\textbf{b}=\begin{pmatrix} e_1\\e_2\\e_3\\O \end{pmatrix}

其中，e1,e2,e3e_1,e_2,e_3 为基向量，OO 为原点。这个坐标系可以完整表述空间里的所有点和向量。

注意，我们不要去考虑e1,e2,e3,Oe_1,e_2,e_3,O 具体等于什么，它就是一个符号。如果你尝试去设e1=(1,0,0)Te_1=(1,0,0)^T 实际上已经在做坐标变换了，变换到i^,j^,k^\hat{i},\hat{j},\hat{k} 的坐标系去了，那么i^,j^,k^\hat{i},\hat{j},\hat{k} 又等于什么呢？这样下去将无止尽。所以我们现在就把它们当作一个符号。
那什么是数字呢？坐标是数字。
考虑空间中的一个点，它一定形如：

P=O+a1e1+a2e2+a3e3=(a1a2a31)(e1e2e3O)=aTb\begin{aligned} P&=O+a_1e_1+a_2e_2+a_3e_3\\ &=\begin{pmatrix}a_1&a_2&a_3&1\end{pmatrix} \begin{pmatrix}e_1\\e_2\\e_3\\O\end{pmatrix}\\ &=\textbf{a}^T\textbf{b} \end{aligned}

其中，那个a\textbf{a} 就是点PP 的坐标，是一个四维的实数向量，特点是第四维为 1。
稍微聪明一点也可以知道，向量也可以用这种方式表示，只不过第四维为 0，因为它不需要加上个原点OO，向量不需要绝对位置，只需要方向和大小。

总结：向量和点都有坐标，坐标是一个实数向量，向量的坐标第四维为 0，点的坐标第四维为 1。
当然，后面为了透视矩阵，我们可能扩展齐次坐标，使得第四维可以为任何数。

处理了其次坐标，我们下面来看仿射变换。
我理解的仿射变换本质上就是在做坐标系间的线性变换 + 平移:
1
相较普通的线性变换，可以看到坐标系原点被平移了。考虑上图两个坐标系b,b′\textbf{b},\textbf{b}' 满足：

b=(e1e2e3O)Tb′=(e1′e2′e3′O′)T\begin{aligned} &\textbf{b}=\begin{pmatrix}e_1&e_2&e_3&O\end{pmatrix}^T\\ &\textbf{b}'=\begin{pmatrix}e_1'&e_2'&e_3'&O'\end{pmatrix}^T\\ \end{aligned}

假设存在坐标系间的变换关系：

{e1′=u1e1+u2e2+u3e3e2′=v1e1+v2e2+v3e3e3′=w1e1+w2e2+w3e3O′=O+p1e1+p2e2+p3e3\begin{cases} e_1'=u_1e_1+u_2e_2+u_3e_3\\ e_2'=v_1e_1+v_2e_2+v_3e_3\\ e_3'=w_1e_1+w_2e_2+w_3e_3\\ O'=O+p_1e_1+p_2e_2+p_3e_3 \end{cases}

即向量e1′,e2′,e3′e_1',e_2',e_3' 和点O′O' 在坐标系b\textbf{b} 中都应该有个坐标，分别为u,v,w,p\textbf{u},\textbf{v},\textbf{w},\textbf{p}。那么，我们可以得到：

(e1′e2′e3′O′)=(u1u2u30v1v2v30w1w2w30p1p2p31)(e1e2e3O)\begin{pmatrix} e_1'\\ e_2'\\ e_3'\\ O' \end{pmatrix}=\begin{pmatrix} u_1&u_2&u_3&0\\ v_1&v_2&v_3&0\\ w_1&w_2&w_3&0\\ p_1&p_2&p_3&1 \end{pmatrix} \begin{pmatrix} e_1\\ e_2\\ e_3\\ O \end{pmatrix}

即b′=Mb\textbf{b}'=M\textbf{b}。这就是仿射变换。仿射变换只有 12 个自由度。
只有在第四维被扩展后，才可能有更多的矩阵形式，但是自由度仍为 12 个。
我们还需要考虑坐标在该变化下的变化。考虑有一个点或一个向量，他在坐标系b\textbf{b} 下的坐标为u\textbf{u}，而在坐标系b′\textbf{b}' 下的坐标为v\textbf{v}，即：

uTb=vTb′\textbf{u}^T\textbf{b}=\textbf{v}^T\textbf{b}'

故有：

uTb=vTb′=vTMb\begin{aligned} \textbf{u}^T\textbf{b}&=\textbf{v}^T\textbf{b}'\\ &=\textbf{v}^TM\textbf{b} \end{aligned}

故有坐标变换：

v=(MT)−1u\textbf{v}=(M^T)^{-1}\textbf{u}

而MTM^T 的最后一行为(0,0,0,1)(0,0,0,1)，也是我们更常见的仿射变化形式。

# 基本的仿射变换：平移 (translation)、旋转 (rotation)、缩放 (scaling)

这部分简写，纯做记录。

  • 平移向量(x,y,z)T(x,y,z)^T

T(x,y,z)=(100x010y001z1)\textbf{T}(x,y,z)=\begin{pmatrix} 1 & 0 & 0 & x\\ 0 & 1 & 0 & y\\ 0 & 0 & 1 & z\\ & & & 1 \end{pmatrix}

  • 沿着向量(x,y,z)T(x,y,z)^T 进行缩放，缩放比例是向量模长

S(x,y,z)=(x0000y0000z01)\textbf{S}(x,y,z)=\begin{pmatrix} x & 0 & 0 & 0\\ 0 & y & 0 & 0\\ 0 & 0 & z & 0\\ & & & 1 \end{pmatrix}

  • 旋转只有绕yy 轴比较特殊，因为只有x×z=−yx\times z=-y，所以旋转方向应取反。
    • 绕xx 轴旋转θ\theta 角度：

    Rx(θ)=(10000cos⁡θ−sin⁡θ00sin⁡θcos⁡θ01)R_x(\theta)=\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & \cos\theta & -\sin\theta & 0\\ 0 & \sin\theta & \cos\theta & 0\\ & & & 1 \end{pmatrix}

    • 绕yy 轴旋转θ\theta 角度：

    Ry(θ)=(cos⁡θ0sin⁡θ00100−sin⁡θ0cos⁡θ01)R_y(\theta)=\begin{pmatrix} \cos\theta & 0 & \sin\theta & 0\\ 0 & 1 & 0 & 0\\ -\sin\theta & 0 & \cos\theta & 0\\ & & & 1 \end{pmatrix}

    • 绕zz 轴旋转θ\theta 角度：

    Rz(θ)=(cos⁡θ−sin⁡θ00sin⁡θcos⁡θ0000101)R_z(\theta)=\begin{pmatrix} \cos\theta & -\sin\theta & 0 & 0\\ \sin\theta & \cos\theta & 0 & 0\\ 0 & 0 & 1 & 0\\ & & & 1 \end{pmatrix}

这部分要讨论几个重点话题：

  1. 一切仿射变化都可以由上述三种变化组合而成。
    首先我们可以不考虑平移（即矩阵最后一列不为 0），因为我们总可以平移到坐标原点，再平移回去。那么我们只需要考虑在普通坐标系下（即 3 维矩阵）时，所有变换都可以表示为旋转和缩放的组合，这句话对应的代数定理就是奇异值分解：
    对于任意一个3×33\times 3 的实数矩阵MM，存在奇异值分解：

    M=UΣVTM=U\Sigma V^T

    其中，UU 和VV 都是正交矩阵 (UTU=UUT=VTV=VVT=IU^TU=UU^T=V^TV=VV^T=I)，Σ\Sigma 是对角矩阵。网上有个图片比较常见：
    2
    它说明正交矩阵就对应了旋转变换，而对角矩阵就对应了缩放变换。然而这并不完全正确。
    对于一个正交矩阵AA (正交矩阵∣A∣2=1|A|^2=1)，若∣A∣=1|A|=1 则它是一个旋转变换，但当∣A∣=−1|A|=-1 时，它还有镜象的作用。
    但是对于缩放变换，当缩放向量中有一维是负数时，它就有镜象的作用。譬如S(1,−1,1)S(1,-1,1) 就是关于xOzxOz 平面镜像。
    故奇异值分解本质严格来说，应该是任何变换都是旋转、镜象和缩放的组合，不过镜象也是缩放中的一种特殊情况。
  2. 绕任意轴旋转矩阵的推导。
    绕一个单位轴n^\hat{n} 旋转θ\theta 角度。我们的基本策略是：将物体和旋转轴一起旋转，然后把旋转轴转到与zz 轴重合后，作用Rz(θ)R_z(\theta)，然后再转回原样。即：

    Rn^(θ)=Rx(−θx)Ry(−θy)Rz(θ)Ry(θy)Rx(θx)R_{\hat{n}}(\theta)=R_x(-\theta_x)R_y(-\theta_y)R_z(\theta)R_y(\theta_y)R_x(\theta_x)

    然后我们尝试计算θx\theta_x，就是把n^\hat{n} 旋转到xOzxOz 平面上，就等于把n^\hat{n} 在yOzyOz 平面上的投影旋转到zz 轴：
    3
    然后θy\theta_y 就是旋转到xOzxOz 平面后，由于绕xx 轴旋转故xx 轴分量大小不变，然后旋转到zz 轴：
    4
    所以有：

    Rx(θx)=(10000nznz2+ny2−nynz2+ny200nynz2+ny2nznz2+ny201)Ry(θy)=(nz2+ny20−nx00100nx0nz2+ny201)R_x(\theta_x)=\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & \frac{n_z}{\sqrt{n_z^2+n_y^2}} & \frac{-n_y}{\sqrt{n_z^2+n_y^2}} & 0\\ 0 & \frac{n_y}{\sqrt{n_z^2+n_y^2}} & \frac{n_z}{\sqrt{n_z^2+n_y^2}} & 0\\ & & & 1 \end{pmatrix}\\ R_y(\theta_y)=\begin{pmatrix} \sqrt{n_z^2+n_y^2} & 0 & -n_x & 0\\ 0 & 1 & 0 & 0\\ n_x & 0 & \sqrt{n_z^2+n_y^2} & 0\\ & & & 1 \end{pmatrix}

    至此推导完毕。
  3. 剪切变换（Shear）
    它不是一种基础变换，但也很重要（在后面投影和透视中）：
    5
    剪切变换的矩阵形式：

    Hx(θ)=(1cotθ00010000101)H_x(\theta)=\begin{pmatrix} 1 & cot\theta & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ & & & 1 \end{pmatrix}\\

# 物体的观察 (Viewing)

我们必须理解清楚一个物体是怎样变换到摄像机胶片上的：

  1. 物体在现实中有一个世界坐标，这个坐标独立于观测者，表示了物体现实中的位置。然后物体自身姿态的平移、旋转变换都应在这个坐标系中首先进行。
  2. 当物体自身不再变化，我们就可以观测它。第一步就是把物体从世界坐标系转化为观察坐标系。即我们摄像机也确定了一个坐标系。这一步的目的是把观测者放到原点，把物体放到坐标轴上：
    6
  3. 根据摄像机的参数 (透视还是正交投影，视角范围多大) 确定观测范围 (view volumn)，把范围外的部分剪裁掉 (clipping)。
  4. 最后把观测范围内的物体投影到胶片 (一个矩形，认为是投影平面) 上。

其中，第 1、2、4 步都可以用矩阵运算描述，第 3 步是三角片剪裁算法，这里不做讨论。第 1、2 步的矩阵对应了 Model-view 矩阵，而第四步则是 Projection 矩阵。

  • 第一步的矩阵可以根据你想要物体是什么姿态来决定。
  • 第二步的矩阵就是在第一部分我们讨论的仿射变化的本质：即坐标系的变换。我们可以把坐标都变换到摄像机坐标系下，当然，摄像机确定的坐标系常用摄像机位置 (点)，摄像机朝向 (向量)，** 摄像机上轴 (向量)** 唯一确定。
    确定了从哪个坐标轴变换到哪个坐标轴，就可以方便地进行坐标变换，这部分对应了 GL 库中的 lookAt 函数。其实这个函数就是实现了一个坐标变换，矩阵形式也较好推导。
    7
  • 第四步的投影矩阵是我们重点讨论对象。可分为平行光投影和透视投影两种。
    其中，两种投影又可以细分为正投影和斜投影。

# 平行光投影

# 正投影

顾名思义，这种投影方式特点为：光束平行，且与投影平面垂直。俯视图为：
8
对于这种投影，我们认为摄像机的观测范围为一个长方体。
其实它的投影矩阵很简单，譬如如果透视平面是z=0z=0，那么矩阵就为：

Morth=(1000010000000001)M_{orth}=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}

就是把zz 分量扔掉。然而在 OpenGL 或 WebGL 中，我们往往认为摄像机有一个观测范围。在正平行光投影下，这个观测范围就是一个长方体：
9
然后实际上在剪裁掉物体长方体外的部分后，其实需要把这个长方体连同里面的物体一起变换到一个边长为 2 的正方体，然后再作用MM 进行正平行投影：
10
这个变换到单位正方体的变换也很简单，就是先平移到原点，然后缩放到单位正方体：

T(−left+right2,−top+bottom2,−far+near2)S(2right−left,2top−bottom,2far−near)T(-\frac{left+right}{2},-\frac{top+bottom}{2},-\frac{far+near}{2})\\ S(\frac{2}{right-left},\frac{2}{top-bottom},\frac{2}{far-near})

所以正投影的最终的矩阵就是：

projectionMatrix=Morth×S×TprojectionMatrix=M_{orth}\times S\times T

# 斜投影

斜投影的特点是：光束平行，但与投影平面不垂直。俯视图为：
11
这样摄像机的观测范围就成了一个平行六面体：
12
它的观测过程实际上就是先做一个剪切变换，再做一个正平行投影（注意是剪切而不是平移！很显然斜投影会导致物体形变）:
13
而这个剪切变换是沿着x,yx,y 轴做的。假设光线和x,yx,y 轴正半轴夹角为θ,ϕ\theta,\phi：
14
则有剪切变换：

H(θ,ϕ)=(10cotθ001cotϕ000101)H(\theta,\phi)=\begin{pmatrix} 1 & 0 & cot\theta & 0 \\ 0 & 1 & cot\phi & 0 \\ 0 & 0 & 1 & 0 \\ & & & 1 \end{pmatrix}

所以斜投影的最终矩阵就是：

projectionMatrix=Morth×S×T×HprojectionMatrix=M_{orth}\times S\times T\times H

# 透视投影

这部分我们只讨论斜投影的情况，因为正投影实际上就是去掉一个剪切变换的HH 矩阵。

  1. 首先，消失点 (COP)，投影平面 (这里也是 viewing plane) 和观测范围如下图所示：
    15
    可见观测范围是一个梯形体。我们首先也还是作一个剪切变化，把z=−farz=-far 和z=−nearz=-near 两矩形的中心剪切到坐标轴上，形式化地：
    16
    我们要把点(left+right2,top+bottom2,−near)(\frac{left+right}{2},\frac{top+bottom}{2},-near) 剪切到(0,0,−near)(0,0,-near)。
    (其中为什么是−near-near 是因为zz 轴其实是垂直屏幕向外，而我们希望观测范围在屏幕内部，所以若near,far>0near,far>0 则−near,−far-near,-far 都在zz 轴负半轴的话就在屏幕内侧了。)
    这个剪切变换的矩阵为：

H(cot−1left+right−2∗near,cot−1top+bottom−2∗near)=(10left+right−2∗near001top+bottom−2∗near000101)H(cot^{-1}\frac{left+right}{-2*near},cot^{-1}\frac{top+bottom}{-2*near})=\begin{pmatrix} 1 & 0 & \frac{left+right}{-2*near} & 0 \\ 0 & 1 & \frac{top+bottom}{-2*near} & 0 \\ 0 & 0 & 1 & 0 \\ & & & 1 \end{pmatrix}

  1. 下一步，我们用缩放变换（沿x,yx,y 轴），把这个梯形体规整成一个每个切面都是一个正方形的梯形体，对任意zz，使得x,yx,y 的范围都是[z,−z][z,-z]（注意是zz 轴负半轴，故z<0z<0）。这个缩放变换的矩阵为：

S(−2∗nearright−left,−2∗neartop−bottom,1)S(\frac{-2*near}{right-left},\frac{-2*near}{top-bottom},1)

  1. 然后最重要的一步来了，我们要有一个 Perspective normalization 的矩阵，它的作用是把规整后的梯形体，仿射变换成一个边长为 2 的正方形：
    17
    为了做到这一点，我们需要扩展 齐次坐标 ，使得其第四维可以不只是 0 或 1，扩展方式为：

(xyzw)=(x/wy/wz/w1)\begin{pmatrix} x \\ y \\ z \\ w \end{pmatrix}=\begin{pmatrix} x/w\\ y/w\\ z/w\\ 1 \end{pmatrix}

即ww 不为 0 时，都可以表示点，且一个点此时就有无穷多个不同的坐标。这样做有什么好处呢？我们考虑如下形式的矩阵，注意它的第四行已经不全为 0，所以作用的是扩展后的齐次坐标：

N=(1000010000αβ00−10)N=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \alpha & \beta\\ 0 & 0 & -1 & 0 \end{pmatrix}

我们看它作用于某个点时：

N(xyz1)=(xyαz+β−z)=(−x/z−y/z−(α+β/z)1)N\begin{pmatrix} x\\ y\\ z\\ 1 \end{pmatrix}=\begin{pmatrix} x\\ y\\ \alpha z+\beta\\ -z \end{pmatrix}= \begin{pmatrix} -x/z\\ -y/z\\ -(\alpha+\beta/z)\\ 1 \end{pmatrix}

为什么我们要用−1-1 乘出来是−z-z 呢？之前提到，因为整个观测范围都在屏幕内测，故都是zz 的负半轴，故zz 值都是负的。所以取−z-z 相当于取了个zz 的绝对值。
然后显然−x/z,−y/z-x/z,-y/z 的范围都是[−1,1][-1,1] 了就，我们希望−(α+β/z)-(\alpha+\beta/z) 的范围也是[−1,1][-1,1]，那么就可以取：

α=−near+farnear−farβ=−2∗near∗farnear−far\alpha=-\frac{near+far}{near-far}\\ \beta=-\frac{2*near*far}{near-far}

最终得到了 Perspective normalization 的矩阵：

N=(1000010000near+farnear−far2∗near∗farnear−far00−10)N=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \frac{near+far}{near-far} & \frac{2*near*far}{near-far}\\ 0 & 0 & -1 & 0 \end{pmatrix}

  1. 最后，我们得到了个边长为 2，中心在原点的正方体。然后我们做一个正平行投影MorthM_{orth} 即可。

故最终的投影矩阵为：

projectionMatrix=Morth×N×S×HprojectionMatrix=M_{orth}\times N\times S\times H