3D Gaussian Splatting 论文与代码笔记

论文: 3D Gaussian Splatting for Real-Time Radiance Field Rendering (SIGGRAPH 2023)

代码仓库: graphdeco-inria/gaussian-splatting

一、3D 高斯的定义（论文 § 4）

每个 3D 高斯本质上是一个带颜色的半透明椭球，由以下参数完全描述：

$s$ 和 $q$ 共同决定协方差矩阵 $\Sigma$，详见下节。

二、协方差矩阵的构建（论文 § 4, Eq.6）

2.1 为什么不直接优化 $\Sigma$ ？

An obvious approach would be to directly optimize the covariance matrix $\Sigma$ to obtain 3D Gaussians that represent the radiance field. However, covariance matrices have physical meaning only when they are positive semi-definite. For our optimization of all our parameters, we use gradient descent that cannot be easily constrained to produce such valid matrices, and update steps and gradients can very easily create invalid covariance matrices. —— 论文 § 4

协方差矩阵必须半正定，梯度下降无法保证这一点。因此将 $\Sigma$ 分解为：

这样无论 $R$、$S$ 取什么值，$\Sigma = (RS)(RS)^T$ 天然半正定。

2.2 对应代码

# gaussian_model.py → setup_functions
def build_covariance_from_scaling_rotation(scaling, scaling_modifier, rotation):
    L = build_scaling_rotation(scaling_modifier * scaling, rotation)  # L = R @ S
    actual_covariance = L @ L.transpose(1, 2)   # Sigma = L @ L^T
    symm = strip_symmetric(actual_covariance)    # 取 6 个独立分量（对称矩阵）
    return symm

# general_utils.py
def build_scaling_rotation(s, r):
    L = torch.zeros((s.shape[0], 3, 3), dtype=torch.float, device="cuda")
    R = build_rotation(r)       # 四元数 → 3x3 旋转矩阵
    L[:,0,0] = s[:,0]           # 对角缩放矩阵 S
    L[:,1,1] = s[:,1]
    L[:,2,2] = s[:,2]
    L = R @ L                   # L = R @ S
    return L

三、从点云初始化高斯（论文 § 5）

SfM 产出的稀疏点云 → 初始化为一组 3D 高斯：

# gaussian_model.py → create_from_pcd
def create_from_pcd(self, pcd, cam_infos, spatial_lr_scale):
    # 位置：直接使用 SfM 点云坐标
    fused_point_cloud = torch.tensor(pcd.points).float().cuda()

    # 颜色：RGB → 0 阶球谐系数，高阶全部初始化为 0
    fused_color = RGB2SH(torch.tensor(pcd.colors).float().cuda())
    features = torch.zeros((N, 3, (self.max_sh_degree + 1) ** 2)).float().cuda()
    features[:, :3, 0] = fused_color

    # 缩放：KNN 找最近邻距离，各向同性初始化
    dist2 = distCUDA2(...)
    scales = torch.log(torch.sqrt(dist2))[..., None].repeat(1, 3)

    # 旋转：单位四元数 [1, 0, 0, 0]（无旋转）
    rots = torch.zeros((N, 4), device="cuda")
    rots[:, 0] = 1

    # 不透明度：初始化为较低值 sigmoid^{-1}(0.1)
    opacities = inverse_sigmoid(0.1 * torch.ones((N, 1)))

四、优化器配置（论文 § 5.1）

# gaussian_model.py → training_setup
l = [
    {'params': [self._xyz],          'lr': position_lr * spatial_lr_scale, "name": "xyz"},
    {'params': [self._features_dc],  'lr': feature_lr,                    "name": "f_dc"},
    {'params': [self._features_rest],'lr': feature_lr / 20.0,             "name": "f_rest"},
    {'params': [self._opacity],      'lr': opacity_lr,                    "name": "opacity"},
    {'params': [self._scaling],      'lr': scaling_lr,                    "name": "scaling"},
    {'params': [self._rotation],     'lr': rotation_lr,                   "name": "rotation"}
]

• 位置 lr 会乘以 spatial_lr_scale（场景大小归一化），并随训练指数衰减
• SH 高阶系数 lr = DC 的 1/20，因为高阶更敏感

