4D Gaussian 数学手册 · 从推导到代码

START WITH SCHOOL MATH

先补 8 个工具，每次只多走一步。

只预设你会四则运算、分数、平方和解简单方程。向量、矩阵、指数、导数都在这里引入。后面的“定义 / 模型选择”不是从某条物理定律必然推出的；“恒等式”可以代数证明；“近似”会写出适用范围。

A · 向量、点积、长度：一串数怎样参加运算？

向量先当作一列数。位置 (3,4) 表示向右 3、向上 4；相加就是每项相加。粗体表示整列，普通字母表示一个数。点积把对应项相乘再相加；勾股定理给出长度。转置 T 只把列改成行，不改变数值。

\mathbf a=\begin{bmatrix}3\\4\end{bmatrix},\quad \mathbf b=\begin{bmatrix}2\\-1\end{bmatrix},\quad \mathbf a+\mathbf b=\begin{bmatrix}5\\3\end{bmatrix},\quad \mathbf a^T\mathbf b=3\cdot2+4(-1)=2

\|\mathbf a\|=\sqrt{3^2+4^2}=5,\qquad \widehat{\mathbf a}=\frac{\mathbf a}{5}=(0.6,0.8)^T

“归一化”在这里指除以长度，结果长度为 1；不是把每个坐标都变成 1。分母为 0 时没有方向，不能归一化。外积则保留所有两两乘积，产出一张表：

\mathbf a\mathbf a^T=\begin{bmatrix}3\cdot3&3\cdot4\\4\cdot3&4\cdot4\end{bmatrix}=\begin{bmatrix}9&12\\12&16\end{bmatrix}

求和符号 ∑ 就是把指定范围内的项全部加起来；连乘符号 ∏ 则把它们全部乘起来。

B · 矩阵：把很多条一次方程装进表格

矩阵乘列向量，就是每一行与这列做一次点积。2×3 表示 2 行 3 列；把 3 个输入变成 2 个输出。矩阵乘法 AB 表示先做 B 再做 A，一般不能交换。

\begin{bmatrix}2&1\\0&3\end{bmatrix}\begin{bmatrix}4\\5\end{bmatrix}=\begin{bmatrix}2\cdot4+1\cdot5\\0\cdot4+3\cdot5\end{bmatrix}=\begin{bmatrix}13\\15\end{bmatrix}

(AB)_{ij}=\sum_k A_{ik}B_{kj},\quad I=\begin{bmatrix}1&0\\0&1\end{bmatrix},\quad A^{-1}A=I,\quad (AB)^T=B^TA^T

I 是原样输出的单位矩阵；逆矩阵是撤销一次变换。不可撤销（例如把某一维压成 0）就没有逆矩阵。旋转保持长度，因此逆旋转等于转置。diag 表示只在对角线上放数，其余为 0。

R^TR=I\ \Longrightarrow\ R^{-1}=R^T,\qquad \operatorname{diag}(2,3)=\begin{bmatrix}2&0\\0&3\end{bmatrix}

C · 指数与对数：从反复相乘到连续增长

2³ 是 2 连乘 3 次；2⁻³ 是它的倒数；2^(1/2) 是 √2。把“份数”继续细分，可以定义任意实数指数。e≈2.71828 是特殊的底数：它的增长率恰好等于自己。ln 是指数的反操作，可用计算器。

e^0=1,\quad e^{a+b}=e^ae^b,\quad e^{-a}=\frac1{e^a},\quad \ln(e^a)=a,\quad \ln(ab)=\ln a+\ln b\ (a,b>0)

指数永远为正。负平方放在指数上，离 0 越远就越接近 0，这就是 Gaussian 的形状。ln 在正数范围单调递增，所以不等式两边取 ln 不会改变方向。

e^{-x^2/2}:\quad x=0\Rightarrow1,\quad x=1\Rightarrow0.60653,\quad x=2\Rightarrow0.13534

为什么 e 的导数还是 e？从 e=lim(1+1/n)ⁿ 引入的实数指数有 lim(h→0)(eʰ−1)/h=1。这个极限是这里新增的基础事实；结合指数乘法规则，就能推出后面使用的导数。

D · 均值、方差与协方差：用“平均乘积”描述形状

平均值就是总和除以个数。方差是“偏离平均值的距离的平方”的平均，不是平均距离。E 读作期望，这里先理解为大量样本的平均；Cov 是协方差。我们描述分布，分母用 N；样本无偏估计的 N−1 是另一个问题。

\mu=\frac1N\sum_{n=1}^N x_n,\quad \operatorname{Var}(X)=\frac1N\sum_n(x_n-\mu)^2,\quad \Sigma_{ij}=\mathbb E[(X_i-\mu_i)(X_j-\mu_j)]

样本 (1,2)、(3,6) 的均值为 (2,4)。偏移为 (−1,−2)、(1,2)，两张外积表相同：

\Sigma=\frac12\left(\begin{bmatrix}1&2\\2&4\end{bmatrix}+\begin{bmatrix}1&2\\2&4\end{bmatrix}\right)=\begin{bmatrix}1&2\\2&4\end{bmatrix}

对角元素是各轴方差；非对角元素表示两个方向是否一起偏大。这里两点共线，所以矩阵不可逆。GS 用三个正尺度和旋转构造有厚度的椭球。“正定”就是任意非零方向上的方差都大于 0。

E · 导数与偏导：从两点之间的斜率开始

汽车走了 Δ路程、用了 Δ时间，比值就是平均速度。把间隔 h 越缩越小，得到某一刻的变化率，记作 df/dx。不是令 h=0 后除以 0，而是观察 h 接近 0 时的比值。

\frac{(x+h)^2-x^2}{h}=\frac{2xh+h^2}{h}=2x+h\ \xrightarrow{h\to0}\ 2x

\frac{(x+h)^3-x^3}{h}=3x^2+3xh+h^2\to3x^2,\qquad \frac{e^{x+h}-e^x}{h}=e^x\frac{e^h-1}{h}\to e^x

偏导数 ∂ 是“先锁住其他输入，只改一个”。函数有多个输出时，把每个输出对每个输入的偏导排成表，叫 Jacobian。函数只有一个输出时，各偏导组成的列叫梯度 ∇。

f(x,y)=x^2+3y:\quad \frac{\partial f}{\partial x}=2x,\quad\frac{\partial f}{\partial y}=3,\quad \nabla f=(2x,3)^T

乘法法则来自 (f+Δf)(g+Δg) 展开，除以 h 后 ΔfΔg/h 趋于 0。用它对 g·(1/g)=1 求导，得到倒数法则；再组合成商法则：

(fg)^{\prime}=f^{\prime}g+fg^{\prime},\qquad (g^{-1})^{\prime}=-\frac{g^{\prime}}{g^2},\qquad \left(\frac fg\right)^{\prime}=\frac{f^{\prime}g-fg^{\prime}}{g^2}

F · 链式法则、反向传播与一步更新

输入多 0.01，让中间量多约 0.03；中间量每多 1 又让输出多约 5。总变化就是 3×5×0.01。这叫链式法则：沿一条路线相乘，多条路线相加。

u=x^2,\quad y=e^u:\qquad \frac{dy}{dx}=\frac{dy}{du}\frac{du}{dx}=e^{x^2}\cdot2x

L=L(a(x),b(x)):\qquad \frac{dL}{dx}=\frac{\partial L}{\partial a}\frac{da}{dx}+\frac{\partial L}{\partial b}\frac{db}{dx}

“反向”是从最终误差向输入传斜率。一个量被使用两次，两条路的贡献都要加。若梯度为正，减小参数能在足够小的步长下减小误差：

L=(\theta-3)^2,\quad L^{\prime}=2(\theta-3),\quad \theta=1,\ \eta=0.1:\quad \theta_{\rm new}=1-0.1(-4)=1.4

η 是学习率。PyTorch 的 backward 负责算梯度；optimizer.step 才修改参数。上面是 SGD，项目教学训练器使用带历史平均和自适应步长的 Adam，两者不是同一个更新式。

G · 配方与求根：条件高斯和椭圆共用的代数

把一次项吸收进平方，叫配方。直接展开右侧可以验证；这也是后面消掉时空交叉项的办法。

x^2+2bx+c=(x+b)^2+c-b^2

解 λ²−sλ+p=0：先把常数移到右边，两边加 (s/2)²，再开平方；平方根有正负两支。

\lambda^2-s\lambda=-p\ \Longrightarrow\ \left(\lambda-\frac s2\right)^2=\frac{s^2-4p}{4}\ \Longrightarrow\ \lambda=\frac{s\pm\sqrt{s^2-4p}}2

矩阵版本也是逐项乘开。不要因为出现粗体就放弃手算：2×2 逆矩阵只需要解两组二元一次方程。

H · 积分与球面：把很多小块加起来

路程可以用很多小段“速度×时间”相加；小段无限细时，和的极限记作积分。积分号右侧的 ds 提醒我们每块的宽度是什么，单位不能丢。

\int_a^b f(s)\,ds=\lim_{N\to\infty}\sum_{j=0}^{N-1}f(a+j\Delta s)\Delta s,\qquad \Delta s=\frac{b-a}{N}

球谐还会用单位球面。方向 v=(x,y,z) 满足 x²+y²+z²=1；经纬角的正弦、余弦可理解为单位圆上点的两个坐标。球面上的“平均”按小面积加权。单位球面积为 4π；由面积元素 sinθ·dθ·dφ 积分得到。

\mathbf v=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),\qquad \int_{S^2}1\,d\Omega=\int_0^{2\pi}\int_0^\pi\sin\theta\,d\theta\,d\phi=4\pi

最后一个积分用 cos 的导数是 −sin：内层等于 [−cosθ] 从 0 到 π 的差，即 2。完整三角函数理论不是 GS 上手的前置门槛；知道它表示方向、能计算基函数，就可以先跑通颜色求和。

00

READ THIS FIRST

先约定符号，再推公式。

基础

空间用三维列向量，时间用独立的标量。先理解共同的投影与合成，再选择具体动态表示。

𝛍 (t)

世界坐标中的 Gaussian 中心，3×1

𝚺 (t)

空间协方差，3×3；单位是长度平方

𝐂

屏幕协方差，2×2；单位是像素平方

t, τ, ρ

归一化时间、时间中心、时间核宽度

α_{i}, T_{i}

像素处 alpha、到达该层前的透射率

𝐑, 𝐉

旋转/线性变换、投影 Jacobian

参数与协方差
相机与投影
椭圆与合成
时间与运动
损失与梯度

原始资料STG · §4.1 ↗Deformation 4D-GS · Wu et al., 2024 ↗Native 4DGS · Yang et al., 2024 ↗

01

THE PRIMITIVE

用协方差描述一个椭球。

必会

中心决定位置，协方差决定形状与方向，不透明度决定该基函数对图像的贡献。

(1.1)三维 Gaussian 核

G (𝐱) = e^{- \frac{1}{2} {(𝐱 - 𝛍)}^{T} 𝚺^{- 1} (𝐱 - 𝛍)}

这里是未归一化的核，不是概率密度。Splat 的 opacity 单独学习，不能额外乘概率密度归一化系数后仍期待相同亮度。

从零推导5 步 · 数字算例 · 源码 · 自测

模型定义 → 坐标变换

为什么一团软椭球会写成“负的二次型的指数”？

用到的小工具A 向量 B 矩阵 C 指数

先只画一条数轴
中心 μ 决定最高点在哪里。尺度 s 决定多远算“远”。先减中心，再除尺度，得到没有单位的距离 z；这样厘米和米不会改变形状。
$z=\frac{x-\mu}{s},\qquad z^2=\frac{(x-\mu)^2}{s^2}$
选择一个平滑的权重
我们希望中心为 1、两侧对称、越远越小。exp(−z²/2) 满足这些要求。选 Gaussian 是建模选择，并非这些要求只允许这一种函数；1/2 是让 s 对应标准差的常用约定。
$g(x)=e^{-z^2/2},\qquad g(\mu)=1,\qquad g(\mu+s)=e^{-1/2}$
把三条独立轴的权重相乘
指数相乘就是指数相加。三维局部坐标 ξ 的每一维都按自己的尺度计量；等权重的位置形成椭球。
$G=\prod_{k=1}^3 e^{-\xi_k^2/(2s_k^2)}=\exp\left[-\frac12\left(\frac{\xi_1^2}{s_1^2}+\frac{\xi_2^2}{s_2^2}+\frac{\xi_3^2}{s_3^2}\right)\right]$
让椭球可以旋转和移动
世界偏移 d=x−μ。用逆旋转 Rᵀ 把它转回椭球自己的坐标，ξ=Rᵀd。把三个平方写成矩阵乘法，再把 ξ 代进去。
$\xi=R^T\mathbf d,\quad G=\exp\left[-\frac12\mathbf d^T R\operatorname{diag}(s_1^{-2},s_2^{-2},s_3^{-2})R^T\mathbf d\right]$
把中间那张表命名为逆协方差
旋转的逆是转置，乘积求逆要倒序。因此括号中恰好是 Σ 的逆。这只是前面逐项平方相加的紧凑写法。
$\Sigma=R\operatorname{diag}(s_1^2,s_2^2,s_3^2)R^T,\qquad G=e^{-\frac12(\mathbf x-\mu)^T\Sigma^{-1}(\mathbf x-\mu)}$

代入数字，亲手核对

取 μ=0、R=I、尺度 (2,1,1)，点 (2,0,0) 刚好离中心一个“长轴尺度”。

r^2=2^2/2^2=1,\qquad G=e^{-1/2}=0.60653066

这一步在代码里在哪里？

GPU 最终评估的是投影后的二维核。local 已经完成“减中心、旋转、除尺度”；dot 就是平方相加。

r² ↔ dot(input.local,input.local)；G ↔ exp(power)

src/splat.wgsl · L130–136 ↗

