3D高斯椭球投影到2D：从CUDA代码到视觉效果的完整解析

3D高斯椭球投影到2D从CUDA代码到视觉效果的完整解析最近在复现和优化一些前沿的神经渲染项目时我反复与一个核心模块打交道如何将三维空间中的高斯椭球高效、准确地“拍扁”到二维图像平面上。这听起来像是一个纯粹的数学变换问题但当你真正深入到代码层面尤其是在CUDA这种并行计算环境中实现时会发现其中充满了工程上的权衡与精妙的设计。这篇文章我想从一个实践者的角度抛开教科书式的推导聊聊从3D高斯椭球的数学表示到forward.cu中那一行行决定最终渲染效果的CUDA代码这中间究竟发生了什么。无论你是正在啃相关论文的研究者还是希望将3D高斯溅射3D Gaussian Splatting技术落地的工程师相信这些结合了代码细节与视觉原理的拆解能带来一些不一样的启发。1. 重新理解3D高斯不止于一个公式提到3D高斯分布我们脑海中首先浮现的可能是那个经典的多元高斯概率密度函数。但在计算机图形学和视觉的语境下我们更倾向于将其视作一个具有空间特性的渲染基元。一个3D高斯椭球由几个核心属性定义中心位置均值 μ、决定其形状和方向的协方差矩阵 Σ以及一个与颜色、透明度相关的外观属性如球谐系数。这个椭球在三维空间中定义了某种“影响范围”离中心越远其影响力或密度按指数衰减。1.1 协方差矩阵的几何直观旋转与缩放为什么是协方差矩阵因为它以一种紧凑的形式同时编码了椭球的朝向、大小和拉伸程度。一个纯对角线的协方差矩阵对应一个轴对齐的椭球。但现实中的物体表面千姿百态我们需要能任意旋转的椭球。这里一个关键的几何事实是任何椭球都可以通过对一个标准球体单位球进行线性变换得到。这个线性变换通常分解为旋转R和缩放S。想象一下你先拿着一个完美的单位球用缩放矩阵S在X, Y, Z三个轴上把它拉长或压扁变成一个轴对齐的椭球。然后再用一个旋转矩阵R对这个轴对齐的椭球进行旋转使其朝向任意方向。这个组合变换矩阵M R * S的协方差矩阵正是我们需要的 Σ不考虑平移平移由均值μ处理。数学上Σ Mᵀ M。这就是为什么在代码中我们总是看到先构造缩放和旋转矩阵然后相乘并转置相乘来得到协方差矩阵。看看computeCov3D这个设备函数的核心片段已做简化示意// 假设 scale 是缩放向量rot 是四元数表示的旋转 glm::mat3 S glm::mat3(1.0f); S[0][0] scale.x; S[1][1] scale.y; S[2][2] scale.z; // 从四元数 rot 计算旋转矩阵 R (代码略) glm::mat3 R quaternionToMatrix(rot); glm::mat3 M R * S; // 注意顺序通常是先缩放后旋转 glm::mat3 Sigma glm::transpose(M) * M; // 这就是3D协方差矩阵注意缩放与旋转矩阵的乘法顺序R * S还是S * R取决于你对坐标系和变换顺序的定义需要与框架的其他部分保持一致。上述代码是一种常见理解方式。1.2 存储优化利用对称性由于协方差矩阵Σ是一个实对称矩阵Σ Σᵀ它只有6个独立元素3x3矩阵的上三角部分。因此在GPU内存极其宝贵的环境下我们绝不会傻傻地存储9个float。在代码中可以看到计算出的Sigma矩阵只将其上三角的6个元素存入数组cov3D[0] Sigma[0][0]; // (0,0) cov3D[1] Sigma[0][1]; // (0,1) cov3D[2] Sigma[0][2]; // (0,2) cov3D[3] Sigma[1][1]; // (1,1) cov3D[4] Sigma[1][2]; // (1,2) cov3D[5] Sigma[2][2]; // (2,2)这种存储方式在后续的所有变换中都必须被牢记因为我们需要在需要时从这6个数重建出完整的3x3矩阵。2. 投影的核心挑战非线性透视与局部线性近似将3D椭球投影到2D图像平面是整个流程中最具技巧性的部分。理想情况下一个3D高斯分布在经过完整的透视投影变换后在图像平面上应该形成一个2D高斯分布。但这里存在一个根本性矛盾透视投影本身是非线性的物体离相机越远成像越小而高斯分布的协方差变换要求变换是线性的才能用矩阵乘法表示。2.1 透视投影的“扭曲”效应让我们回顾一下针孔相机模型。一个世界点P (X, Y, Z)投影到图像坐标(u, v)的公式是u f_x * (X / Z) c_x v f_y * (Y / Z) c_y这里f_x, f_y是焦距c_x, c_y是光心。关键点在于(X/Z)和(Y/Z)这个除法操作它引入了非线性。这意味着即使一个3D高斯椭球在空间中是“均匀”的经过投影后它在图像上的形状也会因为其自身深度Z值的变化而发生复杂的扭曲而不仅仅是简单的缩放和旋转。2.2 雅可比矩阵局部线性化的钥匙解决这个问题的经典方法来自计算机图形学中的可微表面溅射思想即使用局部线性近似。我们不对整个椭球进行非线性变换而是只考虑椭球中心点P_c处的变换行为。在高斯分布的有效支撑集即其影响显著的区域内这个区域相对于整个场景通常很小。因此我们可以用中心点处的一阶泰勒展开即雅可比矩阵J来近似整个局部区域的变换。雅可比矩阵J描述了从相机坐标系下的3D坐标(X_c, Y_c, Z_c)到2D归一化设备坐标(u, v)的映射在点P_c处的局部线性近似。对于上述透视投影公式其雅可比矩阵是一个2x3的矩阵因为输出是2维的uv输入是3维的XYZ。在代码computeCov2D函数中我们能看到它的计算// t 是高斯中心在相机坐标系下的坐标 (t.x, t.y, t.z) glm::mat3 J glm::mat3( focal_x / t.z, 0.0f, -(focal_x * t.x) / (t.z * t.z), 0.0f, focal_y / t.z, -(focal_y * t.y) / (t.z * t.z), 0, 0, 0 // 第三行通常为0因为深度Z不直接影响颜色/alpha );这个矩阵的物理意义非常直观J[0][0] focal_x / t.z: X坐标变化对u坐标的影响与焦距正比与深度反比。J[0][2] -(focal_x * t.x) / (t.z * t.z): Z坐标深度变化对u坐标的影响。如果点不在光轴上t.x ! 0深度变化会导致图像位置横向移动这就是透视效应。第二行同理对应v坐标。有了这个局部线性近似J我们就可以将3D协方差矩阵Σ_3d在相机坐标系下变换到2D图像空间的协方差矩阵Σ_2dΣ_2d J * Σ_3d * Jᵀ这公式保证了变换后的2D分布仍然是一个高斯分布其形状由原3D椭球在图像平面上的“投影椭圆”决定。3. 深入forward.cuCUDA实现的关键步骤与优化理论清晰后我们进入实战环节剖析forward.cu或类似前向投影核函数中的关键实现。这里的目标是并行处理成千上万个高斯椭球为后续的光栅化阶段准备2D参数。3.1 完整的投影流水线一个完整的投影过程在GPU核函数中通常按以下步骤进行每个线程负责一个或一组高斯数据加载从全局内存读取单个高斯的属性位置、旋转四元数、缩放、不透明度、球谐系数等。视锥体裁剪快速判断该高斯是否在相机视锥体内。如果完全不可见则提前终止该线程的处理节省算力。代码中常用tan_fovx和tan_fovy来计算边界。const float limx 1.3f * tan_fovx; const float limy 1.3f * tan_fovy; const float txtz t.x / t.z; const float tytz t.y / t.z; // 将坐标限制在视锥体范围内避免极端值导致数值不稳定 t.x min(limx, max(-limx, txtz)) * t.z; t.y min(limy, max(-limy, tytz)) * t.z;计算3D协方差矩阵调用类似computeCov3D的函数利用缩放和旋转参数计算出世界坐标系下的3D协方差矩阵 Σ_world。坐标系变换将世界坐标系下的3D协方差矩阵 Σ_world 变换到相机坐标系下 Σ_camera。这需要用到视图矩阵viewmatrix的旋转部分因为协方差矩阵是二次型平移不影响形状。设视图矩阵的旋转部分为W则Σ_camera W * Σ_world * Wᵀ。计算投影雅可比矩阵基于高斯中心在相机空间的位置t计算雅可比矩阵J。计算2D协方差矩阵应用公式Σ_2d J * Σ_camera * Jᵀ得到图像空间归一化设备坐标NDC下的2x2协方差矩阵。添加低通滤波与稳定性处理这是一个非常重要的工程trick。理论上投影出的椭圆可能非常小小于一个像素这会导致光栅化时出现锯齿或数值不稳定。因此通常会给2D协方差矩阵的对角线元素加上一个小的常数如0.3相当于强制每个高斯在图像上至少有一个像素的“铺开”范围。cov[0][0] 0.3f; // 增加x方向方差 cov[1][1] 0.3f; // 增加y方向方差准备光栅化数据将计算好的2D椭圆参数中心uv坐标、2x2协方差矩阵、深度t.z、颜色、不透明度等写入到全局内存中的指定数组供后续的排序和光栅化核函数使用。3.2 CUDA编程的特定考量在CUDA中实现上述流程有几点需要特别注意内存访问模式高斯属性在全局内存中的存储应尽量满足合并访问。例如将所有高斯的位置(x,y,z)连续存储旋转四元数(qw,qx,qy,qz)连续存储这样线程束Warp内的线程可以一次性读取连续的内存块效率最高。计算精度与速度在computeCov2D这类函数中涉及大量矩阵乘法和除法。虽然使用double精度更高但考虑到渲染的容错性和性能float通常是更优的选择。像glm这样的数学库在设备端提供了优化过的单精度矩阵运算。线程束内的分支视锥体裁剪会导致线程分支。虽然现代GPU对线程束内发散的分支处理能力有所提升但最好还是尽量减少。一种策略是进行两遍处理第一遍轻量级筛选标记出可见高斯第二遍再对可见高斯进行详细计算。与光栅化的衔接前向投影核函数的输出通常是几个大的数组分别存储所有可见高斯的2D中心、协方差、颜色、深度等。这些数组会被传递给一个基于瓦片Tile-Based的光栅化核函数该函数将屏幕划分为小块每个线程块处理一个瓦片只加载与该瓦片相关的高斯进行计算极大减少不必要的内存访问和计算。4. 从2D椭圆到像素颜色光栅化的最后一环投影得到了2D高斯椭圆的参数但如何高效地将其转换为屏幕上每个像素的颜色这正是3D高斯溅射技术得名的原因——它不像传统三角形光栅化那样处理边而是将每个高斯视为一个带有透明度的“色块”通过溅射的方式贡献到像素上。4.1 可微的Alpha混合对于屏幕上的一个像素p其最终颜色C_p是所有按深度排序后覆盖该像素的高斯椭圆贡献的加权和C_p Σ (c_i * α_i * ∏ (1 - α_j)) i ji这里c_i是第i个高斯根据其球谐系数和视角计算出的颜色。α_i是第i个高斯在该像素点处的不透明度贡献它由2D高斯分布函数的值、高斯自身的不透明度参数以及像素点到椭圆中心的距离共同决定。∏ (1 - α_j)是前面所有高斯带来的透明度累积。这个过程必须是可微的因为整个3D高斯溅射管道通常用于神经渲染的训练中需要通过梯度下降来优化高斯的所有属性位置、颜色、形状等。4.2 高效的范围查询与近似直接为每个像素计算所有高斯的影响是不现实的。实践中采用两种主要优化边界框与层次结构每个2D高斯椭圆都有一个轴对齐的边界框AABB。在光栅化前先计算每个椭圆的屏幕空间AABB。在瓦片化光栅化中可以快速判断哪些椭圆与当前屏幕瓦片相交只加载相交的高斯进行计算。近似评估严格计算2D高斯分布在每个像素点的精确值开销很大。通常采用近似方法例如将椭圆离散化为几个采样点或者使用预计算的查找表。更常见的是利用2D协方差矩阵的逆即精度矩阵来快速计算马氏距离作为评估高斯衰减的依据。下面是一个简化的伪代码展示在像素着色器中如何计算单个高斯的贡献// 输入像素坐标p高斯中心mu_2d2D协方差矩阵Sigma_2d的逆Sigma_inv vec2 delta p - mu_2d; // 计算马氏距离的平方二次型 float mahalanobis_dist2 dot(delta, Sigma_inv * delta); // 使用高斯核函数计算权值exp(-0.5 * dist2) 是标准高斯形式 float alpha_weight opacity * exp(-0.5 * mahalanobis_dist2); // 将权值限制在0-1之间并乘以前序透明度 float effective_alpha 1.0 - exp(-alpha_weight); // 近似保证可微 vec3 color_contrib color * effective_alpha; // 然后进行标准的从前到后的Alpha混合提示在实际实现中exp函数比较昂贵有时会用更简单的衰减函数来近似只要保证函数是光滑可微的即可。4.3 性能与质量的权衡表在实现光栅化时我们常常面临一系列选择下表对比了不同策略的考量策略优点缺点适用场景逐像素精确计算质量最高理论最准确计算量巨大性能瓶颈离线渲染最终帧输出预计算衰减纹理运行时速度快只需一次纹理查找占用额外内存精度受纹理分辨率限制移动端或实时性要求高的场景基于瓦片的并行光栅化充分利用GPU并行性内存访问高效实现复杂度高负载可能不均衡现代GPU上的标准实践用于训练和实时预览低通滤波加常数避免像素级闪烁增强稳定性会使非常小的高斯变模糊损失细节几乎必用是工程上的稳定器在我自己的项目中初期为了快速验证我采用了最简单的逐像素计算结果训练速度慢得令人绝望。后来切换到基于瓦片的实现并加入了严谨的视锥体裁剪和边界框测试性能提升了数十倍。那个给2D协方差矩阵对角线加0.3的操作看似不起眼却实实在在地解决了训练初期因高斯初始化位置随机而导致的画面闪烁问题。