【三维重建】CamP:针对NeRF的摄像机预处理

题目：CamP: Camera Preconditioning for Neural Radiance Fields
来源：KEUNHONG PARK等, Google Research, USA
文章：https://arxiv.org/pdf/2308.10902.pdf
项目：https://github.com/jonbarron/camp_zipnerf

---

---

摘要

  优化后的NeRF（神经辐射场）可以实现物体和大规模场景的高保真三维场景重建，却需要准确的相机参数作为输入（不准确的相机参数会导致模糊的渲染）。相机参数通常估计使用SfM（运动结构）作为预处理步骤（总是存在估计偏差）。先前的工作提出了与NeRF共同优化摄像机参数，但这些方法在具有挑战性的设置中容易出现局部最小值。在这项工作中，我们分析了不同的相机参数如何影响这个联合优化问题，并观察到标准参数对于小扰动表现出很大的幅度差异，这可能导致病态优化问题。我们建议 使用代理问题来计算白化变换（一种统计技术），消除相机参数之间的相关性并将其影响归一化 ，并建议在联合优化过程中使用该变换作为相机参数的preconditioner(预调节器)
  预处理相机优化，在Mip-NeRF 360数据集上显著改善重建质量，相比不优化相机的Zip-NeRF，减少错误率（RMSE）67%，和29%的使用SCNeRF的先进的联合优化方法。方法易于实现，不会显著增加运行时间。
  