fn fs_main(input: VertexOutput) -> @location(0) vec4<f32> {
  let power = -0.5 * dot(input.local, input.local);
  let alpha = min(0.99, input.color_opacity.a * exp(power));
  if (alpha < 0.0039215686) {
    discard;
  }
  return vec4<f32>(input.color_opacity.rgb * alpha, alpha);

停下来想一想：点改成 (4,0,0)，权重是多少？

标准化距离变成 2，平方为 4；exp(−2)≈0.135335。

带着这条公式，查看完整像素算例 →

(1.2)存储参数 → 有效参数

\begin{matrix} s_{k} & = & e^{ℓ_{k}} \\ o & = & sigmoid (β) = \frac{1}{1 + e^{- β}} \end{matrix}

log-scale 保证理论尺度为正；sigmoid 将基准 opacity 映射到 0–1。实际实现仍需防止溢出与下溢。

从零推导4 步 · 数字算例 · 源码 · 自测

参数化选择 + 解方程

训练得到任意实数，怎么保证尺度和 opacity 合法？

用到的小工具C 指数

给尺度一个总为正的出口
原始参数 ℓ 可以是负数；exp(ℓ) 总是正数。训练优化 ℓ，渲染再解码成尺度 s。
$s=e^\ell>0,\qquad \ell=\ln s$
把 opacity 变成胜算
先假设 0<o<1。比值 o/(1−o) 从 0 到无穷大；对它取 ln 后能覆盖所有实数，把这个数叫 logit β。
$\beta=\ln\frac{o}{1-o}\quad\Longrightarrow\quad e^\beta=\frac{o}{1-o}$
把 o 单独解出来
两边乘 1−o，展开，再把含 o 的项移到同一侧，最后除以系数。
$e^\beta-e^\beta o=o\ \Longrightarrow\ e^\beta=o(1+e^\beta)\ \Longrightarrow\ o=\frac{e^\beta}{1+e^\beta}=\frac1{1+e^{-\beta}}$
从文件到 GPU，只激活一次
PLY 的 scale_*、opacity 保存原始值。CPU 为可见性判断解码；GPU 从同一份原始值解码用于绘制。它们是两条并行读取路线，不是把 CPU 已激活的值再激活一次。
$\ell=\ln 2\Rightarrow s=2,\qquad \beta=0\Rightarrow o=\frac12$

代入数字，亲手核对

保存 0.8 的 opacity 时应写 logit，而不是直接写 0.8。

\beta=\ln(0.8/0.2)=\ln4\approx1.386294,\qquad \operatorname{sigmoid}(\ln4)=0.8

这一步在代码里在哪里？

r5.x 是 opacity logit；r5.yzw 是三个 log-scale。

β ↔ splat.r5.x；ℓ₁..₃ ↔ splat.r5.yzw

src/splat.wgsl · L65–73 ↗

  let temporal_scale = max(exp(splat.r1.x), 0.000001);
  let temporal_opacity = exp(-pow(dt / temporal_scale, 2.0));
  let opacity = (1.0 / (1.0 + exp(-splat.r5.x))) * temporal_opacity;
  let scales = exp(splat.r5.yzw);
  let quaternion = normalize(splat.r6 + dt * splat.r7);
  let w = quaternion.x;
  let x = quaternion.y;
  let y = quaternion.z;
  let z = quaternion.w;

教学训练器直接导出 opacity logit 和 log-scale；RGB 在导出前先 sigmoid，和 opacity 不同。

opacity ↔ rows[:,20]；log-scale ↔ rows[:,21:24]

scripts/train_toy_stg.py · L98–106 ↗

    rows = np.zeros((len(teacher[0]),32), dtype="<f4")
    rows[:, :3] = fitted[0].detach().cpu().numpy()
    rows[:, 3] = 0.5
    rows[:, 4] = math.log(0.55)
    rows[:, 8:11] = fitted[1].detach().cpu().numpy()
    rows[:, 17:20] = fitted[2].sigmoid().detach().cpu().numpy()
    rows[:, 20] = fitted[3].detach().cpu().numpy()
    rows[:, 21:24] = math.log(0.12)
    rows[:, 24] = 1  # quaternion wxyz

停下来想一想：scale_0=0 表示零尺寸吗？

不是。exp(0)=1；它表示长度尺度 1。

带着这条公式，查看完整像素算例 →

(1.3)协方差的安全参数化

\begin{matrix} 𝚺 & = & 𝐑 diag ({s_{1}}^{2}, {s_{2}}^{2}, {s_{3}}^{2}) 𝐑^{T} \\ 𝚺 & = & \sum_{k = 1}^{3} 𝐚_{k} {𝐚_{k}}^{T} \\ 𝐚_{k} & = & s_{k} 𝐫_{k} \end{matrix}

rₖ 是旋转矩阵的第 k 列。正尺度时 Σ 正定；轴向量外积求和与矩阵乘法等价，适合直接写进 shader。

从零推导4 步 · 数字算例 · 源码 · 自测

代数恒等式

怎样从三条有方向的轴得到一张协方差表？

用到的小工具B 矩阵 D 方差

从一个没有偏好的小球开始
设局部随机向量 Z 的均值为 0，三维方差均为 1，互不相关。它的协方差就是 I。这里随机变量只是帮助描述形状，渲染时不需要随机采样。
$\mathbb E[Z]=0,\qquad \mathbb E[ZZ^T]=I$
先拉伸，再旋转，最后平移
S 是三个尺度组成的对角矩阵。X−μ=RSZ。把偏移和自己的转置相乘，再求平均；固定矩阵可移到平均之外。
$X=\mu+RSZ,\quad \Sigma=\mathbb E[(RSZ)(RSZ)^T]=RS\,\mathbb E[ZZ^T]\,S^TR^T=RS^2R^T$
把矩阵的列写出来
A=RS 的第 k 列就是旋转的第 k 列乘尺度。矩阵乘法 AAᵀ 的每个元素是各列对应分量的乘积之和。
$A=[\mathbf a_1\ \mathbf a_2\ \mathbf a_3],\quad \mathbf a_k=s_k\mathbf r_k,\quad \Sigma_{ij}=\sum_{k=1}^3 a_{k,i}a_{k,j}$
得到轴外积求和，也解释正定
任意方向 v 的二次型是三个平方的和；三个正尺度对应可逆的 A，非零 v 不可能三个平方都为 0。
$\Sigma=\sum_k\mathbf a_k\mathbf a_k^T,\qquad \mathbf v^T\Sigma\mathbf v=\sum_k(\mathbf a_k^T\mathbf v)^2>0\quad(\mathbf v\ne0)$

代入数字，亲手核对

不旋转、尺度 (2,1,3)，三条轴分别为 (2,0,0)、(0,1,0)、(0,0,3)。

\Sigma=\operatorname{diag}(4,1,9)

这一步在代码里在哪里？

axis0/1/2 是这三条缩放后的轴。投影后仍按同样的外积规则累加，只是向量从三维变成二维。

a₁..₃ ↔ axis0/1/2

src/splat.wgsl · L77–93 ↗

  let axis0 = scales.x * vec3<f32>(
    1.0 - 2.0 * (y * y + z * z),
    2.0 * (x * y + w * z),
    2.0 * (x * z - w * y)
  );
  let axis1 = scales.y * vec3<f32>(
    2.0 * (x * y - w * z),
    1.0 - 2.0 * (x * x + z * z),
    2.0 * (y * z + w * x)
  );
  let axis2 = scales.z * vec3<f32>(
    2.0 * (x * z + w * y),
    2.0 * (y * z - w * x),
    1.0 - 2.0 * (x * x + y * y)
  );

  let center_clip = camera.view_proj * vec4<f32>(center, 1.0);

cov_a 是所有 x 分量平方之和；cov_b 是 xy 乘积之和；cov_c 是所有 y 分量平方之和。

C₁₁,C₁₂,C₂₂ ↔ cov_a,cov_b,cov_c

src/splat.wgsl · L94–101 ↗

  let d0 = projected_axis(center_clip, axis0);
  let d1 = projected_axis(center_clip, axis1);
  let d2 = projected_axis(center_clip, axis2);

  // Project J W Sigma W^T J^T by summing the projected covariance axes.
  let cov_a = dot(vec3<f32>(d0.x, d1.x, d2.x), vec3<f32>(d0.x, d1.x, d2.x)) + 0.3;
  let cov_b = d0.x * d0.y + d1.x * d1.y + d2.x * d2.y;
  let cov_c = dot(vec3<f32>(d0.y, d1.y, d2.y), vec3<f32>(d0.y, d1.y, d2.y)) + 0.3;

停下来想一想：整体位置平移 10 米会改变 Σ 吗？

不会。X 和 μ 一起平移，X−μ 不变。

带着这条公式，查看完整像素算例 →

(1.4)单位四元数与旋转矩阵

𝐑 (𝐪) = [\begin{matrix} 1 - 2 (y^{2} + z^{2}) & 2 (x y - w z) & 2 (x z + w y) \\ 2 (x y + w z) & 1 - 2 (x^{2} + z^{2}) & 2 (y z - w x) \\ 2 (x z - w y) & 2 (y z + w x) & 1 - 2 (x^{2} + y^{2}) \end{matrix}]

这里采用 q=(w,x,y,z)、右手系、主动旋转与列向量，且 ‖q‖=1。内存中的行/列布局和参考实现可能采用转置约定；最终应比较轴向量和 Σ，而非只比较名为 R 的数组。

从零推导5 步 · 数字算例 · 源码 · 自测

旋转的代数表示

为什么四个数能变成 shader 里的九个旋转项？

用到的小工具A 向量 B 矩阵

先约定四元数如何相乘
四元数记为 (w,v)，其中 v=(x,y,z)。下面的规则是定义。叉积 a×b 的三个分量依次为 (a₂b₃−a₃b₂,a₃b₁−a₁b₃,a₁b₂−a₂b₁)，可直接代数计算。
$(a,\mathbf u)(b,\mathbf v)=(ab-\mathbf u^T\mathbf v,\ a\mathbf v+b\mathbf u+\mathbf u\times\mathbf v)$
用单位四元数表示旋转
令 w²+x²+y²+z²=1；逆四元数就是共轭 (w,−v)。把待旋转的点 p 装成 (0,p)，用 q(0,p)q⁻¹ 旋转。角度 θ、单位轴 n 对应 q=(cos(θ/2),n sin(θ/2))；这一步定义了它与几何旋转的联系。
$q(0,\mathbf p)=(-\mathbf v^T\mathbf p,\ w\mathbf p+\mathbf v\times\mathbf p),\qquad q^{-1}=(w,-\mathbf v)$
再乘一次，把向量部分展开
用第一步的乘法规则，得到下式。第二行用叉积逐分量可验证的恒等式 (v×p)×v=||v||²p−(v·p)v；每一步只是乘法与合并同类项。
$\begin{aligned}\mathbf p^{\prime}&=(\mathbf v^T\mathbf p)\mathbf v+w^2\mathbf p+w(\mathbf v\times\mathbf p)-w(\mathbf p\times\mathbf v)-(\mathbf v\times\mathbf p)\times\mathbf v\\&=(w^2-\|\mathbf v\|^2)\mathbf p+2\mathbf v(\mathbf v^T\mathbf p)+2w(\mathbf v\times\mathbf p)\end{aligned}$
把三个坐标轴分别代进去
p=e₁ 时 v×e₁=(0,z,−y)；p=e₂ 时为 (−z,0,x)；p=e₃ 时为 (y,−x,0)。这三次计算的结果就是旋转矩阵的三列。
$\begin{aligned}Re_1&=(w^2+x^2-y^2-z^2,\ 2xy+2wz,\ 2xz-2wy)^T\\Re_2&=(2xy-2wz,\ w^2-x^2+y^2-z^2,\ 2yz+2wx)^T\\Re_3&=(2xz+2wy,\ 2yz-2wx,\ w^2-x^2-y^2+z^2)^T\end{aligned}$
用单位长度替换对角项
例如 w²+x²=1−y²−z²，因此第一个对角元素变成 1−2(y²+z²)。另外两个完全同理，得到上面的九项公式。
$w^2+x^2-y^2-z^2=1-2(y^2+z^2)$

代入数字，亲手核对

绕 z 轴转 90°：q=(√2/2,0,0,√2/2)。第一条轴向上，第二条轴向左。

R=\begin{bmatrix}0&-1&0\\1&0&0\\0&0&1\end{bmatrix},\qquad R(1,0,0)^T=(0,1,0)^T

这一步在代码里在哪里？

按本页的列向量约定，axis0 是 s₁Re₁，axis1 是 s₂Re₂，axis2 是 s₃Re₃。直接比较这三条轴，可避开 GLM 参考代码对矩阵命名的转置差异。

q=(w,x,y,z) ↔ quaternion.xyzw；注意 quaternion.x 存 w

src/splat.wgsl · L77–93 ↗

  let axis0 = scales.x * vec3<f32>(
    1.0 - 2.0 * (y * y + z * z),
    2.0 * (x * y + w * z),
    2.0 * (x * z - w * y)
  );
  let axis1 = scales.y * vec3<f32>(
    2.0 * (x * y - w * z),
    1.0 - 2.0 * (x * x + z * z),
    2.0 * (y * z + w * x)
  );
  let axis2 = scales.z * vec3<f32>(
    2.0 * (x * z + w * y),
    2.0 * (y * z - w * x),
    1.0 - 2.0 * (x * x + y * y)
  );

  let center_clip = camera.view_proj * vec4<f32>(center, 1.0);

停下来想一想：q=(1,0,0,0) 得到什么？

三个轴原样输出，R=I。

带着这条公式，查看完整像素算例 →

原始资料3DGS · Kerbl et al., 2023 ↗本次部署源码 · splat.wgsl ↗

02

WORLD → CAMERA → PIXELS

把三维不确定性投影到屏幕。

必推

本节推导采用 OpenCV 风格相机轴：x 向右、y 向下、z 向前。与 wgpu 的视图和裁剪空间衔接时，需要显式处理轴向和投影约定。

(2.1)世界 → 相机

\begin{matrix} 𝐱_{c} & = & 𝐑_{cw} 𝐱_{w} + 𝐭_{cw} \\ 𝐂_{w} & = & - {𝐑_{cw}}^{T} 𝐭_{cw} \end{matrix}

下标 cw 明确表示 world→camera。相机中心不是平移项 t；camera→world 位姿需要先求逆。

从零推导4 步 · 数字算例 · 源码 · 自测

换坐标 + 解一次方程

为什么相机位置不是文件中的平移 t？

用到的小工具B 矩阵

先把相机搬到原点
相机在世界中的位置为 Cw。一个点相对相机的位移是 xw−Cw，而不是 xw+Cw。
$\mathbf d_w=\mathbf x_w-\mathbf C_w$
再把世界方向转到相机方向
Rcw 把世界向量表示成相机坐标。展开括号，把不含 xw 的项叫 tcw。下标 cw 表示从 world 到 camera。
$\mathbf x_c=R_{cw}(\mathbf x_w-\mathbf C_w)=R_{cw}\mathbf x_w-R_{cw}\mathbf C_w,\qquad \mathbf t_{cw}=-R_{cw}\mathbf C_w$
反解相机中心
两边左乘 Rcwᵀ，利用 RᵀR=I；或令相机中心在相机坐标中等于 0，都能得到同一结果。
$R_{cw}^T\mathbf t_{cw}=-\mathbf C_w\quad\Longrightarrow\quad\mathbf C_w=-R_{cw}^T\mathbf t_{cw}$
放入 GPU 使用的四维表格
点的末尾写 1，平移会生效；方向末尾写 0，平移不会生效。WebGPU 宿主使用相机前方 −Z、上方 +Y；本页手算使用前方 +Z、下方 +Y。两者的相机坐标可用 B=diag(1,−1,−1) 转换。
$V=\begin{bmatrix}R&t\\0&1\end{bmatrix},\qquad V\begin{bmatrix}x\\1\end{bmatrix}=\begin{bmatrix}Rx+t\\1\end{bmatrix},\qquad B=\operatorname{diag}(1,-1,-1)$

代入数字，亲手核对

无旋转，相机在 (2,0,0)，世界点在 (3,1,5)。

\mathbf t=(-2,0,0)^T,\qquad \mathbf x_c=(3,1,5)^T+(-2,0,0)^T=(1,1,5)^T

这一步在代码里在哪里？

宿主传入列主序的 view 和 projection。view_projection=projection*view；不能直接把 camera→world pose 当 view。

R,t ↔ view 的左上 3×3 与平移列

src/stg_pass.rs · L247–250 ↗

        let view = Mat4::from_cols_array_2d(&view_columns);
        let projection = Mat4::from_cols_array_2d(&projection_columns);
        let view_projection = projection * view;
        let [width, height] = viewport.map(|v| v as f32);

教学训练器的相机旋转恒为 I，所以这里只减 cameras。这个简化不能替代任意相机外参。

Cw ↔ cameras；xc ↔ xyz

scripts/train_toy_stg.py · L27–35 ↗

    dt = times[:, None, None] - 0.5
    xyz = mean[None] + motion[None] * dt - cameras[:, None, :]
    x, y, z = xyz.unbind(-1)
    uv = torch.stack((focal*x/z + width/2, focal*y/z + height/2), -1)
    zero = torch.zeros_like(z)
    jacobian = torch.stack((focal/z, zero, -focal*x/z.square(),
                            zero, focal/z, -focal*y/z.square()), -1).reshape(*z.shape, 2, 3)
    covariance = 0.12**2 * (jacobian @ jacobian.transpose(-1, -2))
    covariance = covariance + torch.eye(2, device=device) * 0.3

停下来想一想：相机向右移动 1，静止物体的相机 x 怎样变化？

无旋转时减小 1，所以屏幕上向左移。

带着这条公式，查看完整像素算例 →

(2.2)针孔投影

π (x, y, z) = [\begin{matrix} \frac{f_{x} x}{z} + c_{x} \\ \frac{f_{y} y}{z} + c_{y} \end{matrix}], z > 0

fx、fy 是像素单位焦距；cx、cy 是主点。此处忽略镜头畸变，标定与缩放后的图像内参必须一致。

从零推导4 步 · 数字算例 · 源码 · 自测

相似三角形

三维点为什么要除以深度 z？

用到的小工具A 向量

从侧面看两只相似三角形
一只三角形的横边为 x、纵边为 z；图像平面距离为 f，投影横边为 x′。相同射线给出相同的边长比。我们把成像平面放在前方，因此不用物理底片的倒像符号。
$\frac{x^{\prime}}f=\frac xz\quad\Longrightarrow\quad x^{\prime}=f\frac xz$
从长度换成像素，再加主点
像素宽为 ax，焦距 f/ax 就是 fx。光轴打到图像上的位置是主点 cx，通常在图像附近中心，但不必恰好等于宽度的一半。
$f_x=\frac f{a_x},\qquad u=f_x\frac xz+c_x,\qquad v=f_y\frac yz+c_y$
用齐次坐标统一写法
先乘内参矩阵 K，得到 (fx x+cx z, fy y+cy z,z)，再除以最后一项。乘法是线性的，最后这次除法是非线性的。
$K\mathbf x_c=\begin{bmatrix}f_x&0&c_x\\0&f_y&c_y\\0&0&1\end{bmatrix}\begin{bmatrix}x\\y\\z\end{bmatrix}=\begin{bmatrix}f_xx+c_xz\\f_yy+c_yz\\z\end{bmatrix}$
对应屏幕像素中心
像素索引是整数，但采样位置通常是整数加 0.5。图像宽高缩放时 fx、fy、cx、cy 也必须按相同坐标约定缩放。这里忽略镜头畸变。
$f_x=f_y=100,\ (c_x,c_y)=(64,48),\ (x,y,z)=(1,2,10)\Rightarrow(u,v)=(74,68)$

代入数字，亲手核对

深度变成 20，x 和 y 不变：偏离主点的距离减半。

(u,v)=(100/20+64,\ 200/20+48)=(69,58)

这一步在代码里在哪里？

uv 的两项正是针孔公式。训练器用像素中心网格，和固定图形管线的采样位置衔接。

fx=fy ↔ focal；cx ↔ width/2；cy ↔ height/2

scripts/train_toy_stg.py · L27–35 ↗

    dt = times[:, None, None] - 0.5
    xyz = mean[None] + motion[None] * dt - cameras[:, None, :]
    x, y, z = xyz.unbind(-1)
    uv = torch.stack((focal*x/z + width/2, focal*y/z + height/2), -1)
    zero = torch.zeros_like(z)
    jacobian = torch.stack((focal/z, zero, -focal*x/z.square(),
                            zero, focal/z, -focal*y/z.square()), -1).reshape(*z.shape, 2, 3)
    covariance = 0.12**2 * (jacobian @ jacobian.transpose(-1, -2))
    covariance = covariance + torch.eye(2, device=device) * 0.3

停下来想一想：fx 增加一倍，主点不变时会怎样？

所有点相对主点的横向距离增加一倍，看起来放大。

带着这条公式，查看完整像素算例 →

(2.3)对投影求偏导

𝐉 = \frac{\partial π}{\partial 𝐱_{c}} = [\begin{matrix} \frac{f_{x}}{z} & 0 & - \frac{f_{x} x}{z^{2}} \\ 0 & \frac{f_{y}}{z} & - \frac{f_{y} y}{z^{2}} \end{matrix}]

从零推导4 步 · 数字算例 · 源码 · 自测

逐项求偏导

把“移动一点点，会挪多少像素”写成一张表。

用到的小工具E 导数 B 矩阵

只改 x，锁住 y 和 z
将 x+h 代入 u，再减原来的 u，除以 h。cx 抵消，答案与 h 无关。v 不含 x，所以第二行这一项为 0。
$\frac{u(x+h,y,z)-u(x,y,z)}h=\frac{f_xh/z}h=\frac{f_x}z,\qquad \frac{\partial v}{\partial x}=0$
只改 y
同理，u 不含 y；v 的变化率为 fy/z。
$\frac{\partial u}{\partial y}=0,\qquad \frac{\partial v}{\partial y}=\frac{f_y}z$
最容易漏掉的是深度
这次分母变了。先通分，再约掉 h，最后取 h 接近 0 的极限。不要把 z 当常数带进这一步。
$\frac{f_xx/(z+h)-f_xx/z}{h}=\frac{f_xx[z-(z+h)]}{hz(z+h)}=-\frac{f_xx}{z(z+h)}\to-\frac{f_xx}{z^2}$
把两组结果排成 2×3
第一行对应 u，第二行对应 v；三列依次对应 x,y,z。小位移 δx 经 J 变成像素位移，这只是中心附近的一阶近似。
$J=\begin{bmatrix}f_x/z&0&-f_xx/z^2\\0&f_y/z&-f_yy/z^2\end{bmatrix},\qquad \delta\mathbf p\approx J\delta\mathbf x_c$

代入数字，亲手核对

f=100、中心 (1,2,10)，只向右移动 0.01：屏幕约向右 0.1 像素。

J=\begin{bmatrix}10&0&-1\\0&10&-2\end{bmatrix},\quad J(0.01,0,0)^T=(0.1,0)^T

这一步在代码里在哪里？

torch.stack 的六项逐项对应这张表，reshape(...,2,3) 是把最后六个数排成两行。

J ↔ jacobian；每个 batch、每个 Gaussian 各有一张表

scripts/train_toy_stg.py · L27–35 ↗

    dt = times[:, None, None] - 0.5
    xyz = mean[None] + motion[None] * dt - cameras[:, None, :]
    x, y, z = xyz.unbind(-1)
    uv = torch.stack((focal*x/z + width/2, focal*y/z + height/2), -1)
    zero = torch.zeros_like(z)
    jacobian = torch.stack((focal/z, zero, -focal*x/z.square(),
                            zero, focal/z, -focal*y/z.square()), -1).reshape(*z.shape, 2, 3)
    covariance = 0.12**2 * (jacobian @ jacobian.transpose(-1, -2))
    covariance = covariance + torch.eye(2, device=device) * 0.3

停下来想一想：位于光轴 x=y=0 时，J 的第三列为什么全是 0？

沿光轴前后移动，投影中心仍是主点；但投影大小仍会变化。

带着这条公式，查看完整像素算例 →

(2.4)先证明线性变换的协方差

\begin{matrix} 𝐘 & = & 𝐀 𝐗 + 𝐛 \\ Cov (𝐘) & = & 𝐀 𝚺 𝐀^{T} \end{matrix}

由 Y−E[Y]=A(X−μ)，再对外积取期望得到。平移消去，因此不影响协方差。

从零推导4 步 · 数字算例 · 源码 · 自测

从平均值定义证明

为什么一个线性变换要在协方差左右各乘一次？

用到的小工具D 方差 B 矩阵

写出新均值
Y=AX+b；平均可以分配到加法，也可以把不变的 A、b 提到外面。有限样本时逐项求和就能验证。
$\mu_Y=\mathbb E[AX+b]=A\mu_X+b$
新样本减去新均值
同一个平移 b 被减掉了。协方差描述相对中心的形状，所以不关心中心平移到哪里。
$Y-\mu_Y=AX+b-(A\mu_X+b)=A(X-\mu_X)$
按定义做外积，再求平均
转置会倒转乘法顺序，所以右侧出现 Aᵀ 而不是 A。把固定矩阵移到期望号外面，内部就是原协方差。
$\begin{aligned}\operatorname{Cov}(Y)&=\mathbb E[(Y-\mu_Y)(Y-\mu_Y)^T]\\&=\mathbb E[A(X-\mu_X)(X-\mu_X)^TA^T]\\&=A\Sigma_XA^T\end{aligned}$
把矩阵语言翻回标量
第 i、j 个输出的共同变化来自所有输入分量对的贡献。下面就是上一行的逐项展开，不是另加的一条定理。
$(\Sigma_Y)_{ij}=\sum_k\sum_l A_{ik}(\Sigma_X)_{kl}A_{jl}$

代入数字，亲手核对

X 的二维协方差为 diag(1,4)，x 放大两倍、y 缩小一半。

A=\operatorname{diag}(2,1/2),\qquad A\Sigma_XA^T=\operatorname{diag}(4,1)

这一步在代码里在哪里？

教学例子 Σ=0.12²I，因此 JΣJᵀ 可简化成 0.12²*(J@Jᵀ)。任意各向异性 Σ 不能这样简化。

A ↔ jacobian；@ ↔ 矩阵乘法；transpose(-1,-2) ↔ 转置最后两个维度

scripts/train_toy_stg.py · L27–35 ↗

    dt = times[:, None, None] - 0.5
    xyz = mean[None] + motion[None] * dt - cameras[:, None, :]
    x, y, z = xyz.unbind(-1)
    uv = torch.stack((focal*x/z + width/2, focal*y/z + height/2), -1)
    zero = torch.zeros_like(z)
    jacobian = torch.stack((focal/z, zero, -focal*x/z.square(),
                            zero, focal/z, -focal*y/z.square()), -1).reshape(*z.shape, 2, 3)
    covariance = 0.12**2 * (jacobian @ jacobian.transpose(-1, -2))
    covariance = covariance + torch.eye(2, device=device) * 0.3

停下来想一想：一维尺度放大 3 倍，方差放大几倍？

9 倍，因为偏移先乘 3，再平方。

带着这条公式，查看完整像素算例 →

(2.5)一阶线性化与屏幕协方差

\begin{matrix} π (𝐗) & \approx & π (𝛍_{c}) + 𝐉 (𝐗 - 𝛍_{c}) \\ 𝚺_{c} & = & 𝐑_{cw} 𝚺 {𝐑_{cw}}^{T} \\ 𝐂 & = & 𝐉 𝚺_{c} 𝐉^{T} + ε 𝐈_{2} \end{matrix}

J 在 Gaussian 中心处求值。ε 是像素方差；当前 shader 使用 0.3。透视是非线性的，所以这是一阶近似，不是任意透视变换都精确保留 Gaussian。

从零推导5 步 · 数字算例 · 源码 · 自测

一阶近似 + 线性恒等式

为什么 GS 屏幕椭圆是 J R Σ Rᵀ Jᵀ？

用到的小工具E 导数 D 方差 B 矩阵

先把中心和形状都转到相机坐标
中心需要旋转加平移，形状只需要旋转。使用上一节的协方差恒等式。
$\mu_c=R_{cw}\mu+t_{cw},\qquad \Sigma_c=R_{cw}\Sigma R_{cw}^T$
用中心附近的切线替代弯曲投影
导数告诉我们很小的变化怎样传播。把三维点 X 写成中心加偏移；对每个输出保留一次变化，丢掉二次及更高阶变化。
$\pi(X)=\pi(\mu_c)+J(X-\mu_c)+O(\|X-\mu_c\|^2)$
平移项不改变方差
在线性近似里，π(μc) 是常数，J 是在中心算好的固定矩阵。因此屏幕形状用 JΣcJᵀ。再加 εI 表示额外的各向同性屏幕方差。
$C\approx J\Sigma_cJ^T+\varepsilon I_2,\qquad \varepsilon=0.3\ {\rm pixel}^2$
把矩阵乘法改写成三个小向量
若 Σ=∑ak akᵀ，令 dk=JRcw ak。分配乘法后，C 等于三条投影轴的外积之和。这正是 shader 的计算路线。
$C=\sum_{k=1}^3\mathbf d_k\mathbf d_k^T+\varepsilon I_2,\quad a=\sum d_{k,x}^2+\varepsilon,\quad b=\sum d_{k,x}d_{k,y},\quad c=\sum d_{k,y}^2+\varepsilon$
检查近似的边界
当 Gaussian 的深度厚度相对中心深度很大，或靠近相机时，J 在整团内变化明显，一阶椭圆近似会变差。加 0.3 是过滤策略，不是透视几何必然出现的常数。
$\mathbf d_1=(10,0)^T,\quad\mathbf d_2=(0,20)^T,\quad\mathbf d_3=(-3,-6)^T$

代入数字，亲手核对

用前面的 J 和 Σc=diag(1,4,9)，分别投影长度为 1、2、3 的轴，再累加。

C=\begin{bmatrix}100+9+0.3&18\\18&400+36+0.3\end{bmatrix}=\begin{bmatrix}109.3&18\\18&436.3\end{bmatrix}

这一步在代码里在哪里？

shader 保留 cov_a、cov_b、cov_c 三个独立元素，没有先分配一个 3×3 矩阵。

dk ↔ d0,d1,d2；ε ↔ +0.3

src/splat.wgsl · L94–101 ↗

  let d0 = projected_axis(center_clip, axis0);
  let d1 = projected_axis(center_clip, axis1);
  let d2 = projected_axis(center_clip, axis2);

  // Project J W Sigma W^T J^T by summing the projected covariance axes.
  let cov_a = dot(vec3<f32>(d0.x, d1.x, d2.x), vec3<f32>(d0.x, d1.x, d2.x)) + 0.3;
  let cov_b = d0.x * d0.y + d1.x * d1.y + d2.x * d2.y;
  let cov_c = dot(vec3<f32>(d0.y, d1.y, d2.y), vec3<f32>(d0.y, d1.y, d2.y)) + 0.3;

停下来想一想：为什么主点 cx、cy 不出现在 C 里？

它们只平移投影中心，求导后为 0，不改变局部形状。

带着这条公式，查看完整像素算例 →

(2.6)与 WGSL 对齐的齐次坐标导数

D (\frac{𝐡_{xy}}{𝐡_{w}}) (δ 𝐡) = \frac{δ 𝐡_{xy} 𝐡_{w} - 𝐡_{xy} δ 𝐡_{w}}{{𝐡_{w}}^{2}}

h=PV[μ,1]ᵀ，δh=PV[aₖ,0]ᵀ。将导数按 viewport 尺寸换成像素，再对三条投影轴做外积求和。这个形式把视图坐标的符号差异包含在 P、V 中。

从零推导5 步 · 数字算例 · 源码 · 自测

分式通分 → 方向导数

WGSL 没写 fx/z，怎样确认它算的仍是同一个 Jacobian？

用到的小工具E 导数 B 矩阵

给中心和方向加上不同的末位
中心 h=PV[μ,1]；轴方向 δh=PV[a,0]。沿轴走一个小量 η 后，clip 坐标就是 h+ηδh，因为矩阵乘法是线性的。
$h(\eta)=h+\eta\delta h,\qquad n_x(\eta)=\frac{h_x+\eta\delta h_x}{h_w+\eta\delta h_w}$
先相减再通分
分子里的 hx hw 抵消；剩下两项都含 η，和差商分母中的 η 约掉。
$\frac{n_x(\eta)-n_x(0)}{\eta}=\frac{(h_x+\eta\delta h_x)h_w-h_x(h_w+\eta\delta h_w)}{\eta h_w(h_w+\eta\delta h_w)}=\frac{\delta h_xh_w-h_x\delta h_w}{h_w(h_w+\eta\delta h_w)}$
把 η 缩到 0，y 同理
这里需要 hw 非零。把 x、y 两个结果拼成向量，就是 shader 的分子与分母。
$D(h_{xy}/h_w)[\delta h]=\frac{\delta h_{xy}h_w-h_{xy}\delta h_w}{h_w^2}$
从 NDC 换到像素
NDC 横向从 −1 到 +1 跨 2 个单位，对应宽度 W，故乘 W/2。WebGPU 屏幕原点在左上，实际向下的像素 y 要乘 −H/2；shader 采用向上的内部屏幕坐标，乘 +H/2。两种坐标用 F 翻转转换，整个椭圆和偏移一起翻转，二次型不变。
$\mathbf d_{\uparrow}=\operatorname{diag}(W/2,H/2)\delta\mathbf n,\quad F=\operatorname{diag}(1,-1),\quad C_{\downarrow}=FC_{\uparrow}F^T,\quad \mathbf d_{\downarrow}=F\mathbf d_{\uparrow}$
顶点回到 clip 坐标时为什么又乘 w
想让透视除法后的 NDC 多 δn，就在 clip.xy 里加 δn·hw；除以 hw 后恰好剩 δn。中心 w 对四个角相同，local 的透视插值在这张 billboard 内也就是线性插值。
$\frac{h_{xy}+\delta n\,h_w}{h_w}=\frac{h_{xy}}{h_w}+\delta n$

代入数字，亲手核对

hx=2、hw=4、δhx=1、δhw=0.5、W=800。

\delta n_x=\frac{1\cdot4-2\cdot0.5}{16}=\frac3{16},\qquad d_x=400\cdot\frac3{16}=75\ {\rm pixel}

这一步在代码里在哪里？

公式的每一项在 projected_axis 中都有同名值。max(w²,1e−8) 是数值保护，中心背面的剔除另由 CPU 完成。

h ↔ center_clip；δh ↔ axis_clip

src/splat.wgsl · L41–49 ↗

fn projected_axis(center_clip: vec4<f32>, axis: vec3<f32>) -> vec2<f32> {
  // Exact first-order derivative of ndc = clip.xy / clip.w along a world-space
  // covariance axis. This is the homogeneous form of the perspective Jacobian.
  let axis_clip = camera.view_proj * vec4<f32>(axis, 0.0);
  let w2 = max(center_clip.w * center_clip.w, 0.00000001);
  let ndc_derivative = (
    axis_clip.xy * center_clip.w - center_clip.xy * axis_clip.w
  ) / w2;
  return ndc_derivative * camera.viewport.xy * 0.5;

pixel_offset 先变成 ndc_offset，再乘 center_clip.w 加回 clip.xy。

δn ↔ ndc_offset

src/splat.wgsl · L116–126 ↗

  let pixel_offset = 3.0 * (
    major_axis * sigma_major * corner.x
    + minor_axis * sigma_minor * corner.y
  );
  let ndc_offset = pixel_offset * vec2<f32>(2.0 * camera.viewport.z, 2.0 * camera.viewport.w);

  var output: VertexOutput;
  output.position = center_clip + vec4<f32>(ndc_offset * center_clip.w, 0.0, 0.0);
  output.local = corner * 3.0;
  output.color_opacity = vec4<f32>(clamp(splat.r4.yzw, vec3<f32>(0.0), vec3<f32>(1.0)), opacity);
  return output;

停下来想一想：轴的齐次末位为什么是 0？

它是方向，不是点。末位为 0 才不会把相机平移重复加进轴向量。

带着这条公式，查看完整像素算例 →

手算检查

fx=fy=100，中心为 (1,2,10)，Σc=diag(1,4,9)：

𝐉 = [\begin{matrix} 10 & 0 & -1 \\ 0 & 10 & -2 \end{matrix}]

𝐂 = [\begin{matrix} 109.3 & 18 \\ 18 & 436.3 \end{matrix}]

非零的交叉项来自离轴投影。即使三维协方差是对角矩阵，投影椭圆也可能发生旋转。

原始资料3DGS · Kerbl et al., 2023 ↗本次部署源码 · splat.wgsl ↗COLMAP · 相机坐标约定 ↗

03

THE SCREEN FOOTPRINT

从协方差到一个像素的权重。

必会

协方差告诉我们形状；逆协方差定义到中心的“加权距离”。

(3.1)2×2 协方差的逆

\begin{matrix} 𝐂 & = & [\begin{matrix} a & b \\ b & c \end{matrix}] \\ 𝐐 & = & 𝐂^{- 1} = \frac{1}{a c - b^{2}} [\begin{matrix} c & - b \\ - b & a \end{matrix}] \end{matrix}

正定要求 a>0、ac−b²>0。接近奇异时，逆矩阵会放大误差；需要与 renderer 一致的滤波和数值下限。

从零推导4 步 · 数字算例 · 源码 · 自测

解两组二元一次方程

不背口诀，也能得到 2×2 逆矩阵。

用到的小工具B 矩阵

把逆矩阵的第一列设成未知数
C 乘逆矩阵的第一列应得到 (1,0)。这就是两条一次方程。
$\begin{bmatrix}a&b\\b&c\end{bmatrix}\begin{bmatrix}u\\v\end{bmatrix}=\begin{bmatrix}1\\0\end{bmatrix}\quad\Longrightarrow\quad au+bv=1,\quad bu+cv=0$
消元解出第一列
第一式乘 c，第二式乘 b，相减得到 (ac−b²)u=c。类似消掉 u 可解出 v。令 D=ac−b²。
$D=ac-b^2,\qquad u=\frac cD,\qquad v=-\frac bD$
解第二列，再拼回去
右侧换成 (0,1)，同样消元得到 (−b/D,a/D)。现在把两列排在一起，就得到 Q。
$Q=\frac1D\begin{bmatrix}c&-b\\-b&a\end{bmatrix},\qquad CQ=\frac1D\begin{bmatrix}ac-b^2&0\\0&ac-b^2\end{bmatrix}=I$
为什么要求正定
把 ax²+2bxy+cy² 配方；要让任何非零方向都为正，需要 a>0 且 c−b²/a>0，也就是 D>0。
$ax^2+2bxy+cy^2=a\left(x+\frac ba y\right)^2+\frac Da y^2$

代入数字，亲手核对

C=[[4,1],[1,2]]，D=8−1=7。

Q=\frac17\begin{bmatrix}2&-1\\-1&4\end{bmatrix}

这一步在代码里在哪里？

PyTorch 使用 linalg.inv 显式求逆；GPU 通过特征轴局部坐标得到相同距离，下一节解释这条连接。

Q ↔ inverse；D ↔ det(C)

scripts/train_toy_stg.py · L36–44 ↗

    ys, xs = torch.meshgrid(torch.arange(height, device=device)+0.5,
                           torch.arange(width, device=device)+0.5, indexing="ij")
    delta = torch.stack((xs, ys), -1)[None, None] - uv[:, :, None, None]
    inverse = torch.linalg.inv(covariance)
    power = -0.5 * torch.einsum("bnhwi,bnij,bnhwj->bnhw", delta, inverse, delta)
    temporal = torch.exp(-((times[:, None]-0.5)/0.55).square())
    alpha = (torch.sigmoid(opacity_logits)[None, :, None, None]
             * temporal[:, :, None, None] * torch.exp(power)).clamp(max=0.99)
    alpha = torch.where(alpha >= 1/255, alpha, 0)

停下来想一想：C=diag(4,9) 的逆是什么？

Q=diag(1/4,1/9)，不是 diag(−4,−9)。

带着这条公式，查看完整像素算例 →

(3.2)Mahalanobis 距离与像素 alpha

\begin{matrix} 𝐝 & = & 𝐩 - 𝛍_{2D} \\ r^{2} & = & 𝐝^{T} 𝐐 𝐝 \\ α (𝐩, t) & = & o (t) e^{- \frac{r^{2}}{2}} \end{matrix}

空间核指数含 1/2。当前 fragment shader 另做 alpha≤0.99 和 alpha<1/255 的丢弃；以下训练推导先忽略这些分段操作。

从零推导5 步 · 数字算例 · 源码 · 自测

椭圆换坐标的等价式

CPU 的 dᵀQd 怎样变成 WGSL 的 dot(local,local)？

用到的小工具A 向量 B 矩阵 C 指数

把像素改写为离中心的偏移
像素位置 p 与投影中心必须用同一个屏幕坐标系。离中心远近不能直接用世界距离。
$\mathbf d=\mathbf p-\mu_{2D},\qquad r^2=\mathbf d^TQ\mathbf d$
先把矩阵乘法完全展开
Q 对称，所以中间两项相同，可以合并为 2Q12 dxdy；不是只剩一项。
$r^2=Q_{11}d_x^2+2Q_{12}d_xd_y+Q_{22}d_y^2$
在椭圆自己的两条主轴上计量
C=U diag(λ₁,λ₂)Uᵀ，U 的列是互相垂直的单位轴。先用 Uᵀ 转到主轴坐标，再除以标准差 √λ；得到 local。
$\mathbf l=\operatorname{diag}(\lambda_1^{-1/2},\lambda_2^{-1/2})U^T\mathbf d,\qquad \mathbf l^T\mathbf l=\mathbf d^TU\operatorname{diag}(\lambda_1^{-1},\lambda_2^{-1})U^T\mathbf d=\mathbf d^TC^{-1}\mathbf d$
让两个三角形替我们插值 local
顶点放在主轴各 ±3σ 的矩形角上，并赋 local=corner×3。光栅化后，每个片元收到自己的连续局部坐标，所以直接 dot(local,local) 就能得到 r²。
$\mathbf d=U\operatorname{diag}(\sqrt{\lambda_1},\sqrt{\lambda_2})\mathbf l,\qquad \alpha=o(t)e^{-\|\mathbf l\|^2/2}$
再把实际分支写清楚
shader 把 alpha 截到 0.99，低于 1/255 时丢弃；矩形外没有片元。只有未触发轴长度限幅时，这个 local 才精确对应原始 C。
$\alpha_{\rm cap}=\min(0.99,\alpha),\qquad \alpha_{\rm cap}<1/255\Rightarrow\text{discard}$

代入数字，亲手核对

C=diag(4,1)，像素偏移 (2,1)，o=0.8。

\mathbf l=(2/2,1/1)=(1,1),\quad r^2=2,\quad \alpha=0.8e^{-1}\approx0.294304

这一步在代码里在哪里？

这里没有 inv(C)，因为逆变换已经隐含在顶点给出的 local 和片元插值中。

local ↔ l；power ↔ −r²/2

src/splat.wgsl · L130–136 ↗

fn fs_main(input: VertexOutput) -> @location(0) vec4<f32> {
  let power = -0.5 * dot(input.local, input.local);
  let alpha = min(0.99, input.color_opacity.a * exp(power));
  if (alpha < 0.0039215686) {
    discard;
  }
  return vec4<f32>(input.color_opacity.rgb * alpha, alpha);

major_axis、minor_axis 构成 U，sigma_major/minor 是 √λ。

corner*3 ↔ 顶点处 l

src/splat.wgsl · L116–126 ↗

  let pixel_offset = 3.0 * (
    major_axis * sigma_major * corner.x
    + minor_axis * sigma_minor * corner.y
  );
  let ndc_offset = pixel_offset * vec2<f32>(2.0 * camera.viewport.z, 2.0 * camera.viewport.w);

  var output: VertexOutput;
  output.position = center_clip + vec4<f32>(ndc_offset * center_clip.w, 0.0, 0.0);
  output.local = corner * 3.0;
  output.color_opacity = vec4<f32>(clamp(splat.r4.yzw, vec3<f32>(0.0), vec3<f32>(1.0)), opacity);
  return output;

停下来想一想：local=(0,0) 时 alpha 一定等于 1 吗？

不一定。空间核为 1，但仍要乘时间权重和基础 opacity，并受 0.99 上限限制。

带着这条公式，查看完整像素算例 →

(3.3)椭圆主轴与标准差

\begin{matrix} λ_{\pm} & = & \frac{a + c \pm \sqrt{{(a - c)}^{2} + 4 b^{2}}}{2} \\ r_{\pm} & = & k \sqrt{λ_{\pm}} \end{matrix}

特征向量是主轴方向，√λ 是对应方向的像素标准差。k=3 给出常用的 3σ 覆盖；当前实现对 sigma 额外设了 192px 上限。

从零推导5 步 · 数字算例 · 源码 · 自测

二次方程 + 换坐标

从 a,b,c 求出椭圆的长短轴。

用到的小工具G 配方 B 矩阵

寻找变换后方向不变的向量
若 Cv=λv，v 的方向不变，只拉伸 λ 倍；这就是特征向量和特征值的定义。把右侧移到左侧。
$(C-\lambda I)\mathbf v=0,\qquad C-\lambda I=\begin{bmatrix}a-\lambda&b\\b&c-\lambda\end{bmatrix}$
什么时候有非零的解
两条方程必须不独立，否则只能 v=0；二维时条件就是行列式为 0。展开后得到普通二次方程。
$(a-\lambda)(c-\lambda)-b^2=0\quad\Longrightarrow\quad\lambda^2-(a+c)\lambda+(ac-b^2)=0$
使用配方法求根
把上一式的 s=a+c、p=ac−b² 代入基础课求根公式。判别式化简为 (a−c)²+4b²。
$\lambda_\pm=\frac{a+c\pm\sqrt{(a+c)^2-4(ac-b^2)}}2=\frac{a+c\pm\sqrt{(a-c)^2+4b^2}}2$
解出方向并归一化
第一行是 (a−λ)vx+bvy=0，取 (vx,vy)=(b,λ−a) 恰好满足。除以长度得到单位轴；另一轴垂直于它。b=0 时按 a、c 的大小选坐标轴。
$\mathbf u_+=\frac{(b,\lambda_+-a)^T}{\sqrt{b^2+(\lambda_+-a)^2}},\qquad \mathbf u_-=(-u_{+,y},u_{+,x})^T$
方差开平方才是长度
在特征轴坐标中二次型是 ξ₁²/λ₁+ξ₂²/λ₂。令它等于 k²，就得到椭圆半轴长度 k√λ。
$\frac{\xi_1^2}{\lambda_1}+\frac{\xi_2^2}{\lambda_2}=k^2\qquad\Rightarrow\qquad r_1=k\sqrt{\lambda_1},\quad r_2=k\sqrt{\lambda_2}$

代入数字，亲手核对

C=[[4,1],[1,2]]；λ±=3±√2。

\lambda_+\approx4.414214,\quad\lambda_-\approx1.585786,\quad (3\sigma_+,3\sigma_-)\approx(6.303009,3.777840)

这一步在代码里在哪里？

discriminant、lambda_major/minor 按这一求根式计算。GPU 还把 λ 下限设为 0.1、σ 上限设为 192 像素，属于实现保护。

λ± ↔ lambda_major/minor；u± ↔ major_axis/minor_axis

src/splat.wgsl · L102–115 ↗

  let discriminant = sqrt(max(0.0, (cov_a - cov_c) * (cov_a - cov_c) + 4.0 * cov_b * cov_b));
  let lambda_major = max(0.1, 0.5 * (cov_a + cov_c + discriminant));
  let lambda_minor = max(0.1, 0.5 * (cov_a + cov_c - discriminant));

  var major_axis = vec2<f32>(1.0, 0.0);
  if (abs(cov_b) > 0.00001) {
    major_axis = normalize(vec2<f32>(cov_b, lambda_major - cov_a));
  } else if (cov_c > cov_a) {
    major_axis = vec2<f32>(0.0, 1.0);
  }
  let minor_axis = vec2<f32>(-major_axis.y, major_axis.x);
  let corner = quad_corner(input.vertex_index);
  let sigma_major = min(sqrt(lambda_major), 192.0);
  let sigma_minor = min(sqrt(lambda_minor), 192.0);

停下来想一想：C=diag(9,4)，3σ 椭圆半轴多长？

9 和 6 像素：先开平方得到 3、2，再乘 3。

带着这条公式，查看完整像素算例 →

(3.4)由 alpha 阈值推覆盖范围

r^{2} \leq 2 \log (\frac{o (t)}{α_{min}})

仅在 o(t)≥αmin 时有非空范围。固定 3σ quad 与阈值支持集不完全相同；“Gaussian 理论无限支持”与有限 GPU footprint 需要区分。

从零推导4 步 · 数字算例 · 源码 · 自测

解指数不等式

alpha 小到可以忽略时，到底能截掉多远？

用到的小工具C 指数

写出保留像素的条件
设阈值 αmin>0，先忽略 clamp。o 是当前时间的峰值 opacity，r² 非负。
$o\,e^{-r^2/2}\ge\alpha_{\mathrm{min}}$
先检查有没有解
核的最大值为 1，因此若 o<αmin，连中心都达不到阈值，整个 Gaussian 可跳过。否则除以正数 o。
$o\ge\alpha_{\mathrm{min}}>0\quad\Longrightarrow\quad e^{-r^2/2}\ge\frac{\alpha_{\mathrm{min}}}o$
取对数，再乘负数
ln 单调递增，不等号不变；乘 −2 时不等号要反向。
$-\frac{r^2}{2}\ge\ln\frac{\alpha_{\mathrm{min}}}o\quad\Longrightarrow\quad r^2\le2\ln\frac o{\alpha_{\mathrm{min}}}$
区分理论支持域和实际矩形
精确阈值边界在 local 空间是半径 rmax 的圆，回到屏幕是椭圆。当前 shader 画的是每个主轴 ±3σ 的矩形，再在片元里按 alpha 丢弃，因此两者并不相等。
$o=0.8,\ \alpha_{\mathrm{min}}=1/255:\quad r_{\mathrm{max}}=\sqrt{2\ln204}\approx3.261303$

代入数字，亲手核对

轴向 local=(3.1,0) 的 alpha 仍超过阈值，但超出当前 ±3 的 quad，没有片元；这显示真实的截断近似。

0.8e^{-3.1^2/2}\approx0.006551>1/255\approx0.003922

这一步在代码里在哪里？

corner×3 确定实际几何覆盖范围；它不是上述阈值公式的计算结果。

±3σ ↔ pixel_offset 中的 3.0

src/splat.wgsl · L116–126 ↗

  let pixel_offset = 3.0 * (
    major_axis * sigma_major * corner.x
    + minor_axis * sigma_minor * corner.y
  );
  let ndc_offset = pixel_offset * vec2<f32>(2.0 * camera.viewport.z, 2.0 * camera.viewport.w);

  var output: VertexOutput;
  output.position = center_clip + vec4<f32>(ndc_offset * center_clip.w, 0.0, 0.0);
  output.local = corner * 3.0;
  output.color_opacity = vec4<f32>(clamp(splat.r4.yzw, vec3<f32>(0.0), vec3<f32>(1.0)), opacity);
  return output;

alpha 低于阈值后 discard，两个判断共同决定最终覆盖。

αmin ↔ 0.0039215686

src/splat.wgsl · L130–136 ↗

fn fs_main(input: VertexOutput) -> @location(0) vec4<f32> {
  let power = -0.5 * dot(input.local, input.local);
  let alpha = min(0.99, input.color_opacity.a * exp(power));
  if (alpha < 0.0039215686) {
    discard;
  }
  return vec4<f32>(input.color_opacity.rgb * alpha, alpha);

停下来想一想：o 正好等于 αmin 时，数学上哪里还能保留？

只有 r=0 的中心；实际像素网格可能没有采样点正好落在中心。

带着这条公式，查看完整像素算例 →

原始资料3DGS · Kerbl et al., 2023 ↗本次部署源码 · splat.wgsl ↗

04

ORDER MATTERS

为什么透明度必须按顺序合成。

必推

下式按近→远排列。颜色与 alpha 的贡献受前面所有层的透射率影响。

(4.1)前向透射率

T_{i} = \prod_{j < i}^{} (1 - α_{j}), T_{1} = 1

从零推导4 步 · 数字算例 · 源码 · 自测

剩余比例的连乘

光到达第 i 层时，还剩下多少？

用到的小工具A 向量

把每层理解为覆盖一部分的薄片
alpha 为 0.2，就有 20% 被这层覆盖，80% 继续往后。先用这个透明合成模型建立直觉，它不要求场景真的是玻璃。
$T_1=1,\qquad T_2=1-\alpha_1$
下一层只作用于上一层剩下的部分
若第二层也覆盖 20%，它覆盖的不是原始光的 20%，而是剩余 80% 的 20%。
$T_3=(1-\alpha_1)(1-\alpha_2)$
重复这一规则
到达第 i 层之前已经穿过 i−1 层；∏ 是把这些剩余比例全部相乘。
$T_{i+1}=T_i(1-\alpha_i),\qquad T_i=\prod_{j<i}(1-\alpha_j)$
解释训练代码为什么先补一个 1
cumprod 得到前缀连乘。需要的是当前层之前的连乘，因此数组前面补 1，再丢掉最后一项。
$[\alpha_1,\alpha_2]=[0.2,0.5]\Rightarrow[T_1,T_2,T_3]=[1,0.8,0.4]$

代入数字，亲手核对

两层 alpha 都是 0.5，第二层得到的光为 0.5，背景得到 0.25。

T_2=0.5,\qquad T_3=0.5\cdot0.5=0.25

这一步在代码里在哪里？

先近后远排列 alpha，沿 Gaussian 维度 cumprod；[:, :-1] 取到当前层之前的透射率。

Ti ↔ transmittance；dim=1 ↔ Gaussian 维

scripts/train_toy_stg.py · L45–50 ↗

    order = z.argsort(dim=1)  # Discrete ordering; no gradient through visibility changes.
    alpha = alpha.gather(1, order[:, :, None, None].expand_as(alpha))
    color = torch.sigmoid(rgb_logits)[order]
    transmittance = torch.cumprod(torch.cat((torch.ones_like(alpha[:, :1]), 1-alpha), 1), 1)[:, :-1]
    return (transmittance[:, :, :, :, None] * alpha[:, :, :, :, None]
            * color[:, :, None, None, :]).sum(1)

停下来想一想：若第一层 alpha=1，后面的 Ti 是多少？

全部为 0；当前 shader 把单层 alpha 限到 0.99，因此实际不会由单层直接得到精确的 0。

带着这条公式，查看完整像素算例 →

(4.2)最终像素颜色，包含背景

𝐂_{pixel} = \sum_{i = 1}^{N} T_{i} α_{i} 𝐜_{i} + T_{N + 1} 𝐜_{bg}

从零推导4 步 · 数字算例 · 源码 · 自测

加权相加 + 望远镜求和

每层与背景怎样一起组成最终像素？

用到的小工具A 向量

先看一层和背景
当前层占 alpha，其余部分透出背景。RGB 三个通道分别做同样的运算。
$\mathbf C=\alpha_1\mathbf c_1+(1-\alpha_1)\mathbf c_{\rm bg}$
把背景换成后面的一整幅合成结果
两层时先用第二层和背景合成 B1，再把 B1 放到第一层后面。展开括号。
$\mathbf C=\alpha_1\mathbf c_1+(1-\alpha_1)\alpha_2\mathbf c_2+(1-\alpha_1)(1-\alpha_2)\mathbf c_{\rm bg}$
推广到 N 层
第 i 层之前必须穿过前 i−1 层，所以颜色权重为 Ti αi；背景的权重是穿过全部层的 TN+1。
$\mathbf C=\sum_{i=1}^N T_i\alpha_i\mathbf c_i+T_{N+1}\mathbf c_{\rm bg}$
用权重和检查有没有漏项
Ti αi=Ti−Ti+1；求和时中间项抵消，只剩 1−TN+1。加上背景正好为 1。
$\sum_i T_i\alpha_i+T_{N+1}=\sum_i(T_i-T_{i+1})+T_{N+1}=T_1=1$

代入数字，亲手核对

红在前、蓝在后，两层 alpha=0.5，黑色背景。

\mathbf C=0.5(1,0,0)+0.25(0,0,1)+0.25(0,0,0)=(0.5,0,0.25)

这一步在代码里在哪里？

代码对 T*alpha*color 沿 Gaussian 维求和，背景为黑，所以省略了数值为 0 的背景项。

C ↔ return；N ↔ 张量 dim=1

scripts/train_toy_stg.py · L45–50 ↗

    order = z.argsort(dim=1)  # Discrete ordering; no gradient through visibility changes.
    alpha = alpha.gather(1, order[:, :, None, None].expand_as(alpha))
    color = torch.sigmoid(rgb_logits)[order]
    transmittance = torch.cumprod(torch.cat((torch.ones_like(alpha[:, :1]), 1-alpha), 1), 1)[:, :-1]
    return (transmittance[:, :, :, :, None] * alpha[:, :, :, :, None]
            * color[:, :, None, None, :]).sum(1)

停下来想一想：为什么权重不一定加起来等于 1？

如果只看 Gaussian 权重，和为 1−剩余透射率；把背景权重加进去才等于 1。

带着这条公式，查看完整像素算例 →

(4.3)远→近的预乘 over 操作

{\tilde{𝐜}}_{out} = α_{s} 𝐜_{s} + (1 - α_{s}) {\tilde{𝐜}}_{dst}

graphics pipeline 使用远→近顺序时，shader 输出 RGB×alpha，再配置 premultiplied blending。重复乘 SrcAlpha 会让颜色额外变暗。

从零推导4 步 · 数字算例 · 源码 · 自测

同一合成式的 GPU 写法

为什么 fragment 输出 RGB×alpha，混合状态却用 One？

用到的小工具A 向量

定义预乘颜色
把颜色与它的覆盖率先乘好，记作 c̃。这样透明黑和普通黑在存储时也有一致的含义。
$\widetilde{\mathbf c}=\alpha\mathbf c$
把后面已经画好的结果当 destination
远处先画，当前较近的一层为 source。它留下自己的预乘颜色，并遮掉 destination 的 alpha 比例。
$\widetilde{\mathbf c}_{\rm out}=\widetilde{\mathbf c}_s+(1-\alpha_s)\widetilde{\mathbf c}_{\rm dst}$
映射到固定混合硬件的两个乘数
硬件计算 src×srcFactor + dst×dstFactor。于是 srcFactor=One，dstFactor=OneMinusSrcAlpha。alpha 通道也按 over 规则累积。
$\alpha_{\rm out}=\alpha_s+(1-\alpha_s)\alpha_{\rm dst}$
与近到远公式核对
先画远处蓝，再画近处红；展开正好得到上一式的红 0.5、蓝 0.25。顺序不同，代数结果一致。
$\widetilde{\mathbf c}_{\rm blue}=0.5(0,0,1),\quad \widetilde{\mathbf c}_{\rm out}=0.5(1,0,0)+0.5\,\widetilde{\mathbf c}_{\rm blue}$

代入数字，亲手核对

若 shader 已输出红色 (0.5,0,0)，却又配置 SrcAlpha=0.5，会错误地再乘一次。

0.5\times0.5=0.25\ne0.5

这一步在代码里在哪里？

返回值 vec4(rgb*alpha,alpha) 明确是预乘。

c̃s ↔ input.color_opacity.rgb*alpha

src/splat.wgsl · L130–136 ↗

fn fs_main(input: VertexOutput) -> @location(0) vec4<f32> {
  let power = -0.5 * dot(input.local, input.local);
  let alpha = min(0.99, input.color_opacity.a * exp(power));
  if (alpha < 0.0039215686) {
    discard;
  }
  return vec4<f32>(input.color_opacity.rgb * alpha, alpha);

PREMULTIPLIED_ALPHA_BLENDING 与这份输出匹配；StgPass 按远到近排序，透明 splat 不写深度。

over ↔ blend state；顺序 ↔ radix depth sort

src/stg_pass.rs · L189–204 ↗

                    blend: Some(wgpu::BlendState::PREMULTIPLIED_ALPHA_BLENDING),
                    write_mask: wgpu::ColorWrites::ALL,
                })],
            }),
            primitive: wgpu::PrimitiveState {
                topology: if gpu_variant >= 2 {
                    wgpu::PrimitiveTopology::TriangleStrip
                } else {
                    wgpu::PrimitiveTopology::TriangleList
                },
                ..Default::default()
            },
            depth_stencil: depth.map(|(format, depth_compare)| wgpu::DepthStencilState {
                format,
                depth_write_enabled: Some(false),
                depth_compare: Some(depth_compare),

停下来想一想：预乘颜色仍要乘 Ti 吗？

逐层 GPU over 已经隐式完成遮挡；若在手写近到远累积器里，则仍需乘到达该层的 Ti。

带着这条公式，查看完整像素算例 →

红在前(0.50, 0, 0.25)

蓝在前(0.25, 0, 0.50)

红、蓝两层 alpha 都为 0.5，背景为黑。交换顺序改变结果，所以普通 atomic-add 不能代替有序透明合成。

(4.4)与体渲染的联系

\begin{matrix} T (s) & = & e^{- \int_{s_{0}}^{s} σ (u) d u} \\ α_{i} & = & 1 - e^{- σ_{i} Δ s_{i}} \end{matrix}

第二式来自分段恒定密度的射线积分，帮助理解 transmittance。3DGS 通常直接学习 opacity 并评估投影核，不能把它说成对每个三维 Gaussian 都做了精确物理射线积分。

从零推导4 步 · 数字算例 · 源码 · 自测

薄层极限，非精确 GS 射线积分

连续密度为什么会产生指数透射率？

用到的小工具H 积分 C 指数

把一段路切成 N 个薄层
密度 σ 表示单位长度上的消光强度；薄层厚 Δs 足够小时，透过比例约为 1−σΔs。把总长度 l 平分。
$T_N=\left(1-\sigma\frac lN\right)^N$
让薄层越来越薄
基础极限 (1+a/N)ᴺ→eᵃ，取 a=−σl。常密度下就得到指数透过率；被覆盖比例是 1−T。
$T=e^{-\sigma l},\qquad \alpha=1-e^{-\sigma l}$
密度沿路径变化时分别连乘
取对数把乘法变加法，利用 ln(1−ε)≈−ε，再把细分求和变成积分。
$\ln T\approx-\sum_j\sigma(s_j)\Delta s\ \longrightarrow\ -\int_{s_0}^s\sigma(u)\,du$
回到离散渲染的每一段
每段内部假定密度固定，就用该段长度 Δsi 得到 alpha。相邻段按剩余比例相乘，因此和前面的合成模型接上。
$T(s)=\exp\left[-\int_{s_0}^s\sigma(u)\,du\right],\qquad \alpha_i=1-e^{-\sigma_i\Delta s_i}$

代入数字，亲手核对

σ=2 /米，一段长 0.5 米，乘积 σΔs 没有单位。

T=e^{-1}\approx0.367879,\qquad \alpha=1-e^{-1}\approx0.632121

这一步在代码里在哪里？

对照源码可见，当前 renderer 直接计算 opacity*投影 Gaussian 核，没有沿每条光线积分 σ。这里只解释透明合成与体渲染的关系。

本式 σ、Δs：当前 STG-Lite shader 没有对应字段

src/splat.wgsl · L130–136 ↗

fn fs_main(input: VertexOutput) -> @location(0) vec4<f32> {
  let power = -0.5 * dot(input.local, input.local);
  let alpha = min(0.99, input.color_opacity.a * exp(power));
  if (alpha < 0.0039215686) {
    discard;
  }
  return vec4<f32>(input.color_opacity.rgb * alpha, alpha);

停下来想一想：若把一段恒密度路程平分两半，总透过率会变吗？

不会。e^(−σl/2)×e^(−σl/2)=e^(−σl)。

带着这条公式，查看完整像素算例 →

原始资料3DGS · Kerbl et al., 2023 ↗NeRF · 体渲染背景 ↗

05

THIS VIEWER · STG-LITE

让位置、旋转和可见性随时间变化。

项目核心

下面采用本项目 PLY/WGSL 的实际参数化。归一化时间 t 不自动等于物理秒。

(5.1)物理秒与模型时间

\begin{matrix} t & = & \frac{s - s_{0}}{D} \\ Δ t & = & t - τ \end{matrix}

D 是有来源的时间尺度。例如 loader 使用 frame/N 时，通常需与 N/fps 对齐，而不是擅自换成最后一帧时间。当前研究模型 PLY 不携带 fps/时长；已追溯官方 50 帧训练配置，Sear 播放器采用 D=50/30 秒。

从零推导4 步 · 数字算例 · 源码 · 自测

单位换算

视频的秒、帧号、模型时间如何对齐？

用到的小工具E 导数

先确定时间起点和尺度
物理时刻 s 减去起点 s0，得到经过的秒数；除以 D 秒，得到无单位的模型时间 t。D 必须来自训练约定。
$t=\frac{s-s_0}{D},\qquad s=s_0+Dt$
再对齐每个 Gaussian 自己的中心时刻
τ 是该 Gaussian 最活跃的模型时间。运动系数围绕 τ 展开，所以实际使用 Δt=t−τ。
$\Delta t=\frac{s-s_0}{D}-\tau$
帧号需要匹配分母约定
若训练器定义 t=j/N，帧率为 fps，则 s−s0=j/fps，与 D=N/fps 对齐；最后一帧 t=(N−1)/N，不等于 1。若定义 j/(N−1)，D 则不同。
$\frac{j/fps}{N/fps}=\frac jN$
沿代码找到时间只被换算一次
当前 host 用已核实的片段时长把秒换算为 t；兼容旧接口时先传入 t×10，Rust 再除以 10。这不表示片段长 10 秒。原生 StgPass 直接接收 normalized_time。
$s=2,\ s_0=0,\ D=10,\ \tau=0.5\Rightarrow t=0.2,\ \Delta t=-0.3$

代入数字，亲手核对

若有 100 帧、25fps，并采用 j/N，第 50 帧的模型时间为 0.5。

D=100/25=4\ {\rm s},\qquad t=(50/25)/4=0.5

这一步在代码里在哪里？

当前 LOOP_SECONDS=10，用 rem_euclid 后除以它；这只描述播放器速度。

s ↔ time_seconds；D ↔ LOOP_SECONDS

src/lib.rs · L543–553 ↗

        let normalized_time = time_seconds.rem_euclid(LOOP_SECONDS) / LOOP_SECONDS;

        let presented_key = (normalized_time, view, projection, viewport, width, height);
        if self.render_on_demand && self.last_presented == Some(presented_key) {
            self.telemetry.prepare_ms = prepare_started.elapsed().as_secs_f64() * 1000.0;
            self.telemetry.sort_ms = 0.0;
            self.telemetry.upload_bytes = 0;
            return Ok(());
        }
        self.inspection_frame = None;
        let variant = 9;

原生接入显式传归一化时间；prepare 检查 0..1。

t ↔ normalized_time → camera.scene.x

src/stg_pass.rs · L247–250 ↗

        let view = Mat4::from_cols_array_2d(&view_columns);
        let projection = Mat4::from_cols_array_2d(&projection_columns);
        let view_projection = projection * view;
        let [width, height] = viewport.map(|v| v as f32);

停下来想一想：同一个模型以两倍播放速度播放，运动系数需要改吗？

不需要；改变输入 t 随真实时间推进的速度即可。若重新定义模型时间坐标本身，则系数必须换算。

带着这条公式，查看完整像素算例 →

(5.2)STG-Lite 的时间核

\begin{matrix} ρ & = & e^{ℓ_{τ}} \\ w (t) & = & e^{- {(\frac{Δ t}{ρ})}^{2}} \\ o (t) & = & sigmoid (β) w (t) \end{matrix}

注意指数没有 1/2。若用标准 Gaussian 的标准差 σt 表示同一核，则 σt=ρ/√2。论文写 exp(−sτ·Δt²)，这里 sτ=ρ⁻²；代码对 ρ 还有最小值保护。

从零推导4 步 · 数字算例 · 源码 · 自测

时间核的模型选择

为什么有的时间 Gaussian 指数没有 1/2？

用到的小工具C 指数

定义到活跃时刻的相对距离
Δt=t−τ；ρ 是宽度。先除以宽度再平方，时间单位约掉。STG-Lite 采用 exp(−距离²) 作为可见性包络。
$r_t=\frac{t-\tau}{\rho},\qquad w(t)=e^{-r_t^2}$
用 log-width 保证宽度为正
PLY 保存 ℓτ，不是 ρ 本身。与空间 log-scale 类似，先 exp 解码。官方论文用另一种系数 sτ 乘在平方前，二者可直接换算。
$\rho=e^{\ell_\tau},\qquad w=e^{-\Delta t^2/\rho^2}=e^{-s_\tau\Delta t^2},\qquad s_\tau=\rho^{-2}$
和标准 Gaussian 的宽度比较
两个指数表示相同曲线时，平方项系数必须相等。注意这里的 ρ 不等于标准差 σt。
$\frac1{\rho^2}=\frac1{2\sigma_t^2}\quad\Longrightarrow\quad\sigma_t=\frac{\rho}{\sqrt2}$
把基础不透明度乘上时间权重
基础 opacity 表示最高可见程度；时间核只把它衰减。在中心时刻 w=1，离开中心一倍 ρ 时 w=e⁻¹。
$o(t)=\operatorname{sigmoid}(\beta)w(t),\quad w(\tau)=1,\quad w(\tau+\rho)=e^{-1}$

代入数字，亲手核对

τ=0.5，ρ=0.25，t=0.75，β=0。

\Delta t=0.25,\qquad w=e^{-1}\approx0.367879,\qquad o(t)=0.5w\approx0.183940

这一步在代码里在哪里？

注意 temporal_opacity 的指数没有 0.5；fragment 的空间 power 才有 −0.5。

ρ ↔ temporal_scale；w ↔ temporal_opacity；o(t) ↔ opacity

src/splat.wgsl · L65–73 ↗

  let temporal_scale = max(exp(splat.r1.x), 0.000001);
  let temporal_opacity = exp(-pow(dt / temporal_scale, 2.0));
  let opacity = (1.0 / (1.0 + exp(-splat.r5.x))) * temporal_opacity;
  let scales = exp(splat.r5.yzw);
  let quaternion = normalize(splat.r6 + dt * splat.r7);
  let w = quaternion.x;
  let x = quaternion.y;
  let y = quaternion.z;
  let z = quaternion.w;

停下来想一想：把 ρ 加倍，在同样 Δt 下指数平方项变成几分之一？

四分之一；时间包络因此变宽。

带着这条公式，查看完整像素算例 →

(5.3)三次运动轨迹

𝛍 (t) = 𝛍_{0} + 𝐦_{1} Δ t + 𝐦_{2} {Δ t}^{2} + 𝐦_{3} {Δ t}^{3}

从零推导4 步 · 数字算例 · 源码 · 自测

多项式模型

把随时间移动的中心拆成可存储的系数。

用到的小工具A 向量

从静止模型开始
如果只保存 μ0，任意时刻位置都相同。这里 μ0 表示 t=τ 时的中心，并非一定是视频第一帧的位置。
$\mu(\tau)=\mu_0$
加上一次运动
每经过一个模型时间单位，增加 m1。三维中每个系数都是 x、y、z 三个数。
$\mu(t)=\mu_0+\mathbf m_1\Delta t$
加上弯曲和更灵活的变化
二次项让速度随时间变，三次项让加速度随时间变。这是有限阶模型选择；不要求真实运动必然是三次曲线。
$\mu(t)=\mu_0+\mathbf m_1\Delta t+\mathbf m_2\Delta t^2+\mathbf m_3\Delta t^3$
逐坐标代数计算
对 x 分量先平方、立方，再乘各自系数，最后相加；y、z 一样。三次轨迹也可用 Horner 顺序少做乘法，但当前 shader 直接算 dt2 方便复用。
$\mu_x(t)=\mu_{0x}+\Delta t\big(m_{1x}+\Delta t(m_{2x}+\Delta t\,m_{3x})\big)$

代入数字，亲手核对

μ0=(1,2,10)，m1=(2,0,0)，m2=(0,4,0)，m3=(0,0,8)，Δt=0.5。

\mu=(1,2,10)+(1,0,0)+(0,1,0)+(0,0,1)=(2,3,11)

这一步在代码里在哪里？

9 个 motion 浮点跨越 vec4 边界，所以 quadratic 和 cubic 要重新拼接。

m1 ↔ r2.xyz；m2 ↔ (r2.w,r3.x,r3.y)；m3 ↔ (r3.z,r3.w,r4.x)

src/splat.wgsl · L55–63 ↗

  let dt = camera.scene.x - splat.r0.w;
  let dt2 = dt * dt;
  let motion_linear = splat.r2.xyz;
  let motion_quadratic = vec3<f32>(splat.r2.w, splat.r3.x, splat.r3.y);
  let motion_cubic = vec3<f32>(splat.r3.z, splat.r3.w, splat.r4.x);
  let center = splat.r0.xyz
    + motion_linear * dt
    + motion_quadratic * dt2
    + motion_cubic * dt2 * dt;

CPU 使用同一轨迹来剔除和排序；GPU 使用它来摆放 quad，两边必须一致。

μ(t) ↔ position / center

src/lib.rs · L84–92 ↗

    fn sample(self, time: f32) -> (Vec3, f32) {
        let dt = time - self.temporal_center;
        let dt2 = dt * dt;
        let position = self.base
            + self.motion_linear * dt
            + self.motion_quadratic * dt2
            + self.motion_cubic * (dt2 * dt);
        let temporal = (-(dt / self.temporal_scale).powi(2)).exp();
        (position, sigmoid(self.opacity_logit) * temporal)

停下来想一想：Δt=0 时 motion 再大，中心是什么？

仍是 μ0，因为所有非零阶项都乘了 0。

带着这条公式，查看完整像素算例 →

(5.4)速度取决于时间单位

\begin{matrix} \frac{d 𝛍}{d t} & = & 𝐦_{1} + 2 𝐦_{2} Δ t + 3 𝐦_{3} {Δ t}^{2} \\ \frac{d 𝛍}{d s} & = & \frac{1}{D} \frac{d 𝛍}{d t} \end{matrix}

m₁ 是相对模型时间的系数；只有空间单位与 D 都明确后，才能解释为米/秒。改变归一化规则时，多项式系数也必须对应变换。

从零推导4 步 · 数字算例 · 源码 · 自测

多项式求导 + 单位换算

运动系数和实际速度、加速度差在哪？

用到的小工具E 导数 F 链式法则

从差商验证平方、立方规则
基础课已经用展开证明：t² 的导数是 2t，t³ 是 3t²。τ 固定时，Δt 对 t 的导数为 1。
$\frac{d\Delta t}{dt}=1,\qquad \frac{d\Delta t^2}{dt}=2\Delta t,\qquad \frac{d\Delta t^3}{dt}=3\Delta t^2$
每一项分别求导后相加
μ0 是常数，导数为 0；向量系数不变，对每个分量独立做相同运算。
$\frac{d\mu}{dt}=\mathbf m_1+2\mathbf m_2\Delta t+3\mathbf m_3\Delta t^2$
再求一次导数
在中心时刻 Δt=0，加速度为 2m2；三阶导数为 6m3。因此若从速度、加速度、jerk 写 Taylor 系数，会出现 1/2 和 1/6。
$\frac{d^2\mu}{dt^2}=2\mathbf m_2+6\mathbf m_3\Delta t,\qquad \frac{d^3\mu}{dt^3}=6\mathbf m_3$
从模型时间换成秒
t=(s−s0)/D，因此 dt/ds=1/D。链式法则告诉我们速度除以 D；加速度除以 D²。
$\frac{d\mu}{ds}=\frac1D\frac{d\mu}{dt},\qquad \frac{d^2\mu}{ds^2}=\frac1{D^2}\frac{d^2\mu}{dt^2}$

代入数字，亲手核对

上一例 Δt=0.5，模型时间速度为 (2,4,6)。若 D=10 秒且空间是米：

\frac{d\mu}{ds}=(0.2,0.4,0.6)\ {\rm m/s}

这一步在代码里在哪里？

renderer 只求位置，没有显式求 velocity。需要运动矢量或诊断速度时，应从这些 motion 字段按本式计算；不要声称当前 pass 已输出速度。

导数来源 ↔ motion_linear/quadratic/cubic

src/splat.wgsl · L55–63 ↗

  let dt = camera.scene.x - splat.r0.w;
  let dt2 = dt * dt;
  let motion_linear = splat.r2.xyz;
  let motion_quadratic = vec3<f32>(splat.r2.w, splat.r3.x, splat.r3.y);
  let motion_cubic = vec3<f32>(splat.r3.z, splat.r3.w, splat.r4.x);
  let center = splat.r0.xyz
    + motion_linear * dt
    + motion_quadratic * dt2
    + motion_cubic * dt2 * dt;

停下来想一想：D 加倍、模型不变，同一归一化时刻的物理速度怎样变？

减半。

带着这条公式，查看完整像素算例 →

(5.5)四元数的一阶时间模型

\begin{matrix} 𝐮 (t) & = & 𝐪_{0} + Δ t 𝛚 \\ 𝐪 (t) & = & \frac{𝐮 (t)}{‖ 𝐮 (t) ‖} \end{matrix}

q₀ 与 ω 均为四维系数，且需 u(t)≠0。这里 ω 不是刚体公式中的三维角速度；该模型也不是一般意义上的 SLERP。

从零推导4 步 · 数字算例 · 源码 · 自测

四维系数模型 + 归一化

旋转随时间变化时，为什么先相加再 normalize？

用到的小工具A 向量

先在四维系数空间里走直线
q0、ω 都有四个数，按分量计算 u=q0+Δtω。这里 ω 是四元数系数的变化，不是三维物理角速度。
$u_j(t)=q_{0j}+\Delta t\,\omega_j,\qquad j\in\{w,x,y,z\}$
直线不保证仍在单位球面上
即使 q0 是单位长度，相加后的 u 一般也不是。例如 (1,0,0,0)+(0,0,0,0.5) 的长度是 √1.25。
$\|u\|=\sqrt{u_w^2+u_x^2+u_y^2+u_z^2}$
除以长度，回到单位四元数
每项都除以同一个数，平方和就变成 1。要求 u≠0；归一化后送进前面的旋转矩阵。
$q=\frac u{\|u\|},\qquad \|q\|^2=\frac{\sum_j u_j^2}{\|u\|^2}=1$
理解角度速度不恒定
归一化把系数直线投到单位球面，因此角度一般非线性变化。这个模型不是 SLERP，也不自动满足刚体角速度微分方程。
$q_0=(1,0,0,0),\ \omega=(0,0,0,1):\quad q(t)=\frac{(1,0,0,\Delta t)}{\sqrt{1+\Delta t^2}}$

代入数字，亲手核对

Δt=0.5 时 q≈(0.894427,0,0,0.447214)，对应绕 z 旋转约 53.13°。

\theta=2\arctan(0.5)\approx53.130^\circ

这一步在代码里在哪里？

normalize(splat.r6 + dt*splat.r7) 正是这两步。loader 的有限值检查并不保证所有时刻 u 的范数非零；资产需满足这个前提。

q0 ↔ r6；ω ↔ r7；q ↔ quaternion

src/splat.wgsl · L65–73 ↗

  let temporal_scale = max(exp(splat.r1.x), 0.000001);
  let temporal_opacity = exp(-pow(dt / temporal_scale, 2.0));
  let opacity = (1.0 / (1.0 + exp(-splat.r5.x))) * temporal_opacity;
  let scales = exp(splat.r5.yzw);
  let quaternion = normalize(splat.r6 + dt * splat.r7);
  let w = quaternion.x;
  let x = quaternion.y;
  let y = quaternion.z;
  let z = quaternion.w;

停下来想一想：ω=0 时旋转会随时间变化吗？

不会；始终是 normalize(q0)，前提 q0≠0。

带着这条公式，查看完整像素算例 →

(5.6)归一化的 Jacobian

\begin{matrix} \frac{\partial 𝐪}{\partial 𝐮} & = & \frac{𝐈_{4} - 𝐪 𝐪^{T}}{‖ 𝐮 ‖} \\ \frac{d 𝐪}{d t} & = & \frac{\partial 𝐪}{\partial 𝐮} 𝛚 \end{matrix}

归一化会去掉沿 q 自身方向的变化分量；分母接近零时梯度不稳定。不要跳过零范数与有限值检查。

从零推导4 步 · 数字算例 · 源码 · 自测

逐分量求导，再组成矩阵

normalize 的梯度为什么会减掉 qqᵀ？

用到的小工具E 导数 F 链式法则 A 向量

先给长度一个短名字
令 r=||u||。因为 r²=∑uj²，两边对某个 uj 求偏导：2r ∂r/∂uj=2uj。
$\frac{\partial r}{\partial u_j}=\frac{u_j}{r}$
对 qi=ui/r 使用商法则
δij 在 i=j 时为 1，否则为 0。这只是 ui 对 uj 的偏导，不是神秘函数。
$\frac{\partial q_i}{\partial u_j}=\frac{\delta_{ij}}r-\frac{u_i}{r^2}\frac{u_j}r=\frac{\delta_{ij}}r-\frac{u_i u_j}{r^3}$
用 qi、qj 替换 ui/r、uj/r
把所有 i,j 的项排成一张 4×4 表；对角的 δij 组成 I，qiqj 组成外积 qqᵀ。
$\frac{\partial q}{\partial u}=\frac{I-qq^T}{r}$
把时间这条链接上
u=q0+(t−τ)ω，所以 du/dt=ω。右乘 ω 得到 q 的时间变化率。若变化正好沿着 q，归一化会把它完全消掉。
$\frac{dq}{dt}=\frac{I-qq^T}{r}\omega,\qquad (I-qq^T)q=q-q(q^Tq)=0$

代入数字，亲手核对

u=(2,0,0,0)，r=2，q=(1,0,0,0)。

\frac{\partial q}{\partial u}=\operatorname{diag}(0,1/2,1/2,1/2)

这一步在代码里在哪里？

GPU 的 normalize 只有前向计算。此导数用于理解训练端自动微分，本仓库 toy trainer 固定旋转，未训练 q 或 ω。

前向 q ↔ quaternion；反向 Jacobian ↔ 教学数值检查

src/splat.wgsl · L65–73 ↗

  let temporal_scale = max(exp(splat.r1.x), 0.000001);
  let temporal_opacity = exp(-pow(dt / temporal_scale, 2.0));
  let opacity = (1.0 / (1.0 + exp(-splat.r5.x))) * temporal_opacity;
  let scales = exp(splat.r5.yzw);
  let quaternion = normalize(splat.r6 + dt * splat.r7);
  let w = quaternion.x;
  let x = quaternion.y;
  let y = quaternion.z;
  let z = quaternion.w;

停下来想一想：只把 u 整体乘一个略大的正数，q 会变吗？

不会。沿径向的变化被 I−qqᵀ 消掉。

带着这条公式，查看完整像素算例 →

动一下，观察公式

同一条轨迹，不同的时间权重

这是二维示意：标记沿三次轨迹移动，亮度随时间核变化；不是完整的训练或 rasterizer。

时间 t 0.50宽度 ρ 0.22

时间权重 w(t) = 1.0000

原始资料STG · §4.1 ↗本次部署源码 · splat.wgsl ↗官方 Lite 参数与激活 ↗

06

TWO OTHER FORMULATIONS

“加上时间”，还有不同的数学方式。

理解差异

下面两种表示用于比较；当前 STG-Lite viewer 没有实现这两套 decoder。

形变网络：由规范空间预测动态参数

(6.1)参数形变的概念接口

\begin{matrix} {𝚫𝐆}_{i} (t) & = & 𝐅_{θ} ({𝐆_{i}}^{0}, t) \\ 𝐆_{i} (t) & = & Decode ({𝐆_{i}}^{0}, {𝚫𝐆}_{i} (t)) \end{matrix}

这是方法族的接口表示。Wu et al. 使用时空特征编码与 deformation decoder；平移、尺度、旋转各 head 的增量域和激活规则必须按实现核对。

从零推导4 步 · 数字算例 · 源码 · 自测

另一类模型的接口定义

如果不存运动多项式，时间怎样影响 Gaussian？

用到的小工具B 矩阵 F 链式法则

先保存一个参考状态
规范 Gaussian G⁰ 可以包含中心、log-scale、四元数、opacity 等。规范空间是共同参考，不一定是实际拍到的第一帧。
$\mathcal G^0=(\mu_0,\ell_0,q_0,\beta_0,\ldots)$
让函数读取位置和时间
Fθ 是带可学习参数 θ 的函数。最简单的网络层只是“矩阵乘输入再加偏移”，再经过非线性。例如 ReLU 逐项取 max(0,x)。实际 Wu 等方法还使用时空特征编码，不能被这个小例子替代。
$h=\operatorname{ReLU}(W_1[\mu_0;t]+b_1),\qquad \Delta\mathcal G=W_2h+b_2$
不同输出头要在自己的参数域解码
一个可选示例：中心直接加位移，log-scale 加增量后 exp，四元数相加后归一化。它只是示意，不断言所有 deformation 仓库都这样处理。
$\mu=\mu_0+\Delta\mu,\quad s=e^{\ell_0+\Delta\ell},\quad q=\operatorname{normalize}(q_0+\Delta q)$
把时间分支收口到同一个三维渲染接口
无论时间来自网络还是多项式，只要得到当前时刻的 μ、Σ、o、颜色，就能接投影、排序、合成。额外工作是存储网络参数、执行推理，并匹配训练时的 decoder 规则。
$\mathcal G(t)=\operatorname{Decode}(\mathcal G^0,F_\theta(\mathcal G^0,t))$

代入数字，亲手核对

只看 log-scale head：s0=2，Δℓ=ln1.5。

s=e^{\ln2+\ln1.5}=2\cdot1.5=3

这一步在代码里在哪里？

单独提供的教学函数演示这三种输出域的解码；没有接入 StgPass，也不是 Wu 4D-GS 模型加载器。

Δμ、Δℓ、Δq ↔ deltaPosition、deltaLogScale、deltaQuaternion

public/math/reference.mjs · L9–17 ↗

export function decodeDeformation(base, deltaPosition, deltaLogScale, deltaQuaternion) {
  // Example choice of heads; not a claim about every deformation implementation.
  return {
    position: base.position.map((v, i) => v + deltaPosition[i]),
    scale: base.logScale.map((v, i) => Math.exp(v + deltaLogScale[i])),
    quaternion: normalize(base.quaternion.map((v, i) => v + deltaQuaternion[i])),
  };
}

停下来想一想：log-scale 加 0.1 对真实尺度是加 0.1 吗？

不是，真实尺度乘 exp(0.1)≈1.10517。

带着这条公式，查看完整像素算例 →

原生 4D Gaussian：联合分布 → 固定时刻的 3D 切片

(6.2)四维联合协方差的分块

𝚺_{4} = [\begin{matrix} 𝐀 & 𝐛 \\ 𝐛^{T} & c \end{matrix}] ≻ 0

z=[xᵀ,t]ᵀ；均值 μ₄=[μxᵀ,μt]ᵀ。A 为 3×3 空间块，b 为 3×1 时空交叉项，c>0 为时间方差。四维旋转不等于一个三维四元数。

从零推导4 步 · 数字算例 · 源码 · 自测

联合变量的定义

原生 4D Gaussian 的第四维究竟存了什么？

用到的小工具D 方差 B 矩阵

把空间和时间装进同一列
前三个量是空间坐标，第四个是时间。均值同样有四个分量。单位可以不同，因此协方差的各块单位也不同。
$\mathbf z=(x,y,z,t)^T,\qquad \mu_4=(\mu_x,\mu_y,\mu_z,\mu_t)^T$
逐项写出平均偏移乘积
4×4 表中的每项都是两个坐标偏移的平均乘积。右上角三项描述各空间坐标与时间的共同变化。
$\Sigma_4=\begin{bmatrix}\sigma_{xx}&\sigma_{xy}&\sigma_{xz}&\sigma_{xt}\\\sigma_{xy}&\sigma_{yy}&\sigma_{yz}&\sigma_{yt}\\\sigma_{xz}&\sigma_{yz}&\sigma_{zz}&\sigma_{zt}\\\sigma_{xt}&\sigma_{yt}&\sigma_{zt}&\sigma_{tt}\end{bmatrix}$
把它分成四块
A 为 3×3 空间协方差；b 为 3×1 时空交叉项；c 为时间方差。bᵀ 是 1×3，因此四块能拼成 4×4。
$\Sigma_4=\begin{bmatrix}A&b\\b^T&c\end{bmatrix},\quad [A]=L^2,\quad[b]=L\,T,\quad[c]=T^2$
保证这是一团有厚度的四维形状
要求联合协方差正定。下一节会证明这等价于 c>0 且 S=A−bbᵀ/c 正定。b 的大小并非任意可取。
$c>0,\qquad S=A-\frac{bb^T}{c}\succ0$

代入数字，亲手核对

只有 x 与时间相关的例子：A=diag(2,1,1)、b=(1,0,0)、c=1。

\Sigma_4=\begin{bmatrix}2&0&0&1\\0&1&0&0\\0&0&1&0\\1&0&0&1\end{bmatrix}

这一步在代码里在哪里？

教学函数显式接收 A,b,c。当前 STG 的 32-float 记录没有这些块，不能从 r0..r7 直接假装读取它们。

A,b,c ↔ 同名教学参数；不属于现有模型格式

public/math/reference.mjs · L18–28 ↗

export function condition4D(meanSpace, meanTime, A, b, c, time) {
  if (!(c > 0)) throw new Error('Time variance must be positive');
  // Caller supplies a positive-definite joint covariance.
  const dt = time - meanTime;
  return {
    mean: meanSpace.map((v, i) => v + b[i] * dt / c),
    covariance: A.map((row, i) => row.map((v, j) => v - b[i] * b[j] / c)),
    timeWeight: Math.exp(-0.5 * dt * dt / c),
  };
}

停下来想一想：b=0 表示什么？

在联合 Gaussian 模型中空间和时间独立；固定时刻不会移动空间中心。

带着这条公式，查看完整像素算例 →

(6.3)条件 Gaussian：均值与 Schur 补

\begin{matrix} 𝛍_{x|t} & = & 𝛍_{x} + 𝐛 c^{- 1} (t - μ_{t}) \\ 𝚺_{x|t} & = & 𝐀 - 𝐛 c^{- 1} 𝐛^{T} \end{matrix}

第一式让空间中心随时间线性移动；第二式是条件协方差，不是直接拿 A。对单个固定的联合 Gaussian，条件协方差不随 t 变化。

从零推导5 步 · 数字算例 · 源码 · 自测

消去相关项 → 条件分布

为什么固定时间后要减去 bbᵀ/c，而不能直接用 A？

用到的小工具D 方差 B 矩阵 G 配方

把两个变量都先减去均值
令 X 表示三维空间偏移，T 表示时间偏移。它们均值为 0，协方差块分别是 A、b、c。
$X=\mathbf x-\mu_x,\quad T=t-\mu_t,\quad \mathbb E[XX^T]=A,\quad\mathbb E[XT]=b,\quad\mathbb E[T^2]=c$
尝试从 X 中减去一部分时间变化
设 U=X−kT，希望残差 U 不再与 T 共同变化。展开平均乘积，令它为 0，就能解出 k，而不用背条件均值公式。
$\mathbb E[UT]=\mathbb E[XT]-k\mathbb E[T^2]=b-kc=0\quad\Longrightarrow\quad k=b/c$
展开残差的协方差
四项分别取平均。两个交叉项都是 bbᵀ/c，最后一项加回来一次，因此净减去一次。
$\begin{aligned}\mathbb E[UU^T]&=\mathbb E[(X-kT)(X-kT)^T]\\&=A-bk^T-kb^T+kck^T\\&=A-\frac{bb^T}{c}-\frac{bb^T}{c}+\frac{bb^T}{c}=A-\frac{bb^T}{c}\end{aligned}$
利用联合 Gaussian 的特殊结构
一般分布中“不相关”不一定等于“独立”。但 Gaussian 的指数只有二次项；下一式的分块分解会证明 U 和 T 的指数分开，所以固定 T 不改变 U 的分布。
$\operatorname{Cov}\begin{pmatrix}U\\T\end{pmatrix}=\begin{bmatrix}S&0\\0&c\end{bmatrix},\qquad S=A-bc^{-1}b^T$
把固定时间的偏移加回中心
X=U+(b/c)T。固定时间后，后半项只是平移；残差仍是均值 0、协方差 S 的 Gaussian。
$\mu_{x\mid t}=\mu_x+\frac bc(t-\mu_t),\qquad \Sigma_{x\mid t}=A-\frac{bb^T}c$

代入数字，亲手核对

标量空间也一样：Var(X)=2，Cov(X,T)=1，Var(T)=1。固定 T=3 后，均值移到 3，方差从 2 降为 1。

\mathbb E[X\mid T=3]=3,\qquad \operatorname{Var}(X\mid T)=2-1=1

这一步在代码里在哪里？

逐项构造 shiftedMean 和 A[i][j]−b[i]*b[j]/c；本函数供手算验证，不改变 STG renderer。

条件均值 ↔ mean；S ↔ covariance

public/math/reference.mjs · L18–28 ↗

export function condition4D(meanSpace, meanTime, A, b, c, time) {
  if (!(c > 0)) throw new Error('Time variance must be positive');
  // Caller supplies a positive-definite joint covariance.
  const dt = time - meanTime;
  return {
    mean: meanSpace.map((v, i) => v + b[i] * dt / c),
    covariance: A.map((row, i) => row.map((v, j) => v - b[i] * b[j] / c)),
    timeWeight: Math.exp(-0.5 * dt * dt / c),
  };
}

停下来想一想：单个固定 4D Gaussian 的条件协方差会随 t 改变吗？

不会，它只由 A,b,c 决定；中心随 t 线性变化，时间权重随 t 变化。

带着这条公式，查看完整像素算例 →

(6.4)空间条件核 × 时间边缘核

G_{4} (𝐱, t) = G_{x|t} (𝐱) e^{- \frac{{(t - μ_{t})}^{2}}{2 c}}

对二次型配方即可得到。这里均指未归一化核；时间边缘项必须保留，决定该 primitive 在该时刻的权重。条件化后再走三维投影与合成。

从零推导5 步 · 数字算例 · 源码 · 自测

矩阵配方的完整证明

空间条件核和时间核为什么能相乘？

用到的小工具G 配方 B 矩阵 C 指数

构造一个可撤销的换坐标
令 k=b/c、S=A−bbᵀ/c；下三角块矩阵 L 把独立残差 U 和时间 T 组合成 X=U+kT。
$L=\begin{bmatrix}I&k\\0&1\end{bmatrix},\quad L^{-1}=\begin{bmatrix}I&-k\\0&1\end{bmatrix},\quad \begin{bmatrix}X\\T\end{bmatrix}=L\begin{bmatrix}U\\T\end{bmatrix}$
亲手乘开验证分块协方差
先乘左边两张表，再乘 Lᵀ；右上是 kc=b，左上是 S+kckᵀ=A。这也证明 c>0、S 正定时联合协方差正定。
$L\begin{bmatrix}S&0\\0&c\end{bmatrix}L^T=\begin{bmatrix}S+kck^T&kc\\ck^T&c\end{bmatrix}=\begin{bmatrix}A&b\\b^T&c\end{bmatrix}$
把逆矩阵放进二次型
乘积求逆要倒序。L⁻¹[X;T]=[X−kT;T]，因此不需要展开整个 4×4 逆矩阵，就能把交叉项消掉。
$\begin{bmatrix}X\\T\end{bmatrix}^{T}\Sigma_4^{-1}\begin{bmatrix}X\\T\end{bmatrix}=(X-kT)^TS^{-1}(X-kT)+\frac{T^2}{c}$
指数中的加法变成核的乘法
乘上 −1/2 后取 exp，用 exp(a+b)=exp(a)exp(b)。第一部分以 μx+kT 为中心，第二部分只依赖时间。
$G_4(\mathbf x,t)=\exp\left[-\frac12(\mathbf x-\mu_{x\mid t})^TS^{-1}(\mathbf x-\mu_{x\mid t})\right]\exp\left[-\frac{(t-\mu_t)^2}{2c}\right]$
如果使用概率密度，常数也能分开
det(L)=1，所以 det(Σ4)=det(S)c。归一化常数同样分成三维和一维两份；本页渲染核没有使用这些常数。
$\det\Sigma_4=c\det S,\qquad (2\pi)^{-2}(\det\Sigma_4)^{-1/2}=(2\pi)^{-3/2}(\det S)^{-1/2}(2\pi c)^{-1/2}$

代入数字，亲手核对

前面的例子取 t=1、空间点 (1,0,0)。它正好位于条件空间中心，因此空间核为 1。

G_4((1,0,0),1)=1\cdot e^{-1/2}\approx0.606531

这一步在代码里在哪里？

timeWeight=exp(−0.5*dt²/c) 必须保留；它和条件中心、条件协方差一起给到后续三维 renderer 才能表达这个联合核。

时间边缘核 ↔ timeWeight

public/math/reference.mjs · L18–28 ↗

export function condition4D(meanSpace, meanTime, A, b, c, time) {
  if (!(c > 0)) throw new Error('Time variance must be positive');
  // Caller supplies a positive-definite joint covariance.
  const dt = time - meanTime;
  return {
    mean: meanSpace.map((v, i) => v + b[i] * dt / c),
    covariance: A.map((row, i) => row.map((v, j) => v - b[i] * b[j] / c)),
    timeWeight: Math.exp(-0.5 * dt * dt / c),
  };
}

停下来想一想：这条时间核与 STG-Lite 的 ρ 怎样对应？

若只匹配时间曲线，ρ=√(2c)；这不意味着两种完整表示等价。

带着这条公式，查看完整像素算例 →

(6.5)一个可手算的条件化例子

𝐀 = diag (2, 1, 1), 𝐛 = {(1, 0, 0)}^{T}, c = 1 \Rightarrow 𝚺_{x|t} = 𝐈_{3}

设 μ₄=0，条件中心是 (t,0,0)，时间权重是 exp(−t²/2)。时空相关性 b 编码了这个简单运动。

从零推导4 步 · 数字算例 · 源码 · 自测

逐项代入验算

用一个可手算的 4×4 例子把条件化走完。

用到的小工具B 矩阵 G 配方

先算外积，不把它误写成点积
b 是三维列，bbᵀ 是 3×3；只有左上角为 1。
$b=(1,0,0)^T,\quad c=1,\quad bb^T/c=\begin{bmatrix}1&0&0\\0&0&0\\0&0&0\end{bmatrix}$
相减得到条件协方差
A 的第一个方差从 2 减到 1，其余保持 1。因此固定时间后的空间核是单位球。
$S=\operatorname{diag}(2,1,1)-\operatorname{diag}(1,0,0)=I_3$
条件中心随时间移动
四维均值取 0，b/c=(1,0,0)。所以中心沿 x 轴以每模型时间单位 1 个世界单位移动。
$\mu_{x\mid t}=(t,0,0)^T,\qquad w_t=e^{-t^2/2}$
直接用逆矩阵再验一次
x,t 对应的 2×2 块 [[2,1],[1,1]] 的逆是 [[1,−1],[−1,2]]；二次型配方后恰好是 (x−t)²+t²。再加 y²、z²。
$2t^2-2xt+x^2+y^2+z^2=(x-t)^2+y^2+z^2+t^2$

代入数字，亲手核对

t=2、x=2、y=z=0，空间在中心、时间偏离 2；两种路线都给同一个结果。

G_4=e^{-\frac12(4)}=e^{-2}\approx0.135335

这一步在代码里在哪里？

可以在下载的教学模块中直接调用这个 A,b,c；数值检查同时比较完整四维二次型与分解后的两个核。

A=diag(2,1,1)；b=[1,0,0]；c=1

public/math/reference.mjs · L18–28 ↗

export function condition4D(meanSpace, meanTime, A, b, c, time) {
  if (!(c > 0)) throw new Error('Time variance must be positive');
  // Caller supplies a positive-definite joint covariance.
  const dt = time - meanTime;
  return {
    mean: meanSpace.map((v, i) => v + b[i] * dt / c),
    covariance: A.map((row, i) => row.map((v, j) => v - b[i] * b[j] / c)),
    timeWeight: Math.exp(-0.5 * dt * dt / c),
  };
}

停下来想一想：若误用 A 当条件协方差，哪个轴会过宽？

x 轴，标准差会从 1 错变成 √2。

带着这条公式，查看完整像素算例 →

原始资料Deformation 4D-GS · Wu et al., 2024 ↗Native 4DGS · Yang et al., 2024 ↗Native 4DGS · 公式证明 ↗

07

APPEARANCE ≠ GEOMETRY

几何决定覆盖，外观决定颜色。

理解

不同 Gaussian 方法存储的颜色信息不同。相同字段名不一定有相同激活方式。

(7.1)静态 3DGS 常用的方向外观

𝐜_{i} (𝐯) = \sum_{ℓ = 0}^{L} \sum_{m = - ℓ}^{ℓ} 𝐚_{i ℓ m} Y_{ℓ m} (𝐯)

Yℓm 为球谐基，a 为 RGB 系数；v 必须是按实现约定的单位观察方向。不要把这条外观函数理解成可重光照的材质。

从零推导7 步 · 数字算例 · 源码 · 自测

基函数定义 → 有限展开

SH 不是魔法颜色：它是“几个方向图案乘系数再相加”。

用到的小工具H 积分 E 导数 B 矩阵

先把相机方向化为单位向量
颜色可以随观察方向变化。假设方向函数只依赖 v；SH 选择一组球面上的标准图案 Y，RGB 系数分别控制每张图案的贡献。有限阶展开是一种近似表示。
$\mathbf v=\frac{\mu-\mathbf C_w}{\|\mu-\mathbf C_w\|},\qquad \mathbf c(\mathbf v)\approx\sum_{\ell=0}^L\sum_{m=-\ell}^{\ell}\mathbf a_{\ell m}Y_{\ell m}(\mathbf v)$
从常数图案推导 DC 的系数
为了让一张图案的“平方总量”为 1，常数 k 必须满足球面面积 4π 乘 k²=1。因此 Y00=1/√(4π)。
$\int_{S^2}k^2\,d\Omega=4\pi k^2=1\quad\Longrightarrow\quad Y_{00}=\frac1{\sqrt{4\pi}}\approx0.282095$
一阶图案来自 x、y、z
球面对三个轴完全对称，且 x²+y²+z²=1，因此三个平方的球面积分都为 4π/3。各乘 √(3/4π) 后平方总量为 1。按 Graphdeco 的实数约定，一阶顺序是 −y,z,−x。
$\int x^2\,d\Omega=\int y^2\,d\Omega=\int z^2\,d\Omega=\frac{4\pi}{3},\quad (Y_{1,-1},Y_{10},Y_{11})=\sqrt{\frac3{4\pi}}(-y,z,-x)$
高阶图案也有明确的生成规则
想继续往上，可用 Legendre 多项式的定义：ℓ! 表示从 1 连乘到 ℓ，0!=1；d^ℓ 表示连续求导 ℓ 次。代入 ℓ=0,1,2 就得到下面三项。
$P_\ell(z)=\frac1{2^\ell\ell!}\frac{d^\ell}{dz^\ell}(z^2-1)^\ell,\qquad P_0=1,\quad P_1=z,\quad P_2=\frac{3z^2-1}{2}$
把纬度图案和经度图案组合
对 m≥0，先求关联多项式；再乘 sin(mφ) 或 cos(mφ) 与归一化常数。这里明确采用含 (−1)^m 的约定，避免源码之间符号不一致。
$\begin{gathered}P_\ell^m(z)=(-1)^m(1-z^2)^{m/2}\frac{d^mP_\ell}{dz^m},\qquad N_{\ell m}=\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}}\\Y_{\ell0}=N_{\ell0}P_\ell(\cos\theta),\quad Y_{\ell m}=\sqrt2N_{\ell m}P_\ell^m(\cos\theta)\cos(m\phi)\ (m>0)\\Y_{\ell,-m}=\sqrt2N_{\ell m}P_\ell^m(\cos\theta)\sin(m\phi)\ (m>0)\end{gathered}$
为什么系数可以逐项求
这些图案按上述归一化满足正交性：同一图案的乘积积分为 1，不同图案为 0。将展开式乘 Ykn 再积分，其余项消失，只剩对应系数。GS 实际通常直接优化系数；无需运行时积分。高阶正交性的完整证明可继续读 NIST DLMF，运行时只要评估多项式。
$\int Y_{\ell m}Y_{kn}\,d\Omega=\delta_{\ell k}\delta_{mn}\quad\Longrightarrow\quad \mathbf a_{kn}=\int\mathbf c(\mathbf v)Y_{kn}(\mathbf v)\,d\Omega$
回到具体文件，先检查有没有 SH
本项目 STG-Lite 的 f_dc_0..2 直接保存 RGB。静态 Graphdeco 3DGS 的 DC 则通常按 C0*a00+0.5 转成颜色。相同名字不是相同格式。
$\text{static 3DGS DC:}\quad \mathbf c=0.28209479\,\mathbf a_{00}+0.5$

代入数字，亲手核对

只用 DC 表示红通道 0.8，在静态 3DGS 约定下系数约 1.06347；当前 Lite 则直接存 0.8。

a_{00}=(0.8-0.5)/0.28209479\approx1.063472

这一步在代码里在哪里？

教学模块实现到一阶：C0*a0 + C1*(−y*a1+z*a2−x*a3)。它是用于检查 SH 系数的独立示例，未接入 Lite。

ℓ≤1 的四项 ↔ coefficients[0..3]

public/math/reference.mjs · L29–38 ↗

export function evaluateSH1(coefficients, direction) {
  // Four RGB coefficients, Graphdeco real-SH signs/order. Raw sum, no +0.5/clamp.
  const [x, y, z] = normalize(direction);
  const c0 = 1 / Math.sqrt(4 * Math.PI);
  const c1 = Math.sqrt(3 / (4 * Math.PI));
  const basis = [c0, -c1 * y, c1 * z, -c1 * x];
  return [0, 1, 2].map(channel =>
    basis.reduce((sum, value, index) => sum + value * coefficients[index][channel], 0));
}

真实 shader 读取 r4.yzw 后 clamp，没有求 SH。

当前 RGB ↔ splat.r4.yzw

src/splat.wgsl · L116–126 ↗

  let pixel_offset = 3.0 * (
    major_axis * sigma_major * corner.x
    + minor_axis * sigma_minor * corner.y
  );
  let ndc_offset = pixel_offset * vec2<f32>(2.0 * camera.viewport.z, 2.0 * camera.viewport.w);

  var output: VertexOutput;
  output.position = center_clip + vec4<f32>(ndc_offset * center_clip.w, 0.0, 0.0);
  output.local = corner * 3.0;
  output.color_opacity = vec4<f32>(clamp(splat.r4.yzw, vec3<f32>(0.0), vec3<f32>(1.0)), opacity);
  return output;

停下来想一想：最高阶 L=3 一共需要多少个 RGB 系数？

每阶 2ℓ+1 个，总数 1+3+5+7=16，也就是 (L+1)²。

带着这条公式，查看完整像素算例 →

SPHERICAL HARMONICS / A VISUAL WORKSHOP

球谐：给每个观察方向配一种颜色。

先想象站在同一个 Gaussian 周围，从上、下、左、右看它。它的形状可以不变，颜色却可以随视角变化。SH 用一组固定的“方向图案”，加上每个 Gaussian 自己学到的混合系数，保存这种变化。

从真实 4D 场景选一颗 Gaussian，查看参数与 SH 实验 →

01 · 图案 Y

像调音台里的固定音色。一个基函数输入方向，输出一个有正负的数；它本身不是 RGB。

02 · 系数 a

像每个音色的音量旋钮。每个图案有 R、G、B 三个系数，允许为负；它们是训练要学的参数。

03 · 求和 c

观察方向决定每个图案当前的数值。数值乘系数，再逐项相加，得到这颗 Gaussian 在这个视角下的颜色。

实验一 · 16 张图案，怎样调出颜色？

先点击一个小图案，再拖动它的 RGB 系数。橙色表示基函数为正，蓝色表示为负，灰色接近零。切换到“球面着色”时，球的半径完全不变；“瓣状图”才用离中心的距离表示绝对值。

图案显示最高阶 L载入一个例子

\ell=0

1

\ell=1

3

\ell=2

5

\ell=3

7

正值 +负值 −每张图各自缩放以看清形状，不能用瓣的大小比较不同基函数的数值。

图示缩放

拖动基函数图或下方任一通道球，五张图会一起旋转。可观察背面、顶部与底部；聚焦图后也可用方向键旋转、+/− 缩放、Home 重置。

01 · 一个方向图案，三个独立系数

左图是选中的一个 Yₖ：橙色为正、蓝色为负，只用来标记数值的符号，与 RGB 通道无关。R、G、B 共用这个方向图案，但分别乘以自己的系数。

Y_{2,0}(\mathbf v)=\sqrt{\frac{5}{16\pi}}(3z^2-1)

\Delta\mathbf c_k(\mathbf v)=Y_k(\mathbf v)\begin{bmatrix}a_{k,R}\\a_{k,G}\\a_{k,B}\end{bmatrix}

在通道球旁编辑这个基的 RGB 系数 ↓

02 · 三个通道分别求和，再组成 RGB

下面每颗球都累加当前最高阶以内的所有基函数。只改球下方的 R 系数：R 球和合成球会改变，G、B 球不变。选中一个基函数只是选择要编辑的系数，不会让其他项停止参与求和。

s_j(\mathbf v)=\sum_{k=0}^{(L+1)^2-1}a_{k,j}Y_k(\mathbf v),\qquad j\in\{R,G,B\}

d_j(\mathbf v)=\operatorname{clip}_{[0,1]}\!\left(\tfrac12+s_j(\mathbf v)\right),\qquad \mathbf c_{\rm display}=(d_R,d_G,d_B)

单通道的显示方式

通道本色模式只亮起对应的分量，例如 R 球显示 (dR,0,0)。灰度模式把同一个值复制到三个分量，例如 (dR,dR,dR)，便于比较。所有通道共用固定 0–1 范围，不单独归一化，也不附加光照；RGB 合成球始终显示完整颜色。

方位角 φ：绕 z 轴−45°极角 θ：从 +z 往下45°

R 通道

G 通道

B 通道

→ RGB 合成

把同一方向的 R、G、B 组成一个颜色

RGB

球下方的读数对应同一个白色方向标记。SH 求和允许为负，加 0.5 是这里采用的 3DGS 约定；官方求值只截断负颜色，本页为屏幕显示再截断到 1。被截断后，系数继续变化也可能暂时看不出颜色变化。

这里旋转的是看图的视角：系数、输入方向 v 和颜色计算结果保持不变。正面、侧面、顶部表示图示相机位于 −Y、+X、+Z 轴；下方的 φ、θ 才会改变 SH 的求值方向。

白点代表当前输入方向；虚线圆代表它在背面。这里的坐标是“相机 → Gaussian”的单位向量，所以相机实际位于相反方向。图中的球与瓣都只是函数图，不是场景中的物体或 Gaussian 椭球。

逐项核对：基函数 × 系数 = RGB 贡献

k / (ℓ,m)	当前 Y 值	RGB 系数 a	RGB 贡献 aY

表格只累加最高阶以内的项。RGB 系数分别改变三个通道，基函数的方向图案对三个通道完全相同。

从初中数学开始，推到能写进代码。

不需要先背积分和导数。下面先用长度、面积、对称性解释它们；每一步再给精确公式。积分符号可以先读作“把整颗球上所有小块的结果加起来”。

01 · 方向：从一根箭头变成球上的一点

先减位置，再把长度变成 1
设 Gaussian 中心是 μ，相机是 C。用中心减相机，得到箭头 d。用勾股定理算长度，再把三个分量都除以长度。距离变短了，指向不变。相机与中心重合时长度为零，代码必须处理这个特殊情况。
$\mathbf d=\boldsymbol\mu-\mathbf C,\quad r=\sqrt{d_x^2+d_y^2+d_z^2},\quad (x,y,z)=\mathbf v=\frac{\mathbf d}{r},\quad x^2+y^2+z^2=1$
$\mathbf d=(2,-2,1)\;\Rightarrow\;r=3\;\Rightarrow\;\mathbf v=\left(\frac23,-\frac23,\frac13\right)$
两个角度只是在描述同一个点
θ 是从北极 +z 往下量的角度。直角三角形给出高度 z = cosθ，水平圈半径是 sinθ。再在水平圆内用 φ 分出 x 和 y。θ 不是从赤道量的纬度；代码中的 sin、cos 使用弧度。
$x=\sin\theta\cos\varphi,\qquad y=\sin\theta\sin\varphi,\qquad z=\cos\theta,\qquad \theta_{\rm rad}=\theta_{\rm deg}\frac{\pi}{180}$

02 · “不重复的图案”：加权求和与正交

先从三格表格的点积说起
两张图案如果总是一起变亮，就有重复信息。把对应位置的值相乘再相加，可以测这种重合。比如 [1,1,1] 和 [−1,0,1] 的乘积和是 0：后者的正负抵消，不能冒充前者的整体亮度。把三格换成整个球面，就是函数的内积。
$[1,1,1]\cdot[-1,0,1]=-1+0+1=0,\qquad\langle f,g\rangle=\int_{S^2}f(\mathbf v)g(\mathbf v)\,d\Omega$
球面上的小格子，面积并不都一样
经纬图中同样宽的一格，在北极附近会缩小。纬圈半径是 sinθ，所以球面小块面积约等于“横向 sinθΔφ × 纵向 Δθ”。把格子越分越细，就得到 sinθ dθ dφ。单位球总面积是 4π。正交归一的意思是：不同图案的乘积和为 0，同一图案的平方和为 1。
$d\Omega=\sin\theta\,d\theta\,d\varphi,\quad\int_{S^2}1\,d\Omega=4\pi,\quad\int_{S^2}Y_iY_j\,d\Omega=\begin{cases}1&i=j\\0&i\ne j\end{cases}$

03 · 零阶与一阶：为什么是 0.282 和 0.489？

零阶：所有方向都一样
把基函数设成常数 k。整球的平方和就是“面积 × k²”。让它等于 1，取正根，便得到 DC 常数。它不是拍脑袋调出来的亮度。
$1=\int_{S^2}k^2\,d\Omega=4\pi k^2\quad\Rightarrow\quad Y_{00}=k=\frac1{\sqrt{4\pi}}\approx0.28209479$
一阶：先拿 x、y、z 当三张图案
球的左右两半完全对称：x 的正负抵消，故它与常数正交；xy 也会在 x 取反时抵消，故 x 与 y 正交。三个方向地位相同，所以 x²、y²、z² 的球面积分相等。它们相加恒等于 1，因此每一个积分都是 4π/3。
$\int x^2d\Omega=\int y^2d\Omega=\int z^2d\Omega,\quad\int(x^2+y^2+z^2)d\Omega=4\pi\;\Rightarrow\;\int z^2d\Omega=\frac{4\pi}{3}$
$\int(kz)^2d\Omega=1\;\Rightarrow\;k=\sqrt{\frac3{4\pi}}\approx0.48860251$
正负号与排列：必须和模型文件一致
正交性并不会规定哪个图案必须取负号。这里采用 Graphdeco 3DGS 的实球谐约定：一阶按 −y、z、−x 排列。若只看形状忽略这些负号，读入同一份系数就会得到错误颜色。方向从 μ−C 改成 C−μ，也会把奇数阶翻转。
$(Y_{1,-1},Y_{1,0},Y_{1,1})=\sqrt{\frac3{4\pi}}(-y,z,-x),\qquad Y_{\ell m}(-\mathbf v)=(-1)^\ell Y_{\ell m}(\mathbf v)$

04 · 二阶：从六个乘法式，得到五张新图案

六个二次项里，有一个其实是旧图案
把 x、y、z 两两相乘，得到 x²、y²、z²、xy、yz、xz 六种。但 x²+y²+z²=1 已经是零阶常数。扣掉这个重复方向，剩下五种独立变化。注意“扣掉一个”是五维空间，不是随便删掉一个多项式。
$xy,\quad yz,\quad xz,\quad x^2-y^2,\quad 3z^2-1$
为什么会出现 3z²−1？
z² 从不为负，会夹带一层“整体亮度”。上一节算过 z² 在球上的平均值是 1/3，所以先减掉 1/3，再乘 3 去分母。这样正负面积抵消，保证不再重复 DC。x²−y² 也因为对称性而平均为零。
$z^2-\frac13\quad\longrightarrow\quad3z^2-1=2z^2-x^2-y^2,\qquad\int(3z^2-1)d\Omega=0$
先算两个面积和：z⁴ 与 x²y²
把球切成薄薄的水平带。在单位球上，每条带的面积是 2π × 高度差 Δz；纬圈变短与斜坡变陡刚好抵消。因此 z⁴ 的球面积分变成 2π 乘上 −1 到 1 的 z⁴ 面积。幂函数的面积规则是：指数加一，再除以新的指数。再展开 (x²+y²+z²)²，利用三轴对称，就能求交叉项。
$\int_{S^2}z^4d\Omega=2\pi\int_{-1}^1z^4dz=2\pi\left[\frac{z^5}{5}\right]_{-1}^1=\frac{4\pi}{5}$
$4\pi=\int(x^2+y^2+z^2)^2d\Omega=3\cdot\frac{4\pi}{5}+6\int x^2y^2d\Omega\;\Rightarrow\;\int x^2y^2d\Omega=\frac{4\pi}{15}$
平方和变成 1，常数就自己出现了
xy 的平方和已经有了。对 3z²−1 和 x²−y²，先展开平方，再代入刚才的面积和。每个图案乘上“平方和的倒数开根号”，就全部归一化。最后按约定给 yz、xz 加负号。
$\int(3z^2-1)^2d\Omega=9\frac{4\pi}{5}-6\frac{4\pi}{3}+4\pi=\frac{16\pi}{5}$
$\int(x^2-y^2)^2d\Omega=2\frac{4\pi}{5}-2\frac{4\pi}{15}=\frac{16\pi}{15}$
$\begin{aligned}Y_{2,-2}&=\sqrt{\frac{15}{4\pi}}xy,&Y_{2,-1}&=-\sqrt{\frac{15}{4\pi}}yz,\\Y_{2,0}&=\sqrt{\frac{5}{16\pi}}(3z^2-1),&Y_{2,1}&=-\sqrt{\frac{15}{4\pi}}xz,\\Y_{2,2}&=\sqrt{\frac{15}{16\pi}}(x^2-y^2).&&\end{aligned}$

05 · 更高阶：Legendre 规则从哪里接上来？

“球谐”里的“谐”：一种特殊的多项式
不是把任意多项式归一化就能叫标准球谐。三维齐次调和多项式限制在单位球上，才给出这一族图案。“齐次”指每项总次数相同；“调和”指三个二阶偏导数相加为零。二阶例子 2z²−x²−y² 的结果是 4−2−2=0。单位球上可把它写成 3z²−1，但检查调和性时应使用前面的三维多项式。
$\nabla^2H=\frac{\partial^2H}{\partial x^2}+\frac{\partial^2H}{\partial y^2}+\frac{\partial^2H}{\partial z^2}=0,\qquad H(\lambda x,\lambda y,\lambda z)=\lambda^\ell H(x,y,z)$
先只沿南北变化：逐项算出 P₂ 与 P₃
Legendre 多项式提供不依赖方位角的那一个图案。下面是生成它的 Rodrigues 规则。d/dz 表示“求斜率”：zⁿ 的斜率是 nzⁿ⁻¹；连续做两次就是二阶导。按规则先展开、再逐项求导，不需要记住答案。
$P_\ell(z)=\frac{1}{2^\ell\ell!}\frac{d^\ell}{dz^\ell}(z^2-1)^\ell$
$P_2(z)=\frac18\frac{d^2}{dz^2}(z^4-2z^2+1)=\frac18\frac{d}{dz}(4z^3-4z)=\frac{12z^2-4}{8}=\frac{3z^2-1}{2}$
$P_3(z)=\frac1{48}\frac{d^3}{dz^3}(z^6-3z^4+3z^2-1)=\frac{120z^3-72z}{48}=\frac{5z^3-3z}{2}$
再加入绕圈变化：sin 与 cos 两个搭档
对每个正整数 m，有一张 cos(mφ) 图和一张 sin(mφ) 图。m=0 只需一张，因为 sin0=0。关联 Legendre 函数把南北方向的形状补齐。这里的 (−1)ᵐ 是 Condon–Shortley 相位，必须和前面的一阶符号一致。阶乘 ! 就是从该整数一路乘到 1，0!=1。
$P_\ell^m(z)=(-1)^m(1-z^2)^{m/2}\frac{d^mP_\ell(z)}{dz^m},\qquad N_{\ell m}=\sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}}$
$\begin{aligned}Y_{\ell,0}&=N_{\ell0}P_\ell(z),&&\\Y_{\ell,m}&=\sqrt2N_{\ell m}P_\ell^m(z)\cos(m\varphi)&& (m>0),\\Y_{\ell,-m}&=\sqrt2N_{\ell m}P_\ell^m(z)\sin(m\varphi)&& (m>0).\end{aligned}$
把通用规则亲手化成代码中的 xy
以 ℓ=2、m=2 为例：P₂ 求两次导得到 3，所以 P₂²=3(1−z²)。再用 sin2φ=2sinφcosφ，以及 x=sinθcosφ、y=sinθsinφ，根号和三角函数会合并成 xy。代码使用这些化简后的多项式，所以渲染时不必再求阶乘、积分或三角函数。
$\begin{aligned}Y_{2,-2}&=\sqrt2\sqrt{\frac5{4\pi}\frac{0!}{4!}}\,3(1-z^2)\sin(2\varphi)\\&=6\sqrt{\frac5{48\pi}}(1-z^2)\sin\varphi\cos\varphi=\sqrt{\frac{15}{4\pi}}xy.\end{aligned}$
为何每阶 2ℓ+1 个，L=3 总共 16 个？
m=0 给一张，m=1 到 ℓ 各给 sin、cos 两张，所以是 1+2ℓ。从零阶加到 L 阶，是连续奇数相加；正好拼成边长 L+1 的正方形。每个图案有三个颜色系数，所以 L=3 是 16 个 RGB 向量、48 个标量。
$1+3+5+\cdots+(2L+1)=(L+1)^2,\qquad k=\ell^2+\ell+m$

