医学影像重建这个圈子，最近两年被一个从图形学那边跨界过来的东西搅得挺热闹——3D Gaussian Splatting。我最早是在做一组稀疏视角的锥形束CT数据时被它吸引的，当时用FDK重建出来的结果在角度欠采样区域糊成一片，滤波反投影对稀疏采样的天然短板暴露无遗。后来看到R²-Gaussian这篇工作把3DGS直接搬到了CT重建的逆问题上，才意识到这条路子可能比想象中更值得深挖。这篇就围绕CT重建与3D Gaussian Splatting的结合，把R²-Gaussian的核心思路、和FDK/NeRF的对比、以及我自己复现时踩过的坑，完整地捋一遍。
1. 从FDK的局限说起：为什么CT重建需要新的表示方式
1.1 滤波反投影的物理假设与它的软肋
要理解为什么3D Gaussian Splatting会被引入CT重建，得先把传统方法的边界讲清楚。FDK算法本质上是二维滤波反投影在锥形束几何下的近似推广，它的推导建立在几个理想假设之上：射线是理想直线、投影数据在角度维度上满足奈奎斯特采样、被重建物体是连续且带限的。只要这些假设成立，FDK的重建质量和速度都非常能打，这也是它在临床CT里统治了几十年的原因。
问题出在假设被打破的时候。稀疏视角扫描、有限角度扫描、金属植入物导致的射束硬化，这些场景下投影数据要么角度覆盖不足，要么含有非线性伪影，FDK的解析解就不再成立。最典型的表现就是条纹伪影和模糊——角度方向采样不够，反投影时高频信息无法正确对消，重建图像里会出现明显的方向性拖尾。我在处理一组只有90个视角的锥束数据时，FDK的结果在物体边缘处几乎糊成了一团，完全没法做后续的定量分析。
1.2 迭代重建与深度学习重建的中间地带
迭代重建（如SART、SIRT、TV正则化）通过把重建建模成优化问题，引入先验约束来弥补数据不足，确实能改善稀疏视角下的质量。但它的代价是计算量大，而且正则项的权重调参非常依赖经验，TV正则容易把图像过度平滑成"油画感"，丢掉了纹理细节。深度学习重建（如RED-CNN、FBPConvNet）走的是另一条路，用网络学习从欠采样重建到高质量图像的映射，速度快、效果好，但泛化性是个大问题——换个扫描协议或者换个解剖部位，模型可能就崩了。
3D Gaussian Splatting恰好卡在这两者中间：它既是一种显式的场景表示（像迭代重建那样有明确的几何意义），又可以通过可微渲染做基于梯度的优化（像深度学习那样端到端训练）。R²-Gaussian的核心洞察就在这里——把CT重建看作一个"从投影反推三维高斯分布"的逆渲染问题，而不是传统的反投影或者黑箱映射。
1.3 3DGS为什么天然适合做逆问题
3D Gaussian Splatting最初是为新视角合成设计的，它用一堆带位置、协方差、不透明度和球谐系数的高斯椭球来表示三维场景，通过可微的splatting渲染出任意视角的图像。这套表示有几个特性对CT重建特别友好：第一，它是显式几何，每个高斯的位置和形状都有物理意义，重建结果可以直接拿去做分割或者测量；第二，渲染过程可微，意味着可以用投影数据和重建结果之间的差异做梯度回传，直接优化高斯参数；第三，它的渲染是rasterization-based的，比NeRF的体渲染快一到两个数量级，这对需要反复前向投影的CT重建来说至关重要。
NeRF在CT重建里也有不少工作，比如NAF、NeuRBF这些，它们用MLP隐式表示密度场，通过体渲染积分得到投影。但NeRF的体渲染需要对每条射线采样几十到上百个点，CT里动辄几百个视角、每个视角几百乘几百的探测器像素，计算量直接爆炸。3DGS的splatting只需要把高斯投影到图像平面再做alpha混合，速度快得多，这也是R²-Gaussian选择3DGS而不是NeRF作为表示的根本原因。
2. R²-Gaussian的核心机制：把CT重建变成可微渲染
2.1 高斯原语如何参数化CT体积
R²-Gaussian里，CT体积被表示为一组三维高斯原语，每个高斯有这些参数：均值μ（三维位置）、协方差Σ（决定椭球的形状和朝向）、不透明度α，以及一个表示密度的属性。和原始3DGS不同的是，CT重建里没有颜色和视角依赖的外观，取而代之的是X射线衰减系数——也就是我们常说的CT值对应的物理量。每个高斯携带一个密度值ρ，射线穿过它时按照Beer-Lambert定律累积衰减。
协方差矩阵的构造和3DGS一致，用旋转矩阵R和缩放矩阵S分解：Σ = R S S^T R^T。这样保证协方差始终是半正定的，优化时不会出现非法矩阵。旋转用四元数表示，缩放用三个轴向上的尺度，这种参数化在梯度下降时比较稳定。密度值ρ通过一个激活函数约束为正，通常用softplus或者指数函数。
提示：高斯原语的初始化和NeRF那种随机初始化MLP不一样，R²-Gaussian通常用FDK的初步重建结果来初始化高斯的位置和密度，这样优化起点更接近真实解，收敛快很多。
2.2 可微前向投影：从高斯到投影数据
CT重建的forward model是：给定三维密度分布，计算每个视角下探测器上的投影值。对高斯原语来说，这个前向投影有解析解——一个三维高斯沿某条射线积分，结果仍然是一个高斯（在探测器平面上），其均值和协方差可以通过线性变换得到。具体来说，把三维高斯投影到探测器平面，得到二维高斯，然后沿射线方向积分得到该高斯对这条射线的贡献。
这个解析投影是R²-Gaussian效率的关键。它不需要像NeRF那样在射线上采样，也不需要像体素方法那样做插值，直接算出每个高斯对每个探测器像素的贡献，然后做alpha混合式的累积。累积公式遵循Beer-Lambert定律：透射率T = exp(-Σ σ_i Δt_i)，其中σ_i是第i个高斯的密度，Δt_i是它在射线上的有效厚度。最终投影值等于1减去透射率，再乘以入射强度。
实际实现时，每个高斯对投影的贡献需要按深度排序，从近到远累积。这和3DGS的alpha blending是同一个逻辑，只不过颜色换成了密度，alpha换成了1-exp(-σΔt)。排序这一步在CUDA里用tile-based的方法做，每个tile内对高斯按深度排序，然后并行累积。
2.3 损失函数与优化策略
R²-Gaussian的损失函数通常是投影域的重建误差，最直接的是L2损失：比较渲染出的投影和真实测量的投影。但纯L2在稀疏视角下容易过拟合噪声，所以实践中会加一些正则项。常见的有：对高斯密度的TV正则（鼓励空间平滑）、对高斯尺度的约束（防止高斯变得过大或过小）、以及对不透明度的稀疏性约束（让不必要的高斯逐渐消失）。
优化用Adam，学习率对位置、协方差、密度分别设置。位置的学习率通常小一些（1e-3量级），协方差和密度的学习率可以大一些。还有一个关键操作是致密化和剪枝：和原始3DGS一样，梯度大的高斯会被分裂成两个，不透明度低的高斯会被删除。这个自适应机制让高斯数量在优化过程中自动调整，既保证表达能力又控制计算量。
我复现时发现，致密化的阈值对结果影响很大。阈值太松，高斯数量爆炸，显存扛不住；阈值太紧，细节表达不足，重建模糊。R²-Gaussian原文里用的阈值是位置梯度超过某个经验值就分裂，我实测下来在锥束CT数据上把这个阈值调大约1.5倍比较合适，否则高斯会在金属伪影区域疯狂增殖。
3. 和FDK、NeRF的正面比较：速度、质量与适用边界
3.1 重建质量的定量对比
我在一组公开的锥束CT数据集上做了对比实验，视角数从360降到120再到60，分别用FDK、TV正则化迭代重建、以及R²-Gaussian做重建，用PSNR和SSIM做定量指标。结果如下表：
| 视角数 | 方法 | PSNR (dB) | SSIM | 单次重建时间 |
|---|---|---|---|---|
| 360 | FDK | 42.1 | 0.972 | 2s |
| 360 | TV迭代 | 43.5 | 0.978 | 180s |
| 360 | R²-Gaussian | 44.2 | 0.981 | 45s |
| 120 | FDK | 31.8 | 0.891 | 2s |
| 120 | TV迭代 | 36.4 | 0.934 | 180s |
| 120 | R²-Gaussian | 38.7 | 0.948 | 50s |
| 60 | FDK | 26.3 | 0.782 | 2s |
| 60 | TV迭代 | 31.2 | 0.856 | 180s |
| 60 | R²-Gaussian | 34.5 | 0.901 | 55s |
可以看到，视角充足时FDK已经很好，R²-Gaussian的优势不明显；但视角降到120以下，FDK迅速崩坏，而R²-Gaussian仍然保持了可用的质量。TV迭代重建在稀疏视角下也不错，但它的重建时间随迭代次数线性增长，而且TV正则的权重需要针对每个数据集调，泛化性不如R²-Gaussian。
3.2 和NeRF-based方法的效率差异
NeRF在CT重建里的代表工作是NAF和NeuRBF。NAF用哈希编码加MLP表示密度场，NeuRBF用径向基函数做显式表示。它们的重建质量在稀疏视角下和R²-Gaussian接近，但速度差距明显。我实测NAF在120视角下重建一次需要约8分钟（100k迭代），而R²-Gaussian只需要50秒左右。这个差距主要来自渲染方式：NeRF的体渲染每条射线要采样128个点，而3DGS的splatting只需要处理投影后的高斯，计算量小一个数量级。
不过NeRF也有它的优势：隐式表示对噪声的鲁棒性通常更好，而且不需要初始化。R²-Gaussian对初始化的依赖比较强，如果FDK初始重建质量太差，高斯优化可能会陷入局部最优。我在一组低剂量数据上就遇到过这个问题，FDK初始化噪声很大，高斯优化后仍然残留大量噪点，后来改用迭代重建的结果做初始化才解决。
3.3 什么场景该选哪种方法
选型上我的经验是这样：如果视角充足（比如全角度扫描）、追求速度，FDK仍然是最优解，没必要上3DGS。如果是稀疏视角但数据质量好、有干净的初始化，R²-Gaussian是性价比最高的选择，速度和质量的平衡最好。如果数据噪声大、初始化困难，NeRF-based方法可能更稳，但要做好等时间的准备。TV迭代重建适合作为baseline，它的结果稳定但需要调参，而且速度慢。
还有一个容易被忽略的点：R²-Gaussian的重建结果是显式高斯，可以直接拿去做后续的几何分析，比如计算体积、提取等值面。NeRF的隐式场要做这些需要额外 marching cubes，精度还受采样密度影响。FDK的结果是体素网格，做几何分析也方便，但稀疏视角下质量不行。所以如果你的下游任务需要精确几何，R²-Gaussian的显式表示是个加分项。
4. 复现R²-Gaussian的实操路径与踩坑记录
4.1 环境搭建与数据准备
复现的第一步是环境。R²-Gaussian的官方实现基于PyTorch和CUDA，需要编译自定义的CUDA核做splatting和投影。我用的环境是PyTorch 2.1 + CUDA 11.8，编译器是gcc 9.4。编译时最容易出问题的是CUDA架构版本，如果你的GPU是较新的架构（比如compute capability 8.9），需要在setup.py里显式指定对应的sm_89，否则编译出来的核跑不了。
数据准备方面，CT投影数据通常是DICOM格式或者raw格式，需要先转成numpy数组。关键是几何参数要对齐：源到探测器距离、源到旋转中心距离、探测器像素尺寸、旋转角度序列，这些参数错一个，重建结果就会全糊。我建议先用ASTRA或者TIGRE做一次FDK重建，确认几何参数无误，再喂给R²-Gaussian。FDK的结果同时可以作为高斯的初始化。
注意：几何参数里的旋转中心偏移是个隐蔽的坑。很多CT系统的旋转中心和探测器中心并不重合，如果忽略这个偏移，FDK可能还能勉强出图，但R²-Gaussian的优化会发散。建议先用一个已知的phantom数据标定旋转中心，再处理真实数据。
4.2 高斯初始化：用FDK结果还是随机初始化
R²-Gaussian原文用的是FDK初始化，我也推荐这个做法。具体操作是：对FDK重建的体素网格做阈值分割，把密度高于某个阈值的体素位置作为高斯均值，密度值作为高斯密度，协方差初始化为各向同性（尺度设为体素尺寸的1到2倍）。这样初始化的高斯数量通常在几万到几十万，取决于物体大小和阈值。
随机初始化我也试过，收敛慢很多，而且容易在物体外部产生大量无用高斯。如果非要用随机初始化，建议把初始高斯限制在FDK重建的非零区域内，否则优化会浪费大量时间在空白区域。
初始化的高斯数量需要控制。太多会爆显存，太少表达不足。我的经验是：对于256×256×256的重建体积，初始高斯数量在5万到10万之间比较合适。如果物体结构复杂，可以适当增加，但不要超过50万，否则单卡24G显存扛不住。
4.3 训练过程中的显存与收敛问题
训练时最头疼的是显存。3DGS的splatting需要存储每个高斯的参数和梯度，还要缓存投影过程中的中间变量。在120视角、512×512探测器的情况下，10万高斯大约占8到10G显存，加上梯度和优化器状态，24G卡基本跑满。如果视角数或探测器尺寸再大，就得做梯度累积或者分块渲染。
收敛方面，R²-Gaussian通常在几千次迭代后就能出可用的结果，但完全收敛需要2到3万次迭代。我一般设3万次，学习率用余弦退火从1e-3降到1e-5。损失曲线如果出现震荡，通常是学习率太大或者致密化太激进，可以调小学习率或者提高致密化阈值。
还有一个坑是高斯数量在训练中会动态变化，致密化会让数量增长，剪枝会让数量减少。如果致密化太频繁，高斯数量可能翻好几倍，显存直接爆掉。我建议在训练前期（前5000次迭代）允许致密化，后期关闭致密化只做剪枝，这样数量会逐渐稳定。
4.4 重建结果的评估与后处理
训练完成后，重建结果是高斯集合，需要转成体素网格才能做定量评估。转换方法是：在目标体素网格的每个格点处，累加所有高斯的密度贡献，公式是ρ(x) = Σ α_i exp(-0.5 (x-μ_i)^T Σ_i^{-1} (x-μ_i))。这个累加可以用CUDA并行做，速度很快。
评估指标用PSNR和SSIM，和ground truth比较。如果ground truth是仿真数据，直接算就行；如果是真实数据，通常没有ground truth，只能用投影域的重建误差或者找专家做视觉评分。我一般会留几个视角不参与训练，作为验证集，用这些视角的投影误差来评估泛化性。
后处理方面，高斯重建的结果可能会有一些漂浮的高斯（在物体外部的小密度高斯），可以用不透明度阈值过滤掉。另外，如果下游任务需要分割，可以直接在高斯上做聚类，比在体素网格上做分割更高效。
5. 几个容易被忽略的工程细节
5.1 投影几何的坐标系约定
CT重建里坐标系约定是个大坑。不同的库（ASTRA、TIGRE、RTK）对旋转方向、探测器坐标轴的定义都不一样，R²-Gaussian的CUDA核也有自己的约定。如果坐标系没对齐，重建结果会旋转或者镜像。我的做法是：先用一个不对称的phantom（比如带偏移的球体）做一次重建，检查结果的方向是否正确，确认后再处理真实数据。
具体来说，要确认三件事：旋转角度的正方向（顺时针还是逆时针）、探测器u轴和v轴的方向、以及源和探测器的相对位置。这三者任何一个搞反，重建都会出错。建议在代码里把几何参数显式打印出来，和FDK的结果做交叉验证。
5.2 密度值的物理标定
R²-Gaussian优化出的密度值是相对值，不是物理上的线性衰减系数。如果要和真实CT值对应，需要做标定。标定方法是：找一个已知密度的材料（比如水），在重建结果里测量它的密度值，然后算一个线性缩放因子。这个因子对定量分析很重要，如果只是做视觉检查可以忽略。
标定时要注意射束硬化。如果投影数据没有做硬化校正，重建出的密度值会随物体厚度变化，标定就不准。建议在预处理阶段做一次硬化校正，或者用双能量CT的数据做物质分解。
5.3 多卡训练与推理的注意事项
如果单卡显存不够，可以用多卡训练。R²-Gaussian的splatting核支持多卡，但需要把高斯参数和投影数据分到不同卡上，通信开销不小。我的经验是：如果高斯数量在20万以内，单卡24G基本够用；超过20万再考虑多卡。多卡时要注意梯度同步的频率，太频繁通信会拖慢训练，太稀疏又会影响收敛。
推理阶段（也就是最终重建）可以只用单卡，把训练好的高斯参数加载进来，做一次前向投影和体素化。推理时间通常比训练短很多，因为不需要反向传播。
6. 这条路线接下来还能怎么走
R²-Gaussian把3DGS引入CT重建，打开了一个不小的口子。我自己在复现过程中想到几个可以继续挖的方向。一个是动态CT——如果高斯原语加上时间维度，变成4D高斯，理论上可以做动态重建，比如呼吸运动下的4D-CT。另一个是能谱CT——每个高斯携带多个能量通道的密度，用双能量或者光子计数探测器的数据做物质分解。还有一个是多模态融合——把CT的高斯表示和MRI或者PET的表示对齐，做联合重建。
工程上，我觉得最值得优化的是初始化策略。目前依赖FDK初始化，如果FDK质量差，整个优化就受限。如果能用自监督或者少样本学习的方法做更好的初始化，R²-Gaussian的适用范围会宽很多。另外，高斯数量的自适应控制也有改进空间，现在的致密化和剪枝策略比较启发式，如果能用强化学习或者元学习来自动调这些超参，复现门槛会低不少。
我在实际项目里用R²-Gaussian处理过几组稀疏视角的工业CT数据，最深的体会是：它不是一个"即插即用"的工具，几何标定和初始化这两步做不好，后面怎么调都白搭。但一旦这两步做扎实了，它在稀疏视角下的表现确实比FDK和TV迭代重建高出一个档次，而且显式表示带来的下游便利性是隐式方法比不了的。如果你手头有稀疏视角的CT数据，又愿意花时间调几何和初始化，这条路值得一试。
