对于3D高斯 (Gaussian Splatting) Forward&Backward的粗浅总结

笔者初学小白，有错误的地方欢迎指正！

Basic Math Knowledge

只看前向传播的话，知道协方差矩阵和特征值是啥就可，看反向传播的话需要知道矩阵求梯度的计算法则，建议读一下参考文献的[1]；

• (1) 协方差矩阵是啥[4]：
  椭球方程（半轴长分别是$a,b,c$）:
  $$\frac{(x-\bar{x})}{a^2}+\frac{(y-\bar{y}{)}^2}{b^2}+\frac{(z-\bar{z}{)}^2}{z^2}=1$$
  但是这椭球的三个轴都是平行于坐标轴的，假设对这个椭球的每个点$(x,y,z-\bar{x},\bar{y},\bar{z})$以$(\bar{x},\bar{y},\bar{z})$为中心都做一步旋转$R$，就可以得到一个任意姿态的三维空间的椭球方程$(x_{new},y_{new},z_{new})$：
  $$V(xyz)=[x-\bar{x},y-\bar{y},z-\bar{z}{]}^T$$
  $$V(xyz)=R^{-1}V(xy{z}_{new})$$
  椭球公式变成：
  $$\frac{(R^{-1}[0,:]\cdot V(xy{z}_{new}){)}^2}{a^2}+\frac{(R^{-1}[1,:]\cdot V(xy{z}_{new}){)}^2}{b^2}+\frac{(R^{-1}[2,:]\cdot V(xy{z}_{new}){)}^2}{c^2}=1$$

  $$\begin{gathered} (\frac{R^{-1}[0,0{]}^2}{a^2}+\frac{R^{-1}[1,0{]}^2}{b^2}+\frac{R^{-1}[2,0{]}^2}{c^2})(x_{new}-\bar{x}{)}^2+ \\ (\frac{R^{-1}[0,1{]}^2}{a^2}+\frac{R^{-1}[1,1{]}^2}{b^2}+\frac{R^{-1}[2,1{]}^2}{c^2})(y_{new}-\bar{y}{)}^2+ \\ (\frac{R^{-1}[0,2{]}^2}{a^2}+\frac{R^{-1}[1,2{]}^2}{b^2}+\frac{R^{-1}[2,2{]}^2}{c^2})(z_{new}-\bar{z}{)}^2+ \\ (...)(x_{new}-\bar{x})(y_{new}-\bar{y})+ \\ (...)(y_{new}-\bar{y})(z_{new}-\bar{z})+ \\ (...)(x_{new}-\bar{x})(z_{new}-\bar{z})=1 \end{gathered}$$

  $$\Rightarrow$$
  $$\left[\begin{array}{ccc}{x}_{new}-\bar{x}& y_{new}-\bar{y}& z_{new}-\bar{z}\end{array}\right]R^{-1}\left[\begin{array}{ccc}\frac{1}{a^2}& & \\ & \frac{1}{b^2}& \\ & & \frac{1}{c^2}\end{array}\right](R^{-1}{)}^T\left[\begin{array}{c}{x}_{new}-\bar{x}\\ y_{new}-\bar{y}\\ z_{new}-\bar{z}\end{array}\right]=1$$
  这个新的矩阵就是$({\Sigma }_{3\times 3}{)}^{-1}$, 特征值的定义是: $|X-\lambda I|=0$
  $$\begin{gathered} |R^{-1}\left[\begin{array}{ccc}\frac{1}{a^2}& & \\ & \frac{1}{b^2}& \\ & & \frac{1}{c^2}\end{array}\right](R^{-1}{)}^T-\lambda I|=0 \\ |\left[\begin{array}{ccc}\frac{1}{a^2}& & \\ & \frac{1}{b^2}& \\ & & \frac{1}{c^2}\end{array}\right]R-\lambda R|=0 \\ \Rightarrow \lambda =\left[\begin{array}{ccc}\frac{1}{a^2}& & \\ & \frac{1}{b^2}& \\ & & \frac{1}{c^2}\end{array}\right] \end{gathered}$$
  二阶椭圆方程也同理，因此在求出2D投影的协方差矩阵之后要求逆然后特征值是两个半轴的长度；

• (2) $\nabla$和 $\Delta$分别是啥：
  假设系统的最终输出是$S_o$, $S_o=f(C),C=g(A,B)$
  俺们可以得出：
  $$\begin{gathered} \nabla C:=\frac{d{S}_o}{dC} \\ \nabla A:=\nabla C\frac{\partial C}{\partial A} \end{gathered}$$
  很好，我学了老半天也没学会矩阵求导🐷，问题不大，下面的推导都是逐元素的，只要会高三数学的求导就可以理解。

Vanilla 3DGS, SIGGRAPH2023

https://github.com/graphdeco-inria/gaussian-splatting

概述

（主函数在diff-gaussian-rasterization/cuda_rasterizer/rasterizer_impl.cu/"int CudaRasterizer::Rasterizer::forward(...)"）

• STEP 1: 分配$(p+255/256)$个block, 每个block有256个thread, 对每个高斯做preprocessCUDA(...);
• STEP 2: 生成buffer并对高斯做排序；
• STEP 3: 分配$num\_tiles$个block,每个block有256个thread，对每个pixel做渲染;

STEP 1