06 · 从颜色求系数，再把系数还原成颜色

正交像一把筛子：只留下想要的那项
先看单个颜色通道 f。假设它能由这些基函数的和表达，等式两边都乘 Yⱼ，再对球面求和。因为不同图案乘积的积分全是 0，只有第 j 项留下，而且它的平方积分是 1。于是得到求系数的公式。若只保留有限阶，得到的是在整个球面上的最小平方误差近似。
$f(\mathbf v)=\sum_k a_kY_k(\mathbf v)\quad\Rightarrow\quad\int fY_jd\Omega=\sum_k a_k\int Y_kY_jd\Omega=a_j$
计算机怎样“把所有方向加起来”？
用 N 个近似等面积分布的方向，每个代表球面面积 4π/N，就把积分变成普通循环。如果用等间隔的经纬网格，各点代表的面积不同，需要 sinθ 权重。下面的重建实验采用等面积高度带和黄金角方向，共 4096 个样本；它是数值积分近似，不是精确符号积分。
$a_k\approx\frac{4\pi}{N}\sum_{n=1}^N f(\mathbf v_n)Y_k(\mathbf v_n),\qquad f(\mathbf v)=c_{\rm target}(\mathbf v)-0.5$
算一个真实颜色：把四项全部代进去
取方向 v=(0,0,1)，使用一阶，红通道系数为 [0.4,0,0.6,0]。只有 DC 和 z 那项不为零。把它们乘系数、相加，再按 Graphdeco 的约定加 0.5。这说明系数既不是最终颜色，也不需要限制在 0 到 1。
$c_R^{\rm raw}=0.4(0.28209479)+0.6(0.48860251)=0.40599942,\qquad c_R=\max(0,c_R^{\rm raw}+0.5)=0.90599942$
官方 CUDA 的这一阶段只截掉负值；高于 1 的值仍可能保留。实验的普通屏幕色块另外截到 0–1，属于显示步骤，不能和前面的 SH 求值混为一谈。

