[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

技术文档

文章目录

Preliminary
- 坐标转换的流程
- 坐标转换在3DGS中的代码实现
- render函数接口调用逻辑
Forward
- 总体流程
- preprocessCUDA
- - 视椎体过滤[in_frustum]
  - 3D协方差矩阵计算[**computeCov3D]**
  - 2D协方差计算[**computeCov2D**]
  - 2D高斯投影计算
  - 计算当前高斯在屏幕空间的覆盖范围
  - SH参数转换RGB值
  - 输出结果
- InclusiveSum和duplicateWithKeys
- render
Backward
Reference

这篇文章作为之前3DGS论文解析文章的补充，之前只是从原理和流程上进行了解析，没有深入到内部的数学原理中，这篇文章会结合其他几篇经典的3DGS论文，综合来推导数学原理和代码首先。主线还是沿着render kernel的处理流程来梳理。

💡因为中间涉及到大量数学公式，文章可能存在笔误，欢迎评论区指正，我会及时修正。在此感谢成文过程中的朋友们的讨论和帮助

Preliminary

首先回顾一下之前的内容，在3DGS算法中我们的任务是把三维高斯渲染成一帧图像，这是一个3D到2D的转换，这其中涉及到了几个坐标变换，从物理世界的3D点到图像屏幕的2D坐标系，这个变换是第一部分内容，第二部分分析代码转换的实现，第三部分介绍render的一些代码调用入口

坐标转换的流程

在之前的文章中我们将所有的Gaussian的位置定义在世界坐标系中，从这个坐标系出发，还需要几轮变换才能完成光栅化的渲染，下面的过程主要参考了OpenGL Transformation文章，坐标变换过程如下：
[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

VertexData中的数据是输入信息中的XYZ位置，ObjectCoordinate表示整个重建场景的世界坐标系

ObjectCoordinates下的点经过ModelViewMatrix的变换转换为相机坐标系，ModelViewMatrix表示Camera到World系的转换，代码变量为viewpoint_camera.world_view_transform；

相机系的点经过ProjectionMatrix转换到ClipCoordinates，ProjectionMatrix可以通过FOV和图像长宽重新计算；将ClipCoordinates下的坐标除 $w$ 可以得到NDC(Normalized Device Coordinates)坐标系。整体的投影矩阵为（这里的推导主要是相似三角形和斜截式直线方程的计算，具体推导过程参考https://www.songho.ca/opengl/gl_projectionmatrix.html）；

视椎体到NDC的映射关系为，视椎体 $x$ 坐标的范围从 $[l ， r]$ 到 $[- 1, 1]$ ， $y$ 坐标从 $[b ， t]$ 到 $[- 1, 1]$ ， $z$ 坐标从 $[- n ， - f]$ 到 $[- 1, 1]$ 。

$\\begin{pmatrix}\\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{-2fn}{f-n}\\\\ 0 & 0 & -1 & 0 \\end{pmatrix}\\tag 1$

如果NDC的 $z$ 映射范围设为 $[0, 1]$ ，那么上面的投影公式变为

$\\mathbf{P}=\\begin{pmatrix}\\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}{f-n}&\\frac{-fn}{f-n}\\\\0&0&-1&0\\end{pmatrix} \\tag2$

由于NDC的初始定义为左手系，如果修正为我们熟悉的右手系右下前RDF，那么投影矩阵为：

$\\mathbf{P}_\\text{RDF} = \\begin{pmatrix}1 & 0 & 0 & 0\\\\ 0 & -1 & 0 & 0\\\\ 0 & 0 & 1 & 0\\\\ 0 & 0 & 0 & 1 \\end{pmatrix}\\mathbf{P}\\begin{pmatrix}1 & 0 & 0 & 0\\\\ 0 & -1 & 0 & 0\\\\ 0 & 0 & -1 & 0\\\\ 0 & 0 & 0 & 1 \\end{pmatrix} = \\begin{pmatrix}\\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}{f-n}&\\frac{-fn}{f-n}\\\\0&0&1&0\\end{pmatrix} \\tag3$

但是大家平时用的相机都是自带内参的，在有内参的情况下公式中的每个变量应该如何计算？假设给定相机内参 $f_x,f_y,c_x,c_y$ 、图像的尺寸 $W, H$ 和视椎体的最近和最远的平面坐标 $n, f$ ，计算在Eye坐标系下 $t b r l$ 的表达式(注意正负)

