3DGS前向渲染(1) 算法原理

高斯基元属性

对于 SfM 获得的点云，用三维的高斯函数（kernel）将它们膨胀成带有位置、形状、颜色和不透明度的椭球。这个所谓的膨胀，其实是用 KNN 取最近邻三个点的平均半径作为初始化高斯椭球的主轴长度，并以点的颜色作为椭球的初始颜色。

这些参数都是可学习的，会在后续的流程中以深度学习的方法进行优化，且为了保证这些参数的取值符合自身的物理意义（主轴长度为正、旋转矩阵为正交矩阵、协方差矩阵为半正定对称矩阵、不透明度在 0-1 之间），会将这些参数通过所谓的激活函数之后再使用。高斯椭球的参数在程序中组织如下，其中 xyz 为位置、scaling 和 roration 为控制形状的主轴长度和旋转四元数、features dc 和 features rest 分别为球谐系数的直流分量（0 次球谐系数）和剩余分量（其它阶次球谐系数）、opacity 为不透明度。

之所以选择高斯椭球作为场景的表达基元而不是立方体或其它形式，是因为高斯函数具有良好的数学性质，一是高斯函数对仿射变换封闭（渲染中需要用到仿射变换），二是高斯函数可微，且三维高斯沿某一个轴积分成二维后仍为高斯函数（渲染实质上就是投影到二维平面），三是高斯椭球本身是各向异性的体素，比各向同性的高斯球表现力更强。

位置和形状

这个所谓的椭球形状，是由高斯函数的数学表达式决定的（高维高斯均假设各向异性，即不同方向梯度不同）

\begin{aligned} G(x) &= \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x - \mu)^2}{2\sigma^2} \right)\\ G(\mathbf{x}) &= \frac{1}{2\pi \sqrt{|\mathbf{\Sigma}|}} \exp\left( -\frac{1}{2} (\mathbf{x} - \mathbf{\mu})^\top \mathbf{\Sigma}^{-1} (\mathbf{x} - \mathbf{\mu}) \right),\\ &\text{其中 } \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix},\quad \mathbf{\mu} = \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix},\quad \mathbf{\Sigma} = \begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{bmatrix}\\ G(\mathbf{x}) &= \frac{1}{(2\pi)^{3/2} \sqrt{|\mathbf{\Sigma}|}} \exp\left( -\frac{1}{2} (\mathbf{x} - \mathbf{\mu})^\top \mathbf{\Sigma}^{-1} (\mathbf{x} - \mathbf{\mu}) \right),\\ &\text{其中 } \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix},\quad \mathbf{\mu} = \begin{bmatrix} \mu_1 \\ \mu_2 \\ \mu_3 \end{bmatrix},\quad \mathbf{\Sigma} = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{bmatrix}\\ \mathbf{\mu}为均值&，决定分布的位置\\ \mathbf{\Sigma}为协方&差矩阵，为半正定的对称矩阵，对角线元素为方差，其余元素为协方差，决定分布的形状 \end{aligned}

先以二维为例，令二维高斯函数为某一常数，提取各高次项系数并整理之后可以得到一个标准的椭圆方程；这一点从几何意义上理解更直观：想象一个平面与钟形曲面相交，截出的曲线就是椭圆。三维高斯函数同理，令其为常数并提取各高次项系数，整理之后就是一个标准的椭球面方程。在几何意义上为了可视化出来，可以想象用颜色而非欧式距离来代表函数值，这样，三维高斯函数就是一个实心椭球，从内到外就是颜色渐变的椭球壳。

高斯椭球的位置和形状分别由三维高斯分布中的均值和协方差矩阵确定，均值确定位置很好理解，下面解释协方差如何确定形状。

任意高斯分布可由标准高斯分布通过仿射变换得到，对高斯函数进行仿射变换满足如下规律（推导过程可见 EWA Volume Splatting，简单来说就是写出高斯函数的表达式，应用仿射变换代入并加以整理便可得到这个规律）