实验二 · 同一束“小亮斑”，低阶 SH 能记住多少？

左边是人为定义的方向颜色函数，不是测量材质。右边用相同的投影系数逐级增加阶数。零阶只能记住平均色；一阶出现宽渐变；二、三阶逐渐聚拢。尖锐亮斑仍会被抹宽，附近还可能出现过冲和振铃。

\mathbf b=\frac{(1,-1,1)}{\sqrt3},\quad h(\mathbf v)=e^{(\mathbf v\cdot\mathbf b-1)/0.14},\quad\mathbf c_{\rm target}=(0.12,0.16,0.24)+(0.84,0.65,0.28)h

图示缩放

拖动任意一个球，五张图会同步旋转，便于从同一角度比较。只改变观察视角，拟合系数与整个球面上的 RMSE 不变。也支持方向键、+/− 和 Home。

目标颜色

L = 0 · 1 RGB

L = 1 · 4 RGB

L = 2 · 9 RGB

L = 3 · 16 RGB

RMSE 使用另外 2048 个方向，比较“加 0.5 后、截断前”的重建 RGB 与目标 RGB。它不是图像 PSNR，也不是训练结果。彩色图只负责展示，显示截断会隐藏负值与过冲。低阶 SH 适合平滑方向变化，不能保证逼真地表示镜面尖峰；增加阶数也不会自动获得可重光照的材质。