五、自适应密度控制（论文 § 5.2）

We observe that both [under-reconstruction and over-reconstruction] have large view-space positional gradients. Intuitively, this is likely because they correspond to regions that are not yet well reconstructed, and the optimization tries to move the Gaussians to correct this. —— 论文 § 5.2

核心思想：屏幕空间位置梯度大 → 该区域重建质量差 → 需要增密。

5.1 梯度累积

# gaussian_model.py
def add_densification_stats(self, viewspace_point_tensor, update_filter):
    # 累积屏幕空间 (x, y) 方向的梯度范数（不含深度 z）
    self.xyz_gradient_accum[update_filter] += torch.norm(
        viewspace_point_tensor.grad[update_filter, :2], dim=-1, keepdim=True)
    self.denom[update_filter] += 1

只看 :2（x, y），因为 densification 关心的是 2D 图像平面上哪里还没覆盖好。

5.2 加密与修剪（主入口）

# gaussian_model.py
def densify_and_prune(self, max_grad, min_opacity, extent, max_screen_size, radii):
    grads = self.xyz_gradient_accum / self.denom  # 平均梯度
    grads[grads.isnan()] = 0.0

    self.densify_and_clone(grads, max_grad, extent)  # 克隆
    self.densify_and_split(grads, max_grad, extent)  # 分裂

    # 剪枝：透明度过低 / 屏幕上过大 / 世界空间过大
    prune_mask = (self.get_opacity < min_opacity).squeeze()
    if max_screen_size:
        big_points_vs = self.max_radii2D > max_screen_size
        big_points_ws = self.get_scaling.max(dim=1).values > 0.1 * extent
        prune_mask = prune_mask | big_points_vs | big_points_ws
    self.prune_points(prune_mask)

5.3 克隆（Clone）—— 梯度大 + 尺寸小（under-reconstruction）

小高斯覆盖不够 → 在同一位置复制一份，增加密度。

def densify_and_clone(self, grads, grad_threshold, scene_extent):
    selected = (grads >= grad_threshold) & (max_scaling <= threshold)
    # 直接复制一份完全相同的高斯
    new_xyz = self._xyz[selected]
    ...

5.4 分裂（Split）—— 梯度大 + 尺寸大（over-reconstruction）

大高斯太粗糙 → 拆成 N=2 个更小的，在原高斯分布内随机采样新位置。

def densify_and_split(self, grads, grad_threshold, scene_extent, N=2):
    selected = (grads >= grad_threshold) & (max_scaling > threshold)
    # 在原高斯的分布范围内采样 2 个新位置
    samples = torch.normal(mean=0, std=self.get_scaling[selected])
    new_xyz = rotate(samples) + self.get_xyz[selected]
    new_scaling = original_scaling / (0.8 * N)  # 尺寸缩小
    # 删除原来的大高斯
    self.prune_points(selected)

六、渲染管线（论文 § 4, § 6）

# gaussian_renderer/__init__.py
def render(viewpoint_camera, pc, pipe, bg_color, ...):

    # 1) 创建屏幕空间点，用于捕获 2D 梯度（→ densification）
    screenspace_points = torch.zeros_like(pc.get_xyz, requires_grad=True, device="cuda") + 0
    screenspace_points.retain_grad()

    # 2) 配置光栅化参数
    raster_settings = GaussianRasterizationSettings(
        image_height=..., image_width=...,
        tanfovx=..., tanfovy=...,
        bg=bg_color,
        viewmatrix=viewpoint_camera.world_view_transform,   # W2C 矩阵
        projmatrix=viewpoint_camera.full_proj_transform,     # 完整投影矩阵
        sh_degree=pc.active_sh_degree,
        campos=viewpoint_camera.camera_center,               # SH 需要视线方向
        ...
    )

    # 3) 准备协方差：默认传 scale + rotation，由 CUDA 内部构建 Sigma
    scales = pc.get_scaling
    rotations = pc.get_rotation

    # 4) 准备颜色：默认传 SH 系数，由 CUDA 内部根据视线方向计算 RGB
    shs = pc.get_features

    # 5) 调用 CUDA 可微光栅化器
    rendered_image, radii, depth_image = rasterizer(
        means3D = pc.get_xyz,
        means2D = screenspace_points,
        shs = shs,
        opacities = pc.get_opacity,
        scales = scales,
        rotations = rotations,
        ...
    )

    # 6) 返回
    return {
        "render": rendered_image.clamp(0, 1),   # 渲染图像
        "viewspace_points": screenspace_points,  # 屏幕空间点（带梯度）
        "visibility_filter": (radii > 0).nonzero(),  # 可见高斯
        "radii": radii,                          # 屏幕半径
        "depth": depth_image,                    # 深度图
    }