\begin{aligned} \mathbf{x} &\sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})\\ \mathbf{y} &= \mathbf{A} \mathbf{x} + \mathbf{b}\\ \mathbf{y} &\sim \mathcal{N}(\mathbf{A} \mathbf{\mu} + \mathbf{b},\ \mathbf{A} \mathbf{\Sigma} \mathbf{A}^\top)\\ \end{aligned}

由于仿射变换本身可以表达为旋转、缩放和平移（ $\mathbf{A} \mathbf{x} + \mathbf{b} = \mathbf{R}\mathbf{S} \mathbf{x} + \mathbf{b}$ ）的组合。标准正态分布的协方差矩阵 $\mathbf{\Sigma}=\mathbf{I}$ ，因此变换之后的协方差矩阵 $\mathbf{\Sigma’}$ 可以写成下式，即用旋转矩阵和缩放矩阵构成协方差矩阵。

\mathbf{\Sigma’} = \mathbf{A}\mathbf{I}\mathbf{A}^\top = \mathbf{R}\mathbf{S}\mathbf{S}^\top\mathbf{R}^\top

之所以要这样分解，是为了在参数学习的过程中自然地满足协方差矩阵为半正定的对称矩阵这一约束（旋转矩阵为正交矩阵，满足转置等于逆的性质，而缩放矩阵为对角矩阵，两者合成为 $\mathbf{R}~ \text{diag}(\sigma_1^2 \, \sigma_2^2 \, \sigma_3^2) ~\mathbf{R}^\top$ ，满足半正定对称矩阵的要求）。旋转矩阵控制了椭球绕三个坐标轴的旋转角度，缩放矩阵控制了椭球三个主轴的长度，因此它们构成的协方差矩阵决定了椭球的形状。注：这一点如果从代数上理解，其实就是协方差矩阵的二次型对应的几何图像，即椭圆。协方差矩阵中的方差项和协方差项就分别控制着椭圆的长轴长度和绕坐标轴的旋转角度。

颜色

高斯椭球的颜色由球谐函数确定，高斯泼溅中以此实现不同视角下观察到不同颜色的效果。球谐函数是球坐标系下，定义在单位球面的一组具有正交完备性和旋转不变性的基函数。

之所以使用球坐标系函数，是为了将颜色与观测角度建立联系。而基函数是函数空间中的概念，基函数在其中的作用与向量空间中基向量的作用相同，即基函数可以通过线性组合张成函数空间中的任意函数（其实函数空间就是一种抽象的向量空间，同样满足向量空间中诸如对加法和数乘封闭等一系列性质）。球谐函数是球坐标系下的基函数，类似于笛卡尔坐标系下傅里叶级数中的 $\cos(nx)$ 和 $\sin(nx)$ 的概念，任何球面坐标系的函数都可以展开为无穷个球谐函数，换句话说，多个球谐函数的线性组合，可以近似表示任意的球坐标函数。

球谐函数本身，是球坐标系下的拉普拉斯方程在分离变量 $r$ 后，角度部分通解的正交项。

\begin{aligned} 球坐标拉普拉斯方程:&\nabla^2 f = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial f}{\partial r} \right) + \frac{1}{r^2 \sin\theta} \frac{\partial}{\partial \theta} \left( \sin\theta \frac{\partial f}{\partial \theta} \right) + \frac{1}{r^2 \sin^2\theta} \frac{\partial^2 f}{\partial \varphi^2}\\ 分离变量后角度部分的方程:& \frac{1}{\sin\theta} \frac{\partial}{\partial\theta} \left( \sin\theta \frac{\partial Y}{\partial\theta} \right) + \frac{1}{\sin^2\theta} \frac{\partial^2 Y}{\partial\varphi^2} + l(l + 1)Y = 0\quad (球函数方程)\\ \end{aligned}

球函数方程可以进一步分离变量，最终得到两个常微分方程，分别求解再回代可以得到球函数方程的解如下。其中关于方位角 $\theta$ 的方程又称为连带勒让德方程，方程的解称为连带勒让德函数 $P_l^m (x)$ ，全名为 $l$ 次 $m$ 阶连带勒让德函数（ $\text{degree} = l, \text{order} = m$ ），在计算机中以递推的方式实现。