最后一公里 · 公式怎样走到 GS 训练和 wgpu？

μ(t), C → v
v → Y₀…Y₁₅
Σ aₖYₖ + 0.5 → RGB
projection → α
Σ Tᵢαᵢ RGBᵢ → pixel

通常先为每个 Gaussian、每个相机求一次视角颜色，再让所有被它覆盖的像素使用同一个 RGB。像素之间变化的是高斯核权重与合成顺序，不必在每个 fragment 里重复算 16 项 SH。对动态场景，用当前时刻 μ(t) 算方向；即使系数没变，运动也可能改变观察方向与颜色。

可以直接检查的运行代码

本页实验直接调用 shBasis、sampleSH 与 projectColor。源码页按本次构建生成，可以从符号一路跟到每个中间量。sampleSH 保留了 raw、shifted、shaderColor、display 四层结果，便于定位“为什么颜色不一样”。

public/math/sh.mjs ↗

const direction = mean.map((v, i) => v - camera[i]);
const { basis, contributions, raw, shaderColor } =
  sampleSH(coefficients, direction, 3);
// For each pixel: RGB += transmittance * alpha * shaderColor;

训练：怎么知道该拧哪一个系数？

先固定一个方向，只动一个红通道系数 a。颜色对它是一条直线，斜率就是该方向的 Y 值。如果 Y=0，该方向对这个系数没有直接颜色梯度；要从其他方向获得信息。下面用半平方误差，让斜率可以直接看清。