6.1 CUDA 光栅化器内部流程（论文 § 6）

代码位于 submodules/diff-gaussian-rasterization/cuda_rasterizer/ 下。

N 个 3D 高斯
  → ① 分 16×16 tile
  → ② 视锥剔除 + 3D→2D 投影 + 计算覆盖 tile 数
  → ③ 为每个 高斯×tile 实例生成排序 key = tile_ID | depth
  → ④ GPU Radix Sort 全局排序
  → ⑤ 找出每个 tile 在排序列表中的起止范围
  → ⑥ 逐 tile 前到后 alpha-blending（共享显存 + 饱和退出）
  → 输出 rendered_image + radii + depth

① 分 tile（16×16）

文件: config.h + rasterizer_impl.cu:236-237

// config.h
#define BLOCK_X 16
#define BLOCK_Y 16

// rasterizer_impl.cu:236-237 — 计算 tile 网格大小
dim3 tile_grid((width + BLOCK_X - 1) / BLOCK_X, (height + BLOCK_Y - 1) / BLOCK_Y, 1);
dim3 block(BLOCK_X, BLOCK_Y, 1);
// 例: 1920×1080 → tile_grid = (120, 68) = 8160 个 tile

② 预处理：视锥剔除 + 投影 + 计算覆盖 tile

文件: forward.cu:150-269，preprocessCUDA 内核，每个高斯球一个线程：

// 视锥剔除（第190行）
if (!in_frustum(idx, orig_points, viewmatrix, projmatrix, prefiltered, p_view))
    return;

// 3D → 2D 齐次投影 + 透视除法（第195-197行）
float4 p_hom = transformPoint4x4(p_orig, projmatrix);
float3 p_proj = { p_hom.x * p_w, p_hom.y * p_w, p_hom.z * p_w };

// 构建 3D 协方差 Σ = (RS)(RS)^T（第208行，对应论文 Eq.6）
computeCov3D(scales[idx], scale_modifier, rotations[idx], cov3Ds + idx * 6);

// 投影到 2D 协方差 Σ' = J W Σ W^T J^T（第213行，对应论文 Eq.5）
float3 cov = computeCov2D(p_orig, focal_x, focal_y, tan_fovx, tan_fovy, cov3D, viewmatrix);

// 用 2D 协方差的特征值算屏幕半径（3σ），确定覆盖的 tile 矩形（第237-243行）
float my_radius = ceil(3.f * sqrt(max(lambda1, lambda2)));
getRect(point_image, my_radius, rect_min, rect_max, grid);

// 记录触及的 tile 数（第268行）
tiles_touched[idx] = (rect_max.y - rect_min.y) * (rect_max.x - rect_min.x);

computeCov3D 内部（forward.cu:114-148）还完成了四元数→旋转矩阵→协方差矩阵的构建：

glm::mat3 M = S * R;                          // S 对角缩放 × R 旋转
glm::mat3 Sigma = glm::transpose(M) * M;      // Σ = M^T M，天然半正定

③ 生成排序 key：tile_ID | depth

文件: rasterizer_impl.cu:70-111，duplicateWithKeys 内核。

每个高斯被复制成覆盖 tile 数那么多份，每份生成一个 64-bit key：

// 高 32 位 = tile ID，低 32 位 = 深度（第102-104行）
uint64_t key = y * grid.x + x;      // tile ID
key <<= 32;                          // 左移到高 32 位
key |= *((uint32_t*)&depths[idx]);   // 低 32 位 = 深度的 bit pattern

排序后效果：先按 tile 分组，同一 tile 内按深度从近到远。