一、前言

  NeRF的一个限制因素是假设摄像机参数是已知的。通常，使用结构从运动（SfM）方法，如COLMAP来恢复像机参数。但在给定稀疏或宽基线视图时经常失败。
  鉴于NeRF是可微的图像形成模型，可以通过将从损失到摄像机参数的反向传播梯度来优化摄像机和场景参数。早期的摄像机优化NeRF方法，假设未知的摄像机pose，和已知的内参intrinsics，并通过将每个摄像机分配给 identity pose 来初始化该联合优化。

  场景和摄像机参数的联合恢复是一个容易出现局部极小值的病态问题。这 可以通过一种从粗到细的方法得到部分改善。对于mlp，BARF [Lin等2021]方法中，位置编码的更高频带can be windowed；对于基于网格的方法（如，即时NGP )，可以降低更高分辨率层的影响。虽然这些从粗到细的方法是有效的，但仍需要初始化时大致精确的相机姿态。

  NeRF中的相机参数化，经常被忽视。例如，BARF表示pose为𝔰𝔢(3)，即三维欧式空间的6D螺旋轴参数化，而SCNeRF使用基于Gram-Schmidt的6D方向表示，加上3D平移。选择哪个坐标系来表示平移也很重要，三维物体姿态估计的工作研究了egocentric(相对于相机）和allocentric的表示（相对于世界）以及焦距和translation之间的相互作用。此外，获得良好的相机内参估计（焦距、主点和失真系数）是精确重建的关键，而SfM系统通常估计它们包括径向和切向失真系数。内参可以使用相同的反向传播技术进行优化，如SCNeRF中。

  考虑到可能的摄像机参数化的大空间，场景重建会影响对摄像机参数的估计，反之亦然，而神经场景表示的优化动态（如NeRF中的mlp）的建模具有挑战性。相反，我们分析了相机参数化对一个更简单的代理设置的影响：相机参数化如何影响它前面的点的投影？尽管很简单，这个代理问题仍然产生深刻的见解；它证明了单纯地实现的相机参数化可能是病态的，其中相机参数化的一个轴上的扰动，对投影的影响明显大于其他轴。

  尽管这些条件问题可以在零星的基础上解决（例如，通过手动为相机参数化的每个维度引入缩放超参数），但仍然可以在参数之间留下虚假的相关性，是导致问题的不良条件的本质。我们从数值方法中获得灵感，并为优化问题中为每个摄像机设计一个preconditioner预调节器，使每个参数对场景中点的投影的影响归一化，并解耦每个参数与所有其他参数的影响。 这个预调节器是一个𝑘×𝑘矩阵（其中𝑘是相机参数的数量），在将相机参数传递到NeRF模型之前，它会应用到相机参数上 。我们称这种技术为相机预处理（CamP）。

  论文实验基于最先进的ZipNeRF（使用了具有多分辨率散列网格的抗混叠采样策略）。数据集采用Mip-NeRF 360（perturbed versions）和 synthetic NeRF scenes。最后，证明CamP也改进了ARKit pose的手机序列重建质量。

二、相关工作

  3D Reconstruction。 从图像中恢复一个场景是计算机视觉中一个长期存在的问题。运动结构（SfM）是一种中心技术，其中稀疏关键点匹配是三角化的，并调整束以最小化重投影指标[schon伯格和Frahm 2016]。最近，基于学习的技术已被用于恢复更准确的表示，如基于渲染代理和生成最终图像的实时卷积神经网络的方法[Meshry等，2019；Thies等，2019]，或混合权重[Hedman等，2018]。神经辐射场（NeRF）[米尔登霍尔等人，2020]用MLP的权重对一个场景进行编码，并使用体积渲染来生成一个被捕获场景的高度真实的新颖视图。NeRFs已扩展到多种任务，如互联网照片收集[马丁-布鲁alla等2021]，动态场景[李等2020；Park等2021]，以及照明估计。

  Camera-optimizing NeRFs。 可区分的渲染方法，如在NeRF中使用的体渲染，允许通过场景表示进行反向传播，以更新照相机参数。虽然这些已经在网格和栅格化的背景下进行了探索[Liu等人，2019年；Loper和Black，2014]，但我们关注的是适用于nerf类模型的方法。NeRF−−[Wang et al. 2021b]表明，通过初始化摄像机到origin，可以估计正面场景的pose。BARF [Lin等人2021]采用从粗到细的重建方案[Hertz等人2021；Park等人2021]，逐步引入更高频率的位置编码。通过使用不同分辨率水平的加权方法，精细方法可以适用于多分辨率网格，如NGP [Muller等人，2022]，[Heo等人，2023；Melas-Kyriazi等人，2023]。GNeRF[孟等人. 2021]提出了一种姿态条件的GAN，用于恢复NeRF。其他人已经探索了更适合的架构来执行姿态优化，如高斯[Chng等人2022]或正弦激活[Xia等人2022]。NoPe-NeRF [Bian et al. 2022]使用单眼深度先验来约束场景和相对姿态估计。关键点匹配或密集对应也可用于约束使用射线到射线到射线对应损失的相对姿态估计[Jeong等人2021；Truong等人2023]。DBARF [Chen和Lee 2023]建议使用低频特征图来指导可推广的nerf的束调整[Wang等人2021a；Yu等人2021]

   虽然许多方法忽略了内参的优化，但这是实现高质量重建的关键。SCNeRF [Jeong等人，2021]建模了残差投影矩阵和残差网格参数[格罗斯伯格和Nayar 2005]，并在次采样像素网格上进行插值[Zhou和Koltun 2014]。

   通过使用NeRFs以在线训练的方式积累场景观察，NeRFs也被适应于实时SLAM设置[Rosinol等人2022；Zhu等人2022]。LocalRF [Meuleman等人2023]不依赖SfM，通过使用增量策略，以及类似于BlockNeRF的细分策略[Tancik等人2022]。值得注意的是，LocalRF对translation 和orientation参数使用不同的学习率，突出了服从相机优化的NeRF挑战。

  Preconditioning。 许多计算机图形学任务依赖于解决大规模线性系统或优化非线性目标，如SfM中的bundle adjustment。线性系统的预条件方法能够更快的收敛，但找到一个任务的最佳预调节器是困难的，因为人们必须权衡预调节器的优点和计算它所需的时间。值得注意的是，预调节器可以使用给定任务的领域知识来设计，例如SfM中的匹配图的可见性结构（ visibility-structure of a matching graph）[Kushal和Agarwal 2012] [Szeliski2022, A.5.2]。

   在随机梯度下降（SGD）的情况下，预处理等价于梯度缩放[De Sa 2019]。最近在NeRF优化的背景下提出了Gradient scaling，以减少靠近摄像机的漂浮物，因为这些区域相对于背景被过采样 [Philip and Deschaintre 2023].

  Camera Parameterizations。 对于表示图像形成模型的摄像机参数化的选择，在预测或优化这些参数的应用程序中可能具有显著的影响。例如，众所周知，一些旋转的参数化在SO (3)流形上不是连续的或可微的（例如，欧拉角）。此外，三维重建任务存在歧义，如测量模糊，误差不变为相似变换，平面场景的深度和焦距之间的歧义，以及径向失真和相机位置之间的复杂歧义[Wu 2014]。这些歧义可能会导致优化过程中的不稳定。最后，许多图像形成模型，如传统的透镜失真模型，都是对更复杂现象的低阶近似，如果有足够的数据，可能会明显欠拟合信号[Schops et al. 2020]。

   已经有人提出了一些工作来解决这些问题。Zhou等人[2019]提出了基于Gram-Schmidt 正交化的6D连续参数化。早期的三维姿态估计工作[Kundu等人，2018；Li等人，2018]发现，高度相关的旋转和平移参数会导致较差的结果，并提出了以对象为中心的旋转和平移表示，使它们去关联。波尼马特金等人[2022]进一步扩展到联合姿态和焦距估计。

  我们的工作，首先通过测量摄像机表示的参数空间对场景中三维点投影的影响来分析摄像机表示的参数空间的结构，然后提出一种摄像机优化的预处理方案，在优化过程中对参数进行归一化和去关联。
we first analyze the structure of the parameter space of camera representations by measuring their effects on the projections of 3D points in the scene, and then propose a preconditioning scheme for camera optimization to normalize and decorrelate the parameters during optimization.

三、准备工作

  

3.1 相机优化 NeRFs

  给出𝑛张图像{𝐼_𝑖}和初始相机参数 𝝓_i⁰，其中𝑖∈{1,……,𝑛}。 目标是获得一个优化好的NeRF参数𝜽^∗ 和优化的相机参数𝝓_𝑖^∗的场景模型𝑀_𝜽，来渲染不同视角图像 。NeRF是通过最小化 数据集𝐷上的损失函数L(𝐷;𝜽;𝝓）来训练的，其中像素𝑗有一个相机索引𝑑_𝑗∈{1，…，𝑛}，一个像素位置p_𝑗，和一个RGB颜色c_𝑗。

  类似于非线性最小二乘方法，我们对> 相机参数 𝝓=𝝓⁰+Δ𝝓 的线性化问题进行了优化，通过优化初始值𝝓⁰的残差Δ𝝓 。将相机参数的维度表示为𝑘，即𝝓_𝑖∈R^k。

  相机参数化定义了 一个投影函数$\Pi$，其中 $\Pi$(x；𝝓):R³→R²将一个3D点x投影到一个2D像素位置p 。对应的逆像素映射到点$\Pi$^-1(p,𝑑;𝝓):R²×R→R³，其中𝑑是像素深度。

3.2 Preconditioning

  预处理1[Szeliski2022，a.5.2]是一种使用预处理器来转换优化问题，以具有更好的优化特性( 如one with a lower condition number )。例如，当优化一个值u∈R^𝑘 以最小化𝑓(u)： R^𝑛→R时，我们可以通过计算v=Pu（其中P∈R^𝑘×𝑘是预处理矩阵）来优化 𝑓(P⁻¹v）。一旦找到一个最优的v^∗，就可以得到u的最优值为 u^∗=P⁻¹v^∗。理想情况下，P⁻¹近似于逆的Hessian，从而得到了牛顿-拉夫森方法。

  NeRF通常使用基于梯度的优化器。我们可以通过优化变量 $\tilde{a}$=Pa的预处理版本，并修改目标为 P⁻¹$\tilde{a}$ 的函数，而不是a。如果SGD是优化器，那么梯度缩放技术（如Philip和Deschaintre [2023]）就相当于预处理。

四、分析相机的参数化

   摄像机参数化的选择对输出重建质量有显著影响 。SfM方法，如COLMAP支持透视图或鱼眼相机和不同的镜头失真模型，包括具有径向和失真参数的模型。当使用真实世界的图像时，建模这些效果是至关重要的

  有两种方法可以考虑相机的外参pose：

(i) 作为一个egocentric 的参数化，建模对相机本身的像机pose，
(ii）作为一个allocentric 参数化[Kundu et al. 2018]，根据世界或对象坐标系为相机建模。

刚性6D位姿参数化的具体选择带来了进一步的复杂性：可以使用四元数、 axis-angle (Rodrigues’ formula [Rodrigues 1816]) 或6D连续位姿表示[Zhou et al. 2019] 来表示方向，或𝔰𝔢(3) 的screw axis表示

  受Kundu等人[2018]；Li等人[2018]的启发，我们 可视化了干扰相机模型的参数对投影图像的影响 。考虑一个使用𝔰𝔢(3)螺杆轴参数化S =（r；v）∈R⁶ 来编码rigid pose的相机，并为焦距增加了一个附加的偏置Δ𝑓。图3a显示了渲染对象的motion trails(运动轨迹)，其中摄像机通过干扰参数化的每个轴相同的大小来改变，结果帧被一起平均。请注意， 某些参数对渲染图像的影响更大，例如，与平移相比，图像对焦距扰动的敏感性较低。当考虑到焦距是像素单位，而平移是世界单位时，原因是明显的。

  当使用梯度下降方法进行优化时，这种灵敏度上的差异是有问题的。如果学习速率大到足以影响焦距，其他参数就会overshoot和振荡。然而，如果学习率较小，焦距将优化得太慢。自适应梯度方法，如Adam，通过将更新的大小归一化到类似的范围来改善这一点，但可能忽略参数之间的差异：焦距的更新应该远远大于位置，但Adam的参数更新都有相似的大小。

五、摄像机预处理

  在优化设置中减轻大小差异的一种常见方法，是重新参数化问题。最简单的例子是将参数乘以一个标量，例如，焦距参数化的𝑓=𝑓0+Δ𝑓可以更改为𝑓=𝑓0+𝑠Δ𝑓，其中𝑠是一个大标量。这一变化等效于用等于恒等式的对角预处理器P⁻¹预处理优化问题，除了对角线上对应于焦距的条目，其值为s.（ 原文：This change is equivalent to preconditioning the optimization problem with diagonal preconditioner P⁻¹ equal to the identity, except for the entry on the diagonal corresponding to the focal length, with a value of s）虽然有效，但每个参数的缩放不能解耦高度相关的系数。我们研究了相机投影函数 $\Pi$(x；𝝓) 对输入参数的敏感性。考虑投影函数的一个增强版本： $\Pi$^m(𝝓): R^𝑘→R^2𝑚，它是在𝑚个点 {x_𝑗}^m _𝑗=1上验证的$\Pi$的拼接。这组点可以被认为是被重建的三维场景内容的代理；我们想测量这些点的二维投影，如何作为每个摄像机参数的函数而变化。

  摄像机参数空间的任意点𝝓⁰处，每个摄像机参数对每个投影点p_𝑗=$\Pi$(x_𝑗; 𝝓）的影响，可以用雅可比矩阵表示：

矩阵Σ_Π = J_Π^⊤J_Π ∈ R^𝑘×𝑘 的第𝑟𝑠个条目等于:

在对角线上，这告诉我们由改变第𝑟个摄像机参数𝝓[𝑟]引起的平均运动幅度，而在对角线上，这告诉我们摄像机参数𝝓[𝑟]和𝝓[𝑠]之间的运动有多密切相关。On the diagonal this tells us the average motion magnitude induced by varying the 𝑟-th camera parameter 𝝓[𝑟], and off the diagonal this tells how closely correlated motion is between camera parameters 𝝓[𝑟] and 𝝓[𝑠].

这也被称为零分量分析（ZCA）白化变换。

  图2提供完整方法概述。图3的可视化过程，注意预调节器如何增加焦距的幅度（第一个参数），同时减少𝑦-平移的幅度（第六个参数）的幅度。

   预调节器仅对于初始估计摄像机参数𝝓⁰ 的代理投影目标是最优的。但在优化过程中，相机估计偏离𝝓⁰，达到次优。我们可以考虑使用更新的最优预处理器动态地重新参数化预处理相机参数，但这将使Adam [Kingma和Ba 2014]等优化器中使用的梯度动量的平均值不准确。

  以上分析假设照相机是相互独立地参数化的。但同一场景img是用相同的相机拍摄的，因此在优化过程中共享内参。我们通过施加一个额外的损失项L_share，以最小化共享的内参方差来处理。

  在实际应用中，我们在协方差矩阵ΣΠ的对角线上增加一个小的衰减参数，为了更好地在取逆矩阵平方根之前适应矩阵

其中，𝜆和𝜇为超参数。请注意，当包括透镜失真参数时，这一点尤其重要，因为非常小的变化可能会极大地影响点投影。

六、实验

6.1 实施细节

  模型参数。实验结合ZipNeRF重建方法，采用从粗到细的策略，逐步增加了更高分辨率的NGP的贡献。通过调整学习率，而不是调整特征来均衡网格的贡献。采用Adam优化器。相机参数，学习率从10⁻³衰减到10⁻⁴。在8个NVIDIA V100GPU上训练25,000步，训练时间约为2小时。预处理：超参数𝜆=10⁻¹和𝜇=10⁻⁸。对于具有共享内参的场景（360，扰动360），将焦距的方差损失加权为10⁻¹，主要点加权为10⁻²，径向失真加权为10⁻²。 对于NeRF-Synthetic数据集，不使用这些损失。

  相机参数化。所有的相机参数化都被实现为残差Δ𝝓_𝑖，并应用于初始相机参数𝝓_𝑖⁰。具体相机参数：

SCNeRF：实验使用SCNeRF的参数化：6D旋转、3D平移（𝑡_𝑥，𝑡_𝑦，𝑡_𝑧）、焦距（𝑓_𝑢，𝑓_𝑣）、主点（𝑢₀，𝑣₀）和两个径向扭曲系数（𝑘1，𝑘2）。我们省略了光线参数，也不使用在SCNeRF中提出的投影射线损失。

SE3：使用BARF中的 𝔰𝔢(3)相机参数化。这对应于BARF的一个改进版本，它已经适应了instant-NGP设置。

SE3+focal+intrinsic：我们用额外的焦距、主点和两个径向透镜畸变参数来增加SE3公式。焦距残差被参数化为一个乘法对数尺度，这样𝑓‘=𝑓expΔ𝑓。

FocalPose+Intrinsics：我们使用FocalPose中提出的 joint pose和焦距参数化（以对象为中心的表示而设计）。与SE3参数化类似，我们添加了principal point（主点）和两个径向透镜失真参数。我们注意到FocalPose被提出用于迭代姿态估计，因此我们将其重新参数化为关于初始估计的残差，而不是之前的估计。

   mip-NeRF 360：数据集有4个室内（分辨率1560×1040）和5个室外（分辨率1228×816）场景，逐场景100~300视图。

   扰动-360：为了评估对相机误差的鲁棒性，分别用N(0,0.005) 和 N(0,0.005) 通过观察相机的视点和位置来扰动相机。用随机尺度 exp N(0,log(1.02))来扰动焦距，并用相同的随机尺度因子exp N(0,log(1.05))来扰动到原点的距离（场景中心）和焦距。最后，所有摄像机的初始失真参数设为零。

   Perturbed-Synthetic：数据集包含8个对象（100个视点训练，200个不同的视点用于计算测试）。用 N(0,0.1) 来干扰相机的观察点和位置。我们用随机尺度因子 expN(0,log(1.2))对焦距进行扰动，用相同的随机尺度因子 expN(0,log(1.1))扰动到原点（场景中心）和焦距。我们对合成数据集禁用失真。请注意，我们对合成数据集使用了更强的扰动，因为对象与无界场景相比表现出更弱的视角失真。

  为了对测试摄像机进行定量评估，我遵循BARF协议，采用光度测试时间优化方法(photometric test-time optimization) 评估摄像机的参数。所有方法均使用带有预处理的 FocalPose+Intrinsics 参数化，并使用Adam优化，100步的学习率为0.001。

6.2 实验效果

  表1表明，相机优化的NeRF比基线有所改善。除了SE3参数化，我们的预处理方法的表现略优于非预处理方法。图1在360数据集上进行了定性比较，展示了在具有挑战性的户外场景的差异。

七.代码

7.1 环境搭建

# 克隆项目到本地
git clone https://github.com/jonbarron/camp_zipnerf.git
cd camp_zipnerf

# 创建环境
conda create --name camp_zipnerf python=3.11
conda activate camp_zipnerf

# Prepare pip（有pip的可以省略）
conda install pip
pip install --upgrade pip

# pip安装依赖库
pip install -r requirements.txt

# 手动安装 pycolmap  (don't use pip's! It's different).
git clone https://github.com/rmbrualla/pycolmap.git ./internal/pycolmap

# 确认所有单元测试无误.
./scripts/run_all_unit_tests.sh

7.2 训练

调用 train.py，参数为数据集路径(要求NeRF格式，带有colmap的初始pose)，和权重路径（可为空文件夹）

  renderings, ray_history = model.apply(
      variables,                                     # params
      key if config.randomized else None,
      rays,                                          # near,far,origins,viewdirs (16384,1,1,3)
      train_frac=train_frac,
      compute_extras=compute_extras,
      zero_glo=False,
  )

1. 射线采样

sdist = stepfun.sample_intervals(  key, sdist,  logits_resample, num_samples,   
        single_jitter=self.single_jitter,  domain=(init_s_near, init_s_far))                                                      # (10,65) 0~1之间

tdist = s_to_t(sdist)                                  # 映射到0~505.84之间

2. 将上步的distance intervals 转为高斯

   gaussians = render.cast_rays( tdist, rays.origins, rays.directions,  
                          rays.radii, self.ray_shape, diag=False )
   

展开：
directions：(10,3)array
([[-0.8701506 , -0.45856315, -0.18043813],
[-0.3565826 , -0.54477096, 0.758995 ],
[-0.36425403, 0.08462194, -0.92744714],
[-0.50708354, 0.00654202, -0.861872 ],
[-0.57497436, 0.8106414 , 0.11074794],
[-0.25367722, -0.70381755, -0.66354257],
[ 0.6983983 , 0.5845935 , -0.41290492],
[ 0.43995845, 0.03484576, 0.89734185],
[ 0.54235345, 0.15224242, -0.8262415 ],
[-0.58228904, 0.34765905, -0.73489624]

radii :(10,1) Array([[5.3690898e-04], [1.6620250e-04], [8.7686157e-04],
[8.1245729e-04], [7.6067407e-04], [2.1810175e-04], [9.7374397e-04],
[1.7704273e-04], [6.6734880e-05], [8.0434268e-04]

t0 = tdist[Ellipsis, :-1]        # (10,64)
t1 = tdist[Ellipsis, 1:]
if ray_shape == 'cone':
  gaussian_fn = conical_frustum_to_gaussian

means, covs = gaussian_fn(directions, t0, t1, radii, diag=False)

	def gaussianize_frustum(t0, t1):
		"""Convert intervals along a conical frustum into means and variances."""
		# A more stable version of Equation 7 from https://arxiv.org/abs/2103.13415.
		s = t0 + t1
		d = t1 - t0
		eps = jnp.finfo(jnp.float32).eps ** 2
		ratio = d**2 / jnp.maximum(eps, 3 * s**2 + d**2)
		t_mean = s * (1 / 2 + ratio)
		t_var = (1 / 12) * d**2 - (1 / 15) * ratio**2 * (12 * s**2 - d**2)
		r_var = (1 / 16) * s**2 + d**2 * (5 / 48 - (1 / 15) * ratio)
		return t_mean, t_var, r_var

r_var *= base_radius**2
mean, cov = lift_gaussian(d, t_mean, t_var, r_var, diag)  # (10,64,3)(10,64,3,3)

	def lift_gaussian(d, t_mean, t_var, r_var, diag):
		 """函数用于将沿着射线定义的高斯分布提升到3D坐标空间。它首先计算了均值mean，"""
		 """然后根据diag参数的不同进行不同的处理：如果diag为False，它会计算完整的"""
		 """协方差矩阵cov，其中包括t_var和r_var的所有协方差元素，以及与射线方向正交的元素."""
		 mean = d[Ellipsis, None, :] * t_mean[Ellipsis, None]
		
		 d_mag_sq = jnp.maximum(1e-10, jnp.sum(d**2, axis=-1, keepdims=True))
		
		 if diag=false：
		   d_outer = d[Ellipsis, :, None] * d[Ellipsis, None, :]
		   eye = jnp.eye(d.shape[-1])
		   null_outer = eye - d[Ellipsis, :, None] * (d / d_mag_sq)[Ellipsis, None, :]
		   t_cov = t_var[Ellipsis, None, None] * d_outer[Ellipsis, None, :, :]
		   xy_cov = r_var[Ellipsis, None, None] * null_outer[Ellipsis, None, :, :]
		   cov = t_cov + xy_cov
		   return mean, cov

3. predict_density

self.pos_basis_t = jnp.array(geopoly.generate_basis(self.basis_shape, self.basis_subdivisions)).T  # (3,21)
lifted_means, lifted_vars = coord.lift_and_diagonalize(means, covs, self.pos_basis_t)

lift_and_diagonalize 接受三个参数 mean、cov 和 basis，分别代表均值向量、协方差矩阵和基础矩阵。该函数的作用是将均值向量和协方差矩阵投影到给定的基础上，并对投影后的协方差矩阵进行对角化处理

 x = []                        # Encode input positions.                                                       
 x.append(
     self.posenc_feature_scale
     * coord.integrated_pos_enc(
         lifted_means,
         lifted_vars,
         self.min_deg_point,
         self.max_deg_point,
     )
 )

x = jnp.concatenate(x, axis=-1)                                 # (504)

inputs = x
# Evaluate network to produce the output density.
for i in range(self.net_depth):                                 # 8次
    x = density_dense_layer(self.net_width)(x)                           # (256)
    x = self.net_activation(x)
  if i % self.skip_layer == 0 and i > 0:                                 # skip = 4
    x = jnp.concatenate([x, inputs], axis=-1)

raw_density = final_density_dense_layer(1)(x)[Ellipsis, 0]
return raw_density, x

以下是具体函数，可忽略

def generate_basis(
    base_shape, angular_tesselation, remove_symmetries=True, eps=1e-4
):
  """Generates a 3D basis by tesselating a geometric polyhedron.

  Args:
    base_shape: string, the name of the starting polyhedron, must be either
      'tetrahedron', 'icosahedron' or 'octahedron'.
    angular_tesselation: int, the number of times to tesselate the polyhedron,
      must be >= 1 (a value of 1 is a no-op to the polyhedron).
    remove_symmetries: bool, if True then remove the symmetric basis columns,
      which is usually a good idea because otherwise projections onto the basis
      will have redundant negative copies of each other.
    eps: float, a small number used to determine symmetries.

  Returns:
    basis: a matrix with shape [3, n].
  """
  if base_shape == 'icosahedron':
    a = (np.sqrt(5) + 1) / 2
    verts = np.array([
        (-1, 0, a),
        (1, 0, a),
        (-1, 0, -a),
        (1, 0, -a),
        (0, a, 1),
        (0, a, -1),
        (0, -a, 1),
        (0, -a, -1),
        (a, 1, 0),
        (-a, 1, 0),
        (a, -1, 0),
        (-a, -1, 0),
    ]) / np.sqrt(a + 2)
    faces = np.array([
        (0, 4, 1),
        (0, 9, 4),
        (9, 5, 4),
        (4, 5, 8),
        (4, 8, 1),
        (8, 10, 1),
        (8, 3, 10),
        (5, 3, 8),
        (5, 2, 3),
        (2, 7, 3),
        (7, 10, 3),
        (7, 6, 10),
        (7, 11, 6),
        (11, 0, 6),
        (0, 1, 6),
        (6, 1, 10),
        (9, 0, 11),
        (9, 11, 2),
        (9, 2, 5),
        (7, 2, 11),
    ])
  elif base_shape == 'octahedron':
    verts = np.array(
        [(0, 0, -1), (0, 0, 1), (0, -1, 0), (0, 1, 0), (-1, 0, 0), (1, 0, 0)]
    )
    corners = np.array(list(itertools.product([-1, 1], repeat=3)))
    pairs = np.argwhere(compute_sq_dist(corners.T, verts.T) == 2)
    faces = np.sort(np.reshape(pairs[:, 1], [3, -1]).T, 1)
  else:
    raise ValueError(f'base_shape {base_shape} not supported')
  verts = tesselate_geodesic(verts, faces, angular_tesselation)

  if remove_symmetries:
    # Remove elements of `verts` that are reflections of each other.
    match = compute_sq_dist(verts.T, -verts.T) < eps
    verts = verts[~np.any(np.triu(match), axis=0), :]

  basis = verts[:, ::-1]
  return basis

4. render

raw_density = raw_density_flat.reshape(means.shape[:-1])          # (10,64)
x = x_flat.reshape(means.shape[:-1] + (x_flat.shape[-1],))        # (10,64,256)
raw_grad_density = raw_grad_density_flat.reshape(means.shape)     # (10,64,3)
normals = -ref_utils.l2_normalize(raw_grad_density)               # (10,64,3)

grad_pred = None
normals_pred = None
normals_to_use = normals

density = self.density_activation(raw_density + self.density_bias)  + -1.0   # (10,64)


rendering = render.volumetric_rendering(
    ray_results['rgb'],                                     # None ,因此该函数不执行
    weights,                                                # (16384,1,1,64)
    tdist,
    bg_rgbs,
    compute_extras,
    extras={
        k: v
        for k, v in ray_results.items()
        if k.startswith('normals') or k in extras_to_render
    },
    percentiles=percentiles,
)

# 5.计算 loss--------------------------------------------------------------------------------------
losses['data'], stats = compute_data_loss(batch, renderings, rays, config)
losses['interlevel'] = interlevel_loss(ray_history, config)

7.0 更新相机参数

主要函数：
latent_params = self.create_params(cameras)
包含以下两项内容，4个键值：
**IntrinsicCameraDelta.create_params(self, cameras),
**SE3WithFocalCameraDelta.create_params(self, cameras),

# 1.生成 principal_point_bias 与 radial_distortion_bias--------------------------------------------------
if self.use_principal_point:
  params['principal_point_bias'] = self.param(
      'principal_point_bias',
      jax.nn.initializers.zeros,
      (*cameras.shape, 2),
  )                                                             # param['principal_point_bias']:(88,2)
if self.use_radial_distortion:
  params['radial_distortion_bias'] = self.param(
      'radial_distortion_bias',
      jax.nn.initializers.zeros,
      (*cameras.shape, self.num_radial_distortion_coeffs),
  )                                                             # param['radial_distortion_bias']:(88,2)

# 2.生成 screw_axis_bias 与 log_focal_scale：------------------------------------------------------------
def create_params(self, cameras: jaxcam.Camera) -> chex.ArrayTree:
  params = {
      **SE3CameraDelta.create_params(self, cameras),            # param['screw_axis_bias']:(88,6)
  }
			def create_params(self, cameras: jaxcam.Camera) -> chex.ArrayTree:
			    return {
			        'screw_axis_bias': self.param(
			            'screw_axis_bias', jax.nn.initializers.zeros, (*cameras.shape, 6))}

  if self.use_log_scales:
    params['log_focal_scale'] = self.param(
        'log_focal_scale', jax.nn.initializers.zeros, cameras.shape   # params['log_focal_scale']:(88)


# 3.用以上的latent_params，跟预计算的hessian矩阵jtj交互---------------------------------------------------
flat_params, _ = jax.vmap(lambda p: jax.flatten_util.ravel_pytree(p)[0])(latent_params)   # (88,11)
# pytree 即latent_params在 JAX 中通常指代一个类似树状结构的数据结构，比如由多个参数组成的字典、列表、元组等。
# 这些参数可能是神经网络的权重、偏置等。因对输入的字典 pytree 中的每个参数进行扁平化操作，将其转换为一维数组。vmap 函数将扁平化操作映射到输入的所有参数上

jtj = self.variable('precondition', 'jtj', init_jtj, (len(cameras), flat_params.shape[1])  # (88,11,11))
      def init_jtj(_):
        return self.compute_approximate_hessian(camera_params, cameras, random.PRNGKey(0))
# 这段代码是用于计算相机参数的近似 Hessian 矩阵的
_compute_approximate_hessian 函数：
接受相机参数 camera_params、点云 points 和相机对象 camera 作为输入。
首先使用 jax.flatten_util.ravel_pytree 函数将相机参数扁平化，以便后续计算。
然后定义了一个内部函数 _project_points，用于计算 3D 到 2D 的投影。在这个函数中，首先将扁平化的相机参数还原成原始结构，然后通过调用 self.transform_camera 函数对相机进行变换，最后使用 jaxcam.project 函数进行投影计算。
调用 jax.jacfwd 函数计算 _project_points 函数的雅可比矩阵，即相对于相机参数的梯度。
对雅可比矩阵进行处理，忽略超出相机视野的点，并计算雅可比矩阵的乘积的转置矩阵。最后取所有点的平均值，得到近似的 Hessian 矩阵。
compute_approximate_hessian 函数：

接受相机参数 camera_params、相机对象列表 cameras 和 PRNGKey rng 作为输入。
首先使用 random.split 函数生成新的 PRNGKey，并将其分割为与相机数量相同的子 PRNGKey。
调用 jax.vmap 函数对 transform_camera 和 create_points 函数进行向量化映射，以处理多个相机和点云。
最后调用 v_compute_jtj 函数（即 _compute_approximate_hessian 函数的向量化版本）来计算每个相机的近似 Hessian 矩阵，并返回结果。
总的来说，这段代码主要是对相机参数进行变换和投影计算，然后计算这些计算的雅可比矩阵，并进一步处理得到近似的 Hessian 矩阵。

camera_params = self._compute_camera_params_from_latent( latent_params, jtj.value)

# 4.更新camera----------------------------------------------------------------------------------------
next_jtj = self.compute_approximate_hessian(camera_params, cameras, key)    # 输入：可学习变量+真实相机参数

        if self.use_log_scales:
		      new_focal_length = camera.focal_length * jnp.exp(camera_params['log_focal_scale'])  # 更新focal

        SE3CameraDelta.transform_camera(camera_params, camera.replace(focal_length=new_focal_length)

			def transform_camera(
			     self, camera_params: chex.ArrayTree, camera: jaxcam.Camera
			 ) -> jaxcam.Camera:
			   # Convert camera to screw axis representation.
			   translation = spin_math.matmul(-camera.orientation, camera.position)   # (3,3)*(3) =(3)
			   transform = rigid_body.rp_to_se3(camera.orientation, translation)      # (4, 4) homogeneous transformation matrix (rotating by R and translating by p)
			   screw_axis = rigid_body.log_se3(transform)                             # A 6-vector encoding a screw axis of motion
			   new_screw_axis = screw_axis + camera_params['screw_axis_bias']         # (6)
			
			   new_transform = rigid_body.exp_se3(new_screw_axis)
			   new_orientation, new_translation = rigid_body.se3_to_rp(new_transform)
			   new_position = spin_math.matmul(-new_orientation.T, new_translation)
			   return camera.replace( orientation=new_orientation,  position=new_position)

        points = jax.vmap(self.create_points)(transformed_cameras, keys)   # (88,5000,3) camera.position + rays * t_dist

  
  

$\sqrt{d}$ $\frac{1}{0.24}$ $\bar{x}$ $\hat{x}$ $\tilde{x}$ $ϵ$
$\varphi$

  
  

---

总结

提示：这里对文章进行总结：