• diff-gaussian-rasterization/cuda_rasterizer/forward.cu/"FORWARD::preprocess"/"preprocessCUDA":
  • 计算${\Sigma }_{2\times 2}$
  • 计算gaussians的投影半径radii, 计算gaussians属于哪个tiletiles_touched
  • 如果用的sh，将sh转成rgb；
  • 记录gaussians的像素坐标points_xy_image;
  • 将$[{\Sigma }_{2\times 2},opacity]$存成conic_opacity；
  • 记录gaussians接触到的tiles的数目tiles_touched，注意这里是数目不是序号;

STEP 2

• duplicateWithKeys和cub::DeviceRadixSort::SortPairs:
  • 创建buffer并记录每个gaussiangaussian_values_unsorted touch到了哪些tile_idsgaussian_keys_unsorted；

  // diff-gaussian-rasterization/cuda_rasterizer/rasterizer_impl.cu
  __global__ void duplicateWithKeys
  // gaussian_keys_unsorted.shape[0] = sum(each_gaussian_touch_tiles_num)
  // gaussian_keys_unsorted的每一项记录的是tile的索引;
  // gaussian_values_unsorted的每一项记录的是gaussian的索引;
  

  • cub::DeviceRadixSort::SortPairs是CUB库的一个内置函数，经过排序后得到①point_list: 排序之后的 tilegaussians 对应的gaussians_id; ② point_list_keys: 排序之后的 tilegaussians 对应的 tile_id;
  • identifyTileRanges: 这里直接按照tile_id切开就行，然后得到的imgState.ranges是dim2, shape=(N_tile, 2)，用这个可以直接索引每个tile下的gaussians_id们;

STEP 3

渲染人们比较熟悉，分成<numtiles_unumtiles_v, 1616>渲染，每个thread是这样的：

• diff-gaussian-rasterization/cuda_rasterizer/forward.cu/"renderCUDA(...)": 这个函数的最外层的for循环是针对tile滴；
  • 拿到当前pixel下的从前到后排好序的gaussians_ids：
    uint2 range = ranges[block.group_index().y * horizontal_blocks + block.group_index().x];
  • 下面是两层循环，第一层循环的循环次数是toDo代表当前block(tile)下高斯数量；内层循环的循环次数是BLOCK_SIZE(256)，代表对当前tile的pixel数；
  • 迭代每个高斯之后要block.sync()即所有当前block下的所有thread等同步；
  • 对于内层循环干的活儿很简单：
    • 计算$G$:

    float power = -0.5f * (con_o.x * d.x * d.x + con_o.z * d.y * d.y) - con_o.y * d.x * d.y;
    

    • 计算$\alpha$并按3通道循环累加颜色（如果$\alpha <1/255$或者后面的累计不透明度T<1e-4 都直接跳过，不再做累加颜色运算）
  • 结尾还有一个循环次数为256的小循环，就是把上面计算的结果赋值给pixel；
    梯度传播图
    这里的球谐是根据sh直接线性计算出来的，我这里的$c$就是计算好的colors_precomp。
    

Forward👉

① ② $(q,s)\Rightarrow {\Sigma }_{3\times 3}$

• 数学公式很简单：
  $$\begin{gathered} q={\left[\begin{array}{cccc}w& x& y& z\end{array}\right]}^T/\sqrt{w^2+x^2+y^2+z^2} \\ R=\left[\begin{array}{ccc}1-2q_2^2-2q_3^2& 2q_1{q}_2-2q_0{q}_3& 2q_1{q}_3+2q_0{q}_2\\ 2q_1{q}_2+2q_0{q}_3& 1-2q_1^2-2q_3^2& 2q_2{q}_3-2q_0{q}_1\\ 2q_1{q}_3-2q_0{q}_2& 2q_2{q}_3+2q_0{q}_1& 1-2q_1^2-2q_2^2\end{array}\right] \\ {\Sigma }_{3\times 3}=RS(RS{)}^T \end{gathered}$$
• 代码里是这样的:

  # scene/gaussian_model.py
  self.rotation_activation = torch.nn.functional.normalize
  
  @property
  def get_rotation(self):
      return self.rotation_activation(self._rotation)
  
  # gaussian_render/__init__.py
  rotations = pc.get_rotation
  // diff-gaussian-rasterization/cuda_rasterizer/forward.cu
  // __global__ void preprocessCUDA(...)
  computeCov3D(scales[idx], scale_modifier, rotations[idx], cov3Ds + idx * 6);
  cov3D = cov3Ds + idx * 6;
  
  // diff-gaussian-rasterization/cuda_rasterizer/forward.cu
  // __device__ void computeCov3D(...)
  // Create scaling matrix
  glm::mat3 S = glm::mat3(1.0f);
  S[0][0] = mod * scale.x;
  S[1][1] = mod * scale.y;
  S[2][2] = mod * scale.z;
  // Compute rotation matrix from quaternion
  glm::mat3 R = glm::mat3(
      1.f - 2.f * (y * y + z * z), 2.f * (x * y - r * z), 2.f * (x * z + r * y),
      2.f * (x * y + r * z), 1.f - 2.f * (x * x + z * z), 2.f * (y * z - r * x),
      2.f * (x * z - r * y), 2.f * (y * z + r * x), 1.f - 2.f * (x * x + y * y)
  );
  glm::mat3 M = S * R;
  // Compute 3D world covariance matrix Sigma
  glm::mat3 Sigma = glm::transpose(M) * M;
  // Covariance is symmetric, only store upper right
  cov3D[0] = Sigma[0][0];
  cov3D[1] = Sigma[0][1];
  cov3D[2] = Sigma[0][2];
  cov3D[3] = Sigma[1][1];
  cov3D[4] = Sigma[1][2];
  cov3D[5] = Sigma[2][2];
  