\begin{aligned} 球函数方程:&Y(\theta, \varphi) = \Theta(\theta) \Phi(\varphi) \\ 解得:&\Phi(\varphi) = \text{e}^{\text{i} m\varphi} ,\, m = 0, \pm 1, \cdots, \pm l &\quad(极角方程的解)\\ &\Theta(\theta) = P_l^m (\cos\theta) ,\, m = 0, \pm 1, \cdots, \pm l &\quad(方位角方程的解)\\ 其中&P_l^m(x) = (-1)^m (1 - x^2)^{\frac{m}{2}} \frac{\text{d}^m}{\text{d}x^m} \left(P_l(x)\right)&\quad(递推公式)\\ &P_l(x) = \frac{1}{2^l \cdot l!} \frac{\text{d}^l}{\text{d}x^l} \left[(x^2 - 1)^l\right]&\quad(通项公式)\\ 球函数方程的解为:&Y(\theta, \varphi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} P_l^m (\cos\theta) \text{e}^{\text{i} m\varphi}, \, m = 0, \pm 1, \pm 2, \cdots \\ \end{aligned}

高斯泼溅中使用的是归一化的球谐函数的分量，归一化即在球函数方程解的基础上增加了归一化系数，分量是指在 $m>0$ 时取 $\text{e}$ 指数的实数部分， $m<0$ 时取 $\text{e}$ 指数的虚数部分，最终得到高斯泼溅中球谐函数的表达式如下。

\begin{aligned} Y_l^m (\theta, \varphi) &= K_l^m Y_{lm} (\theta, \varphi)\\ &=\begin{cases} \sqrt{2} K_l^m \cos(m\varphi) P_l^m(\cos\theta) & m > 0 \\ \sqrt{2} K_l^m \sin(-m\varphi) P_l^{-m}(\cos\theta) & m < 0 \\ K_l^0 P_l^0(\cos\theta) & m = 0 \end{cases}\\ 其中K_l^m &= \sqrt{ \frac{2l + 1}{4\pi} \frac{(l - |m|)!}{(l + |m|)!} }\\ \end{aligned}

解释完了球谐函数本身，再聊聊正交完备性和旋转不变性。正交完备性包括正交归一性和完备性两个性质。函数空间的正交归一性定义为不同函数在定义域上的内积为 0，函数自身的内积为 1。完备性是指函数展开为基函数表示时，当次数趋于无穷时，级数和平均收敛于原函数，这保证了展开结果会趋近于被展开的函数。而旋转不变性是指，当原函数发生旋转时，只需对生成原函数的线性组合系数（也称为广义傅里叶系数）进行相应的变换，就可以保证变换后的系数能等价还原出新函数，而无需对新函数进行展开并计算广义傅里叶系数（生成系数的计算开销很大）。这具体表现为，当光源发生旋转时，只要同步计算出变换后的广义傅里叶系数，就能保证画面的光照效果不会抖动跳变。

对于给定的球谐函数次数，球谐函数和对应的球谐系数的数量是确定的，用于表示颜色的原函数可以按如下式子展开(以高斯泼溅中支持的最高次数三次为例)，其中球谐函数线性组合的系数就称为球谐系数。由原函数生成线性组合系数的过程称为投影，而用系数和基函数组合成原函数这个逆过程称为重建。

\begin{aligned} f &= \sum_{l} \sum_{m=-l}^{l} \mathbf{c}_{l}^{m} y_{l}^{m}(\theta, \phi) \\ &= \mathbf{c}_{0}^{0} y_{0}^{0}+ \\ &\quad \mathbf{c}_{1}^{-1} y_{1}^{-1}+\mathbf{c}_{1}^{0} y_{1}^{0}+\mathbf{c}_{1}^{1} y_{1}^{1}+ \\ &\quad \mathbf{c}_{2}^{-2} y_{2}^{-2}+\mathbf{c}_{2}^{-1} y_{2}^{-1}+\mathbf{c}_{2}^{0} y_{2}^{0}+\mathbf{c}_{2}^{1} y_{2}^{1}+\mathbf{c}_{2}^{2} y_{2}^{2}+ \\ &\quad \mathbf{c}_{3}^{-3} y_{3}^{-3}+\mathbf{c}_{3}^{-2} y_{3}^{-2}+\mathbf{c}_{3}^{-1} y_{3}^{-1}+\mathbf{c}_{3}^{0} y_{3}^{0}+\mathbf{c}_{3}^{1} y_{3}^{1}+\mathbf{c}_{3}^{2} y_{3}^{2}+\mathbf{c}_{3}^{3} y_{3}^{3} \end{aligned}