④ GPU Radix Sort（基数排序）

文件: rasterizer_impl.cu:306-311，使用 NVIDIA CUB 库：

cub::DeviceRadixSort::SortPairs(
    binningState.list_sorting_space,
    binningState.sorting_size,
    binningState.point_list_keys_unsorted, binningState.point_list_keys,   // key: tile|depth
    binningState.point_list_unsorted, binningState.point_list,             // value: gaussian ID
    num_rendered, 0, 32 + bit);   // 只排必要的 bit 数，跳过高位零（性能优化）

一次全局排序代替逐像素排序——这是论文强调的关键性能优化。

⑤ 确定每个 tile 的高斯范围

文件: rasterizer_impl.cu:116-138，identifyTileRanges 内核。

排序后同一 tile 的高斯连续排列，遍历一遍找出边界：

uint32_t currtile = key >> 32;            // 从 key 提取 tile ID
if (currtile != prevtile) {               // tile 边界
    ranges[prevtile].y = idx;             // 上一个 tile 结束
    ranges[currtile].x = idx;             // 当前 tile 开始
}

⑥ 逐 tile α-blending（前到后）

文件: forward.cu:274-397，renderCUDA 内核——光栅化器的核心。

每个 tile 一个 thread block（256 线程），每个线程负责一个像素：

// 加载该 tile 的高斯范围（第305行）
uint2 range = ranges[block.group_index().y * horizontal_blocks + block.group_index().x];

float T = 1.0f;    // 初始透射率

// 分批协作加载到共享显存（第310-312行）
__shared__ int collected_id[BLOCK_SIZE];
__shared__ float2 collected_xy[BLOCK_SIZE];
__shared__ float4 collected_conic_opacity[BLOCK_SIZE];

// 对每个高斯，计算该像素处的高斯值（第352行，对应论文 Eq.2）
float power = -0.5f * (con_o.x * d.x * d.x + con_o.z * d.y * d.y) - con_o.y * d.x * d.y;
float alpha = min(0.99f, con_o.w * exp(power));

// 前到后 alpha-blending（第371-372行，对应论文 Eq.3）
C[ch] += features[...] * alpha * T;   // C += c_i × α_i × T_i
T = T * (1 - alpha);                  // 更新透射率

// α 饱和 → 该像素停止（第364-367行）
if (test_T < 0.0001f) { done = true; continue; }

// 整个 tile 所有像素都 done → thread block 提前退出（第326-328行）
int num_done = __syncthreads_count(done);
if (num_done == BLOCK_SIZE) break;

最后写出结果，保存 final_T（反向 pass 要用）：

final_T[pix_id] = T;                                    // 残余透射率 → 反向 pass 用
n_contrib[pix_id] = last_contributor;                    // 最后贡献者 ID → 反向 pass 用
out_color[...] = C[ch] + T * bg_color[ch];              // 残余透射率 × 背景色

反向 pass

文件: backward.cu，BACKWARD::render()。

复用前向 pass 的排序数组和 tile 范围，从后往前遍历高斯列表。不存储逐步的中间不透明度，而是利用前向 pass 保存的 final_T 和每个点的 α 反推中间系数。

总结：论文流程 → 代码映射

七、训练主循环（论文 § 5, train.py）

7.1 训练时间线

iter:   0       500           3000          7000        15000          30000
        |        |             |             |            |              |
        v        v             v             v            v              v
     warm-up   开始           重置          评估+保存     停止           评估+保存
     不densify  densify       opacity       (test)      densify        训练结束
               (每100轮)     (每3000轮)                 (此后只优化)

SH:    0     →1(1k)     →2(2k)     →3(3k)      之后保持 3 阶

7.2 初始化阶段

训练正式循环开始之前，需要完成模型、场景、优化器等的初始化：

# train.py:48-64
first_iter = 0
tb_writer = prepare_output_and_logger(dataset)           # 创建输出目录 + TensorBoard
gaussians = GaussianModel(dataset.sh_degree, opt.optimizer_type)  # 创建高斯模型（默认 sh_degree=3）
scene = Scene(dataset, gaussians)                        # 加载场景：读取点云/图片，初始化高斯球
gaussians.training_setup(opt)                            # 配置 Adam 优化器，为6组参数分别设置 lr
if checkpoint:                                           # 支持从 checkpoint 恢复训练
    (model_params, first_iter) = torch.load(checkpoint)
    gaussians.restore(model_params, opt)