③ $(JW,{\Sigma }_{3\times 3})\Rightarrow {\Sigma }_{2\times 2}$

• 数学公式很简单：
  $$\begin{gathered} {\Sigma }_{2\times 2}=JW{\Sigma }_{3\times 3}(JW{)}^T \\ \left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\end{array}\right]\left[\begin{array}{c}\frac{x}{z}\\ \frac{y}{z}\\ 1\end{array}\right] \end{gathered}$$
• 代码里是这样的：

  // diff-gaussian-rasterization/cuda_rasterizer/forward.cu
  // __global__ void preprocessCUDA(...)
  float3 p_orig = { orig_points[3 * idx], orig_points[3 * idx + 1], orig_points[3 * idx + 2] };
  float3 cov = computeCov2D(p_orig, focal_x, focal_y, tan_fovx, tan_fovy, cov3D, viewmatrix);
  
  // ↓
  
  // diff-gaussian-rasterization/cuda_rasterizer/forward.cu
  //  __device__ float3 computeCov2D(...)
  float3 t = transformPoint4x3(mean, viewmatrix);
  // 区别一
  // tan_fovx = H/2 / fv 感觉就是 cv/fv
  const float limx = 1.3f * tan_fovx;
  const float limy = 1.3f * tan_fovy;
  const float txtz = t.x / t.z;
  const float tytz = t.y / t.z;
  t.x = min(limx, max(-limx, txtz)) * t.z;
  t.y = min(limy, max(-limy, tytz)) * t.z;
  glm::mat3 J = glm::mat3(
      focal_x / t.z, 0.0f, -(focal_x * t.x) / (t.z * t.z),
      0.0f, focal_y / t.z, -(focal_y * t.y) / (t.z * t.z),
      0, 0, 0);
  glm::mat3 W = glm::mat3(
      viewmatrix[0], viewmatrix[4], viewmatrix[8],
      viewmatrix[1], viewmatrix[5], viewmatrix[9],
      viewmatrix[2], viewmatrix[6], viewmatrix[10]);
  glm::mat3 T = W * J;
  glm::mat3 Vrk = glm::mat3(
      cov3D[0], cov3D[1], cov3D[2],
      cov3D[1], cov3D[3], cov3D[4],
      cov3D[2], cov3D[4], cov3D[5]);
  glm::mat3 cov = glm::transpose(T) * glm::transpose(Vrk) * T;
  // 区别二
  // Apply low-pass filter: every Gaussian should be at least
  // one pixel wide/high. Discard 3rd row and column.
  cov[0][0] += 0.3f;
  cov[1][1] += 0.3f;
  return { float(cov[0][0]), float(cov[0][1]), float(cov[1][1]) };
  

• 这里和公式有两处小区别：
  • 区别一：这里写的min,max很迷惑，写成if else就懂了

    limx = 1.3f * cv/fv;
    const float txtz = t.x / t.z;
    if txtz < -limx:
        txtz = -limx
    if txtz > limx:
        txtz = limx
    t.x *= txtz
  

  • 区别二：让高斯投影大一些，保证每个高斯的投影半径至少大于1个pixel。（感觉像是抑制低频采样高频的情况）

④⑤ $({\Sigma }_{2\times 2},u{v}_{pix},u{v}_{proj})\Rightarrow G$ & $(G,o)\Rightarrow \alpha$

• 数学公式是这样的：
  $$\begin{gathered} G=exp(-\frac{\left[\begin{array}{cc}u-\bar{u}& v-\bar{v}\end{array}\right]{\Sigma }_{2\times 2}^{-1}\left[\begin{array}{c}u-\bar{u}\\ v-\bar{v}\end{array}\right]}{2}) \\ \alpha =G\cdot o \end{gathered}$$
• 代码是这样的：

  // diff-gaussian-rasterization/cuda_rasterizer/forward.cu
  // __global__ void __launch_bounds__(BLOCK_X * BLOCK_Y) renderCUDA(...)
  // 对于循环
  for (int i = 0; i < rounds; i++, toDo -= BLOCK_SIZE)
      // 对于每个pixel循环
      for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++)
          ...
          float2 xy = collected_xy[j];
          float2 d = { xy.x - pixf.x, xy.y - pixf.y };
          float4 con_o = collected_conic_opacity[j];
          float power = -0.5f * (con_o.x * d.x * d.x + con_o.z * d.y * d.y) - con_o.y * d.x * d.y;
          ...
          float alpha = min(0.99f, con_o.w * exp(power));
  

⑥ $(\alpha ,c)\Rightarrow rgb$

• 数学公式是这样的:
  $$rgb=\sum _{i=0}^N(\prod _{j=0}^i(1-{\alpha }_j)){\alpha }_i$$