c_R^{\rm raw}=\sum_k a_{kR}Y_k,\quad\frac{\partial c_R^{\rm raw}}{\partial a_{kR}}=Y_k,\quad E=\frac12(c_R-c_R^*)^2,\quad\frac{\partial E}{\partial a_{kR}}=(c_R-c_R^*)Y_k

例如 Y=0.4886，预测红色 0.7，目标 0.8，梯度为 −0.04886。学习率 0.1 时，a 要增加 0.004886，颜色约增加 0.002387。预测向目标靠近，说明正负号合理。这是单通道、单方向、未触发 clamp 的局部例子。

a^{\rm new}=a-\eta\frac{\partial E}{\partial a},\qquad\Delta c_R=Y\Delta a=0.4886\cdot0.004886\approx0.002387

真正的 GS 训练从渲染像素的损失反传，链条还要乘本层合成权重 Tᵢαᵢ，再把所有像素、相机的贡献加起来。下式假设只改 SH 系数、alpha 和顺序固定且颜色没有被下界截断。若 shifted 小于 0，官方 clamp 分支的颜色梯度为 0；恰好等于 0 的导数按实现约定。

\frac{\partial\mathcal L}{\partial a_{ikR}}=\sum_{\text{views},\mathbf p}\frac{\partial\mathcal L}{\partial C_R(\mathbf p)}\,T_i(\mathbf p)\alpha_i(\mathbf p)\,Y_k(\mathbf v_i)

球面积分的投影只是本页的教学拟合方法。真实 GS 通常通过图像损失联合优化位置、尺度、不透明度和 SH 系数，训练视角也不一定均匀覆盖球面。因此不能把本页的积分循环当成完整 GS 训练器。

存储与接口：Rust / WGSL 接入时核对什么？

L=3 的 16 个 RGB 系数在紧密排列的 float32 文件中占 16×3×4=192 字节。但 WGSL 的 array<vec3<f32>,16> 元素步长是 16 字节，占 256 字节。最容易核对的上传格式是 16 个 vec4，每个 w 填 0，Rust 对应 [[f32;4];16]。这是系数块，不含位置和其他参数；混合结构仍要单独核对每个字段的 offset。

// Rust: upload layout for the coefficient block
#[repr(C)]
struct ShBlock { coefficients: [[f32; 4]; 16] }

// WGSL: sketch of the accumulation, once basis[] is computed
// coeff: array<vec4<f32>, 16>; basis: array<f32, 16>;
var rgb = vec3<f32>(0.5);
for (var k = 0u; k < count; k += 1u) {
    rgb += coeff[k].xyz * basis[k];
}
rgb = max(rgb, vec3<f32>(0.0));

还要核对：阶数是否表示“最高阶”还是“阶数数量”、RGB 通道排列、DC 与其余项的文件布局、基函数符号、方向朝向、是否已经加过 0.5、颜色空间与曝光。上面是接口示意，不是已经接入本项目的 SH shader。

查表：零阶到三阶的完整 16 项

以下均假设 x²+y²+z²=1，使用本页和 Graphdeco 相同的实球谐符号与索引。三阶公式里的 5z²−1 可以用 4z²−x²−y² 替换，5z²−3 可以用 2z²−3x²−3y² 替换；单位球上完全相等。

k	(ℓ,m)	Yℓm(x,y,z)
`0`	$(0,0)$	$\frac{1}{\sqrt{4\pi}}$
`1`	$(1,-1)$	$-\sqrt{\frac{3}{4\pi}}y$
`2`	$(1,0)$	$\sqrt{\frac{3}{4\pi}}z$
`3`	$(1,1)$	$-\sqrt{\frac{3}{4\pi}}x$
`4`	$(2,-2)$	$\sqrt{\frac{15}{4\pi}}xy$
`5`	$(2,-1)$	$-\sqrt{\frac{15}{4\pi}}yz$
`6`	$(2,0)$	$\sqrt{\frac{5}{16\pi}}(3z^2-1)$
`7`	$(2,1)$	$-\sqrt{\frac{15}{4\pi}}xz$
`8`	$(2,2)$	$\sqrt{\frac{15}{16\pi}}(x^2-y^2)$
`9`	$(3,-3)$	$-\sqrt{\frac{35}{32\pi}}y(3x^2-y^2)$
`10`	$(3,-2)$	$\sqrt{\frac{105}{4\pi}}xyz$
`11`	$(3,-1)$	$-\sqrt{\frac{21}{32\pi}}y(5z^2-1)$
`12`	$(3,0)$	$\sqrt{\frac{7}{16\pi}}z(5z^2-3)$
`13`	$(3,1)$	$-\sqrt{\frac{21}{32\pi}}x(5z^2-1)$
`14`	$(3,2)$	$\sqrt{\frac{105}{16\pi}}z(x^2-y^2)$
`15`	$(3,3)$	$-\sqrt{\frac{35}{32\pi}}x(x^2-3y^2)$

面试自测：一句话解释 SH 在 GS 里做什么？

对每个 Gaussian，SH 用一组固定的球面方向基和学习到的 RGB 系数表示随视角变化的颜色；渲染时将相机到中心的方向归一化，求基函数值并加权求和，再进入投影与 alpha 合成。它不直接改变 Gaussian 的位置或协方差，也不等价于可重光照的 BRDF。

推导与实现的一手依据NIST · Spherical harmonics ↗NIST · Legendre / Rodrigues ↗Graphdeco · SH conventions ↗Graphdeco · CUDA forward ↗W3C · WGSL alignment ↗

(7.2)先合成特征，再解码

\begin{matrix} 𝐅 (𝐩, t) & = & \sum_{i}^{} T_{i} α_{i} 𝐟_{i} (t) \\ 𝐂 & = & 𝐃_{θ} (𝐅, view inputs) \end{matrix}

这是 STG 完整版特征渲染的概念式；特征组成与 decoder 输入要跟源码核对。当前 Lite 使用直接 RGB，没有在运行时调用这个 MLP。

从零推导4 步 · 数字算例 · 源码 · 自测

颜色合成推广到特征

为什么有的 STG 先画特征，再通过网络得到颜色？

用到的小工具A 向量 B 矩阵

把三个颜色数字换成 D 个特征数字
alpha 仍由几何与时间决定，前缀透射率也相同。只把每层的 RGB 列换成长一点的特征列。
$\mathbf f_i(t)\in\mathbb R^D,\qquad F_d(\mathbf p,t)=\sum_iT_i\alpha_i f_{i,d}(t)$
每个通道独立做同样的加权和
堆起来就是 F=∑Tiαifi。若背景特征不为 0，也要像颜色一样加上剩余透射率乘背景特征。这里概念式默认背景特征为 0。
$\mathbf F=\sum_iT_i\alpha_i\mathbf f_i$
用 decoder 把特征和方向变成颜色
示意网络先拼接 F 与观察方向，矩阵乘法加偏移，再做非线性，最后输出三维。真实 STG 的特征分组和输入应按论文与 checkpoint 配套实现。
$\mathbf h=\operatorname{ReLU}(W_1[\mathbf F;\mathbf v]+b_1),\qquad \mathbf C=D_\theta(\mathbf F,\mathbf v)$
解释为什么不能随意交换先后顺序
非线性函数一般不满足 D(a+b)=D(a)+D(b)。所以“先解码每个 Gaussian 再合成颜色”不一定等于“先合成特征再解码”。
$D(x)=x^2:\quad D(0.5\cdot1+0.5\cdot3)=4\ne0.5D(1)+0.5D(3)=5$

代入数字，亲手核对

两层特征 (1,2)、(3,0)，alpha 都为 0.5、近到远，黑背景特征。

\mathbf F=0.5(1,2)+0.25(3,0)=(1.25,1)

这一步在代码里在哪里？

教学函数复用透明合成计算任意 D 维特征；示例 decoder 的权重只是手写练习，不是 STG 训练权重。

F ↔ compositeFeatures 返回值

public/math/reference.mjs · L39–49 ↗