bg_color = [1, 1, 1] if dataset.white_background else [0, 0, 0]
background = torch.tensor(bg_color, dtype=torch.float32, device="cuda")

iter_start = torch.cuda.Event(enable_timing = True)      # CUDA 计时器，用于统计每轮耗时
iter_end = torch.cuda.Event(enable_timing = True)

use_sparse_adam = opt.optimizer_type == "sparse_adam" and SPARSE_ADAM_AVAILABLE
depth_l1_weight = get_expon_lr_func(opt.depth_l1_weight_init, opt.depth_l1_weight_final, max_steps=opt.iterations)  # 深度正则权重的衰减函数

viewpoint_stack = scene.getTrainCameras().copy()          # 训练相机列表（用于无放回采样）
viewpoint_indices = list(range(len(viewpoint_stack)))

7.3 每轮迭代做什么

# train.py:73
for iteration in range(first_iter, opt.iterations + 1):   # 默认 30000 轮

① 更新学习率 + 提升 SH 阶数

# train.py:91-95
gaussians.update_learning_rate(iteration)   # 位置 lr 按指数衰减 + 曝光 lr 衰减

if iteration % 1000 == 0:
    gaussians.oneupSHdegree()               # SH 阶数: 0→1→2→3，每 1000 轮提升一阶

每 1000 轮提升一阶 SH，是一种课程学习——先学均匀颜色，再学视角相关的高光。

② 随机选一个训练视角

# train.py:97-103
if not viewpoint_stack:                                      # 所有视角用完 → 重新装填（新 epoch）
    viewpoint_stack = scene.getTrainCameras().copy()
    viewpoint_indices = list(range(len(viewpoint_stack)))
rand_idx = randint(0, len(viewpoint_indices) - 1)
viewpoint_cam = viewpoint_stack.pop(rand_idx)                # 无放回随机采样：每个 epoch 内每个视角恰好用一次
vind = viewpoint_indices.pop(rand_idx)

• 每轮迭代用 1 个视角
• 每个 epoch 用完 所有视角（假设有 200 张训练图，那 200 轮迭代 = 1 个 epoch）
• 总共 30000 轮，如果有 200 张图，就跑了 150 个 epoch

③ 渲染

# train.py:109-116
bg = torch.rand((3), device="cuda") if opt.random_background else background  # 可选随机背景色（数据增强）

render_pkg = render(viewpoint_cam, gaussians, pipe, bg,
                    use_trained_exp=dataset.train_test_exp, separate_sh=SPARSE_ADAM_AVAILABLE)
image, viewspace_point_tensor, visibility_filter, radii = \
    render_pkg["render"],              \  # 渲染图像 (3, H, W) → 用于算 loss
    render_pkg["viewspace_points"],    \  # 屏幕空间点（带梯度）→ 用于 densification
    render_pkg["visibility_filter"],   \  # 可见高斯掩码 → 只对可见点做梯度累积
    render_pkg["radii"]                   # 屏幕半径 → 用于剪枝判断

if viewpoint_cam.alpha_mask is not None:                    # 如果有 alpha mask，遮掉无效区域
    alpha_mask = viewpoint_cam.alpha_mask.cuda()
    image *= alpha_mask

④ 计算损失（论文 Eq.7）

# train.py:119-140
gt_image = viewpoint_cam.original_image.cuda()               # 真值图像
Ll1 = l1_loss(image, gt_image)                               # L1 损失：逐像素对齐颜色
if FUSED_SSIM_AVAILABLE:
    ssim_value = fused_ssim(image.unsqueeze(0), gt_image.unsqueeze(0))  # 融合版 SSIM（更快）
else:
    ssim_value = ssim(image, gt_image)

loss = (1.0 - opt.lambda_dssim) * Ll1 + opt.lambda_dssim * (1.0 - ssim_value)  # 默认 lambda=0.2