• 代码是这样的:

  // diff-gaussian-rasterization/cuda_rasterizer/forward.cu
  // __global__ void __launch_bounds__(BLOCK_X * BLOCK_Y) renderCUDA(...)
  // 对于循环
  for (int i = 0; i < rounds; i++, toDo -= BLOCK_SIZE)
      // 对于每个pixel循环
      for (int j = 0; !done && j < min(BLOCK_SIZE, toDo); j++)
          ...
          float test_T = T * (1 - alpha);
          // Eq. (3) from 3D Gaussian splatting paper.
          for (int ch = 0; ch < CHANNELS; ch++)
              C[ch] += features[collected_id[j] * CHANNELS + ch] * alpha * T;
  
          T = test_T
  
  // __global__ void __launch_bounds__(BLOCK_X * BLOCK_Y) renderCUDA(...)
  // All threads that treat valid pixel write out their final
  // rendering data to the frame and auxiliary buffers.
  if (inside)
  {
      final_T[pix_id] = T;
      n_contrib[pix_id] = last_contributor;
      for (int ch = 0; ch < CHANNELS; ch++)
          out_color[ch * H * W + pix_id] = C[ch] + T * bg_color[ch];
  }
  

Backward👈

• 反向传播的任务也很简单，根据$\nabla rgb\Rightarrow (\nabla q,\nabla s,\nabla o,\nabla c)$，反传的主函数在diff-gaussian-rasterization/cuda_rasterizer/rasterizer_impl.cu/"RasterizeGaussiansBackwardCUDA(...)"/"CudaRasterizer::Rasterizer::backward(...)", 和前向传播正好相反，先BACKWARD::render(...)然后再BACKWARD::preprocess(...)；