export function compositeFeatures(alphas, features, background) {
  // Near to far; any feature dimension, including ordinary RGB.
  const result = background.map(() => 0);
  let transmittance = 1;
  for (let i = 0; i < alphas.length; i++) {
    for (let d = 0; d < result.length; d++) {
      result[d] += transmittance * alphas[i] * features[i][d];
    }
    transmittance *= 1 - alphas[i];
  }
  return result.map((v, d) => v + transmittance * background[d]);

Lite 的特殊情况直接合成三通道 RGB，没有 feature render target 或 MLP 调用。

D=3 且直接颜色 ↔ 当前 color_opacity.rgb

src/splat.wgsl · L130–136 ↗

fn fs_main(input: VertexOutput) -> @location(0) vec4<f32> {
  let power = -0.5 * dot(input.local, input.local);
  let alpha = min(0.99, input.color_opacity.a * exp(power));
  if (alpha < 0.0039215686) {
    discard;
  }
  return vec4<f32>(input.color_opacity.rgb * alpha, alpha);

停下来想一想：只有线性 decoder 就能随意交换合成和解码吗？

纯线性且背景处理一致时可以；有偏置时还要考虑权重之和，不能省略背景后直接搬过去。

带着这条公式，查看完整像素算例 →

原始资料3DGS · Kerbl et al., 2023 ↗STG · §4.1 ↗本次部署源码 · splat.wgsl ↗

08

PIXEL ERROR → PARAMETERS

渲染是前向过程，训练把误差传回来。

必推

先固定点集合、排序与可见性分支，理解局部梯度。完整训练还包含分段裁剪、离散的密度控制与具体实现的梯度路径。

(8.1)常见的 photometric 目标

ℒ = (1 - λ) \frac{{‖ \hat{I} - I ‖}_{1}}{M} + λ (1 - SSIM (\hat{I}, I)) + ℛ (θ)

M 为参与 L1 平均的像素通道数；R 为方法特定的正则，可能为零。这是常见目标结构，不是所有 4DGS 共享的一套固定 loss。训练相机、时间采样与 mask 必须记录。

从零推导6 步 · 数字算例 · 源码 · 自测

误差度量的选择

渲染得到图像之后，训练到底在减少哪个数字？

用到的小工具D 方差 E 导数 F 链式法则

把预测和目标逐通道相减
把所有参与训练的像素通道排成 M 个数。L1 取误差绝对值再平均；MSE 先平方再平均。二者都非负，完全一致时为 0，但它们不是同一个损失。
$e_j=\widehat I_j-I_j,\qquad L_1=\frac1M\sum_j|e_j|,\qquad L_{\rm MSE}=\frac1M\sum_j e_j^2$
为什么还可能看局部结构
只比较单个像素会忽略局部亮度和对比度关系。SSIM 先在窗口内求均值、方差与协方差，再构造相似度；下面给出一种常见加权窗口版本，权重 wn 非负且和为 1。
$\begin{gathered}\mu_a=\sum_nw_na_n,\quad \sigma_a^2=\sum_nw_n(a_n-\mu_a)^2,\quad \sigma_{ab}=\sum_nw_n(a_n-\mu_a)(b_n-\mu_b)\\\operatorname{SSIM}(a,b)=\frac{(2\mu_a\mu_b+C_1)(2\sigma_{ab}+C_2)}{(\mu_a^2+\mu_b^2+C_1)(\sigma_a^2+\sigma_b^2+C_2)}\end{gathered}$
说明 SSIM 是定义，不是假装能唯一推出来
亮度项利用 2ab=a²+b²−(a−b)²：均值差越小，比例越接近 1。结构项利用两个中心化信号的平方差展开得到类似性质。C1、C2 是防止小分母的正常数；窗口和边界处理由实现决定。
$\frac{2\mu_a\mu_b+C_1}{\mu_a^2+\mu_b^2+C_1}=1-\frac{(\mu_a-\mu_b)^2}{\mu_a^2+\mu_b^2+C_1}$
把相似度转成要最小化的损失
1−SSIM 在相似度最高时最小。λ 是两种误差的权重，R(θ) 是额外正则，例如限制运动大小；不需要时为 0。这是常见训练目标的结构，不是所有 4DGS 的统一公式。
$L=(1-\lambda)L_1+\lambda(1-\operatorname{SSIM})+\mathcal R(\theta),\qquad \mathcal R=\gamma\sum_i\|\mathbf m_i\|^2\ \text{(example)}$
当前 toy trainer 使用更容易核对的 MSE
平方的导数是 2e，目标图像视为常量；除以 M 得到每个预测通道的起始梯度。L1 在非零误差处则只取符号，零点使用选定的次梯度。
$\frac{\partial L_{\rm MSE}}{\partial\widehat I_j}=\frac{2(\widehat I_j-I_j)}M,\qquad \frac{\partial L_1}{\partial\widehat I_j}=\frac{\operatorname{sign}(e_j)}M\quad(e_j\ne0)$
把误差传回去，再更新
zero_grad 清除上一轮梯度；backward 按链式法则累计；step 使用优化器修改参数。教学页面的一步演示用 SGD 便于手算，实际 train_toy_stg.py 用 Adam。
$\theta_{\rm new}=\theta-\eta\nabla_\theta L\qquad\text{(SGD teaching step)}$

代入数字，亲手核对

预测 (0.6,0.2,0.1)，目标 (0.5,0.3,0.1)，三通道平均。

L_{\rm MSE}=(0.1^2+(-0.1)^2+0)/3=0.0066667,\quad \nabla_{\widehat I}L=(0.0666667,-0.0666667,0)

这一步在代码里在哪里？

实际代码 square().mean() 是 MSE，没有 SSIM 或运动正则。Adam 的 fitted 列表决定哪些变量会被训练。

L ↔ loss；∇L ↔ .grad；参数更新 ↔ optimizer.step()

scripts/train_toy_stg.py · L80–89 ↗

    optimizer = torch.optim.Adam(fitted, lr=0.025)
    history = []
    started = time.perf_counter()
    for step in range(args.steps):
        optimizer.zero_grad(set_to_none=True)
        loss = (render(*fitted, times, cameras)-truth).square().mean()
        if not torch.isfinite(loss):
            raise RuntimeError("Non-finite training loss")
        loss.backward()
        optimizer.step()

停下来想一想：为什么 backward 后画面可能完全没变？

backward 只算梯度；尚未 optimizer.step，或者变量不在优化器参数列表里。

带着这条公式，查看完整像素算例 →

(8.2)Gaussian 核对投影均值的梯度

\nabla_{𝛍_{2D}} α = α 𝐐 𝐝

推导：d=p−μ，∂(−½dᵀQd)/∂μ=Qd，再乘指数本身。这里固定协方差，忽略 clamp/cutoff。

从零推导5 步 · 数字算例 · 源码 · 自测

指数链式法则

把 Gaussian 中心挪向一个像素，alpha 会怎样变？

用到的小工具E 导数 F 链式法则 B 矩阵

把指数部分命名为 E
固定 C、o，令 d=p−μ。Q 对称，先按 x、y 分量把 E 展开，避免直接背矩阵求导规则。
$E=-\frac12(Q_{11}d_x^2+2Q_{12}d_xd_y+Q_{22}d_y^2),\qquad \alpha=o\,e^E$
对中心 x 分量求导
dx 对 μx 的导数是 −1，dy 不变。平方导数的 2 与前面的 1/2 抵消，两个负号又抵消。
$\frac{\partial E}{\partial\mu_x}=-\frac12\left(2Q_{11}d_x(-1)+2Q_{12}d_y(-1)\right)=Q_{11}d_x+Q_{12}d_y$
对中心 y 分量同样计算
合起来正好就是矩阵 Q 乘偏移 d。
$\frac{\partial E}{\partial\mu_y}=Q_{12}d_x+Q_{22}d_y,\qquad \nabla_\mu E=Qd$
外层指数再乘一次自己
exp 的导数还是 exp，o 视为常量，因此外层导数 o exp(E) 恰好等于 alpha。
$\nabla_\mu\alpha=\frac{\partial\alpha}{\partial E}\nabla_\mu E=\alpha Qd$
注意这只是局部偏导
中心的真实三维位置还会改变投影大小、方向外观和排序。此处锁住 C，下一条完整链式公式再把其他路径加回来。
$\delta\alpha\approx(\alpha Qd)^T\delta\mu$

代入数字，亲手核对

C=diag(4,1)、d=(2,1)、o=0.8，alpha≈0.294304。

Qd=(0.5,1)^T,\quad \nabla_\mu\alpha\approx(0.147152,0.294304)^T

这一步在代码里在哪里？

delta=pixel−uv、inverse=inv(covariance)、exp(power) 组成这条路径；PyTorch 自动微分负责同样的计算。

μ ↔ uv；d ↔ delta；Q ↔ inverse

scripts/train_toy_stg.py · L36–44 ↗

    ys, xs = torch.meshgrid(torch.arange(height, device=device)+0.5,
                           torch.arange(width, device=device)+0.5, indexing="ij")
    delta = torch.stack((xs, ys), -1)[None, None] - uv[:, :, None, None]
    inverse = torch.linalg.inv(covariance)
    power = -0.5 * torch.einsum("bnhwi,bnij,bnhwj->bnhw", delta, inverse, delta)
    temporal = torch.exp(-((times[:, None]-0.5)/0.55).square())
    alpha = (torch.sigmoid(opacity_logits)[None, :, None, None]
             * temporal[:, :, None, None] * torch.exp(power)).clamp(max=0.99)
    alpha = torch.where(alpha >= 1/255, alpha, 0)

停下来想一想：像素恰好在中心，均值偏导是多少？

0，因为 d=0；峰值附近向任意方向小挪动的影响从二阶开始。

带着这条公式，查看完整像素算例 →

(8.3)对屏幕协方差的梯度

\frac{\partial α}{\partial 𝐂} = \frac{α}{2} 𝐐 𝐝 𝐝^{T} 𝐐

利用 d(C⁻¹)=−C⁻¹(dC)C⁻¹。这里固定投影中心且 C 对称；若只存独立的 (a,b,c)，交叉项 b 的导数要合并两个对称元素的贡献。

从零推导5 步 · 数字算例 · 源码 · 自测

逆矩阵微分，拆成代数步骤

怎样知道把椭圆变宽会让某个像素变亮多少？

用到的小工具B 矩阵 E 导数 F 链式法则

从逆矩阵的定义出发
Q=C⁻¹，所以 CQ=I。让 C 增加很小的 δC，Q 增加 δQ，乘开后减去 CQ=I。
$(C+\delta C)(Q+\delta Q)=I\quad\Longrightarrow\quad C\delta Q+\delta C Q+\delta C\delta Q=0$
只保留一次小量
δCδQ 是二次小量，在导数极限中消失。左乘 Q 得到 δQ=−QδCQ。这是逆矩阵求导的来源。
$C\delta Q+\delta C Q\approx0\quad\Longrightarrow\quad\delta Q\approx-Q\,\delta C\,Q$
把它代入指数的变化
固定 d，E=−dᵀQd/2。两个负号抵消；令 v=Qd，则结果只剩 vᵀδCv/2。
$\delta E=-\frac12d^T\delta Qd\approx\frac12d^TQ\,\delta C\,Qd=\frac12v^T\delta C v$
逐元素读出系数
把 vᵀδCv 展开为 ∑ij vi vj δCij。每个 δCij 前面的系数就是 E 对那个矩阵元素的偏导；再乘指数导数 alpha。
$\delta\alpha\approx\frac\alpha2\sum_{ij}v_iv_j\delta C_{ij}\quad\Longrightarrow\quad\frac{\partial\alpha}{\partial C}=\frac\alpha2vv^T=\frac\alpha2Qdd^TQ$
实际只存一个 b 时要加两项
C12 和 C21 同时等于 b；b 增加会同时改变两个元素。因此 db 的导数是两项之和，不能只读矩阵右上角。
$C=\begin{bmatrix}a&b\\b&c\end{bmatrix}:\quad \frac{\partial\alpha}{\partial(a,b,c)}=\left(\frac\alpha2v_x^2,\ \alpha v_xv_y,\ \frac\alpha2v_y^2\right)$

代入数字，亲手核对

继续用 alpha≈0.294304、v=(0.5,1)。

\frac{\partial\alpha}{\partial(a,b,c)}\approx(0.036788,0.147152,0.147152)

这一步在代码里在哪里？

torch.linalg.inv 的反向传播与这条逆矩阵规则一致；WGSL 仅负责前向，没有对应训练 backward。

C ↔ covariance；∂α/∂C ↔ autograd 中间梯度

scripts/train_toy_stg.py · L36–44 ↗

    ys, xs = torch.meshgrid(torch.arange(height, device=device)+0.5,
                           torch.arange(width, device=device)+0.5, indexing="ij")
    delta = torch.stack((xs, ys), -1)[None, None] - uv[:, :, None, None]
    inverse = torch.linalg.inv(covariance)
    power = -0.5 * torch.einsum("bnhwi,bnij,bnhwj->bnhw", delta, inverse, delta)
    temporal = torch.exp(-((times[:, None]-0.5)/0.55).square())
    alpha = (torch.sigmoid(opacity_logits)[None, :, None, None]
             * temporal[:, :, None, None] * torch.exp(power)).clamp(max=0.99)
    alpha = torch.where(alpha >= 1/255, alpha, 0)

停下来想一想：像素位于中心时，只变宽度会改变未归一化核的峰值吗？

不会。d=0，协方差梯度为 0；它仍等于 o。

带着这条公式，查看完整像素算例 →

(8.4)对 opacity logit 的梯度

\frac{\partial α}{\partial β} = α (1 - sigmoid (β))

从零推导4 步 · 数字算例 · 源码 · 自测

sigmoid 的导数

为什么训练 opacity 用 logit，而梯度里出现 1−o？

用到的小工具E 导数 F 链式法则 C 指数

把 sigmoid 写成倒数
令 o=(1+e⁻ᵝ)⁻¹。外层倒数求导给负号，内层 e⁻ᵝ 又给负号，两者抵消。
$\frac{do}{d\beta}=-(1+e^{-\beta})^{-2}(-e^{-\beta})=\frac{e^{-\beta}}{(1+e^{-\beta})^2}$
把结果改写成 o 的乘积
一个分母给 o，剩下的比例等于 1−o。
$1-o=\frac{e^{-\beta}}{1+e^{-\beta}},\qquad \frac{do}{d\beta}=o(1-o)$
乘上固定的时间核和空间核
令 K=w(t)G(p)，与 β 无关；alpha=oK。求导得到 K o(1−o)，把 Ko 换回 alpha。
$\frac{\partial\alpha}{\partial\beta}=K\,o(1-o)=\alpha(1-o)$
再接像素损失的上游梯度
β 改变 alpha，alpha 改变像素，像素改变 loss。每个受影响像素的贡献都要累加。
$\frac{\partial L}{\partial\beta_i}=\sum_p\frac{\partial L}{\partial\alpha_i(p)}\,\alpha_i(p)(1-o_i)$

代入数字，亲手核对

β=0，o=0.5，K=0.6，所以 alpha=0.3。

\frac{\partial\alpha}{\partial\beta}=0.3(1-0.5)=0.15

这一步在代码里在哪里？

opacity_logits 参与 torch.sigmoid；导出时写回原始 logit，GPU 再解码。

β ↔ opacity_logits；o ↔ torch.sigmoid(...)

scripts/train_toy_stg.py · L36–44 ↗

    ys, xs = torch.meshgrid(torch.arange(height, device=device)+0.5,
                           torch.arange(width, device=device)+0.5, indexing="ij")
    delta = torch.stack((xs, ys), -1)[None, None] - uv[:, :, None, None]
    inverse = torch.linalg.inv(covariance)
    power = -0.5 * torch.einsum("bnhwi,bnij,bnhwj->bnhw", delta, inverse, delta)
    temporal = torch.exp(-((times[:, None]-0.5)/0.55).square())
    alpha = (torch.sigmoid(opacity_logits)[None, :, None, None]
             * temporal[:, :, None, None] * torch.exp(power)).clamp(max=0.99)
    alpha = torch.where(alpha >= 1/255, alpha, 0)

停下来想一想：为什么 β 非常大时难以继续训练 opacity？

sigmoid 接近 1，o(1−o) 接近 0，梯度饱和。

带着这条公式，查看完整像素算例 →

(8.5)位置梯度不只有投影均值这一条路

\nabla_{𝛍_{c}} ℒ = 𝐉^{T} \nabla_{𝛍_{2D}} ℒ + {(\frac{\partial vec (𝐂)}{\partial 𝛍_{c}})}^{T} \nabla_{vec (𝐂)} ℒ

本式固定外观与 opacity，只展开投影均值和协方差两条路径。J 随三维位置变化，C=JΣcJᵀ 也随位置变化；只乘 Jᵀ 会漏掉后者。SH/decoder 若使用位置相关的观察方向，还需补外观路径。转回世界坐标的梯度要乘对应旋转的转置。

从零推导6 步 · 数字算例 · 源码 · 自测

多条路径的链式法则

训练位置时，为什么只乘 Jᵀ 会漏梯度？

用到的小工具F 链式法则 B 矩阵

先画出同一位置的两条出口
三维中心同时决定屏幕均值 μ2D=π(μc) 和 Jacobian J(μc)；后者又决定 C。暂时固定颜色、opacity 和世界协方差。
$\mu_c\to\mu_{2D}\to\alpha\to L,\qquad \mu_c\to J\to C\to\alpha\to L$
先对一个三维坐标分量展开
比如 x：loss 随两个均值分量和四个协方差元素变化；每条路径都是“上游梯度×局部斜率”，最后相加。若只存 a,b,c，则按三个独立量写等价求和。
$\frac{\partial L}{\partial x}=\sum_{j=1}^2\frac{\partial L}{\partial\mu_{2D,j}}\frac{\partial\mu_{2D,j}}{\partial x}+\sum_{a,b=1}^2\frac{\partial L}{\partial C_{ab}}\frac{\partial C_{ab}}{\partial x}$
把 x,y,z 三条式子排成列
vec 只是按固定顺序把矩阵元素排成向量；一旦两边使用相同顺序，就能把逐项求和写成矩阵乘法。转置负责把输出梯度送回输入。
$\nabla_{\mu_c}L=J^T\nabla_{\mu_{2D}}L+\left(\frac{\partial\operatorname{vec}C}{\partial\mu_c}\right)^T\nabla_{\operatorname{vec}C}L$
C 的变化本身也有左右两条路
C=JΣcJᵀ+εI；δJ 可以出现在左边，也可以出现在右边。Σc 暂时固定，ε 是常数。
$\delta C=(\delta J)\Sigma_cJ^T+J\Sigma_c(\delta J)^T$
直接给一个均值路径为零的例子
球位于光轴上，x=y=0，均匀世界尺度 s。前后移动不会改变屏幕中心，但投影方差 f²s²/z² 会改变。因此只保留第一项会把深度梯度误算成 0。
$C=(f^2s^2/z^2+\varepsilon)I_2,\qquad \frac{\partial C}{\partial z}=-\frac{2f^2s^2}{z^3}I_2$
其余可训练参数也沿自己的边连接
回到世界位置需左乘 Rcwᵀ。若训练 log-scale，Σ=∑akakᵀ 且 ak=e^ℓk rk，对 ℓk 求导得到 2akakᵀ；若训练旋转，还需接四元数矩阵及 normalize 的导数。
$\nabla_{\mu_w}L=R_{cw}^T\nabla_{\mu_c}L,\qquad \frac{\partial\Sigma}{\partial\ell_k}=2a_ka_k^T$

代入数字，亲手核对

f=100、s=0.1、z=10，投影方差=1.3；深度变化率为 −0.2。

C=1.3I_2,\qquad \partial C/\partial z=-0.2I_2

这一步在代码里在哪里？

xyz 同时进入 uv 和 jacobian，这就是两条计算图路径。不要 detach jacobian 后仍声称训练的是完整位置梯度。

μc ↔ xyz；μ2D ↔ uv；C ↔ covariance

scripts/train_toy_stg.py · L27–35 ↗

    dt = times[:, None, None] - 0.5
    xyz = mean[None] + motion[None] * dt - cameras[:, None, :]
    x, y, z = xyz.unbind(-1)
    uv = torch.stack((focal*x/z + width/2, focal*y/z + height/2), -1)
    zero = torch.zeros_like(z)
    jacobian = torch.stack((focal/z, zero, -focal*x/z.square(),
                            zero, focal/z, -focal*y/z.square()), -1).reshape(*z.shape, 2, 3)
    covariance = 0.12**2 * (jacobian @ jacobian.transpose(-1, -2))
    covariance = covariance + torch.eye(2, device=device) * 0.3

停下来想一想：在哪段代码最容易漏掉第二条路径？

预先缓存或 detach J、C 的训练实现；当前 toy 的 J 由可导 xyz 重新计算。

带着这条公式，查看完整像素算例 →

(8.6)合成对 alpha 与颜色的梯度

\begin{matrix} \frac{\partial 𝐂}{\partial α_{i}} & = & T_{i} (𝐜_{i} - 𝐁_{i}) \\ \frac{\partial 𝐂}{\partial 𝐜_{i}} & = & T_{i} α_{i} 𝐈_{3} \end{matrix}

Bᵢ 是从第 i 层后面开始、按自身透射率合成的尾部颜色，包含背景。增加前景 alpha 会增加本层贡献，同时遮住后面颜色。此式假定排序固定。

从零推导4 步 · 数字算例 · 源码 · 自测

分离前缀、本层和尾部

让前景更不透明，也会挡住后景；梯度怎样同时记住两件事？

用到的小工具F 链式法则 A 向量

固定顺序，只观察第 i 层
Ai 表示前面各层已经贡献的颜色；Ti 是到达本层前的透射率；Bi 是从本层之后开始独立合成的尾部颜色，包含背景。
$\mathbf C=\mathbf A_i+T_i\left[\alpha_i\mathbf c_i+(1-\alpha_i)\mathbf B_i\right]$
锁住其他层，对 alpha 求导
Ai、Ti、ci、Bi 都视为常量。第一项贡献 ci，第二项贡献 −Bi。
$\frac{\partial\mathbf C}{\partial\alpha_i}=T_i(\mathbf c_i-\mathbf B_i)$
锁住 alpha，对当前颜色求导
RGB 各通道独立相乘，当前颜色的每一通道只改变输出的同一通道，因此 Jacobian 是系数乘 I3。
$\frac{\partial\mathbf C}{\partial\mathbf c_i}=T_i\alpha_iI_3$
接上标量 loss
像素 loss 的 RGB 梯度是三维列；对 alpha 的梯度是它与颜色差做点积。这样颜色差的正负才正确参与参数更新。
$\frac{\partial L}{\partial\alpha_i}=(\nabla_{\mathbf C}L)^TT_i(\mathbf c_i-\mathbf B_i),\qquad \nabla_{\mathbf c_i}L=T_i\alpha_i\nabla_{\mathbf C}L$

代入数字，亲手核对

前景红、后景蓝 alpha=0.5、黑背景；第一层 T1=1，B1=(0,0,0.5)。

\frac{\partial\mathbf C}{\partial\alpha_1}=(1,0,-0.5)^T

这一步在代码里在哪里？

autograd 穿过 alpha 的直接颜色贡献和 cumprod 的遮挡贡献；手写 backward 时两者都要保留。

Ai/Bi ↔ 合成链前缀/后缀；Ti ↔ transmittance

scripts/train_toy_stg.py · L45–50 ↗

    order = z.argsort(dim=1)  # Discrete ordering; no gradient through visibility changes.
    alpha = alpha.gather(1, order[:, :, None, None].expand_as(alpha))
    color = torch.sigmoid(rgb_logits)[order]
    transmittance = torch.cumprod(torch.cat((torch.ones_like(alpha[:, :1]), 1-alpha), 1), 1)[:, :-1]
    return (transmittance[:, :, :, :, None] * alpha[:, :, :, :, None]
            * color[:, :, None, None, :]).sum(1)

停下来想一想：前景颜色和后方合成颜色完全一样时，alpha 的颜色梯度是多少？

0；增加遮挡不会改变这一个像素的颜色。

带着这条公式，查看完整像素算例 →

(8.7)时间核与运动系数的导数

\begin{matrix} \frac{d w}{d t} & = & - \frac{2 Δ t}{ρ^{2}} w \\ \frac{\partial w}{\partial ℓ_{τ}} & = & \frac{2 {Δ t}^{2}}{ρ^{2}} w \\ \frac{\partial 𝛍}{\partial 𝐦_{k}} & = & {Δ t}^{k} 𝐈_{3} \end{matrix}

这些是数学函数的局部导数。官方 STG 训练在部分时间路径上使用 detach；“纸面上可导”不等于实现对所有变量都传播该梯度。

从零推导5 步 · 数字算例 · 源码 · 自测

时间链的逐项求导

时间、宽度、运动系数分别怎样接收梯度？

用到的小工具E 导数 F 链式法则 C 指数

给时间核的指数命名
E=−Δt²/ρ²；固定 τ、ρ 时，Δt 对 t 的导数是 1。先求 E 的导数，再乘外层 w。
$\frac{dE}{dt}=-\frac{2\Delta t}{\rho^2},\qquad \frac{dw}{dt}=w\frac{dE}{dt}=-\frac{2\Delta t}{\rho^2}w$
对 log-width 求导，不是直接对宽度求导
把 ρ=e^ℓτ 代入，E=−Δt²e^(−2ℓτ)。指数求导带出 −2，与前面负号抵消。
$\frac{\partial E}{\partial\ell_\tau}=2\Delta t^2e^{-2\ell_\tau}=\frac{2\Delta t^2}{\rho^2},\qquad \frac{\partial w}{\partial\ell_\tau}=\frac{2\Delta t^2}{\rho^2}w$
时间中心的符号相反
Δt=t−τ，对 τ 的偏导为 −1；所以只看时间核这条路，τ 的导数与 t 相反。若 τ 同时进入运动，还必须累加运动路径。
$\frac{\partial w}{\partial\tau}=\frac{2\Delta t}{\rho^2}w$
运动系数只出现在一次乘法里
第 k 阶系数每个分量乘同一个 Δt^k；只影响同一空间分量，因此 Jacobian 是这个标量乘 I3。
$\frac{\partial\mu}{\partial\mathbf m_k}=\Delta t^kI_3,\qquad \nabla_{\mathbf m_k}L=\Delta t^k\nabla_\mu L$
时间变化同时走多条路
一条改 opacity，一条改位置，还有旋转和可能的特征时间分支。数学导数存在不代表实现一定传递；detach 会有意切断某条边。
$\frac{dL}{dt}=\frac{\partial L}{\partial w}\frac{dw}{dt}+(\nabla_\mu L)^T\frac{d\mu}{dt}+(\nabla_q L)^T\frac{dq}{dt}+\cdots$

代入数字，亲手核对

Δt=ρ=0.25 时 w=e⁻¹。

\frac{dw}{dt}=-8e^{-1}\approx-2.943036,\qquad \frac{\partial w}{\partial\ell_\tau}=2e^{-1}\approx0.735759

这一步在代码里在哪里？

GPU 前向分别计算 dt、center、quaternion、temporal_opacity；当前 toy 只训练线性 motion，没有训练 τ、ρ 或 ω。

τ ↔ r0.w；ℓτ ↔ r1.x；mk ↔ motion 字段

src/splat.wgsl · L55–63 ↗

  let dt = camera.scene.x - splat.r0.w;
  let dt2 = dt * dt;
  let motion_linear = splat.r2.xyz;
  let motion_quadratic = vec3<f32>(splat.r2.w, splat.r3.x, splat.r3.y);
  let motion_cubic = vec3<f32>(splat.r3.z, splat.r3.w, splat.r4.x);
  let center = splat.r0.xyz
    + motion_linear * dt
    + motion_quadratic * dt2
    + motion_cubic * dt2 * dt;

停下来想一想：Δt=0 时，一次 motion 参数从这一时刻的图像得到什么位置梯度？

0，因为 ∂μ/∂m1=Δt I=0；需要其他时刻来约束运动。

带着这条公式，查看完整像素算例 →

(8.8)用有限差分检验一个梯度

\frac{\partial ℒ}{\partial θ_{k}} \approx \frac{ℒ (θ + h 𝐞_{k}) - ℒ (θ - h 𝐞_{k})}{2 h}

选 float64 的小例子，避开深度交换、near plane、cutoff 和 clamp 边界；比较多个 h 的误差趋势。单一 h 恰好接近不能证明整个 backward 正确。

从零推导5 步 · 数字算例 · 源码 · 自测

用两次前向验证反向

怎样不用相信手推或 autograd，就独立检查一个梯度？

用到的小工具E 导数 F 链式法则

一次只改一个参数
ek 表示只有第 k 项为 1 的向量。固定图像、相机、时间、其他参数，分别把 θk 增加 h 和减少 h。
$L_+=L(\theta+h e_k),\qquad L_-=L(\theta-h e_k)$
用中间对称的两点求斜率
两次参数相差 2h，所以分母是 2h。它近似中心处的变化率。
$g_{\rm FD}=\frac{L_+-L_-}{2h}\approx\frac{\partial L}{\partial\theta_k}$
解释为什么对称差分更准确
平滑函数局部展开的二次项在相减时抵消；除以 2h 后剩下的主要误差是 h² 级别。O(h²) 表示误差随 h² 缩小的量级，不是精确常数。
$L(\theta\pm h)=L(\theta)\pm hL^{\prime}+\frac{h^2}{2}L^{\prime\prime}\pm\frac{h^3}{6}L^{\prime\prime\prime}+\cdots\quad\Rightarrow\quad g_{\rm FD}=L^{\prime}+O(h^2)$
不要把 h 越小越好当口诀
过大的 h 有截断误差，过小的 h 会把相近浮点数相减，舍入误差被放大。用 float64，比较 1e−4、1e−5、1e−6 多组，而不是只挑一个看起来好的结果。
$E_{\rm abs}=|g_{\rm analytic}-g_{\rm FD}|,\qquad E_{\rm scaled}=\frac{E_{\rm abs}}{\max(1,|g_{\rm analytic}|,|g_{\rm FD}|)}$
确保两次前向处在同一个光滑分支
深度换序、clip/cutoff 边界、alpha 上限、零范数会使局部光滑假设失效。不要通过放宽容差隐藏这类变化；先选稳定例子，再单独测试分支。
$L(\theta)=(\theta-3)^2,\quad\theta=1:\quad g_{\rm FD}=\frac{(-2+h)^2-(-2-h)^2}{2h}=-4$

代入数字，亲手核对

把解析梯度用于一次小更新后，再跑前向。下降只能作为额外检查，不能替代逐参数有限差分。

\theta_{\rm new}=1-0.1(-4)=1.4,\qquad L:4\to2.56

这一步在代码里在哪里？

现有数值检查分别验证均值、协方差独立项、logit、时间和归一化导数。新增像素实验还检查包含 C 路径的完整位置梯度。

gFD ↔ finiteJacobian；解析值 ↔ expected

scripts/check-math.mjs · L40–49 ↗

function finiteJacobian(fn, x, step) {
  return tr(x.map((_, k) => {
    const a = [...x], b = [...x]; a[k] += step; b[k] -= step;
    return fn(a).map((v, j) => (v - fn(b)[j]) / (2 * step));
  }));
}
function finiteCheck(label, fn, x, expected) {
  for (const step of [1e-4, 1e-5, 1e-6]) close(`${label} h=${step}`, finiteJacobian(fn, x, step), expected, 2e-6);
}

停下来想一想：analytic=0、finite difference≈0 时只算相对误差有什么问题？

分母接近 0，比例会失真；需要绝对误差与带尺度下限的误差一起看。

带着这条公式，查看完整像素算例 →

原始资料3DGS · Kerbl et al., 2023 ↗STG · §4.1 ↗CMU · 可微 GS 编程练习 ↗

09

FROM EQUATIONS TO CODE

会解释，也要能检查。

实践

建议先完成六个小练习，再沿同一顺序走读 renderer。

投影。从 u=fx·x/z+cx 推出 J 的第一行，说出单位。
协方差。从期望定义证明 Cov(AX+b)=AΣAᵀ，解释为什么透视需要近似。
混合。手算红蓝两层各 0.5 alpha，解释前后顺序的差别。
时间。把归一化轨迹改成秒，检查速度系数与时间核宽度。
条件化。完成第 6 节的 Schur 补例子，说清它与 STG 三次轨迹的差异。
梯度。固定排序，对 μ、C 或 opacity 做一次 finite difference。

阶段	当前仓库位置	核对什么
参数读入	src/lib.rs · parse_stg_ply	字段、log-scale、logit、wxyz
时间求值	ResearchSplat::sample / vs_main	CPU 与 GPU 的同一运动和时间核
投影椭圆	projected_axis / vs_main	齐次导数、轴外积、像素单位
像素合成	fs_main / blend state	指数、cutoff、预乘、排序
原生宿主	src/stg_pass.rs	矩阵、viewport、时间契约
教学训练	scripts/train_toy_stg.py	前向、loss、backward、导出

原始资料本次部署源码 · splat.wgsl ↗项目源码 ↗完整 Field Guide ↗

ONE RECORD → ONE PIXEL → ONE UPDATE

从文件里的 32 个数，走到一次参数更新。

WATCH A TRAINING ITERATION

把训练慢放：从三帧图像到一次参数更新

这是在浏览器里实际运行的双精度 CPU 教学训练器：2 颗 STG-Lite Gaussian、3 个时间、1 台相机。目标图固定，只优化 A 的位置 μ₀x、运动 m₁x 和不透明度 β；没有密度控制或完整场景重建。

1. 读取训练样本
2. 把参数变成当前时刻的 Gaussian
3. 投影到屏幕
4. 覆盖、排序、合成像素
5. 计算误差
6. 反向传播并累积梯度
7. 准备一次 SGD 更新
8. 应用更新，再渲染核对

读取训练样本

目标是固定的三帧合成图，时间分别为 0.25、0.5、0.75。每帧取 12×10 个像素，三帧一起组成这次更新的 batch。迭代次数增加不会改变目标图。

N=3\times12\times10=360,\qquad (t_n,\mathbf p_n,\mathbf y_n)_{n=1}^{N}

这一步的实际 JavaScript 源码

public/math/training.mjs · L10–19 ↗

export function createTrainingData() {
  const truth=initialRecords();
  truth[0][0]=1.12; truth[0][8]=0.8; truth[0][20]=0.2;
  const {width,height,left,top,times}=trainingSpec;
  return times.map(time=>({time,samples:Array.from({length:width*height},(_,i)=>{
    const pixel=[left+i%width+0.5,top+Math.floor(i/width)+0.5];
    return {pixel,target:forward(truth,time,pixel).color};
  })}));
}

把参数变成当前时刻的 Gaussian

同一套参数生成三个时刻的形状。μ₀x 是基准位置，m₁x 是水平方向的一阶运动系数；β 经 sigmoid 变成基础不透明度。本实验只学习这三个参数，其他字段和 Gaussian B 固定。

\Delta t=t-\tau,\quad\boldsymbol\mu(t)=\boldsymbol\mu_0+\mathbf m_1\Delta t+\mathbf m_2\Delta t^2+\mathbf m_3\Delta t^3,\quad o=\sigma(\beta)

这一步的实际 JavaScript 源码

public/math/trace.mjs · L42–58 ↗

export function evaluateRecord(row, time, pixel) {
  const dt = time - row[3], rho = Math.max(Math.exp(row[4]), 1e-6);
  const temporal = Math.exp(-((dt / rho) ** 2));
  const mean = row.slice(0, 3).map((v, i) =>
    v + row[8 + i] * dt + row[11 + i] * dt ** 2 + row[14 + i] * dt ** 3);
  const u = row.slice(24, 28).map((v, i) => v + dt * row[28 + i]);
  const norm = Math.hypot(...u);
  if (!(norm > 0)) throw new Error('Undefined quaternion');
  const q = u.map(v => v / norm), [w, x, y, z] = q;
  const scales = row.slice(21, 24).map(Math.exp);
  // Columns of the active, right-handed rotation, exactly as splat.wgsl axes.
  const axes = [
    [1 - 2 * (y*y + z*z), 2 * (x*y + w*z), 2 * (x*z - w*y)],
    [2 * (x*y - w*z), 1 - 2 * (x*x + z*z), 2 * (y*z + w*x)],
    [2 * (x*z + w*y), 2 * (y*z - w*x), 1 - 2 * (x*x + y*y)],
  ].map((axis, i) => scaled(axis, scales[i]));
  const sigma = axes.map(outer).reduce((a, b) => a.map((row, i) => plus(row, b[i])));

投影到屏幕

位置经相机变换得到屏幕中心；协方差经 Jacobian 变成屏幕椭圆。中间图现在只画椭圆轮廓。这里的相机固定，学习者改变的是 Gaussian 参数。

\mathbf u=\pi(\boldsymbol\mu(t)),\qquad C=J\Sigma J^T+0.3I

这一步的实际 JavaScript 源码

public/math/trace.mjs · L59–67 ↗

  const [mx, my, mz] = mean, f = camera.focal;
  if (!(mz > 0)) throw new Error('Teaching camera requires positive depth');
  const uv = [f * mx / mz + camera.cx, f * my / mz + camera.cy];
  const J = [[f / mz, 0, -f * mx / mz ** 2], [0, f / mz, -f * my / mz ** 2]];
  const projectedAxes = axes.map(axis => matrixVector(J, axis));
  const covariance = multiply(multiply(J, sigma), transpose(J))
    .map((row, i) => row.map((v, j) => v + (i === j ? 0.3 : 0)));
  const [[a, b], [, c]] = covariance, determinant = a*c - b*b;
  const inverse = [[c / determinant, -b / determinant], [-b / determinant, a / determinant]];

覆盖、排序、合成像素

先算原始 alpha，再设 0.99 上限；quad 外或 alpha 小于 1/255 的贡献归零。按近到远累积颜色，每层乘以前面剩余的透射率，最后加背景。预测图由当前参数真正计算出来。

\widetilde\alpha_i=o_iw_i(t)e^{-\frac12\mathbf d_i^TC_i^{-1}\mathbf d_i},\quad T_i=\prod_{j<i}(1-\alpha_j),\quad\widehat{\mathbf y}=\sum_iT_i\alpha_i\mathbf c_i+T_{\rm bg}\mathbf c_{\rm bg}

这一步的实际 JavaScript 源码

public/math/trace.mjs · L94–109 ↗

export function forward(records, time = defaults.time, pixel = defaults.pixel, options = {}) {
  const expected = options.target ?? target;
  const splats = records.map((row, id) => ({ id, ...evaluateRecord(row, time, pixel) }))
    .sort((a, b) => a.mean[2] - b.mean[2]);
  let remaining = 1, color = [0, 0, 0];
  for (const splat of splats) {
    splat.transmittance = remaining;
    splat.contribution = scaled(splat.color, remaining * splat.alpha);
    color = plus(color, splat.contribution);
    remaining *= 1 - splat.alpha;
  }
  color = plus(color, scaled(background, remaining));
  const error = color.map((v, i) => v - expected[i]);
  return { splats, remaining, color, target: [...expected], loss: dot(error, error) / 3 };
}

计算误差

逐像素、逐 RGB 通道比较预测和目标，再对全部 360 个像素取平均。右图亮的地方误差较大。点击任一图可以追踪一个像素，但不会改变训练 batch。

\mathcal L=\frac1{3N}\sum_{n=1}^{N}\sum_{j\in\{R,G,B\}}(\widehat y_{n,j}-y_{n,j})^2

这一步的实际 JavaScript 源码

public/math/trace.mjs · L106–107 ↗

  const error = color.map((v, i) => v - expected[i]);
  return { splats, remaining, color, target: [...expected], loss: dot(error, error) / 3 };

反向传播并累积梯度

从误差倒着走：RGB → alpha → 位置或不透明度。位置还会通过投影 Jacobian 影响椭圆，因此需要两条路径。运动系数的梯度再乘 Δt；将所有像素、三个时间的贡献平均，得到下表的 batch 梯度。

\frac{\partial\mathcal L_n}{\partial\theta}=\frac{\partial\mathcal L_n}{\partial\alpha_A}\frac{\partial\alpha_A}{\partial\theta},\quad\frac{\partial\mathcal L_n}{\partial m_{1x}}=\Delta t\frac{\partial\mathcal L_n}{\partial\mu_{0x}}

这一步的实际 JavaScript 源码

public/math/trace.mjs · L110–141 ↗

export function gradient(records, parameter, time = defaults.time, pixel = defaults.pixel, options = {}) {
  // Differentiate record A, while holding the visibility/sort branches fixed.
  // Supported fields: base mean x/y/z or opacity logit. Scales/rotation remain fixed.
  if (![0, 1, 2, 20].includes(parameter)) throw new Error('Unsupported teaching parameter');
  const result = forward(records, time, pixel, options);
  const index = result.splats.findIndex(s => s.id === 0), s = result.splats[index];
  let tail = [...background];
  for (let j = result.splats.length - 1; j > index; j--) {
    const layer = result.splats[j];
    tail = plus(scaled(layer.color, layer.alpha), scaled(tail, 1-layer.alpha));
  }
  const colorSlope = scaled(s.color.map((v, i) => v - tail[i]), s.transmittance);
  const lossColorSlope = result.color.map((v, i) => 2 * (v-result.target[i]) / 3);
  const lossAlphaSlope = dot(lossColorSlope, colorSlope);
  let meanPath = 0, covariancePath = 0, opacityPath = 0;
  if (s.branch === '光滑区间') {
    if (parameter === 20) opacityPath = s.alpha * (1-s.baseOpacity);
    else {
      const [x, y, z] = s.mean, f = camera.focal;
      const dJ = parameter === 0 ? [[0, 0, -f/z**2], [0, 0, 0]] :
        parameter === 1 ? [[0, 0, 0], [0, 0, -f/z**2]] :
        [[-f/z**2, 0, 2*f*x/z**3], [0, -f/z**2, 2*f*y/z**3]];
      const left = multiply(multiply(dJ, s.sigma), transpose(s.J));
      const right = multiply(multiply(s.J, s.sigma), transpose(dJ));
      const dC = left.map((row, i) => plus(row, right[i]));
      meanPath = s.alpha * dot(s.qd, s.J.map(row => row[parameter]));
      covariancePath = s.alpha / 2 * dot(s.qd, matrixVector(dC, s.qd));
    }
  }
  const alphaSlope = meanPath + covariancePath + opacityPath;
  const analytic = lossAlphaSlope * alphaSlope;
  const finite = (options.finite === false ? [] : [1e-4, 1e-5, 1e-6]).map(h => {

准备一次 SGD 更新

此刻梯度已经算完，参数还没有修改。三个参数使用同一份旧参数求得的梯度，同时执行 θ−ηg。学习率越大，移动越远，但不保证损失下降。

\theta_j^{\rm next}=\theta_j-\eta g_j,\qquad g_j=\frac1N\sum_n\frac{\partial\mathcal L_n}{\partial\theta_j}

这一步的实际 JavaScript 源码

public/math/training.mjs · L42–50 ↗

export function trainingIteration(records, data, rate=1) {
  if(!Number.isFinite(rate)||rate<=0)throw new RangeError('Learning rate must be positive');
  const beforeRecords=copyRecords(records),before=evaluateTraining(beforeRecords,data,true);
  const afterRecords=copyRecords(beforeRecords);
  // All gradients come from the same pre-update batch; commit simultaneously.
  trainingParameters.forEach((field,j)=>{afterRecords[0][field]-=rate*before.gradients[j];});
  const after=evaluateTraining(afterRecords,data);
  if(!Number.isFinite(after.loss)||!afterRecords.flat().every(Number.isFinite))throw new Error('Non-finite training state');
  return {beforeRecords,before,afterRecords,after,rate};

应用更新，再渲染核对

参数在进入这一步时只更新一次。重新渲染同一批三帧，比较前后损失，向曲线追加一个实测点。上一小步会撤销这次更新，便于核对；下一轮重新计算梯度，不会复用旧梯度。

\mathcal L_{\rm before}=\mathcal L(\theta),\qquad\mathcal L_{\rm after}=\mathcal L(\theta-\eta g)

这一步的实际 JavaScript 源码

public/math/training.mjs · L42–50 ↗

export function trainingIteration(records, data, rate=1) {
  if(!Number.isFinite(rate)||rate<=0)throw new RangeError('Learning rate must be positive');
  const beforeRecords=copyRecords(records),before=evaluateTraining(beforeRecords,data,true);
  const afterRecords=copyRecords(beforeRecords);
  // All gradients come from the same pre-update batch; commit simultaneously.
  trainingParameters.forEach((field,j)=>{afterRecords[0][field]-=rate*before.gradients[j];});
  const after=evaluateTraining(afterRecords,data);
  if(!Number.isFinite(after.loss)||!afterRecords.flat().every(Number.isFinite))throw new Error('Non-finite training state');
  return {beforeRecords,before,afterRecords,after,rate};

查看哪一帧（全部三帧都参与训练）追踪像素 0–119

场景时间 t 选择画面；训练轮数记录参数更新了几次。

固定目标

当前预测

逐像素误差

三图共享白色像素标记。误差图显示每像素 RGB 均方根误差，固定 0–0.25 映射为暗到亮，超过 0.25 饱和；下方数值始终是未缩放的 MSE。

更新前 batch MSE

更新后 batch MSE

已完成参数更新0

参数	更新前 θ	batch 梯度 g	改变量 −ηg	候选 θ − ηg

跟踪这一个像素

实际 batch 损失（初始值 + 每次已完成更新）

与 PyTorch、真实训练和引擎的对应关系

步骤 1–5 对应取 batch、render 和计算 loss；步骤 6 对应 loss.backward()；步骤 7–8 对应 optimizer.step() 和下一次前向。本页每轮重新累积梯度，等价于先清零；位置梯度包含中心与协方差两条路径。这里使用解析梯度与 SGD，下面链接的 PyTorch toy trainer 使用自动微分与 Adam。

可见性、alpha 阈值和排序有离散分支，梯度按当前分支计算；跨过边界或学习率过大时损失可能上升。单相机、少量参数的合成拟合不能证明真实多视角重建能力。RGB 在本 STG-Lite 格式中是直接颜色，当前实验不训练 SH。

PyTorch · train_toy_stg.py ↗ · JavaScript · training.mjs ↗

展开原有单像素实验与完整数值表

这是可运行的双精度 CPU 教学算例。两颗 STG-Lite Gaussian、一个像素、已知相机、固定 RGB 目标。每次改动都会重新求时间、椭圆、覆盖、合成与梯度；它不是浏览器里的完整 GPU renderer。屏幕坐标向右、向下，深度向前。

模型时间 t 0.50像素横坐标 u 74.5对 A 的哪个参数求导？

相机 fx=fy=100，主点 (64,48)，画幅 128×96；像素纵坐标固定 68.5。世界到相机旋转为 I、平移为 0。背景 (0.04,0.05,0.08)，目标 (0.45,0.18,0.22)。更改时间和像素时目标保持固定，方便观察梯度；这不是一套多帧训练标签。

预测像素

目标像素(0.45, 0.18, 0.22)

RGB MSE对 3 个通道取平均

正在计算初始例子…

下载会把当前参数写成 binary_little_endian、32×float32 的真实 STG-Lite 记录；计算过程用 float64，导出时会有 float32 舍入。PLY 不含相机和时间单位，使用原生渲染器时仍需提供上面的相机约定。

跟着一条数据记录检查内存布局与张量形状

浮点索引 / 字节	PLY 字段	数学量	WGSL
0–2 / 0–11	x,y,z	τ 时刻的 μ₀	r0.xyz
3 / 12–15	trbf_center	τ	r0.w
4 / 16–19	trbf_scale	ln ρ	r1.x
5–7 / 20–31	nx,ny,nz	此模型布局占位，渲染不使用	r1.yzw
8–16 / 32–67	motion_0..8	m₁,m₂,m₃，各 3 个数	r2,r3,r4.x 跨行重组
17–19 / 68–79	f_dc_0..2	已激活的直接 RGB	r4.yzw
20 / 80–83	opacity	β，未激活 logit	r5.x
21–23 / 84–95	scale_0..2	ln s₁,ln s₂,ln s₃	r5.yzw
24–27 / 96–111	rot_0..3	q₀，wxyz	r6
28–31 / 112–127	omega_0..3	四维系数 ω	r7

一个记录 128 字节，恰好对应 8 个 vec4。CPU 只读取排序需要的中心、运动和 opacity；GPU storage buffer 保留整条记录，instance index 决定当前画哪条。

教学 PyTorch：B=相机/时间组合数，N=Gaussian 数。xyz: [B,N,3]；uv: [B,N,2]；J: [B,N,2,3]；C: [B,N,2,2]；delta: [B,N,H,W,2]；alpha/T: [B,N,H,W]；color: [B,N,3]；输出图像: [B,H,W,3]。None 增加可广播维度，sum(1) 沿 Gaussian 维合成。

从教学算例继续到真正的训练 → 导出 → Rust/wgpu

用 python scripts/train_toy_stg.py --steps 400 --device cuda 跑合成多相机训练；无 CUDA 时使用 --device cpu。需要 torch 与 numpy。它学习 mean、linear motion、RGB logits、opacity logits，固定尺度、旋转和时间核。
loss.backward() 把 RGB 误差经合成、Gaussian、相机和时间路径传给这四组参数。输出在 artifacts/toy-training/，其中 PLY 的 RGB 已 sigmoid，opacity 仍为 logit。
原生宿主读取 PLY，调用 StgPass::prepare 提供 V、P、viewport 和归一化 t；CPU 求值后从远到近排序，把索引传给 GPU。
每实例画 6 个顶点组成的 quad；vertex 算当前椭圆；fragment 求权重；固定混合单元做预乘 over。当前引擎没有反向传播 kernel，训练在 PyTorch 侧完成。
把真实拍摄数据接进来，还需相机标定/位姿、统一时间、初始化、训练采样、密度控制与评估。这个五点合成训练器没有包含这些步骤；完整路线见 Field Guide。

实际优化器 Adam 把梯度 g 的历史平均 m 和平方平均 v 结合，使用每个参数自己的步长。它是优化算法的设计选择，不是 Gaussian 渲染公式的必然结论；本页按钮用 SGD 以便每个数字都能手算。

m_k=\beta_1m_{k-1}+(1-\beta_1)g_k,\quad v_k=\beta_2v_{k-1}+(1-\beta_2)g_k^2

\widehat m_k=\frac{m_k}{1-\beta_1^k},\quad \widehat v_k=\frac{v_k}{1-\beta_2^k},\quad \theta_k=\theta_{k-1}-\eta\frac{\widehat m_k}{\sqrt{\widehat v_k}+\epsilon}

平方、开方、除法逐参数进行。若梯度恒为 g、m₀=0，展开递推是 (1−β₁)(1+β₁+…+β₁^(k−1))g=(1−β₁^k)g，所以除以 1−β₁^k 消掉从零开始的偏差；v 同理。这里 β₁、β₂ 是 Adam 的衰减系数，与 opacity logit β 无关。

下载并独立运行：像素前向与解析梯度 · 其他方法的教学函数 · 完整 toy trainer · 源码指纹。

PRIMARY SOURCES

与论文、与源码一起读。

推导采用统一符号；模型选择、数学恒等式和实现近似分别标明。其他方法的教学函数与当前 STG-Lite 渲染器明确区分。

3DGS · Kerbl et al., 2023 ↗STG · §4.1 ↗Deformation 4D-GS · Wu et al., 2024 ↗Native 4DGS · Yang et al., 2024 ↗STG 官方实现与 Lite 分支 ↗COLMAP 相机与文件约定 ↗CMU Learning for 3D Vision · Assignment 4 ↗3DV 2024 · 3DGS 作者教程 ↗MathML Core · 本页公式排版标准 ↗NIST DLMF · 球谐定义与正交性 ↗NIST DLMF · Legendre 多项式 ↗3DGS 官方 SH 符号与系数顺序 ↗KaTeX · 构建时排版 ↗

主公式使用 MathML；逐步推导在构建时排版为 KaTeX HTML + 可访问 MathML。字体和样式均由本站提供，无外部 CDN。源码快照的 SHA-256 可在指纹清单核对。

每一步，都能算出来。

先补 8 个工具，每次只多走一步。

先约定符号，再推公式。

用协方差描述一个椭球。

先只画一条数轴

选择一个平滑的权重

把三条独立轴的权重相乘

让椭球可以旋转和移动

把中间那张表命名为逆协方差

这一步在代码里在哪里？

给尺度一个总为正的出口

把 opacity 变成胜算

把 o 单独解出来

从文件到 GPU，只激活一次

这一步在代码里在哪里？

从一个没有偏好的小球开始

先拉伸，再旋转，最后平移

把矩阵的列写出来

得到轴外积求和，也解释正定

这一步在代码里在哪里？

先约定四元数如何相乘

用单位四元数表示旋转

再乘一次，把向量部分展开

把三个坐标轴分别代进去

用单位长度替换对角项

这一步在代码里在哪里？

把三维不确定性投影到屏幕。

先把相机搬到原点

再把世界方向转到相机方向

反解相机中心

放入 GPU 使用的四维表格

这一步在代码里在哪里？

从侧面看两只相似三角形

从长度换成像素，再加主点

用齐次坐标统一写法

对应屏幕像素中心

这一步在代码里在哪里？

只改 x，锁住 y 和 z

只改 y

最容易漏掉的是深度

把两组结果排成 2×3

这一步在代码里在哪里？

写出新均值

新样本减去新均值

按定义做外积，再求平均

把矩阵语言翻回标量

这一步在代码里在哪里？

先把中心和形状都转到相机坐标

用中心附近的切线替代弯曲投影

平移项不改变方差

把矩阵乘法改写成三个小向量

检查近似的边界

这一步在代码里在哪里？

给中心和方向加上不同的末位

先相减再通分

把 η 缩到 0，y 同理

从 NDC 换到像素

顶点回到 clip 坐标时为什么又乘 w

这一步在代码里在哪里？

从协方差到一个像素的权重。

把逆矩阵的第一列设成未知数

消元解出第一列

解第二列，再拼回去

为什么要求正定

这一步在代码里在哪里？

把像素改写为离中心的偏移

先把矩阵乘法完全展开

在椭圆自己的两条主轴上计量

让两个三角形替我们插值 local

再把实际分支写清楚

这一步在代码里在哪里？

寻找变换后方向不变的向量

什么时候有非零的解

使用配方法求根

解出方向并归一化

方差开平方才是长度

这一步在代码里在哪里？

写出保留像素的条件

先检查有没有解

取对数，再乘负数

每一步，
都能算出来。