与函数的傅里叶级数类似，球谐函数的次数 $l$ 越高，实现的颜色表达就越细致。在三阶的情况下，按 RGB 三通道计算，完整的球谐系数存储为一个 $3\times 16$ 的矩阵。

高斯渲染管线

高斯椭球是三维场景结构的表示，需要一定的步骤将其转化到二维图像平面。从具体操作上来说，其实是将高斯椭球的所有参数转化到二维平面，再离散到像素上。位置和协方差使用不同的处理管线，颜色和不透明度不涉及维度的概念，因此也不需要所谓的转化，颜色只需按给定的观察方位和高斯的球谐系数分别计算颜色的 rgb 分量即可，不透明度在这一步无需处理。注：不透明度用于后续合成图像时的 $\alpha$ -blending，为了叙述方便，将这部分内容划分到后面的“并行渲染管线”章节介绍，故不在这里涉及不透明度的处理。

位置处理管线

对于位置参数，采用世界坐标系—>相机坐标系—>归一化设备坐标系—>像素坐标系的流程。

由于高斯椭球中的位置表示在 world frame（世界坐标系）下，而图像是针对于某个具体的相机视角的，因此需要先将位置的表示方式从 world frame 转化到 camera frame 下，这一步转换称为观测变换，简称 W2C（world to camera）。

转换是坐标系之间的变换，实质是平移和旋转的合成，而为了将两个操作合为一个矩阵，将空间位置用齐次向量 $[x, y, z, 1]^\top$ 表示。根据从 COLMAP 中得到的表示平移的向量 $\text{tvec} = [t_x, t_y, t_z]^\top$ ，和表示旋转的四元数 $\text{qvec} = \{q_0, q_1, q_2, q_3\}$ 得到的旋转矩阵 $\mathbf{R}=\begin{bmatrix} r_{xx} & r_{xy} & r_{xz} \\ r_{yx} & r_{yy} & r_{yz} \\ r_{zx} & r_{zy} & r_{zz} \end{bmatrix}$ (四元数转旋转矩阵的过程有固定的计算步骤，此处省略；另外，高斯泼溅中所有的四元数转旋转矩阵的计算均使用右手系下的公式)，可以按如下方式构造 W2C 的变换矩阵 $\text{view matrix}$ （完整的 W2C 矩阵推导过程可见参考资料，此处只给出由 tvec 和 qvec 构造矩阵的过程）。不难发现，W2C 变换矩阵其实就是相机的外参矩阵

\text{view matrix} = \begin{bmatrix} r_{xx} & r_{xy} & r_{xz} & t_x \\ r_{yx} & r_{yy} & r_{yz} & t_y \\ r_{zx} & r_{zy} & r_{zz} & t_z \\ 0 & 0 & 0 & 1 \end{bmatrix}\\

不妨假设位置在世界坐标系下表示为 $[x_w, y_w, z_w, 1]^\top$ ，在相机坐标系下表示为 $[x_c, y_c, z_c, 1]^\top$ ，则通过以下坐标变换可以得到相机坐标系下的位置表示

\begin{bmatrix} x_c\\ y_c\\ z_c\\ 1 \end{bmatrix} = \begin{bmatrix} r_{xx} & r_{xy} & r_{xz} & t_x \\ r_{yx} & r_{yy} & r_{yz} & t_y \\ r_{zx} & r_{zy} & r_{zz} & t_z \\ 0 & 0 & 0 & 1 \end{bmatrix}\\ \begin{bmatrix} x_w\\ y_w\\ z_w\\ 1 \end{bmatrix}