⑥⑤④ $\nabla rgb\Rightarrow (\nabla o,\nabla c,\nabla u{v}_{proj},\nabla (({\sum }_{2\times 2}{)}^{-1})$

• 数学公式是这样的：

  • $$\nabla rgb\Rightarrow \nabla c$$
    $$\begin{gathered} rgb+=T_i{\alpha }_i{c}_i \\ \nabla c+=\nabla rgb{T}_i{\alpha }_i \end{gathered}$$
    $$({\Sigma }_{2\times 2}{)}^{-1}=\left[\begin{array}{cc}\alpha & \beta \\ \beta & \gamma \end{array}\right]$$

  • $$\nabla rgb\Rightarrow (\nabla o,\nabla u{v}_{proj},\nabla ({\Sigma }_{2\times 2}{)}^{-1})$$
    $$\begin{gathered} \nabla {\alpha }_i=\nabla rgb{T}_i{c}_i \\ G=exp(-\frac{\left[\begin{array}{cc}u-\bar{u}& v-\bar{v}\end{array}\right]{\Sigma }_{2\times 2}^{-1}\left[\begin{array}{c}u-\bar{u}\\ v-\bar{v}\end{array}\right]}{2})=exp[\frac{(u-\bar{u}{)}^2\alpha +2(u-\bar{u})(v-\bar{v})\beta +(v-\bar{v}{)}^2\gamma }{-2}] \\ \alpha =G\cdot o \end{gathered}$$
    $$\nabla o=\nabla \alpha \cdot G$$

$$\nabla G=\nabla \alpha \cdot o$$
$$\begin{gathered} \nabla \bar{u}=\nabla G\cdot G\cdot [(u-\bar{u})\alpha +(v-\bar{v})\beta ] \\ \nabla \bar{v}=\nabla G\cdot G\cdot [(v-\bar{v})\gamma +(u-\bar{u})\beta ] \end{gathered}$$

$$\begin{gathered} \nabla \alpha =\nabla G\cdot G\cdot \frac{(u-\bar{u}{)}^2}{-2} \\ \nabla \beta =\nabla G\cdot G\cdot \frac{(u-\bar{u})(v-\bar{v})}{-1} \\ \nabla \gamma =\nabla G\cdot G\cdot \frac{(v-\bar{v}{)}^2}{-2} \end{gathered}$$

• 代码是这样的：
  • $\nabla rgb\Rightarrow \nabla c$

    // diff-gaussian-rasterization/cuda_rasterizer/backward.cu
    // renderCUDA(...)
    const float dchannel_dcolor = alpha * T;
    const float dL_dchannel = dL_dpixel[ch];
    
    atomicAdd(&(dL_dcolors[global_id * C + ch]), dchannel_dcolor * dL_dchannel);
      - $$\nabla rgb \Rightarrow (\nabla o, \nabla uv_{proj}, \nabla (\Sigma_{2 \times 2})^{-1} )$$
        - `gdx`是$$G \cdot (v-\bar{v})$$
    // diff-gaussian-rasterization/cuda_rasterizer/backward.cu
    // renderCUDA(...)
    const float2 d = { xy.x - pixf.x, xy.y - pixf.y };
    const float gdx = G * d.x;
    const float gdy = G * d.y;
    
    const float ddelx_dx = 0.5 * W;
    const float dL_dG = con_o.w * dL_dalpha;
    const float dG_ddelx = -gdx * con_o.x - gdy * con_o.y;
    atomicAdd(&dL_dmean2D[global_id].x, dL_dG * dG_ddelx * ddelx_dx);
    
    atomicAdd(&dL_dconic2D[global_id].x, -0.5f * gdx * d.x * dL_dG);
    
    // \nabla opacity.
    atomicAdd(&(dL_dopacity[global_id]), G * dL_dalpha);
    

③②① $(\nabla u{v}_{proj},\nabla (({\Sigma }_{2\times 2}{)}^{-1})\Rightarrow (\nabla s,\nabla q,\nabla xyz)$

• 数学公式是这样的：

  • $\nabla u{v}_{proj}\Rightarrow \nabla xyz$
    $$\begin{gathered} \nabla x=\nabla u\frac{\partial u}{\partial y}+\nabla v\frac{\partial v}{\partial x} \\ 记w2c为(R,t): \\ u=fu\cdot \frac{y_{cam}}{z_{cam}}+cu=fu\cdot \frac{R_{21}x+R_{22}y+R_{23}z+t_2}{R_{31}x+R_{32}y+R_{33}z+t_3}+cu \\ v=fu\cdot \frac{y_{cam}}{z_{cam}}+cu=fv\cdot \frac{R_{11}x+R_{12}y+R_{13}z+t_2}{R_{31}x+R_{32}y+R_{33}z+t_3}+cv \\ \nabla x=\nabla u\cdot \frac{R_{21}{z}_{cam}-R_{31}{y}_{cam}}{z_{cam}^2}+\nabla v\cdot \frac{R_{11}{z}_{cam}-R_{31}{x}_{cam}}{z_{cam}^2} \\ （\nabla y,\nabla z一样，此处省略） \end{gathered}$$

  • $$\nabla ({\Sigma }_{2\times 2}{)}^{-1}\Rightarrow \nabla {\Sigma }_{3\times 3}$$
    $$\begin{gathered} M:=RS \\ \nabla ({\Sigma }_{3\times 3})=\nabla ({\Sigma }_{2\times 2}^{-1})\cdot \frac{\partial [(M{\Sigma }_{3\times 3}{M}^T{)}^{-1}]}{\partial {\Sigma }_{3\times 3}} \end{gathered}$$

很好我是废物，这种复杂矩阵求导不会算，那就拆开写变量哈；

$$\begin{gathered} {\Sigma }_{2\times 2}=\left[\begin{array}{cc}a& b\\ b& c\end{array}\right] \\ ({\Sigma }_{2\times 2}{)}^{-1}=\left[\begin{array}{cc}c& -b\\ -b& a\end{array}\right]/(ac-b^2):=\left[\begin{array}{cc}\alpha & \beta \\ \beta & \gamma \end{array}\right] \\ \{\begin{array}{c}\nabla a=\nabla \alpha \cdot \frac{-c^2}{(ac-b^2{)}^2}+\nabla \beta \cdot \frac{2bc}{(ac-b^2{)}^2}+\nabla \gamma \cdot \frac{-b^2}{(ac-b^2{)}^2}\\ \nabla b=\nabla \alpha \cdot \frac{2bc}{(ac-b^2{)}^2}+\nabla \beta \cdot \frac{-ac-b^2}{(ac-b^2{)}^2}+\nabla \gamma \cdot \frac{2ab}{(ac-b^2{)}^2}\\ \nabla c=\nabla \alpha \cdot \frac{-b^2}{(ac-b^2{)}^2}+\nabla \beta \cdot \frac{ab}{(ac-b^2{)}^2}+\nabla \gamma \cdot \frac{-a^2}{(ac-b^2{)}^2}\end{array} \end{gathered}$$

$$\begin{gathered} {\Sigma }_{3\times 3}:=\left[\begin{array}{ccc}{\Sigma }_{3,3}[0]& {\Sigma }_{3,3}[1]& {\Sigma }_{3,3}[3]\\ {\Sigma }_{3,3}[1]& {\Sigma }_{3,3}[2]& {\Sigma }_{3,3}[4]\\ {\Sigma }_{3,3}[3]& {\Sigma }_{3,3}[4]& {\Sigma }_{3,3}[5]\end{array}\right] \\ T:=WJ=\left[\begin{array}{cc}T[0,0]& T[0,1]\\ T[1,0]& T[1,1]\\ T[2,0]& T[2,2]\end{array}\right] \\ {\Sigma }_{2\times 2}=T^T{\Sigma }_{3\times 3}T\Rightarrow \\ \{\begin{array}{c}\alpha =T[0,0{]}^2{\Sigma }_{3,3}[0]+T[1,0{]}^2{\Sigma }_{3,3}[1]+T[2,0{]}^2{\Sigma }_{3\times 3}[3]\\ +T[0,0{]}^2{\Sigma }_{3,3}[1]+T[1,0{]}^2{\Sigma }_{3,3}[2]+T[2,0{]}^2{\Sigma }_{3,3}[4]\\ +T[0,0{]}^2{\Sigma }_{3,3}[3]+T[1,0{]}^2{\Sigma }_{3,3}[4]+T[2,0{]}^2{\Sigma }_{3,3}[5]\\ \beta =T[0,0]T[0,1]{\Sigma }_{3\times 3}[0]+T[1,0]T[0,1]{\Sigma }_{3\times 3}[1]+T[2,0]T[0,1]{\Sigma }_{3\times 3}[3]\\ +T[0,0]T[1,1]{\Sigma }_{3\times 3}[1]+T[1,0]T[1,1]{\Sigma }_{3\times 3}[2]+T[2,0]T[2,1]{\Sigma }_{3\times 3}[4]\\ +T[0,0]T[2,1]{\Sigma }_{3\times 3}[3]+T[1,0]T[2,1]{\Sigma }_{3\times 3}[4]+T[2,0]T[2,1]{\Sigma }_{3\times 3}[5]\\ \gamma =T[0,1{]}^2{\Sigma }_{3,3}[0]+T[1,1{]}^2{\Sigma }_{3,3}[1]+T[2,1{]}^2{\Sigma }_{3\times 3}[3]\\ +T[0,1{]}^2{\Sigma }_{3,3}[1]+T[1,1{]}^2{\Sigma }_{3,3}[2]+T[2,1{]}^2{\Sigma }_{3,3}[4]\\ +T[0,1{]}^2{\Sigma }_{3,3}[3]+T[1,1{]}^2{\Sigma }_{3,3}[4]+T[2,1{]}^2{\Sigma }_{3,3}[5]\end{array} \end{gathered}$$

$$\begin{gathered} \nabla {\Sigma }_{3\times 3}[0]=\nabla \alpha \cdot \frac{\partial \alpha }{\partial {\Sigma }_{3\times 3}[0]}+\nabla \beta \cdot \frac{\partial \beta }{\partial {\Sigma }_{3\times 3}[0]}+\nabla \gamma \cdot \frac{\partial \gamma }{\partial {\Sigma }_{3\times 3}[0]} \\ \Rightarrow \end{gathered}$$
$$\begin{gathered} \nabla {\Sigma }_{3\times 3}[0]=[T[00{]}^2,T[00]T[01],T[01{]}^2]\cdot [\nabla \alpha ,\nabla \beta ,\nabla \gamma ] \\ \nabla {\Sigma }_{3\times 3}[1]=[T[00{]}^2+T[10{]}^2,T[10]T[01]+T[00]T[11],T[01{]}^2+T[11{]}^2]\cdot [\nabla \alpha ,\nabla \beta ,\nabla \gamma ] \end{gathered}$$
剩下的一样的方式自己推吧，此处省略；

• $\nabla {\Sigma }_{3\times 3}\Rightarrow (\nabla R,\nabla S)\Rightarrow (\nabla q,\nabla s)$

  • $\nabla {\Sigma }_{3\times 3}\Rightarrow (\nabla R,\nabla S)$
    $$\left[\begin{array}{ccc}{\Sigma }_{3\times 3}[0]& {\Sigma }_{3\times 3}[1]& {\Sigma }_{3\times 3}[3]\\ {\Sigma }_{3\times 3}[1]& {\Sigma }_{3\times 3}[2]& {\Sigma }_{3\times 3}[4]\\ {\Sigma }_{3\times 3}[3]& {\Sigma }_{3\times 3}[4]& {\Sigma }_{3\times 3}[5]\end{array}\right]=RS(RS{)}^T$$

    $$\Rightarrow$$

    $$\begin{gathered} {\Sigma }_{3\times 3}[0]=R_{00}^2{S}_0^2+R_{01}^2{S}_1^2+R_{02}^2{S}_2^2 \\ {\Sigma }_{3\times 3}[1]=R_{00}{R}_{10}{S}_0^2+R_{01}{R}_{11}{S}_1^2+R_{02}{R}_{12}{S}_2^2 \\ {\Sigma }_{3\times 3}[2]=R_{10}^2{S}_0^2+R_{11}^2{S}_1^2+R_{12}^2{S}_2^2 \\ {\Sigma }_{3\times 3}[3]=R_{00}{R}_{20}{S}_0^2+R_{01}{R}_{21}{S}_1^2+R_{02}{R}_{22}{S}_2^2 \\ {\Sigma }_{3\times 3}[4]=R_{10}{R}_{20}{S}_0^2+R_{11}{R}_{21}{S}_1^2+R_{12}{R}_{22}{S}_2^2 \\ {\Sigma }_{3\times 3}[5]=R_{20}^2{S}_0^2+R_{21}^2{S}_1^2+R_{22}^2{S}_2^2 \end{gathered}$$
    $$\Rightarrow$$

  • $\nabla {\sigma }_{3\times 3}:=[\nabla {\Sigma }_{3\times 3}[0],...,\nabla {\Sigma }_{3\times 3}[5]]$

    $$\begin{gathered} \nabla R_{00}=[2S_0^2{R}_{00},R_{10}{S}_0^2,0,R_{20}{S}_0^2,0,0]\cdot \nabla {\sigma }_{3\times 3} \\ \nabla R_{01}=[2S_1^2{R}_{01},R_{11}{S}_1^2,0,R_{21}{S}_1^2,0,0]\cdot \nabla {\sigma }_{3\times 3} \\ \nabla R_{02}=[2S_2^2{R}_{02},R_{12}{S}_2^2,0,R_{22}{S}_2^2,0,0]\cdot \nabla {\sigma }_{3\times 3} \\ \nabla R_{10}=[0,R_{00}{S}_0^2,0,2R_{10}{S}_0^2,R_{20}{S}_0^2,0]\cdot \nabla {\sigma }_{3\times 3} \end{gathered}$$
    剩下的省略，一样的推导方式；
    $$\Rightarrow$$
    $$\begin{gathered} \nabla S_0=[2R_{00}^2{S}_0,2R_{00}{R}_{10}{S}_0,2R_{10}^2{S}_0,2R_{00}{R}_{20}{S}_0,2R_{10}{R}_{20}{S}_0,2R_{20}^2{S}_0]\cdot \nabla {\sigma }_{3\times 3} \\ \nabla S_1=[2R_{01}^2{S}_1,2R_{01}{R}_{11}{S}_1,2R_{11}^2{S}_1,2R_{01}{R}_{21}{S}_1,2R_{11}{R}_{21}{S}_1,2R_{21}^2{S}_1]\cdot \nabla {\sigma }_{3\times 3} \\ \nabla S_2=[2R_{02}^2{S}_2,2R_{02}{R}_{12}{S}_2,2R_{12}^2{S}_2,2R_{02}{R}_{22}{S}_2,2R_{12}{R}_{22}{S}_2,2R_{22}^2{S}_2]\cdot \nabla {\sigma }_{3\times 3} \end{gathered}$$

  • $(\nabla R,\nabla S)\Rightarrow (\nabla q,\nabla s)$
    $$\begin{gathered} q={\left[\begin{array}{cccc}w& x& y& z\end{array}\right]}^T/\sqrt{w^2+x^2+y^2+z^2} \\ R=\left[\begin{array}{ccc}1-2q_2^2-2q_3^2& 2q_1{q}_2-2q_0{q}_3& 2q_1{q}_3+2q_0{q}_2\\ 2q_1{q}_2+2q_0{q}_3& 1-2q_1^2-2q_3^2& 2q_2{q}_3-2q_0{q}_1\\ 2q_1{q}_3-2q_0{q}_2& 2q_2{q}_3+2q_0{q}_1& 1-2q_1^2-2q_2^2\end{array}\right] \\ \Rightarrow \\ \nabla q_0=-2q_3\nabla R_{01}+2q_2\nabla R_{02}+2q_3\nabla R_{10}-2q_2\nabla R_{20} \end{gathered}$$
    …
    $$\Rightarrow$$
    $$\begin{gathered} \nabla w=\nabla q_0\frac{\partial q_0}{\partial w}+\nabla q_1\frac{\partial q_1}{\partial w}+\nabla q_2\frac{\partial q_2}{\partial w}+\nabla q_3\frac{\partial q_3}{\partial w} \\ =\nabla q_0\frac{x^2+y^2+z^2}{(w^2+x^2+y^2+z^2{)}^{1.5})}+\nabla q_1\frac{-x^2}{(w^2+x^2+y^2+z^2{)}^{1.5})}+\nabla q_2\frac{-y^2}{(w^2+x^2+y^2+z^2{)}^{1.5})}+\nabla q_3\frac{-z^2}{(w^2+x^2+y^2+z^2{)}^{1.5})} \end{gathered}$$