# 深度正则（可选）：如果有可靠的单目深度先验，加入逆深度 L1 约束
Ll1depth_pure = 0.0
if depth_l1_weight(iteration) > 0 and viewpoint_cam.depth_reliable:
    invDepth = render_pkg["depth"]
    mono_invdepth = viewpoint_cam.invdepthmap.cuda()
    depth_mask = viewpoint_cam.depth_mask.cuda()
    Ll1depth_pure = torch.abs((invDepth  - mono_invdepth) * depth_mask).mean()
    Ll1depth = depth_l1_weight(iteration) * Ll1depth_pure    # 权重随训练衰减
    loss += Ll1depth

• 只用 L1：像素级对齐了，但图像可能模糊、缺乏纹理细节（L1 对均匀误差不敏感）
• 只用 SSIM：结构好看了，但颜色可能整体偏移（SSIM 对全局亮度偏移不太敏感）
• 两者结合：L1 管”颜色准不准”，D-SSIM 管”看起来像不像”

⑤ 反向传播

# train.py:142
loss.backward()
# 梯度回传到所有高斯参数（_xyz, _scaling, _rotation, _opacity, _features_dc/rest）
# 同时 viewspace_point_tensor.grad 被填充 → 用于 densification 的屏幕空间梯度

⑥ 自适应密度控制（仅 iter 500~15000）

# train.py:163-174
if iteration < opt.densify_until_iter:                       # 默认 15000 轮后停止 densify
    # 记录每个可见高斯在屏幕上的最大半径（用于后续 prune 大高斯）
    gaussians.max_radii2D[visibility_filter] = torch.max(
        gaussians.max_radii2D[visibility_filter], radii[visibility_filter])
    gaussians.add_densification_stats(viewspace_point_tensor, visibility_filter)  # 累积屏幕空间梯度

    if iteration > opt.densify_from_iter and iteration % opt.densification_interval == 0:
        # 默认: iter>500 且每 100 轮执行一次 clone + split + prune
        size_threshold = 20 if iteration > opt.opacity_reset_interval else None
        gaussians.densify_and_prune(opt.densify_grad_threshold, 0.005,
                                    scene.cameras_extent, size_threshold, radii)

    if iteration % opt.opacity_reset_interval == 0 \
       or (dataset.white_background and iteration == opt.densify_from_iter):
        gaussians.reset_opacity()            # 每 3000 轮重置 opacity → 让优化器重新筛选有用的高斯

后半段（15000轮之后）停止 densify，只做参数精细优化，不再增删高斯球。

⑦ 优化器更新

# train.py:177-186
if iteration < opt.iterations:
    gaussians.exposure_optimizer.step()                      # 更新曝光补偿参数
    gaussians.exposure_optimizer.zero_grad(set_to_none = True)
    if use_sparse_adam:
        visible = radii > 0
        gaussians.optimizer.step(visible, radii.shape[0])    # 稀疏 Adam：只更新可见高斯（更快）
        gaussians.optimizer.zero_grad(set_to_none = True)
    else:
        gaussians.optimizer.step()                           # 标准 Adam：更新所有高斯参数
        gaussians.optimizer.zero_grad(set_to_none = True)

⑧ 保存与评估

# train.py:157-161  保存模型（默认在 iter 7000 和 30000）
training_report(tb_writer, iteration, Ll1, loss, l1_loss,
                iter_start.elapsed_time(iter_end), testing_iterations,
                scene, render, ...)                          # 写 TensorBoard 日志 + 在测试集上评估 PSNR
if (iteration in saving_iterations):
    print("\n[ITER {}] Saving Gaussians".format(iteration))
    scene.save(iteration)                                    # 保存到 point_cloud/iteration_N/point_cloud.ply

# train.py:188-190  保存 checkpoint（可从此恢复训练）
if (iteration in checkpoint_iterations):
    print("\n[ITER {}] Saving Checkpoint".format(iteration))
    torch.save((gaussians.capture(), iteration),
               scene.model_path + "/chkpnt" + str(iteration) + ".pth")

7.4 单轮流程总结

更新学习率 → 提升SH阶 → 选视角 → 渲染 → 算 loss → 反向传播 → densify(前半段) → Adam 更新 → 保存/评估