darshanmakwana412/tsplat

来源：https://github.com/darshanmakwana412/tsplat

tsplat

在终端中运行高斯溅射渲染，仅需CPU，甚至可通过SSH使用

demo screenshot

tsplat 使用 Unicode 半块字符或任何支持的图形协议，直接在终端中渲染 3D 高斯溅射场景（基于 https://repo-sam.inria.fr/fungraph/3d-gaussian-splatting/）。用 Rust 编写，目前仅支持 CPU，无需 GPU 或显示服务器，甚至可通过 SSH 工作。

安装

需要 Rust。

cargo install --git https://github.com/darshanmakwana412/tsplat

或者克隆并在本地构建：

git clone https://github.com/darshanmakwana412/tsplat
cd tsplat
cargo build --release

快速开始

你需要一个 INRIA 3DGS 格式的 .ply 场景。预训练的 garden 场景（https://repo-sam.inria.fr/fungraph/3d-gaussian-splatting/）是一个很好的初次测试：

tsplat path/to/scene.ply
# 提高或移除默认的 20 万 splat 上限
tsplat path/to/scene.ply --max-splats 500000
tsplat path/to/scene.ply --no-cap
# 当场景看起来褪色时
tsplat path/to/scene.ply --raw-opacity

注意：如果你的终端不支持任何图形协议，将回退到半块字符渲染，可能看起来像 Minecraft 风格。

控制键

基准测试

cargo run --release --bin bench_forward
cargo run --release --bin bench_forward -- --threads 1,2,4,8 --splats 200000
cargo run --release --bin bench_forward -- --ply path/to/scene.ply

高斯溅射教程

我最近发现了一些相关的笔记，决定在这个周末将其数字化，在这个过程中我用 Rust 重新实现了前向光栅化过程，并决定写一篇向所有人解释高斯溅射的教程。

• 什么是高斯溅射？
• 前向传递流水线
• 第一步：投影溅射
  • 1.1：构建 3D 协方差矩阵
  • 1.2：变换到视图空间
  • 1.3：投影到 2D
• 第二步：计算边界框
• 第三步：深度排序
• 第四步：图块分箱
• 第五步：Alpha 混合
  • 图块并行混合
• 优化
  • 快速近似 exp
  • 行级提前退出
  • 临时缓冲区复用

什么是高斯溅射？

一个 3D 高斯溅射是空间中的一个有向椭球体，携带颜色和不透明度。你可以把它想象成一个模糊的彩色斑点。一个场景由成千上万个这样的斑点组成，当你从某个视角观察它们时，它们重叠并混合形成最终图像。

每个高斯溅射我们用以下属性表示：

pub struct Splat {
    pub pos: Vec3,      // 世界空间中的中心位置
    pub scale: Vec3,    // 沿每个局部轴的大小
    pub rot: Quat,      // 用单位四元数表示的方向
    pub color: Vec3,    // RGB 颜色（已从球谐函数解码）
    pub opacity: f32,   // 该斑点的不透明度，范围为 [0, 1]
}

尺度以对数空间存储，不透明度以 logits 存储，颜色以球谐系数存储，四元数归一化为单位长度，以确保值处于各自的有效范围内。

单个高斯溅射的球谐函数 (SH) 系数只是定义在单位球面上的颜色函数的频域表示。那么为什么用球谐函数呢？因为在现实世界中，表面的颜色取决于观察方向。SH 系数可以紧凑地编码这种与视点相关的外观。SH 函数按频带组织（就像音乐中的八度），频带越高，系数越多，因此能捕捉更精细的细节。INRIA 3DGS 格式最多存储到第 3 频带（每个 RGB 溅射 48 个系数）。

要解码 band-0，band-0 SH 基函数为 Y_0^0 = \frac{1}{2\sqrt{\pi}} \approx 0.282。从 SH 系数到 RGB 的转换公式为：

\text{color} = \text{clamp}\left(0.5 + C_0 \cdot f_{dc},\ 0,\ 1\right)

其中 C_0 = Y_0^0，f_{dc} 是文件中 3 分量的直流系数。

pub const SH_C0: f32 = 0.28209479177387814;
pub fn sh_band0_to_rgb(f_dc: Vec3) -> Vec3 {
    (Vec3::splat(0.5) + SH_C0 * f_dc).clamp(Vec3::ZERO, Vec3::ONE)
}