值得注意的是，从 COLMAP 中得到的 qvec 和 tvec 是 W2C 的，即 qvec 表示世界系的基底在相机系中的姿态（且已经做了归一化），tvec 表示世界系的原点在相机系中的位置。以坐标变换来理解，设世界系为旧基，世界系下坐标为 $\mathbf{x}=(x_1, x_2, \cdots, x_n)^\top$ ；相机系为新基，相机系下坐标为 $\mathbf{y}=(y_1, y_2, \cdots, y_n)^\top$ ；过渡矩阵 $\mathbf{P}$ 为用世界系表示相机系的变换矩阵，那么 qvec 对应的旋转矩阵和 tvec 对应的平移向量构成的变换矩阵是用相机系表示世界系的变换矩阵，即 $\mathbf{P}^{-1}$ 。根据坐标变换公式 $\mathbf{y} = \mathbf{P}^{-1} \mathbf{x}$ ，因此 $\mathbf{P}^{-1}$ 记为 W2C 矩阵，即世界坐标系转换到相机坐标系的变换矩阵；而如果想获得相机在世界坐标系下的位置，则应该从过渡矩阵 $\mathbf{P}$ 取。另外，由于 CUDA 中矩阵的内存布局为列主序；而 torch 和 numpy 中矩阵的布局均为行主序，因此需要传入 CUDA 部分参与运算的 view matrix 和 projection matrix 在代码中均进行了转置。

转化到相机坐标系之后，还需要投影到二维平面，这一步骤称为投影变换，具体来说是透视投影变换。投影变换的具体操作，是用内参构造投影矩阵（注意不是使用内参矩阵），再用这个矩阵将坐标从相机坐标系变换到 NDC（Normalized Device Coordinate，归一化坐标系）；另外，此处的投影变换指的是计算机图形学范畴的概念，是将相机坐标系转化到 NDC 的步骤；后文关于协方差处理中的投影变换，则是指计算机视觉中的概念，即将相机坐标系转化到图像平面上的步骤。所谓的 NDC，在不同的图形学 API 中定义不同，OpenGL 和 glm 中定义为 $[-1, 1]^3$ 的右手系立方体空间，DirectX 中定义为 $[-1, 1]^2, z \in [0, 1]$ 的左手系长方体空间，不同的 NDC 定义，推导出的投影矩阵是有差别的。高斯泼溅中使用 DirectX 的定义，将 znear 平面映射到 z=0 平面，zfar 平面映射到 z=1 平面。

在这样的规定下，结合由视锥范围计算出的 $t、b、l、r$ ，记近平面和远平面的 z 值分别为 n 和 f，推导出投影矩阵 $\text{projection matrix}$ 如下。推导过程可见参考资料，此处省略；另外，投影矩阵的推导方式有两种，一种是按照 GAMES 101 中的方式先将视锥非线性压缩为长方体再正交投影到 NDC，另一种是直接推导投影矩阵；两种方式的推导结果在[0,2]和[1, 2]两个元素上相差一个负号，相当于将变换之后的图像旋转 180°，对高斯椭球的参数学习没有影响；此处列出的为前一种推导方式的结果)

\text{projection matrix} = \begin{bmatrix} \displaystyle \frac{2n}{\,r-l\,} & 0 & \displaystyle -\frac{r+l}{\,r-l\,} & 0 \\[8pt] 0 & \displaystyle \frac{2n}{\,t-b\,} & \displaystyle -\frac{t+b}{\,t-b\,} & 0 \\[8pt] 0 & 0 & \displaystyle \frac{f}{\,n-f\,} & \displaystyle -\frac{nf}{\,n-f\,} \\[8pt] 0 & 0 & 1 & 0 \end{bmatrix}