• 代码是这样的：

  • backward.cu/preprocess_CUDA(...): $\nabla u{v}_{proj}\Rightarrow \nabla xyz$

    // diff-gaussian-rasterization/cuda_rasterizer/backward.cu
    // preprocessCUDA(...)
    float3 m = means[idx];
    // Taking care of gradients from the screenspace points
    float4 m_hom = transformPoint4x4(m, proj);
    float m_w = 1.0f / (m_hom.w + 0.0000001f);
    // 
    glm::vec3 dL_dmean;
    float mul1 = (proj[0] * m.x + proj[4] * m.y + proj[8] * m.z + proj[12]) * m_w * m_w;
    float mul2 s= (proj[1] * m.x + proj[5] * m.y + proj[9] * m.z + proj[13]) * m_w * m_w;
    dL_dmean.x = (proj[0] * m_w - proj[3] * mul1) * dL_dmean2D[idx].x + (proj[1] * m_w - proj[3] * mul2) * dL_dmean2D[idx].y;
    dL_dmean.y = (proj[4] * m_w - proj[7] * mul1) * dL_dmean2D[idx].x + (proj[5] * m_w - proj[7] * mul2) * dL_dmean2D[idx].y;
    dL_dmean.z = (proj[8] * m_w - proj[11] * mul1) * dL_dmean2D[idx].x + (proj[9] * m_w - proj[11] * mul2) * dL_dmean2D[idx].y;
      - `backward.cu/computeCov2DCUDA`: $$\nabla (\Sigma_{2\times 2})^{-1} \Rightarrow \nabla \Sigma_{3\times 3}$$
    // diff-gaussian-rasterization/cuda_rasterizer/backward.cu
    // computeCov2DCUDA(...)
    
    // \nabla(\alpha, \alpha\beta, \alpha \gammma) → \nabla (a, b, c) 
    float a = cov2D[0][0] += 0.3f;
    float b = cov2D[0][1];
    float c = cov2D[1][1] += 0.3f;
    float denom = a * c - b * b;
    float dL_da = 0, dL_db = 0, dL_dc = 0;
    float denom2inv = 1.0f / ((denom * denom) + 0.0000001f);
    dL_da = denom2inv * (-c * c * dL_dconic.x + 2 * b * c * dL_dconic.y + (denom - a * c) * dL_dconic.z);
    dL_dc = denom2inv * (-a * a * dL_dconic.z + 2 * a * b * dL_dconic.y + (denom - a * c) * dL_dconic.x);
    dL_db = denom2inv * 2 * (b * c * dL_dconic.x - (denom + 2 * b * b) * dL_dconic.y + a * b * dL_dconic.z);
    
    // \nabla(a, b, c) → \nabla(Σ_{3,3})
    dL_dcov[6 * idx + 0] = (T[0][0] * T[0][0] * dL_da + T[0][0] * T[1][0] * dL_db + T[1][0] * T[1][0] * dL_dc);
    dL_dcov[6 * idx + 3] = (T[0][1] * T[0][1] * dL_da + T[0][1] * T[1][1] * dL_db + T[1][1] * T[1][1] * dL_dc);
    dL_dcov[6 * idx + 5] = (T[0][2] * T[0][2] * dL_da + T[0][2] * T[1][2] * dL_db + T[1][2] * T[1][2] * dL_dc);
    
    dL_dcov[6 * idx + 1] = 2 * T[0][0] * T[0][1] * dL_da + (T[0][0] * T[1][1] + T[0][1] * T[1][0]) * dL_db + 2 * T[1][0] * T[1][1] * dL_dc;
    dL_dcov[6 * idx + 2] = 2 * T[0][0] * T[0][2] * dL_da + (T[0][0] * T[1][2] + T[0][2] * T[1][0]) * dL_db + 2 * T[1][0] * T[1][2] * dL_dc;
    dL_dcov[6 * idx + 4] = 2 * T[0][2] * T[0][1] * dL_da + (T[0][1] * T[1][2] + T[0][2] * T[1][1]) * dL_db + 2 * T[1][1] * T[1][2] * dL_dc;
    

  • backward.cu/computeCov3D: $$\nabla {\Sigma }_{3\times 3}\Rightarrow (\nabla R,\nabla S)\Rightarrow (\nabla q,\nabla s)$$

    • dL_dMt就是:
      $$\left[\begin{array}{ccc}2{\Sigma }_{3\times 3}[0]& {\Sigma }_{3\times 3}[1]& {\Sigma }_{3\times 3}[2]\\ {\Sigma }_{3\times 3}[1]& 2{\Sigma }_{3\times 3}[2]& {\Sigma }_{3\times 3}[3]\\ {\Sigma }_{3\times 3}[3]& {\Sigma }_{3\times 3}[4]& 2{\Sigma }_{3\times 3}[5]\end{array}\right]R^{T}S$$

    按照代码：
    $$\begin{gathered} \nabla S_0=[R_{00},R_{10},R_{20}]\cdot \\ [2{\Sigma }_{3\times 3}[0]R_{00}{S}_0+{\Sigma }_{3\times 3}[1]R_{01}{S}_0+{\Sigma }_{3\times 3}[3]R_{02}{S}_0, \\ {\Sigma }_{3\times 3}[1]R_{00}{S}_0+2{\Sigma }_{3\times 3}[2]R_{01}{S}_0+{\Sigma }_{3\times 3}[4]R_{02}{S}_0, \\ {\Sigma }_{3\times 3}[3]R_{00}{S}_0+{\Sigma }_{3\times 3}[4]R_{01}{S}_0+2{\Sigma }_{3\times 3}[5]R_{02}{S}_0] \\ =[2R_{00}^2{S}_0,R_{01}{R}_{10}{S}_0+R_{00}{R}_{11}{S}_0,2R_{01}{R}_{10}{S}_0,R_{00}{R}_{02}{S}_0+R_{00}{R}_{20}{S}_0,...]\cdot \nabla {\sigma }_{3\times 3} \end{gathered}$$

    // diff-gaussian-rasterization/cuda_rasterizer/backward.cu
    // computeCov3D(...)
    
    glm::mat3 dL_dSigma = glm::mat3(
        dL_dcov3D[0], 0.5f * dL_dcov3D[1], 0.5f * dL_dcov3D[2],
        0.5f * dL_dcov3D[1], dL_dcov3D[3], 0.5f * dL_dcov3D[4],
        0.5f * dL_dcov3D[2], 0.5f * dL_dcov3D[4], dL_dcov3D[5]
    );
    glm::mat3 dL_dM = 2.0f * M * dL_dSigma;
    glm::mat3 Rt = glm::transpose(R);
    glm::mat3 dL_dMt = glm::transpose(dL_dM);
    
    
    // \nabla q
    glm::vec4 dL_dq;
    dL_dq.x = 2 * z * (dL_dMt[0][1] - dL_dMt[1][0]) + 2 * y * (dL_dMt[2][0] - dL_dMt[0][2]) + 2 * x * (dL_dMt[1][2] - dL_dMt[2][1]);
    dL_dq.y = 2 * y * (dL_dMt[1][0] + dL_dMt[0][1]) + 2 * z * (dL_dMt[2][0] + dL_dMt[0][2]) + 2 * r * (dL_dMt[1][2] - dL_dMt[2][1]) - 4 * x * (dL_dMt[2][2] + dL_dMt[1][1]);
    dL_dq.z = 2 * x * (dL_dMt[1][0] + dL_dMt[0][1]) + 2 * r * (dL_dMt[2][0] - dL_dMt[0][2]) + 2 * z * (dL_dMt[1][2] + dL_dMt[2][1]) - 4 * y * (dL_dMt[2][2] + dL_dMt[0][0]);
    dL_dq.w = 2 * r * (dL_dMt[0][1] - dL_dMt[1][0]) + 2 * x * (dL_dMt[2][0] + dL_dMt[0][2]) + 2 * y * (dL_dMt[1][2] + dL_dMt[2][1]) - 4 * z * (dL_dMt[1][1] + dL_dMt[0][0]);
    
    
    // \nabla S.
    glm::vec3* dL_dscale = dL_dscales + idx;
    dL_dscale->x = glm::dot(Rt[0], dL_dMt[0]);
    dL_dscale->y = glm::dot(Rt[1], dL_dMt[1]);
    dL_dscale->z = glm::dot(Rt[2], dL_dMt[2]);
    

Reference

• [1] 矩阵求导：https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf
• [2] https://dyn4j.org/2010/01/sat/
• [3] https://github.com/joeyan/gaussian_splatting/blob/main/MATH.md
• [4] 俺把带注释的cuda代码存放到了: https://vzia2ov9de.feishu.cn/file/L5RFbdLEpotN3ExPCSWcHo9En6y?from=from_copylink