前向传递流水线

前向传递将一组 3D 高斯 + 一个摄像机转换为 2D 图像。以下是渲染流水线的概述：

pipeline

第一步：投影溅射

1.1：构建 3D 协方差矩阵

对于每个溅射，给定原始的 (尺度, 旋转) 对，我们需要构建一个 3D 协方差矩阵 \Sigma，描述高斯在空间中的形状和方向。公式为：

\Sigma = R \cdot S \cdot S^T \cdot R^T

其中 R 是从四元数得到的 3×3 旋转矩阵，S 是尺度对角矩阵。令 M = R·S，则简化为：

\Sigma = M \cdot M^T

let r_mat = Mat3::from_quat(s.rot);
let s_mat = Mat3::from_diagonal(s.scale);
let m = r_mat * s_mat;
let cov3d = m * m.transpose();

注意：为什么我们要这样分解协方差？
协方差矩阵只有在半正定时才有物理意义。梯度下降很难约束生成有效的矩阵。通过将协方差表示为 M \cdot M^T，它被保证是半正定的——A^T A 形式的矩阵总是半正定的。这是一种重新参数化技巧：我们分别优化无约束的 尺度 和 旋转，由此推导出的协方差始终有效。

这个矩阵实际上长什么样？对于尺度为 (0.1, 0.05, 0.02)、旋转为单位矩阵的溅射：

\Sigma = \begin{pmatrix} 0.01 & 0 & 0 \\ 0 & 0.0025 & 0 \\ 0 & 0 & 0.0004 \end{pmatrix} \quad \begin{aligned} &= \text{diag}(0.1^2,\; 0.05^2,\; 0.02^2) \end{aligned}

在单位旋转下，它只是对角线上尺度的平方，是一个轴对齐的椭球体。

1.2：变换到视图空间

我们刚刚计算的 3D 协方差位于世界空间。要将其投影到摄像机的图像平面上，首先需要将其旋转到视图空间——即摄像机位于原点、朝向 −z 的坐标系。对于溅射中心，这只需用 4×4 视图矩阵进行矩阵-向量乘法：

let p_view4 = view * Vec4::new(s.pos.x, s.pos.y, s.pos.z, 1.0);
let p_view = Vec3::new(p_view4.x, p_view4.y, p_view4.z);
if p_view.z > -znear || p_view.z < -zfar { return None; }
let zc = -p_view.z;

注意 zc = -p_view.z。我们的视图空间是右手系，摄像机朝向 −z，因此摄像机前方的点 z 为负。我们使用 zc（前方为正）作为深度用于排序和投影。

对于协方差，我们用视图矩阵的 3×3 部分 W 旋转它：

\Sigma_{view} = W \cdot \Sigma \cdot W^T

let w_mat = Mat3::from_mat4(view);
let w_mat_t = w_mat.transpose();
let cov3d_view = w_mat * cov3d * w_mat_t;

这仅仅是协方差矩阵的标准基变换公式。椭球体的形状不会因此改变——实际上我们只是用摄像机坐标系重新表达它。

1.3：投影到 2D

现在我们在视图空间中有了一个 3D 高斯，需要将其投影到 2D 图像平面上。投影是透视的，这意味着 3D 高斯不会投影为精确的 2D 高斯，因为透视是一种非线性变换。但我们可以用投影函数的雅可比矩阵对其进行局部线性化，结果足够接近。

投影函数将视图空间中的 3D 点 (x, y, z) 映射到像素坐标 (u, v)：

u = f_x \cdot \frac{x}{z_c} + c_x v = f_y \cdot \frac{y}{z_c} + c_y

其中 f_x, f_y 是焦距，c_x, c_y 是主点（图像中心）。为简化计算，我们也可以假设 f_x = f_y。

在溅射中心处计算的投影雅可比矩阵 J 为：

J = \begin{bmatrix} \frac{f_x}{z_c} & 0 & \frac{f_x \cdot x_v}{z_c^2} \\ 0 & \frac{f_y}{z_c} & \frac{f_y \cdot y_v}{z_c^2} \\ 0 & 0 & 0 \end{bmatrix}