注：内参矩阵是计算机视觉的概念，直接将三维空间转化到像素平面，用于只关心真实的像素投影的情况；而计算机图形学中一般会先转化到 NDC 坐标系，主要是为了兼容不同的屏幕分辨率，降低渲染管线的耦合度。后续图形硬件对空间进行裁剪做可见性测试以及硬件栅格化都可以针对这个统一的空间进行，最后再单独映射到屏幕分辨率，避免了与分辨率产生耦合。高斯泼溅中使用了这套渲染管线，因此使用投影矩阵而不是内参矩阵。

不妨假设位置在相机坐标系下表示为 $[x_w, y_w, z_w, 1]^\top$ ，在 NDC 坐标系下表示为 $[x_n, y_n, z_n, w_n]^\top$ ，则通过以下坐标变换可以得到相机坐标系下的位置表示

\begin{bmatrix} x_n\\ y_n\\ z_n\\ w_n \end{bmatrix} = \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} & -\frac{nf}{n-f} \\ 0 & 0 & 1 & 0 \end{bmatrix}\\ \begin{bmatrix} x_c\\ y_c\\ z_c\\ 1 \end{bmatrix}

对于齐次坐标，所有分量同时除以同一个数得到的结果与原坐标等价；为了保证 w 分量保持为 1，对结果同时除以 $w_n$ ，即 $[x_n, y_n, z_n, w_n]^\top = [\frac{x_n}{w_n}, \frac{y_n}{w_n}, \frac{z_n}{w_n}, 1]^\top$ ，这一步称为透视除法。

由于后续只关心屏幕空间二维投影，故只取结果的 x 和 y 分量，相当于转化到图像坐标系。此时已经将三维椭球转化到了 $[-1, 1]^2$ 的平面范围内，而图像尺寸是 $[w, h]$ ，因此还需要拉伸回这个尺寸；此外，图像坐标系的原点在图像正中心，目标的像素坐标系在图像的左上角，因此相差一个平移。这两步合起来将 $[-1, 1]^2$ 的平面坐标系转化到像素坐标系，称为视口变换。不妨假设位置在 $[-1, 1]^2$ 平面范围内表示为 $(x_n, y_n)$ ，在像素坐标系下表示为 $(x_p, y_p)$ ，则视口变化可以表示为

\begin{aligned} x_p = \frac{w}{2}\cdot x_n + \frac{w}{2}\\ y_p = \frac{h}{2}\cdot y_n + \frac{h}{2}\\ \end{aligned}

但考虑到下一步的像素离散操作，需要从不考虑面积的点转化到有面积的像素上，所以这一步在代码实现中其实将 NDC 中的点映射到了像素中心而不是边缘，避免了边界处偏移半像素的问题。对于索引为 $(i, j)$ 的像素，像素中心为 $(i+0.5, j+0.5)$ ，因此在上述转化的基础上只需偏移 0.5 个单位即可。另外需要注意的是，像素索引范围等于像素坐标或者说像素数量 - 1。以一维的四个像素为例说明这一点，浅蓝色代表像素索引，浅绿色代表像素坐标，同时也是像素范围。

高斯泼溅的代码中使用向左偏移 0.5 个单位的写法，将 $(-1, 1)$ 映射到 $(-0.5, \text{width}-0.5)$ ，而不是 $(0.5, \text{width}+0.5)$ 。映射到 $(-0.5, \text{width}-0.5)$ 之后再转为 int 截去小数部分(截取小数部分的操作对正数相当于向下取整，对负数则是向上取整)，这样正好使得索引控制在 $[0, \text{width}-1] \times [0, \text{height}-1]$ 。反之如果映射到 $(0.5, \text{width}+0.5)$ ，就需要编写额外的逻辑来去掉多出来的像素。

协方差处理管线

对于形状参数，采用世界坐标系—>相机坐标系—>图像平面坐标系的流程，具体来说是将协方差矩阵进行层层传递的过程。

首先根据旋转四元数转化得到的旋转矩阵 $\mathbf{R}$ 和三个主轴的缩放 $\sigma_i$ 构造世界坐标系下的三维协方差矩阵 $\mathbf{V}_k'' = \mathbf{R}~ \cdot\text{diag}(\sigma_1^2 \, \sigma_2^2 \, \sigma_3^2)\cdot ~\mathbf{R}^\top$

