简介
三维cryoET重建的保真度进一步受到采集过程中物理扰动的影响。这些扰动以各种形式表现出来，例如连续采集之间的样本漂移，导致连续投影未对准，或者由于未散射的电子而导致二维投影中的局部变形。
传统的冷冻电子断层扫描工作流程需要对齐倾斜序列图像，然后应用标准重建算法，例如迭代最小二乘法（ILS）、代数重建技术（ART）、同步迭代重建技术（SIRT）或更常见的滤波反投影（FBP）。
冷冻电子断层扫描流程的第一步是收集三维样本在不同倾斜角度的投影图像。 玻璃化样品只能承受有限剂量的电子而不会受到严重损坏。 这导致信噪比 (SNR) 非常低。 通过直接电子探测器以高帧速率记录每个视角的多个短曝光快照[7]。 尽管每部电影的帧具有极低的信噪比，但可以对齐和平均，补偿给定样本方向上的样本漂移[18]。 然后可以使用曝光加权对对齐的帧进行平均以创建投影图像[19]。 来自所有不同视角的投影图像的集合构成了所谓的倾斜系列。
此外，对比度传递函数（CTF）定义了影响每个观察图像的频率相关衰减和相位反转。CTF是一种与散焦相关的震荡函数，由透射电子显微镜中的光学元件以及样品在焦平面上的不完美放置引起。尽管目标散焦是由显微镜操作员定义的，但是需要通过将CTF函数拟合到图像的功率谱来测量倾斜序列中每个图像的有效散焦。本文中，用CTFFIND4或Warp等现有软件实现此目的。
cryoET重建流程中遇到的主要挑战是电子束引起的连续倾斜之间的大量样品漂移和变形。将图像配准为倾斜序列并校正失真的过程成为倾斜序列对齐。流行方法基本依赖于复杂的局部跟踪和连续图像之间的特征关联。本文重点关注的是没有基准标记的更一般的情况，即原位冷冻电子断层扫描的情况。
AreTomo 是最近流行的解决cryoET对齐问题的软件，通过估计位移、面内旋转和局部变形来对齐，然后再应用FBP算法获得tomogram。
与本文方法类似的[25]联合估计断层图像和倾斜系列位移变形，但是依赖于基于网格的体积表示，并使用原始算法仅估计位移。本文感兴趣的是估计更复杂的局部变形。
另外由于cryoET信噪比极低，去噪是处理重建体积的另一个关键步骤。
隐式表示：神经辐射场（NeRF）模型可以用于表示大型3D场景，而且还用于解决涉及偏微分方程的高维逆问题。它们使用隐式神经网络定义3D对象，然后一致地合成与观察结果相对应的2D场景。和基于网格的体积表示相比，NeRF模型通常表现出较少的参数数量，这使得它们适合表示高分辨率细节，同时保持相对容易训练。本文采用了最近的拓展，引入了多分辨率空间编码，显著提高了训练效率。这在数据量可能非常大的cryoET中很有用。正如STEM所证明的，NeRF模型在表示生物体机方面表现出了良好的前景，并可以和变形估计相结合，如Nerfies。
方法
Image formation model
设
ρ
:
R
3
→
R
+
\rho : \mathbb{R}^3 \rightarrow \mathbb {R}^+
ρ:R3→R+ 是有界、非负的体积密度函数。冷冻电子断层正向模型可以描述为对由大小为
N
×
N
N\times N
N×N 的
M
M
M个离散图像组成的带噪声的倾斜系列的观察，即：
y
m
=
H
m
D
(
γ
m
∗
)
Π
R
t
i
l
t
(
θ
m
)
ρ
+
η
m
\mathbf{y}_m=H_mD(\gamma^*_m)\Pi R_{tilt}(\theta_m)\rho+\eta_m
ym=HmD(γm∗)ΠRtilt(θm)ρ+ηm，其中
m
=
1
,
.
.
.
,
M
m=1,...,M
m=1,...,M。
式中：
- R t i l t ( θ m ) R_{tilt}(\theta_m) Rtilt(θm)表示绕第二轴旋转（倾斜） θ m \theta_m θm角度，使得对于所有 x ∈ R 3 x\in\mathbb{R}^3 x∈R3， R t i l t ( θ m ) ρ ( x ) = ρ ( R t i l t ( α ) x ) R_{tilt}(\theta_m)\rho(x)=\rho(R_{tilt}(\alpha)x) Rtilt(θm)ρ(x)=ρ(Rtilt(α)x)。
这句话描述了旋转操作 R tilt ( θ m ) R_{\text{tilt}}(\theta_m) Rtilt(θm) 对体积密度函数 ρ \rho ρ 的作用。具体来说， R tilt ( θ m ) R_{\text{tilt}}(\theta_m) Rtilt(θm) 是绕第二轴旋转（倾斜） θ m \theta_m θm 角度的旋转操作。对于所有的 x ∈ R 3 x \in \mathbb{R}^3 x∈R3， R tilt ( θ m ) R_{\text{tilt}}(\theta_m) Rtilt(θm) 将 x x x 旋转到一个新的位置，并且在新位置上的密度值等于在原位置上旋转后的密度值。
-
沿第三轴的积分投影 Π \Pi Π定义为：
( Π u ) ( x 1 , x 2 ) = ∫ − ∞ ∞ u ( x 1 , x 2 , x 3 ) d x 3 (\Pi u)(x_1,x_2) = \int_{-\infty}^{\infty} u(x_1,x_2,x_3)dx_3 (Πu)(x1,x2)=∫−∞∞u(x1,x2,x3)dx3 -
D ( γ m ∗ ) D(\gamma^*_m) D(γm∗) 表示影响倾斜序列的第m个变形， H m H_m Hm将离散化建模为 N × N N\times N N×N像素的有限集合，并且 η m \eta_m ηm是独立的加性随机噪声。这种cryoET图像形成过程如图2所示。在cryoET中可以合理地假设CTF和观察方向已知。为了简化，本文忽略CTF的影响。因此需要估计的是体积密度 ρ \rho ρ和基于M个观察 y m \mathbf{y}_m ym的变形 D ( γ m ∗ ) D(\gamma^*_m) D(γm∗)。
总之，密度的准确重建受到三个主要挑战：
- 由于角度 θ m \theta_m θm未覆盖整个范围导致的楔形缺失；
- 未知的变形 D ( γ m ) D(\gamma^m) D(γm)；
- 极低的SNR。
Deformation model
机械平台漂移对应的是全局变形，即倾斜序列的二维图像的移位和面内旋转。电子束造成的样品局部变形，不能通过刚性变换来建模。
本文将变形算子
D
D
D分解成三个基本算子：位移S，面内旋转R和连续局部变形
ϕ
\phi
ϕ。
D ( γ ) = Φ ( C ) R ( α ) S ( τ ) D(\gamma)=\Phi(C)R(\alpha)S(\tau) D(γ)=Φ(C)R(α)S(τ)
需要估计参数 γ = { C , α , τ } \gamma = \{C, \alpha, \tau\} γ={C,α,τ} 。 令 u : R 2 → R u : \mathbb{R}^2 \rightarrow \mathbb{R} u:R2→R 表示一个二维连续图像，变形模型的三个分量的作用如下：
a) The shift operator
S
(
τ
)
S(\tau)
S(τ): 由平移向量
τ
∈
R
2
\tau \in \mathbb{R}^2
τ∈R2 参数化：
S
(
τ
)
u
:
x
→
u
(
x
−
τ
)
S(\tau)u : x \rightarrow u(\mathbf{x} - \tau)
S(τ)u:x→u(x−τ).
b) The in-plane rotation operator R ( α ) R(\alpha) R(α): 由旋转角度参数化 α ∈ [ 0 , 2 π ] \alpha \in [0, 2\pi] α∈[0,2π]:
R ( α ) u : x → u ( R α x ) R(\alpha)u : \mathbf{x} \rightarrow u(R_{\alpha}\mathbf{x} ) R(α)u:x→u(Rαx)
其中 R α = ( cos ( α ) − sin ( α ) sin ( α ) cos ( α ) ) R_{\alpha} = \begin{pmatrix} \cos(\alpha) & -\sin(\alpha) \\ \sin(\alpha) & \cos(\alpha) \end{pmatrix} Rα=(cos(α)sin(α)−sin(α)cos(α)) 是 2D rotation matrix.
c) Local continuous deformation: 定义了一个连续变形操作 Φ ( C ) \Phi(C) Φ(C) 作用于投影图像：
Φ ( C ) y : x → y ( x + ϕ ( C ) ( x ) ) \Phi(C)y : \mathbf{x} \rightarrow y(\mathbf{x} + \phi(C)(\mathbf{x} )) Φ(C)y:x→y(x+ϕ(C)(x))
其中 Φ ( C ) : R 2 → R 2 \Phi(C) : R^2\rightarrow R^2 Φ(C):R2→R2 是一个定义视场任何给定位置处局部位移的连续函数。
本文假设 Φ ( C ) \Phi(C) Φ(C) 是通过离散位移场的插值来定义的。我们在固定控制点定义 N 1 × N 2 N_1\times N_2 N1×N2位移向量。这些位移向量的幅度和方向是定义函数 ϕ \phi ϕ的参数C。一共有 2 × N 1 × N 2 2\times N_1 \times N_2 2×N1×N2 个参数。使用双线性或者双三次插值来获得任意连续坐标处的位移（图6）。一种可能的替代方案是使用隐式神经网络来模拟变形场。
注意：
- 本文将变形建模为作用在倾斜系列的二维投影上的运算符。和在体积空间定义变形模型相比，简化了估计过程，但仍具有足够的表现力，可以在实际数据重建中近似体积变形。
- 为了尽可能通用，我们不假设参数系数{ γ m {\gamma_m} γm}之间存在任何时间或者顺序关系。强制加入这样的关系可以在复杂情境（例如高噪声情况）中对参数的估计进行正则化。可以使用特定设计的正则化项来纳入这样的先验假设。
- 根据经验性地观察，将三类变形分别建模时性能要好得多，类似的观察已经在Nerfies中提出过，但是和其不同的是本文使用单独的参数化来捕获从一个图像到下一个图像的局部变形，因为不假设它们之间存在连续的变形。
Volume density
所有变形算子自然地被定义为作用在体积密度的坐标空间上。因此拥有一种允许快速实现这些变形的表示形式至关重要。
因此，我们通过一个由参数 ψ \psi ψ参数化的隐式神经表示 V ( ψ ) : R 3 → R V(\psi) : \mathbb{R}^3 \rightarrow \mathbb{R} V(ψ):R3→R 来表示体积 ρ \rho ρ。
使用 instant-NGP中描述的多分辨率哈希编码作为位置编码。它定义了不同分辨率的多个特征网格，并使用这些特征的连续插值作为多层感知器(MLP)的输入。插值的特征图包含了查询位置的丰富表示。
这种表示进一步促进了在计算前向成像模型中沿射线进行密度积分的数值近似；对于固定的基于网格的表示，通常必须借助插值。
参数估计
联合估计参数
γ
m
=
(
C
m
,
α
m
,
τ
m
)
m
=
1
M
{γ_m =(C_m,α_m,τ_m)}^M_{m=1}
γm=(Cm,αm,τm)m=1M 和
ψ
ψ
ψ，通过最小化观测数据和根据我们的重建计算出的倾斜序列之间的平均绝对误差，并用估计变形进行校正：
在成像深度学习中，使用 ℓ1-范数代替统计最优（对于独立同分布高斯噪声）ℓ2-范数，从经验上看，它会有更快的收敛和更好的重建质量。
注意：
4. 公式（6）允许在像素空间中进行小批量优化，这与[16]中提出的公式相反。此功能对于处理冷冻电子断层扫描实验中收集的大尺寸倾斜系列至关重要。
本文中假设变形算子是“小”的（接近恒等式），这在冷冻电子断层扫描中确实是这样。该先验是通过最小化变形算子的参数大小来强制执行的：
使用mini-batch梯度下降来最小化(8)。
本文没有对体积密度施加任何显式的正则化，如总变差。作者观察到这是不必要的，因为隐式神经表示的有利归纳偏差会产生抑制噪声和增强细节的重建结果。类似的观察也在其他地方得到过，并从理论上和神经网络架构的频谱特性联系起来。事实上，本文发现隐式正则化甚至部分恢复了缺失楔形信息（图7、图12）。
实验
仿真数据集
我们使用 SHREC 2021 数据集，该数据集由以 1 nm/体素采样的断层图像组成，总大小为 512×512×512 体素。 样本板仅占据沿 z 方向的中央 180 个体素。 因此，我们将体积裁剪为 512 × 512 × 180 的大小，仅保留 z 指数在 256 − 90 和 256 + 90 之间的中心切片。每个断层图像由随机方向的各种大小和形状的不同模拟蛋白质组成。把这个当成实验中的初始层析图。
随后，应用由方程(1)描述的正向模型生成一个倾斜系列，其中包含61个投影图像，观察角度在-60到+60度之间，步长为2度。变形参数独立随机采样，振幅较小，最大移位为12像素，最大平面内旋转为0.01度，局部变形最多引起2像素的位移。真实的局部变形使用模型(5)获得；为了避免逆问题，变形定义在5×5的控制点网格上，而我们的算法使用更细的10×10控制点网格，这些控制点网格不包括真实控制点作为子集。我们扰动倾斜系列以添加全局和局部变形，并添加高斯噪声以获得最终信噪比为10 dB（干净信号y和扰动信号y’之间的信噪比）。
图3展示了倾斜序列中的一些图像。
用FSC评估重建质量。壳半径 r 处两个离散体积 ρ,ρ′ 之间的 FSC 定义为
FSC 以频率相关的方式测量估计体积和原始体积之间的相关性，提供有关可以恢复的细节数量的信息。
此外，我们还报告了估计值与真实体积之间的归一化相关系数 (CC)。 相关系数的优点是用单个数字概括重建质量，从而更容易比较不同的实验设置。 两个离散体积 ρ,ρ′ 之间的相关系数定义为：
请注意，FSC 和 CC 是两个不同的指标，应理解为互补的。 估计的体积有时可能会遭受全局偏移，这对重建的质量没有影响，但如果不考虑的话，可能会恶化 FSC 或 CC 度量的结果。 因此，在报告 FSC 或 CC 之前，我们会系统地将估计量与实际情况进行调整。
和其他常见方法对比
将ICE-TIDE和AreTomo、Etomo这两个流行的对齐算法进行比较。
我们还报告了使用 FBP 直接应用于原始（变形和噪声）投影（FBP w/ def）的重建体积。 这对应于忽略变形，因此应该导致最差的重建。 最后，我们在没有变形但存在噪声的情况下使用 FBP 报告重建体积（FBP w/o def）。 这代表了我们希望在不补偿缺失楔块的情况下实现的最佳估计。
- 仅依靠获得的预测：
通过解决问题(8)得到的体积是通过对隐式神经网络V(ψ)进行采样得到的。实践中有时希望仅使用倾斜系列观测到的重建，而不使用深度神经网络的输出。幸运的是，我们的方法可以用于保留仅估计的变形，这些变形已被证明是非常准确的，以创建一个对齐的倾斜系列版本。因此，可以使用标准的重建算法如FBP或SART来重建体积。我们将这个结果报告为FBP + ICETIDE，并展示它的质量几乎与使用FBP在真实未变形的噪声倾斜系列上获得的重建相同（这在现实中是不可用的）。当应用倾斜系列对一些去噪方法需要的“奇/偶”版本时，这是特别有用的。
重建一个单独体积
我们研究了每种方法在 SHREC 数据集模型 0上的性能。我们对倾斜序列进行降级，使得无噪声和变形投影的 SNR 达到 10 dB。 将估计体积的质量与图 4 中的 FSC 曲线进行比较。我们观察到 AreTomo 和 Etomo 表现不佳。 相比之下，ICE-TIDE 表现出很强的性能，优于最佳 FBP 重建（FBP w/o def）。 证实了隐式神经网络的归纳偏差充当了有效的正则化器，与冷冻电子发射体积的空间统计完美匹配：ICE-TIDE 产生的估计体积恢复了一些缺失的楔形，尽管它没有接受过这方面的明确训练。
图 7 中显示了估计体积的中心切片的傅里叶变换。观察到对于所有 FBP 重建，仅由于空间离散误差和噪声而填充了缺失的楔形。 相反，由于 ICE-TIDE 估计的体积不受缺失楔形的限制，因此它可以在那里假设任意值； 同样，推断值结构良好，代表自然体积密度。
在存在噪声的情况下，重要的是尽早停止（8）的迭代解析，以免过度拟合噪声。 在实验中我们观察到过拟合状态在训练损失稳定后很长时间才出现，这提供了方便的早停信号。 从图 5 中重建体积的正交 XYZ 切片观察到，与 FBP 重建相比ICE-TIDE 中的条纹伪影明显减少 。
我们注意到，使用 FBP 重建的体积（仅由 ICE-TIDE 估计变形参数）紧密遵循 FSC 曲线 FBP w/o def。 在图 5 中观察到，其中非常小的结构的密度与从没有变形的投影获得的密度相似。
在图 6 中直观地比较了估计局部变形和真实局部变形的位移矢量场。我们观察到 AreTomo 错误地估计了真实变形，而 ICE-TIDE 能够准确捕获大多数局部变形。
Time & Memory
ICE-TIDE 的一项关键优势是它可以执行批量采样，然后以内存换取时间，而无需对倾斜系列数据进行下采样。 因此，它可以在低资源硬件上运行，对于固定的重建精度，计算时间相应更长。 在表 I 中，我们比较了 AreTomo 和 Etomo 的时间和内存使用情况。
ICE-TIDE 在 CPU Intel i9 和 GPU NVIDIA GeForce RTX 3090 的计算机上，不到 20 分钟即可完成重建。
抗噪声鲁棒性
评估倾斜序列上不同噪声水平的重建质量。 图 8 中报告了 SNR 范围为 -20dB 至 30dB 的相关系数。 虽然高 SNR 数据实际上并不可行，但它可以让我们了解 ICE-TIDE 理论上可以捕获的细节水平（如果可以通过任何方式衰减噪声）。 在中等 SNR (> 10dB) 下，使用 ICE-TIDE 获得的相关性稳定在 0.9 左右，而在纯 FBP 方法或 AreTomo 中它仍然显着增加。
训练轮次
多个体积的比较
在SHREC2021的10个体积上进行了实验，本文方法表现都不错。
真实数据集实验
使用了T. kivui的断层图，原始数据集可在EMPIAR-11058上获得。由于实验获取的倾斜系列没有GT重建，我们将我们模型的重建结果与使用多个软件包组合得到的结果进行比较。
原始倾斜系列已通过eTomo进行了补丁跟踪对齐，然后在IMOD中使用了加权反投影获得了第一个断层图，然后使用cryo-CARE对从奇偶原始帧重建的倾斜系列对进行去噪。最终的重建结果可以在EMD-15056中找到。图11显示了重建体积的正交切片。我们观察到，ICE-TIDE不仅在重建质量上达到了与人工组合多个软件包并进行人工监督相媲美的水平，而且在分辨率方面要好得多。
这表明ICE-TIDE确实能够同时进行变形估计、去噪和重建步骤，甚至在某种程度上恢复了缺失的楔形和稀疏的角度采样。正如在对合成数据进行的实验中观察到的那样，我们注意到ICE-TIDE以一致的方式填补了缺失楔形中的傅立叶系数，参见图12。图13中显示了ICE-TIDE估计的局部变形以及倾斜系列中3个连续视图的5×5控制点网格。
讨论和限制性
ICE-TIDE的成功基于对冷冻电子断层扫描的采集流程进行了仔细的建模，包括各种样本变形。值得注意的是，ICE-TIDE在合成测量上产生了高质量的重建结果，并且在实际测量中也产生了高质量的重建结果，尽管为了简单起见并没有对CTF的影响进行建模，并且考虑了一个简单的噪声模型。
调整ICE-TIDE以补偿CTF的影响将需要一些原创性思维。一个直接的扩展(1)，假设CTF在实空间中具有c×c像素的支持，将使计算成本增加c2倍——这是不可接受的减速，因为应该取c相对较大。对于低电子区域的采集，噪声最好由依赖于信号强度的随机变量进行建模，例如泊松分布。然而，这种建模不匹配似乎并没有显著影响重建的质量，正如我们的实验所观察到的那样。然而，通过简单地改变数据保真度项(6)，可以实现更准确的噪声模型，参见[46]，[47]。
考虑低电子制度开启了将倾斜系列对齐和重建过程直接扩展到直接读取剂量分数数据的可能性，即在短暴露时间内获得的每个视角的原始电影帧，目前在单独的阶段事先对齐和平均。请注意，这也将按比例增加所需的计算资源，与分数的数量成正比。
虽然我们没有明确地强加先验知识约束以估计缺失的楔形信息，但我们观察到ICE-TIDE的估计可以以相当准确的方式检索这些信息。这是由于通过隐式神经网络对体积密度进行参数化，其归纳偏差使其成为逼真体积的有效先验。换句话说，我们选择的表示方式允许网络以合理的方式在观察到的傅立叶系数之间进行插值。还可以推断出有关缺失楔形的信息，例如，假设同等可能地观察到感兴趣体积的每个补丁或该补丁的任何旋转，从物理观点来看，这种假设是有意义的，因为体积中的许多特征，如蛋白质，实际上几乎同等可能地出现在任何方向。这样的先验可以作为附加的正则化项（8）中加以考虑。这与IsoNet和DeepDeWedge中提出的想法有关。
结论
提出了ICE-TIDE，用于解决在存在影响噪声投影的物理变形的情况下进行冷冻电子断层重建的问题。ICE-TIDE利用了一种隐式神经网络架构，能够以非凡的保真度学习产生观察到的2D投影的底层3D体积。有趣的是，我们的方法不仅可以提供准确的倾斜系列对准，还可以直接产生断层图。此外，由于冷冻电子断层成像中采用的离散且范围受限的倾斜方案，网络学习恢复了未直接记录在倾斜系列中的信息，因此在相当大程度上补偿了缺失的楔形区域。
局部和全局变形可能具有不同的时间分辨率。我们忽略了这一点，并假设所有变形具有相同的时间分辨率。我们的模型更具表现力，可以以高度保真度恢复被测样本的图像。对变形的时间性质的更精确了解可以进一步作为先验知识纳入到我们的方法中。
展示了模型的未知参数可以通过对数据保真度损失进行梯度下降来简单地估计。虽然这里专注于冷冻电子断层成像的挑战性问题，但唯一的要求是对正演算子有精确的了解。因此，更一般的反问题可以用这种形式来解决。