这个矩阵的结构非常稀疏，9 个条目中只有 4 个非零。我们可以直接做完整的 JCJ^T，进行两次 3×3 矩阵乘法（约 54 次标量乘法）。但由于 J 的第三行全为零，我们只需要 J Cov3D_view J^T 的左上角 2\times 2 部分。同时 J 的前两行各有只有两个非零条目。因此，我们可以通过手动展开乘积，只用约 20 次标量乘法计算 2D 协方差：

let c = &cov3d_view;
let inv_zc = 1.0 / zc;
let inv_zc2 = inv_zc * inv_zc;
let j00 = fx * inv_zc;
let j02 = fx * xv * inv_zc2;
let j11 = fy * inv_zc;
let j12 = fy * yv * inv_zc2;

// Row 0 of J * C: [j00*c00 + j02*c20, j00*c01 + j02*c21, j00*c02 + j02*c22]
let t0x = j00 * c.x_axis.x + j02 * c.z_axis.x;
let t0y = j00 * c.y_axis.x + j02 * c.z_axis.y;
let t0z = j00 * c.x_axis.z + j02 * c.z_axis.z;

// Row 1 of J * C: [j11*c10 + j12*c20, j11*c11 + j12*c21, j11*c12 + j12*c22]
let t1y = j11 * c.y_axis.y + j12 * c.z_axis.y;
let t1z = j11 * c.y_axis.z + j12 * c.z_axis.z;

// 2D cov = (J*C) * J^T, top-left 2x2:
let cov2d_00 = t0x * j00 + t0z * j02 + eps2d;
let cov2d_01 = t0y * j11 + t0z * j12;
let cov2d_11 = t1y * j11 + t1z * j12 + eps2d;

注意对角线上的 eps2d。这是一个小的扩张（默认 0.3），用于数值稳定性，它确保 2D 协方差是严格正定的（而不仅仅是半正定），这意味着它总是可逆的。

注意：eps2d 技巧为什么有效
根据构造，2D 协方差 JCJ^T 只是半正定的（A^T A 形式）。但我们稍后需要对其求逆（用于在每个像素处计算高斯值）。奇异矩阵不可逆。在对角线上添加 eps2d 相当于给矩阵加上 \lambda I。对于任意向量 x：

x^T \cdot (A^T A + \lambda I) \cdot x = \|Ax\|^2 + \lambda \|x\|^2 > 0

对于任意非零 x 这是严格为正的，正是正定矩阵的定义，可逆，且所有特征值严格为正。

接下来我们求 2\times 2 协方差的逆。对于 2\times 2 矩阵，逆有封闭形式：

\begin{bmatrix} a & b \\ b & d \end{bmatrix}^{-1} = \frac{1}{ad - b^2} \begin{bmatrix} d & -b \\ -b & a \end{bmatrix}

let det = cov2d_00 * cov2d_11 - cov2d_01 * cov2d_01;
if det <= 0.0 { return None; }
let inv_det = 1.0 / det;
let cov2d_inv = Mat2::from_cols(
    Vec2::new(cov2d_11 * inv_det, -cov2d_01 * inv_det),
    Vec2::new(-cov2d_01 * inv_det, cov2d_00 * inv_det),
);

屏幕位置是标准的透视除法：

let sx = fx * xv * inv_zc + cx;
let sy = fy * yv * inv_zc + cy;

至此，对于每个溅射我们有了：屏幕位置 (sx, sy)，深度 zc，以及逆 2D 协方差矩阵 cov2d_inv。这就是在屏幕上任何像素处计算该高斯值所需的一切。

第二步：计算边界框

我们不需要在屏幕的每个像素处计算溅射，因为 2D 高斯具有无限支撑，它永远不会真正达到零。但我们可以计算一个边界框，包围高斯具有可见效果的区域，只计算该框内的像素。

原始的 3DGS 代码计算 2D 协方差的两个特征值 \lambda_1、\lambda_2（椭圆两个主轴上的方差），取 r = 3\sqrt{max(\lambda_1, \lambda_2)}（3σ 规则，覆盖 99.7% 的高斯），然后用该半径的圆作为边界框。这很简单但很浪费。当高斯被拉长（一个特征值远大于另一个）时，边界圆会包含大量空区域。

bounding box

我们可以通过观察以下两点来创建更紧凑的边界框：