$\\left\\{\\begin{matrix} t = & n*c_{y}/f_{y}\\\\ b = & -n*(H-c_{y})/f_{y} \\\\ r = & n*(W-c_{x})/f_{x} \\\\ l = & -n*c_{x}/f_{x} \\end{matrix}\\right. \\tag4$

带入公式3可知

$\\mathbf{P}_\\text{RDF} = \\left(\\begin{array}{cccc} \\frac{2f_{x}}{W} & 0 & -1+\\frac{2c_{x}}{W} & 0\\\\ 0 & \\frac{2f_{y}}{H} & -1+\\frac{2c_{y}}{H} & 0\\\\ 0 & 0 & \\frac{f}{f-n} & -\\frac{fn}{f-n}\\\\ 0 & 0 & 1 & 0 \\end{array}\\right) \\tag5$

假设Eye相机系下齐次坐标点 $[X, Y, Z, 1]$ ，经过 $\\mathbf{P}_\\text{RDF}$ 投影矩阵变换，除 $w_{clip}=Z$ ，得到NDC下坐标为

$\\left\\{ \\begin{aligned}x_{ndc} & =\\frac{2f_{x}X}{WZ}-1+\\frac{2c_{x}}{W}\\\\ y_{ndc} & =\\frac{2f_{y}Y}{HZ}-1+\\frac{2c_{y}}{H}\\\\ z_{ndc} & =\\frac{f(Z-n)}{Z(f-n)} \\end{aligned} \\right. \\tag{6}$

最后经过ViewportTransform，转到WindowsCoordinates，得到图像坐标系的 $u ， v$

[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

经过推导可得NDC到WindowCoordinate转换公式如下，公式中 $\\mathbf{W}=W,\\mathbf{H}=H$ ，即为图像分辨率， $\\mathbf{F}=1， \\mathbf{N}=0$ ， $\\mathbf{X}$ 轴起始坐标 $\\mathbf{X}=0$ ， $\\mathbf{Y}$ 轴起始坐标 $\\mathbf{Y}=0$ ，这是我们刚才设定好的，

$\\left\\{\\begin{aligned}x_w &= \\frac{\\mathbf{W}}{2}x_n+(\\mathbf{X}+\\frac{\\mathbf{W}}{2}) \\\\y_w &= \\frac{\\mathbf{H}}{2}y_n+(\\mathbf{Y}+\\frac{\\mathbf{H}}{2}) \\\\z_w &= \\frac{\\mathbf{F}-\\mathbf{N}}{2}z_n+\\frac{\\mathbf{F}+\\mathbf{N}}{2}\\end{aligned}\\right. \\tag{7}$

将公式6的值带入公式7，得到我们平时常用的内参投影矩阵，可以认为Projection Matrix+Viewport Transform就是相机内参矩阵

$\\left\\{ \\begin{aligned}u & =\\frac{(x_{ndc}+1)W}{2}=\\frac{f_{x}X}{Z}+c_{x}\\\\ v & =\\frac{(y_{ndc}+1)H}{2}=\\frac{f_{y}Y}{Z}+c_{y} \\end{aligned} \\right. \\tag{8}$

坐标转换在3DGS中的代码实现

在Projection Matrix之前的步骤，都是正常的三维坐标系变换，Projection Matrix的计算我们可以参考StreetGaussian中的实现，getProjectionMatrix是3DGS的实现，getProjectionMatrixK是加入相机内参之后的计算方式。

render函数接口调用逻辑

GaussianRasterizationSettings：输入配置项，包括宽高，视场角，C2W矩阵，内参矩阵等参数

GaussianRasterizer：光栅器初始化

渲染步骤是GaussianRasterizer的forward函数，从这里开始就要进入cuda的kernel部分，具体跳转要参考ext.cpp文件中Python和Cuda的函数映射关系；

CUDA的程序在diff-gaussian-rasterization仓库中；

rasterize_points.cu文件中的RasterizeGaussiansCUDA函数是Python调用的入口；在输入一系列参数之后，调用rasterizer_impl.cu文件中的CudaRasterizer::Rasterizer::**forward()**函数；

Forward

首先明确每个Gaussian的属性，我们从总体流程入手；然后详细解释Cuda代码中的主要函数的作用

Attribute Variable Name N Params Valid Range Activation Mean Position xyz 3 (-inf, inf) none Color rgb 3 [0, 1] none or sigmoid Opacity opacity 1 [0, 1] sigmoid Scale scale 3 (0, inf) exponential Spherical Harmonic Coefficients sh 0, 9, 24, 45 none

总体流程

CudaRasterizer::Rasterizer::forward()流程如下：

从FOV反算焦距；
构建了几何状态GeometryState，就是将变量对应到显存，这里是一串指针内存地址的操作，就是提前开辟好空间；
根据长宽构建tile网格，后面的计算都是基于tile进行的；为什么要设计Tile进行计算，很多文章中也提到了，就是为了并行加速，3DGS中设计的BLOCK为16。
构建图像状态ImageState，映射到显存；
前处理函数FORWARD::preprocess，映射到preprocessCUDA函数；这一步功能包括对所有Gaussian进行视锥剔除，NDC投影，协方差计算，半径估计，通过球谐函数计算RGB，将所有结果将更新到GeometryState中；
通过InclusiveSum计算每个tile受到影响的高斯数量累积值；
构建binning状态BinningState，同时构建[Tile序号|Depth]对作为key
duplicateWithKeys函数用来构造键值对，value值是Gaussian的ID
cub::DeviceRadixSort::SortPairs函数对刚才的键值对排序，根据key值，先根据Tile，Tile相等时根据depth排序；
identifyTileRanges函数确定每个Tile的在排序之后的list中的范围，也就是每个Tile读取哪些Gaussian
最后调用render函数，映射为forward.cu中的renderCUDA 函数

[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

图片是来自于FlashGS论文，可以很直观的表示3DGS的渲染计算过程

并行处理流程优化也有人做过尝试，在BalanceGS中让BLOCK和Tile尺寸适配数据集的图像尺寸可以提高处理速度，同时像素层面并行也可以改为高斯层面的并行。如果想了解更底层GPU并行可以参考这篇论文

preprocessCUDA

视椎体过滤[in_frustum]

正如之前的推导，这里用到了viewmatrix和projmatrix 两个变量。viewmatrix是world2eye的转换，将Gaussian转到相机系；projmatrix 是全量投影矩阵，viewmatrix @ projection ，projection是我们之前推导的 $\\mathbf{P}_{\\text{RDF}}$

但是在实际代码中，只用了viewmatrix 在Eye Coordinate下进行判断；（从注释代码的情况看，似乎有打算在NDC空间下判断）

__forceinline__ __device__ bool in_frustum(int idx,const float* orig_points,const float* viewmatrix,const float* projmatrix,bool prefiltered,float3& p_view){float3 p_orig = { orig_points[3 * idx], orig_points[3 * idx + 1], orig_points[3 * idx + 2] };// Bring points to screen spacefloat4 p_hom = transformPoint4x4(p_orig, projmatrix);float p_w = 1.0f / (p_hom.w + 0.0000001f);float3 p_proj = { p_hom.x * p_w, p_hom.y * p_w, p_hom.z * p_w };p_view = transformPoint4x3(p_orig, viewmatrix);if (p_view.z <= 0.2f)// || ((p_proj.x  1.3 || p_proj.y  1.3))){if (prefiltered){printf(\"Point is filtered although prefiltered is set. This shouldn\'t happen!\");__trap();}return false;}return true;}

3D协方差矩阵计算[computeCov3D]

计算公式为公式9，因为协方差矩阵是一个半正定矩阵，所以通过保存cov3d矩阵上三角6个元素值，来保证对称性；但是对称矩阵并不代表是半正定矩阵，根据判定条件如果实对称矩阵A的所有特征值都大于0，可以判断A是一个半正定矩阵。因此在论文中作者通过旋转矩阵R和尺度值S构建一个协方差矩阵；（具体推导可以看Ref1的Optimizing 3D Gaussians部分）

$\\mathbf{{\\Sigma}_3}= \\mathbf{RS}\\mathbf{S}^T\\mathbf{R}^T = \\mathbf{RS}(\\mathbf{SR})^T \\tag{9}$

在4D-Gaussian中，其实这个协方差变成了一个4维方阵， $\\mathbf{S}$ 矩阵可以在对角线增加时间 $t$ 变量， $\\mathbf{R}$ 矩阵通过两个四元数进行双边乘法，可以得到最终的结果（此处有点复杂，可以参考论文和Ref5的乘法）。通过在 $\\mathbf{S}$ 和 $\\mathbf{R}$ 矩阵引入时间 $t$ 变量，可以推导条件高斯矩阵；

假设有4D高斯分布

$\\begin{bmatrix} \\mathbf{x} \\\\ \\mathbf{y} \\end{bmatrix} \\sim \\mathcal{N} \\left( \\mu = \\begin{bmatrix} \\boldsymbol{\\mu}_x \\\\ \\boldsymbol{\\mu}_y \\end{bmatrix}, \\Sigma =\\begin{bmatrix} \\Sigma_{xx} & \\Sigma_{xy} \\\\ \\Sigma_{yx} & \\Sigma_{yy} \\end{bmatrix} \\right) \\tag{10}$

我们希望在给定时间 $t_0$ 条件下，计算空间上的条件分布；

多元高斯边缘分布的定理，已知多元高斯分布为

$\\mathbf{x} \\mid \\mathbf{y} = \\mathbf{y}_0 \\sim \\mathcal{N}(\\boldsymbol{\\mu}_{x|y}, \\Sigma_{x|y}) \\tag{11}$

该条件分布的条件均值为：

$\\boldsymbol{\\mu}_{x|y} = \\boldsymbol{\\mu}_x + \\Sigma_{xy} \\Sigma_{yy}^{-1} (\\mathbf{y}_0 - \\boldsymbol{\\mu}_y) \\tag{12}$

条件协方差为

$\\Sigma_{x|y} = \\Sigma_{xx} - \\Sigma_{xy} \\Sigma_{yy}^{-1} \\Sigma_{yx}$

借助这个定理，可以计算出在 $y=t_0$ 时刻的条件高斯分布，代码实现参考4DGaussianSplatting仓库（作者数学功力有限，就不放详细推导过程了，如果有兴趣可以自己推导一下，利用到二次型配方法和Schur公式来完成~）

2D协方差计算[computeCov2D]

矩阵偏导数需要借助雅克比矩阵进行，2D协方差的计算公式为13，其中 $\\mathbf{J}$ 为雅克比矩阵， $\\mathbf{W}$ 表示viewmatrix的旋转部分也就是之前推导中提到的ProjectionMatrix；

$\\mathbf{{\\Sigma}}_2=\\mathbf{J}\\mathbf{W}\\mathbf{\\Sigma}_3\\mathbf{W}^T\\mathbf{J}^T \\tag{13}$

根据EWA Splating论文，雅克比矩阵 $\\mathbf{J}$ 是像素坐标对Gaussian中心点的偏导数，如果是pinhole投影，计算公式如下：

$\\mathbf{J} = \\begin{bmatrix} f_x / Z_c & 0 & -f_x X_c/Z_c^2 \\\\\\ 0 & f_y/Z_c & -f_y Y_c/Z_c^2\\end{bmatrix} \\tag{14}$

如果是Fisheye投影，可以参考FisheyeGS，推导步骤如下：

常见鱼眼相机是等距Equidistant投影，3D点 $\\mathbf{P} = (x, y, z)$ 到像素点 $(u, v)$ 的计算方式为

$\\frac{x}{r} \\cdot \\theta,\\quad v = \\frac{y}{r} \\cdot \\theta \\tag{15}$

其中向量长度 $\\sqrt{x^2 + y^2}$ ，光轴夹角 $\\theta = \\arctan\\left( \\frac{r}{z} \\right)$

我们要计算的雅克比矩阵

$\\mathbf{J} = \\frac{d(u, v)}{d(x, y, z)} = \\begin{bmatrix}\\frac{\\partial u}{\\partial x} & \\frac{\\partial u}{\\partial y} & \\frac{\\partial u}{\\partial z} \\\\\\frac{\\partial v}{\\partial x} & \\frac{\\partial v}{\\partial y} & \\frac{\\partial v}{\\partial z}\\end{bmatrix} \\tag{16}$

那么 $(u, v)$ 对 $(x, y, z)$ 的偏导数分别为：

$\\begin{aligned}\\frac{\\partial u}{\\partial x} & = \\frac{\\partial}{\\partial x} \\left( \\frac{x}{r} \\cdot \\theta \\right) = \\theta \\cdot \\frac{\\partial}{\\partial x} \\left( \\frac{x}{r} \\right) + \\frac{x}{r} \\cdot \\frac{\\partial \\theta}{\\partial x} \\\\\\frac{\\partial u}{\\partial y} & = -\\frac{x y}{r^3} \\cdot \\theta + \\frac{x}{r} \\cdot \\frac{\\partial \\theta}{\\partial y} \\\\\\frac{\\partial u}{\\partial z} & = \\frac{x}{r} \\cdot \\frac{\\partial \\theta}{\\partial z}\\end{aligned}$

中间变量 $\\theta$ 的偏导数

$\\begin{aligned}\\frac{\\partial \\theta}{\\partial x} & = \\frac{1}{1 + (r/z)^2} \\cdot \\frac{1}{z} \\cdot \\frac{x}{r} = \\frac{x z}{r (r^2 + z^2)} \\\\\\frac{\\partial \\theta}{\\partial y} &= \\frac{y z}{r (r^2 + z^2)}\\\\\\frac{\\partial \\theta}{\\partial z} &= -\\frac{r}{r^2 + z^2}\\end{aligned}$

以及 $\\frac{\\partial}{\\partial x} \\left( \\frac{x}{r} \\right) = \\frac{r^2 - x^2}{r^3}$ 。

结合内参的焦距 $f_x, f_y$ ，我们得到如下的公式，为了方便书写和对应代码，其中 $\\frac{z}{r^2 (r^2 + z^2)}$ ， $\\frac{\\theta}{r^3}$ ，具体代码实现可以参考FisheyeGS的实现

$\\mathbf{J}_{\\text{fisheye}} =\\begin{bmatrix}f_x (x^2 a + y^2 b) & f_x xy (a - b) & -f_x x / (r^2 + z^2) \\\\f_y xy (a - b) & f_y (y^2 a + x^2 b) & -f_y y / (r^2 + z^2)\\end{bmatrix} \\tag{17}$

2D高斯投影计算

从3DGS投影到2D平面是公式18，这个公式在最后的渲染阶段会计算，现在我们要做的就是提前计算这个公式中用到的相关变量2D协方差 $\\Sigma_{2}$ ，Gaussian的2D中心点 $\\boldsymbol{\\mu}$

$\\exp\\left( - \\frac{1}{2}(\\mathbf{x} - \\boldsymbol{\\mu})^T \\Sigma_2^{-1}(\\mathbf{x} - \\boldsymbol{\\mu})\\right),\\\\ \\mathbf{x} = [u,v]^T,\\,\\boldsymbol{\\mu} = [\\bar{u},\\bar{v}]^T \\tag{18}$

其中 $\\mathbf{{\\Sigma}}_2^{-1}$ 的计算可以通过如下公式计算，但是2Dcov是半正定矩阵所以 $b = c$ ，所以代码中简化为3个数字。

$\\Sigma_{2}=\\begin{bmatrix}a & b\\\\ c & d \\end{bmatrix},\\quad\\Sigma_{2}^{-1}=\\frac{1}{ad-bc}\\begin{bmatrix}d & -b\\\\ -c & a \\end{bmatrix} \\tag{19}$

如果是采用Mip-Splating的2DFilter的去毛刺方案，这里要在最外层加一个系数，系数的计算公式在Mip-Splatting论文解析文章中有写，这里不再赘述。因为2D高斯投影到图像上面还是一个椭圆，因此使用一个窗函数对保留最主要的部分，剔除概率密度较低的区域达到去除毛刺的效果；

计算当前高斯在屏幕空间的覆盖范围

代码注释中解释这里是通过计算cov2D的特征值来实现，具体计算过程就是通过对角线元素计算矩阵的迹，通过欧矩阵的迹得到特征值，两个特征值代表椭圆的长轴和短轴； $\\det$ 是协方差矩阵的行列式。（具体推导计算过程可以看参考文献1的Tile Based Rasterization部分，推导过程也是比较常用的矩阵性质）

$\\lambda_{1,2}=\\frac{\\text{Tr}}{2}\\pm \\sqrt{ \\left(\\frac{\\text{Tr}}{2}\\right)^2 - \\det } \\tag{20}$

三维高斯在二维的投影是一个椭圆，这里通过 $3\\sigma$ 来计算最大覆盖半径，得到最终影响范围也就是和tile是否相交，在代码中为了简化采用了正方形；

利用公式8将NDC坐标投影到像素系，得到当前高斯在屏幕上的中心位置；最后利用getRect获取这个高斯圆在tile网络上面的包围盒（如果没有，则跳过）；

// Compute extent in screen space (by finding eigenvalues of// 2D covariance matrix). Use extent to compute a bounding rectangle// of screen-space tiles that this Gaussian overlaps with. Quit if// rectangle covers 0 tiles. float mid = 0.5f * (cov.x + cov.z);float lambda1 = mid + sqrt(max(0.1f, mid * mid - det));float lambda2 = mid - sqrt(max(0.1f, mid * mid - det));float my_radius = ceil(3.f * sqrt(max(lambda1, lambda2)));//坐标投影float2 point_image = { ndc2Pix(p_proj.x, W), ndc2Pix(p_proj.y, H) };// 范围判断uint2 rect_min, rect_max;getRect(point_image, my_radius, rect_min, rect_max, grid);if ((rect_max.x - rect_min.x) * (rect_max.y - rect_min.y) == 0)return;

SH参数转换RGB值

普通的3DGS颜色就按照3DGS论文解析文章中提到的球谐函数转换公式计算。

但是我们想引入时间维度怎么做呢？4DGaussianSplatting论文给出了解释，通过时间的傅里叶调制来使得球谐函数随着时间变化；构造和时间相关的傅里叶余弦基函数对SH参数进行调制，公式表达如下； $T$ 表示时间周期， $t$ 表示当前Gaussian时间到当前时间戳的距离（这也是一个优化量）

$C(\\mathbf{d}, t) = \\sum_i a_i Y_i(\\mathbf{d}) + b_i Y_i(\\mathbf{d}) \\cos\\left( \\frac{2\\pi t}{T} \\right) + c_i Y_i(\\mathbf{d}) \\cos\\left( \\frac{4\\pi t}{T} \\right) + \\dots \\tag{21}$

代码实现参考4DGaussianSplatting

输出结果

最终preprocessCUDA这个函数输出了相机系Z作为深度depths，当前Gaussian的在2维平面的最大半径radii，像素坐标位置points_xy_image，二维协方差的逆矩阵和不透明度conic_opacity，当前Gaussian的tile面积tiles_touched 。

InclusiveSum和duplicateWithKeys

结束前处理之后，就是对tile和Gaussian进行分组；
[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

InclusiveSum这个函数是计算前缀和的，根据代码注释可以得到这一点；这里解释具体含义，以代码注释为例tiles_touched = [2, 3, 0, 2, 1] ，经过前缀和函数得到point_offsets = [2, 5, 5, 7, 8] 。它表示：第 0 个高斯影响 tile0和1；第 1 个高斯影响 tile2~4；第 2 个不影响；第 3 个影响 tile5~6；第 4 个影响 tile7。

现在我们已知了每个Gaussian的影响范围，也就是能影响哪些tile，duplicateWithKeys 这个函数作用就是将每个Gaussian影响的tile复制多份，比如Gaussian0影响了tile0和1，那么就会复制为两份。

point_list_unsorted = [0, 0]point_list_keys_unsorted = [tileKey(0), tileKey(1)]

在得到所有tile|depth的键值对之后，再进行排序，可以参考上图中的c。最后在排序后的列表中计算得到每个tile的起点和终点Gaussian索引，用于后续的多tile并行计算；

render

好的，在经过了一系列的前置计算，数据准备之后，终于来到了render光栅化的部分。也就是forward.cu中的renderCUDA 函数

根据3DGS论文中的Eq2和代码实现，修改其形式为

$\\sum_{i=1}^{N} \\mathbf{c}_{i}\\alpha_{i}T_{i}+\\mathbf{c}_{bg}T_{N+1}\\\\\\text{where}\\;\\; T_{i} = \\prod_{j=1}^{i-1}(1 - \\alpha_{j}),\\;\\alpha_{i} = o_i g_{i}(u, v) \\tag{22}$

其中 $\\mathbf{c}_{i}$ 和 $o_i$ 为第 $i^{th}$ 个高斯球的颜色和不透明度， $g_{i}(u, v)$ 为第 $i^{th}$ 个高斯球在像素坐标 $u, v$ 处的概率。 $T_i$ 可以通过 $T_i = T_{i-1}(1-\\alpha_{i-1})$ 迭代计算得到，另外也可以写成下面这种形式：

$T_i = 1 - \\sum_{j=1}^{j-1}\\alpha_jT_j \\tag{23}$

函数首先定位了当前tile，涉及到的像素索引pix_id，需要处理的Gaussian范围range ，收集到当前tile用到的像素坐标，GaussID，2D协方差和不透明度；

然后开始计算自然常数 $e$ 的幂power ，设当前坐标和gauss中心坐标的向量为 $(\\mathbf{x} - \\boldsymbol{\\mu})^T=[d_x , d_y]$

$\\begin{aligned} \\text{power} & = - \\frac{1}{2}(\\mathbf{x} - \\boldsymbol{\\mu})^T \\Sigma_2^{-1}(\\mathbf{x} - \\boldsymbol{\\mu}) \\\\ & = - \\frac{1}{2} \\begin{bmatrix} d_x & d_y \\end{bmatrix} \\begin{bmatrix} a & b\\\\ c & d \\end{bmatrix} \\begin{bmatrix} d_x \\\\ d_y \\end{bmatrix} \\\\ & = -0.5(a\\cdot {d}_x^2+d\\cdot {d}_y^2)-b\\cdot {d}_x \\cdot {d}_y \\end{aligned} \\tag{24}$

带入公式12和公式15，我们可以得到当前Gaussian的不透明度贡献值alpha；由于自然对数的幂函数是一个单调递增函数，所以距离Gaussian中心越远的坐标位置，alpha越小，颜色越浅，乘上其余Gaussian的 $T$ 值，对当前像素颜色贡献就越小；

然后根据公式16计算累积透光率 $T_i$ ，最后alpha和T都会作用在输出结果的每一个通道上（这里不单是颜色结果，还包括深度等其他的信息的处理，代码中逆深度invdepth也是用这种方式进行渲染）；计算出所有颜色之后根据公式15，对每个像素位置的加权上背景颜色，得到了最后的渲染结果；

Backward

一个巨大的数学计算工程Orz，留坑，有时间会开一篇文章写这个求导过程；

Reference

公式推导参考文章：gaussian_splatting/MATH.md at main · joeyan/gaussian_splatting
3DGS：[2308.04079] 3D Gaussian Splatting for Real-Time Radiance Field Rendering
A Survey on 3D Gaussian Splatting：[2401.03890] A Survey on 3D Gaussian Splatting
4D Gaussian Splatting：[2409.04751] Fisheye-GS: Lightweight and Extensible Gaussian Splatting Module for Fisheye Cameras
Fisheye-GS：[2409.04751] Fisheye-GS: Lightweight and Extensible Gaussian Splatting Module for Fisheye Cameras
Quaternion kinematics for the error-state Kalman filter：[2311.16493] Mip-Splatting: Alias-free 3D Gaussian Splatting
Mip-Splatting：[2311.16493] Mip-Splatting: Alias-free 3D Gaussian Splatting
FlashGS：[2408.07967] FlashGS: Efficient 3D Gaussian Splatting for Large-scale and High-resolution Rendering
Balanced 3DGS：[2412.17378] Balanced 3DGS: Gaussian-wise Parallelism Rendering with Fine-Grained Tiling

如果你喜欢我的文章欢迎点赞、关注，同时非常欢迎一起探讨交流，谢谢支持；

[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

文章目录

Preliminary

坐标转换的流程

坐标转换在3DGS中的代码实现

render函数接口调用逻辑

Forward

总体流程

preprocessCUDA

视椎体过滤[in_frustum]

3D协方差矩阵计算[computeCov3D]

2D协方差计算[computeCov2D]

2D高斯投影计算

计算当前高斯在屏幕空间的覆盖范围

SH参数转换RGB值

输出结果

InclusiveSum和duplicateWithKeys

render

Backward

Reference

公告

DeepSeek全套部署资料免费下载

免费可商用字体批量下载

标签

[万字长文]3DGS中forward过程的公式推导及代码分析_3d gs 公式推导

文章目录

Preliminary

坐标转换的流程

坐标转换在3DGS中的代码实现

render函数接口调用逻辑

Forward

总体流程

preprocessCUDA

视椎体过滤[in_frustum]

3D协方差矩阵计算[computeCov3D]

2D协方差计算[computeCov2D]

2D高斯投影计算

计算当前高斯在屏幕空间的覆盖范围

SH参数转换RGB值

输出结果

InclusiveSum和duplicateWithKeys

render

Backward

Reference

相关问题

公告

DeepSeek全套部署资料免费下载

免费可商用字体批量下载

标签