世界坐标系到相机坐标系的转化为观测变换，本质上是平移+旋转构成仿射变换。高斯协方差矩阵在仿射变换下的规律如下(推导过程见 EWA Volume Splatting，此处直接给出结论)，将位置 $\mathbf{t}_k$ 处、形状为协方差矩阵 $\mathbf{V}_k''$ 的高斯 kernel 记为

r_k''(t) = G_{\mathbf{V}_k''}(\mathbf{t}-\mathbf{t}_k) = \frac{1}{2|\mathbf{V}_k''|^{-\frac{1}{2}}}\text{exp}\left( \left(\mathbf{t}-\mathbf{t}_k\right)^\top \mathbf{V}_k''^{-1} \left(\mathbf{t}-\mathbf{t}_k\right) \right)

设仿射变换为 $\mathbf{u} = \varphi(\mathbf{t}) = \mathbf{W}\mathbf{t} + \mathbf{d}$ ，变换前后的高斯 kernel 满足如下等式，可见变换之后协方差矩阵传递为 $\mathbf{W}\mathbf{V}_k'\mathbf{W}^\top$

G_{\mathbf{V}_k''}(\varphi^{-1}(\mathbf{u})-\mathbf{t}_k) = \frac{1}{|\mathbf{W}^{-1}|}G_{\mathbf{V}_k'}(\mathbf{u}-\mathbf{u}_k) = r_k'(u), 其中\mathbf{V}_k' = \mathbf{W}\mathbf{V}_k''\mathbf{W}^\top\\

为了处理非线性的投影变换，高斯泼溅中采用了 EWA Volume 中的方法，先将相机空间变换到投影平面为图像平面的光线空间，再用局部仿射变换来近似投影变换。为了理解推导过程，需要先介绍这个所谓的光线空间。光线空间是用于参数化所有可能光线的非线性空间，以下的介绍针对 EWA Volume Splatting 中的光线空间，与计算机图形学或计算机视觉中一般意义上的光线空间有区别。给定投影中心和投影平面，从投影中心发出光线穿过投影平面，令光线空间中，点的坐标为 $\mathbf{x} = [x_0, x_1, x_2]^\top$ 。这三个坐标分量的含义如下：穿过投影中心和投影平面上坐标为 $(x_0, x_1)$ 的点的光线称为观测光线； $x_2$ 表示从投影中心出发，沿观测光线方向、移动 $x_2$ 的欧式距离到达该点。光线空间可以看作视锥的参数化描述；光线由像素位置决定，点的位置由光线方向和距离决定。

EWA Volume Splatting 指出，通过将相机空间变换到这样的光线空间，可以将非线性的投影变换在局部近似为仿射变换这样的线性变换(这里的线性变换是指变换将直线映射为直线，不是线性代数意义上的线性变换)。在相机空间下，直接进行透视投影需要除以深度，引入了非线性，雅可比矩阵对深度变化敏感，不适合用局部线性近似的方式来处理；转化到光线空间后，由于深度信息在光线空间的方向向量表示下被归一化掉了，深度在光线空间表现为线性缩放，可以采用局部线性近似处理。之所以是局部近似，是因为近似的实质是泰勒一阶展开，只对光线穿过的点是精确的，而对于整个投影的对象，并不能在全局采用这种近似方法。这种局部一阶近似提供了一种高效执行透视投影的方式，可实现与 OpenGL 中投影变换相同的功能。

高斯泼溅中将投影平面设为图像平面，该平面到投影中心的距离为焦距 $f$ ，由于要将空间坐标映射到图像平面，故需要用水平像素尺寸下的表示 $f_x$ 和竖直像素尺寸下的表示 $f_y$ 。设空间点的在相机空间下的坐标表示为 $\mathbf{u} = [u_0, u_1, u_2]^\top$ ，在光线空间下的坐标表示为 $\mathbf{x} = [x_0, x_1, x_2]^\top$ ，两者映射关系记为 $\mathbf{x} = m(\mathbf{u})$ ，令 $l = \vert\vert(u_0, u_1, u_2)^\top\vert\vert$

\begin{aligned} \begin{bmatrix} x_0\\ x_1\\ x_2 \end{bmatrix} = m(\mathbf{u}) = \begin{bmatrix} f_x \frac{u_0}{u_2}\\ f_y \frac{u_1}{u_2}\\ l \end{bmatrix} \end{aligned}

通过这样的映射，即可将相机空间投影到图像平面，代替了直接进行透视投影的操作；而这一步映射本身也是非线性的，用局部仿射变换来近似表示。局部仿射变换定义为 $m(\mathbf{u})$ 在 $\mathbf{u}_k$ 处的一阶泰勒展开，根据雅可比矩阵的定义，计算向量对向量的导数可以求出雅可比矩阵

\begin{aligned} \mathbf{x} = m_{\mathbf{u}_k}(\mathbf{u}) &= \mathbf{x}_k + \mathbf{J}_{\mathbf{u}_k}(\mathbf{u} - \mathbf{u}_k), 其中 \mathbf{J}_{\mathbf{u}_k} = \left. \frac{\partial m(\mathbf{u})}{\partial \mathbf{u}} \right|_{\mathbf{u} = \mathbf{u}_k}, \mathbf{x}_k = m(\mathbf{u}_k)\\ \mathbf{J}_{\mathbf{u}_k} &= \begin{bmatrix} \displaystyle \frac{\partial m_1(\mathbf{u})}{\partial u_1} & \displaystyle \frac{\partial m_1(\mathbf{u})}{\partial u_2} & \displaystyle \frac{\partial m_1(\mathbf{u})}{\partial u_3} \\[8pt] \displaystyle \frac{\partial m_2(\mathbf{u})}{\partial u_1} & \displaystyle \frac{\partial m_2(\mathbf{u})}{\partial u_2} & \displaystyle \frac{\partial m_2(\mathbf{u})}{\partial u_3} \\[8pt] \displaystyle \frac{\partial m_3(\mathbf{u})}{\partial u_1} & \displaystyle \frac{\partial m_3(\mathbf{u})}{\partial u_2} & \displaystyle \frac{\partial m_3(\mathbf{u})}{\partial u_3} \end{bmatrix} = \begin{bmatrix} \displaystyle \frac{f_x}{u_2} & \displaystyle 0 & \displaystyle -\frac{f_x u_0}{u_2^2} \\[8pt] \displaystyle 0 & \displaystyle \frac{f_y}{u_2} & \displaystyle -\frac{f_y u_1}{u_2^2} \\[8pt] \displaystyle \frac{u_0}{l} & \displaystyle \frac{u_1}{l} & \displaystyle \frac{u_2}{l} \\[8pt] \end{bmatrix} \end{aligned}

根据前面的高斯协方差矩阵在仿射变换下的规律，应用局部仿射变换前后的高斯 kernel 满足如下等式，变换之后协方差矩阵传递为 $\mathbf{J}\mathbf{W}\mathbf{V}_k''\mathbf{W}^\top\mathbf{J}^\top$

G_{\mathbf{V}_k'}(m^{-1}(\mathbf{x})-\mathbf{u}_k) = \frac{1}{|\mathbf{J}^{-1}|}G_{\mathbf{V}_k}(\mathbf{x}-\mathbf{x}_k) = r_k(x), 其中\mathbf{V}_k = \mathbf{J}\mathbf{V}_k'\mathbf{J}^\top = \mathbf{J}\mathbf{W}\mathbf{V}_k''\mathbf{W}^\top\mathbf{J}^\top\\

由于后续处理只关心二维高斯椭圆的协方差矩阵，因此代码实现时会将雅可比矩阵的第三行置零，最终结果矩阵左上角 $2\times 2$ 的分块矩阵作为二维协方差，记为 $\mathbf{\Sigma}$ ，至此便得到了高斯椭球渲染到图像平面的 footprint。后续计算中会使用 $\mathbf{\Sigma}^{-1}$ ，因此在代码中 preprocess 部分会进行求逆，并将结果记为为 $\text{conic}$ 。