1. 沿每个轴的范围是 k\sqrt{\Sigma_{ii}}，其中 \Sigma_{ii} 是该轴对应的 2D 协方差对角线元素。对于一个拉长的椭圆，短轴的范围远小于长轴的范围，因此边界框更紧凑。
2. 经典的 3\sigma 规则是保守的。一个微弱溅射（低不透明度）不需要 3\sigma，因为它的贡献很快就会降到可见阈值以下。截止距离 k 可以在溅射的不透明度低于某个阈值 \tau 时计算出来。基于此，每个像素的值为

\alpha = \text{opacity} \cdot \exp(-\tfrac{1}{2} \cdot \mathbf{d}^T \Sigma^{-1} \mathbf{d})

我们希望 \alpha \geq \tau，这可以整理为：

\mathbf{d}^T \Sigma^{-1} \mathbf{d} \leq 2 \ln\left(\frac{\text{opacity}}{\tau}\right) = k^2

对于一个接近不透明的溅射（不透明度 = 1），k^2 = 2ln(255) = 11.1，接近 3\sigma 的值 9。对于一个微弱溅射（不透明度 = 0.1），k^2 = 2ln(255) = 11.1，因此框和计算量会大幅减小。

if s.opacity <= alpha_threshold { return None; }
let k2 = (2.0 * (s.opacity / alpha_threshold).ln()).min(max_k2);
if !(k2 > 0.0) { return None; }
let rx_f = (k2 * cov2d_00).sqrt();
let ry_f = (k2 * cov2d_11).sqrt();
if !rx_f.is_finite() || !ry_f.is_finite() || rx_f < 1.0 || ry_f < 1.0 { return None; }
let x0 = (sx - rx_f).floor() as i32;
let y0 = (sy - ry_f).floor() as i32;
let x1 = (sx + rx_f).ceil() as i32;
let y1 = (sy + ry_f).ceil() as i32;
// 裁剪到帧缓冲。
let x0 = x0.max(0);
let y0 = y0.max(0);
let x1 = x1.min(w_i - 1);
let y1 = y1.min(h_i - 1);
if x0 > x1 || y0 > y1 { return None; }

另外，由于每个溅射相互独立，投影计算是令人尴尬的并行。我们使用 rayon 的 par_iter().filter_map() 在所有核心上并行投影所有溅射：

pub fn project(
    splats: &[Splat],
    camera: &OrbitCamera,
    params: &RenderParams,
    pool: &Option,
) -> Vec {
    // ... 预计算视图矩阵、内参等（每帧一次）...
    let do_project = || {
        splats
            .par_iter()
            .filter_map(|s| {
                // ... 上述所有数学运算，返回 Some(Projected) 或 None ...
            })
            .collect()
    };
    match pool.as_ref() {
        Some(p) => p.install(do_project),
        None => do_project(),
    }
}

现在我们可以将所有内容打包到一个 Projected 结构中：

pub struct Projected {
    pub screen: Vec2,      // 像素中心 (sx, sy)
    pub depth: f32,        // zc（前方为正）
    pub cov2d_inv: Mat2,   // 逆 2D 协方差
    pub bbox: [i32; 4],    // 包含边界：x0, y0, x1, y1
    pub color: Vec3,       // RGB
    pub opacity: f32,      // [0, 1]
}

第三步：深度排序

对于 Alpha 混合，我们使用反向画家算法（https://en.wikipedia.org/wiki/Painter’s_algorithm）。这意味着在混合前需要对投影的溅射按深度排序。顺序很重要：如果一个近处的溅射遮挡了远处的溅射，它必须首先合成，以便在混合过程中具有更大的权重。

由于所有深度均为正（我们已经剔除了摄像机后面的点），因此 f32 的位模式在重新解释为 u32 时保持了浮点顺序。这使我们能够使用 depth.to_bits() 作为排序键进行实际排序。

我们使用简单的 2 遍 16 位基数排序。对于我们的输入规模（~10万–20万个元素），这比基于比较的排序（如 Rust 的 sort_unstable_by_key）更快，并且运行时间为 O(n)，常数小。2 遍 16 位比 4 遍 8 位效果更好：更少的遍数意味着更少的数据遍历，并且 65536 个条目的直方图（每个 256KB）可以舒适地放在我笔记本的 L2 缓存中。在 20 万溅射下，速度持续快 2 倍。