Fast deep learning based reconstruction for limited angle tomography

基于快速深度学习的有限角度断层成像重建

Abstract

计算机断层扫描的一个主要挑战是从不完整的数据中重建物体。对于这些问题，一个日益流行的解决方案是将深度学习模型整合到重建算法中。本文介绍了一种将傅里叶神经算子(FNO)集成到滤波后反投影(FBP)重建方法中的新方法，得到了FNO反投影(FNO-BP)网络。我们采用矩条件进行正弦图外推，以帮助模型从有限的数据中减轻伪影。值得注意的是，我们的深度学习架构保持了与经典滤波后投影(FBP)重建相当的运行时，确保了推理和训练期间的快速性能。我们在Helsinki Tomography Challenge 2022的背景下评估了我们的重建方法，并将其与常规FBP方法进行了比较。

1 Background, motivation, and specific contributions

1.1 层析成像原理

断层扫描是一种通过从物体周围的方向反复探测物体的粒子/波来对物体内部结构进行非侵入性成像的技术。这些粒子/波在与(内部结构)物体相互作用后被检测到。因此，观测到的数据，从此称为正弦图，代表了人们试图确定的内部结构的嘈杂的间接观测。因此层析成像包括一种重建方法，这是一种从正弦图中恢复内部结构的计算方法。

上面的一个众所周知的例子是X射线计算机断层扫描(CT)，用X射线光子探测物体。对象的性质取决于用例，例如，在医学成像中，它是一个人，而在无损检测中，它通常是一个制造对象。内部结构由(空间变化的)线性衰减系数表示，通常可视化为图像。数据(sinogram)由X射线强度测量组成，沿着穿过物体的一组特定的线。采集几何是指这些线(线的流形)排列的数学规范。

1.2 CT重建方法的演变

从正弦图中恢复图像的重建方法是CT成像的关键部分。许多重建方法是在连续介质设置中导出的，其中试图恢复的图像由二维或三维空间上的实值函数表示，而正弦图表示直线上的实值函数的样本。此外，在CT中，人们通常采用比尔-朗伯定律来模拟X射线光子如何与组织/物质相互作用。将这种简化的物理模型与适当的测量数据预处理相结合，可以将重建任务重新表述为对采集几何结构给出的线流形进行射线变换的反演。

接下来，我们描述了基于上述重构任务重构的CT重构方法的演变。从数学的角度来看，重构的第一步早在1917年就已经开始了，当时Johan Radon发表了Radon变换的解析反演公式[29]，给出了平面环境下射线变换的反演公式。1938年，Fritz John导出了一个类似的一般维射线变换显式反演公式[13]。从那时起，数学界的许多工作都致力于刻画这种积分变换的注入性、范围和稳定性，并推导出它们的显式反演公式，参见[19]对这一数学发展的概述，参见[35]对历史发展的广泛描述。

射线变换的精确反演公式(当它被限制为来自特定采集几何形状的线时)允许人们从表示无噪声连续体数据的正弦图中精确地恢复图像。然而，反演射线变换是一个不适定问题，这意味着精确的反演将不可避免地涉及放大正弦图中小的测量/采样误差的步骤。因此，精确反演公式的实际适用性有限，但它们是解析(重建)方法的基础。这些方法使用傅立叶分析和采样理论的原理来稳定本质不稳定的反演，例如，通过恢复图像的带限制部分，而滤除高频成分，如噪声。这种方法是由Bracewell和Riddle首创的，他们在1967年发表了平面环境下平行线的滤波反投影(滤波后投影)方法[7]。从那时起，FBP的许多变体被开发出来，主要是为了使其适应其他采集几何形状[34,22,9,8,15,25,14]。另一类密切相关的解析方法是傅里叶方法，它代表了投影切片定理的直接实现(它将射线变换与傅里叶变换联系起来)。

除了解析方法之外，还有另外三种用于射线变换反求的重建算法。一种是基于模型的迭代方法，如代数重构技术、同时迭代重构技术、最小二乘共轭梯度法等。这里，射线变换通过迭代方案近似反转，该方案在极限下被设计为收敛到与正弦图最大一致的图像。典型的设置是使用最小化数据保真度项的迭代方案，就像在计算最小二乘解时一样。然而，迭代在收敛之前很久就停止了(提前停止)，以避免由于不适定性而产生的不稳定性，参见[24，第5.3节]了解更多细节。另一类是基于模型的变分方法，如quadratic Tikhonov或全变分正则化。这里，我们最小化一个结合了数据保真度术语和正则化器的目标。数据保真度确保了数据的一致性(如在基于模型的迭代方法中)，而正则化器则会惩罚不稳定性，参见[33,6]对这些方法的广泛调查。最后一类重建方法是基于深度学习的方法。首先将重建方法看作贝叶斯推理中的统计估计量。估计器是一种重构方法(不是重构)，然后从训练数据中学习。可以访问的训练数据的类型以及选择如何设置学习问题决定了要学习哪个估计器。这里的深度神经网络架构是可能估计器的参数化。通常使用领域适应的架构来编码一些底层物理，因为这些架构比一般的深度神经网络架构具有更好的泛化属性，参见[4，第5节]了解更多关于这些问题的信息。

1.3 临床CT重建面临的挑战

CT重建的总体目标是有一种重建方法，可以及时产生“高质量”的3D图像，并且不需要使用超级计算资源的内存占用。

对于全剂量3D临床CT，基于FBP的分析方法在平衡图像质量与重建速度和内存占用方面提供了一种折衷方法，但仍难以超越[27]。然而，分析方法的一个缺点是它们需要专门针对数据采集几何进行定制，例如，如果不诉诸近似，很难利用所有测量数据。当CT数据采样不足(稀疏视图CT)或采样不均匀(如有限角度CT)时，从分析方法获得的图像质量就会成为一个问题。分析方法中使用的稳定策略对于全剂量CT来说是足够的，但在这些情况下，稳定策略的强度不足以提供有用的重建。这在有限角度CT中尤其明显，因为CT扫描仪中的x射线源不能绕物体移动。这就产生了一个正弦图，它的射线变换值只在有限的方向范围内的直线上，这种设置已知会显著增加重建问题的不适定性。分析方法的最后一个缺点是它们本质上是确定性的，因此它们不适合考虑噪声中的统计特性。当数据噪声很大(低剂量CT)时，这是一个问题。

为了解决分析方法的这些缺陷，促进了基于模型的迭代/变分方法的发展，其中一些方法已经进入临床实践[36,30]。与解析方法相比，迭代/变分方法提供了结合更复杂(和更强大)稳定策略的可能性。其次，与分析方法不同，它们也不是针对特定的采集几何形状量身定制的。此外，对于x射线光子与组织/物质的相互作用，他们可以结合一个比分析方法中使用的比尔-朗伯定律提供的更精确的物理模型。最后，迭代/变分方法还可以考虑数据中噪声的统计特性，这在低剂量CT中是必不可少的。另一方面，这些基于模型的方法有很长的运行时间，在时间紧迫的用例中变得令人望而却步，比如创伤中的CT成像。

基于深度学习的方法形成了重建的最后一类。最先进的方法依赖于领域适应物理的深度神经网络架构[3]，在图像质量方面优于基于模型的方法，同时具有更接近分析方法的重建运行时[4,36,30]。通过巧妙的工程设计，人们可以限制它们的内存占用，从而在高端游戏图形处理单元(GPU)硬件上也在3D设置中训练这些适应领域的架构[31]。然而，这并没有解决训练所花费的时间。这里的主要瓶颈是训练涉及计算大量的前投影(对sinogram中所有直线上的射线变换进行评估)和相应的反投影。

这给我们带来了本文所考虑的挑战，即开发一种基于深度神经网络的重建方法，该方法提供比分析方法(在有限的数据设置中)更好的重建质量，同时具有与分析方法相当的运行时和内存占用。这一要求限制了深度神经网络的重建，例如，通过展开一些迭代方案[3,1,4]获得的非常成功的重建架构涉及几个正向和反向投影，因此这些都是不合格的。事实上，拥有与分析方法相当的运行时和内存占用意味着使用一种架构，这种架构避免了前投影，并合并了一个(或几个)反投影。

1.4 具体贡献

其中一个贡献是扇形波束有限角层析成像中的正弦图外推，该外推基于射线变换的距离条件。另一个贡献是引入了新的FNO反投影(FNO-BP)网络架构，该架构专门适用于层析成像重建。它避免了前投影，并结合了一个反投影，因此它在推理期间的运行时间和内存占用与FBP相当。同时，在HTC-2022有限角度层析成像数据上测试，其重建质量优于FBP[20]。

2 Computationally efficient reconstruction methods

我们在这里提供了一个简单的概述计算高效的重建方法。这是指从正弦图中恢复图像的方法，它只涉及计算一个(或几个)反向投影操作。此外，人们可以完全避免计算正向预测。这与基于模型的迭代和变分方法以及大多数最先进的基于深度学习的方法形成对比，后者依赖于包含前向预测的领域适应架构[1,4]。

2.1 比尔-朗伯定律

当用x射线光子照射物体时，光子将与物体的空间变化的线性衰减系数相互作用。这是物体的“内部结构”，人们可以用X射线感知，它在数学上用函数$f:{\mathbb{R}}^d\to \mathbb{R}$表示，如果所有X射线都在一个固定的平面2D横截面上，则d = 2，否则d = 3。

辐射输运方程为物体和X射线光子之间的相互作用提供了一个精确的模型，X射线光子沿着穿过物体的线$\ell$传播[5]。简化这个模型(假设单色X射线并忽略散射事件)得到比尔-朗伯定律：
$$\frac{I_{\mathrm{out}}(\ell )}{I_{\mathrm{in}}(\ell )}=e^{{\int }_{\ell }f(x)d\mu (x)}⟹g(\ell ):=-\log \left(\frac{I_{\mathrm{out}}(\ell )}{I_{\mathrm{in}}(\ell )}\right)={\int }_{\ell }f(x)d\mu (x)\mathrm{for}\ell \in M.\text{\hspace{1em}(1)}$$
其中，$M$是穿过物体的直线的集合，$\mu$是直线上的1D勒贝格量，$\ell \in M$， $I_{in}(\ell )$是x射线光子在与物体相互作用前沿着$\ell$运动的强度，$I_{out}$是与物体相互作用后X射线光子的相应强度。一般来说，在层析数据采集过程中，沿着有限多条线测量$I_{out}(\ell )$(透射率)和相应的入射强度$I_{in}(\ell )$。数据采集几何指定集合$M$以及$M$中对应于测量数据的实际采样线。从数据中计算出的项$g(\ell )$称为沿$\ell$的投影(或吸收)。一个正弦图是函数$g:M\to \mathbb{R}$，即当$\ell \in M$时，所有投影的集合是变化的。

2.2 积分几何的概念

描述本文所考虑的深度神经网络结构需要引入积分几何中的一些概念。

我们首先定义限定在给定的一组直线$M$上的射线变换为：
$$\mathcal{P}(f)(\ell ):={\int }_{\ell }f(x)d\mu (x)\text{for a line}\ell \in M.\text{\hspace{1em}(2)}$$
射线变换在数学上是一个线性算子，它将一个实值函数$f:{\mathbb{R}}^d\to \mathbb{R}$映射到$M$中直线上的一个实值函数。只要f是可测量的并且有紧支持，它就是定义良好的(事实上，它足以要求$f$比任何多项式衰减得更快，参见[24，第2.1节])。接下来，比尔-朗伯定律允许我们将CT中的重建任务看作是对限制在$M$内的射线变换进行反演。要了解这一点，请注意CT重建的目的是从(1)中的正弦图$g:M\to \mathbb{R}$中恢复图像$f:{\mathbb{R}}^d\to \mathbb{R}$。但是$g(\ell )=\mathcal{P}(f)(\ell )$，因此重建在数学上与对限制在$M$内的$\mathcal{P}$进行反演是一样的。

在大多数重建方法中，相应的反投影是一个关键的组成部分。当作用于函数$g:M\to \mathbb{R}$时，定义为：
$${\mathcal{P}}^{\ast }(g)(x):={\int }_{M_x}g(\ell )d\sigma (\ell )\text{for}x\in {\mathbb{R}}^d.\text{\hspace{1em}(3)}$$
这里，$M_x\subset M$表示穿过$x\in {\mathbb{R}}^d$的线 $\ell \in M$的子集，$\sigma$是在$M_x$上对应的(Haar)测度(见[19，第3.3.4节])。我们看到${\mathcal{P}}^{\ast }$是一个线性算子，它将直线上的函数$g:M\to \mathbb{R}$映射到${\mathbb{R}}^d$上的实值函数。现在可以证明${\mathcal{P}}^{\ast }$是射线变换$\mathcal{P}$的形式伴随(见[19，定理3.29])。

我们最后提到Radon变换，它在${\mathbb{R}}^d$中的超平面(而不是直线)上对$f:{\mathbb{R}}^d\to \mathbb{R}$进行积分。由于直线和超平面在${\mathbb{R}}^2$中重合，我们得到射线变换和Radon变换在平面设置中重合(当d = 2时)。

2.3 滤波反投影(FBP)

FBP方法基于以下恒等式[24，定理2.13]，它将函数$f$与其在M上的正弦图$g:=\mathcal{P}(f)$联系起来：
$$f\ast K={\mathcal{P}}^{\ast }(k⊛g)\mathrm{where}K:={\mathcal{P}}^{\ast }(k)\text{\hspace{1em}(4)}$$
这里，左边的$\ast$是${\mathbb{R}}^d$上的函数之间的d维卷积，而右边的$\ast$是直线上$\ell \in M$的函数之间的(d−1)维卷积。对于后者，在直线$\ell$处的卷积由超平面上通过与$\ell$正交的原点的“共”(d−1)维卷积给出[24，定理(2.30)]。

请注意，正弦图$g=\mathcal{P}(f)$是无噪声连续体数据，后者意味着假设射线变换在$M$中的所有直线上都给定，函数$K:M\to \mathbb{R}$称为重构核。它决定了一个寻求恢复的特征。当k可以选择时，FBP方法理论上是精确的，因此$K=\delta$，因为这对应于从连续体无噪声数据中精确地恢复$f$。然而，$\mathcal{P}$的反演是一个不适定问题，因此在这种重构核的选择下计算右侧将涉及不稳定步骤。FBP的思想是通过选择重构核$K$使$K\approx \delta$来稳定它。

明确的重建核存在于平行和扇形梁几何平面设置。这同样适用于3D中的平行光束几何形状。对于三维的锥型几何形状，情况更为复杂。自Katsevich[15]的工作以来，在任何曲线(圆形，螺旋，…)上的源设置中也存在显式重构核。． 然而，这些理论上精确的FBP方法并没有利用所有的数据，这反过来又使它们对噪声过于敏感。由于这个原因，已经开发了近似FBP类型的方法，如Feldkamp-Davis-Kress (FDK)用于圆形源的锥束几何形状，FDK的变体用于其他锥束采集几何形状。

平面设置中的一个特殊挑战是当$M$是不完全的，这意味着${\mathbb{R}}^2$中的直线上有一个开放集不在M中，在这种设置中，没有选择重建核的好方法。在未在开子集上给出的正弦图上应用标准FBP重建意味着隐含地假设正弦图在该缺失区域为零。这几乎总是与测量数据不一致，从而导致非常差的重建。Sinogram extrapolation (Sinogram外推法)指的是在缺失区域中预测数据的方法，以便用非零的数据填充。有一些正弦图外推方法，但是当它们与标准FBP结合时，这些方法通常不能给出足够的精度，因此所得的重建仍然很差。为了改善这种情况，必须将正弦图外推与改进的重构算子结合起来。

2.4 平行和扇形波束坐标

在实际情况中，$M$中无法访问完整的连续体数据，而是根据采集几何形状确定的方案对数据进行采样。因此，我们需要一个适合于给定几何的$M$的参数化。最直接的参数化是平面平行光束几何，其中一条线的参数化是由它的方向和它到原点的距离的角度。在这样的几何中，射线变换及其伴随函数显式地表示为：
P ( f ) ( θ , s ) = ∫ − ∞ ∞ f ( s cos ⁡ ( θ ) − t sin ⁡ ( θ ) , s sin ⁡ ( θ ) + t cos ⁡ ( θ ) ) d t P ∗ ( g ) ( x , y ) = 1 π ∫ 0 π g ( θ , cos ⁡ ( θ ) x + sin ⁡ ( θ ) y ) d θ .

$$\begin{align} \mathcal{P}(f)(\theta,s) &=\int_{-\infty}^{\infty}f\big(s\cos(\theta)-t\sin(\theta),s\sin(\theta)+t\cos(\theta)\big)dt \tag{5} \\ \mathcal{P}^{*}(g)(x,y) &=\frac{1}{\pi}\int_{0}^{\pi}g\big(\theta,\cos(\theta)x+\sin(\theta)y\big)d\theta. \tag{6} \end{align}$$ P(f)(θ,s)P∗(g)(x,y)​=∫−∞∞​f(scos(θ)−tsin(θ),ssin(θ)+tcos(θ))dt=π1​∫0π​g(θ,cos(θ)x+sin(θ)y)dθ.​(5)(6)​
这个参数化决定了线的长度为穿过点(s cos(θ)，s sin(θ))垂直于方向(cos(θ)，sin(θ))的线。

而平面平行光束几何是最简单的工作，它不是最常见的。在收集数据集时使用的一种更典型的采样方案[20,21]是扇形波束几何，具有平面探测器和半径R > 0的圆形轨迹上的源。在扇形波束采集几何中，一个坐标确定线源在圆轨道上的位置，另一个坐标确定线从该点开始的方向。在这个参数化中，射线变换可以写成[11]
P ( f ) ( β , u ) = ∫ − ∞ ∞ f ( R cos ⁡ ( β ) ( 1 − t ) + u t sin ⁡ ( β ) , R sin ⁡ ( β ) ( 1 − t ) − u t cos ⁡ ( β ) ) u 2 + R 2 d t P ∗ ( g ) ( x , y ) = 1 2 π ∫ 0 2 π g ( β , R x sin ⁡ β − y cos ⁡ β R − x cos ⁡ β − y sin ⁡ β ) R ( x sin ⁡ β − y cos ⁡ β ) 2 + ( R − x cos ⁡ β − y sin ⁡ β ) 2 ( R − x cos ⁡ β − y sin ⁡ β ) 2 d β .

$$\begin{align} \mathcal{P}(f)(\beta,u)&=\int_{-\infty}^{\infty}f\big(R\cos(\beta)(1-t)+ut\sin(\beta),R\sin(\beta)(1-t)-ut\cos(\beta)\big)\sqrt{u^{2}+R^{2}}dt \tag{7} \\ \mathcal{P}^{*}(g)(x,y)&=\frac{1}{2\pi}\int_{0}^{2\pi}g\Big(\beta,R\frac{x\sin\beta-y\cos\beta}{R-x\cos\beta-y\sin\beta}\Big)\frac{R\sqrt{(x\sin\beta-y\cos\beta)^{2}+(R-x\cos\beta-y\sin\beta)^{2}}}{(R-x\cos\beta-y\sin\beta)^{2}}d\beta. \tag{8} \end{align}$$ P(f)(β,u)P∗(g)(x,y)​=∫−∞∞​f(Rcos(β)(1−t)+utsin(β),Rsin(β)(1−t)−utcos(β))u2+R2 ​dt=2π1​∫02π​g(β,RR−xcosβ−ysinβxsinβ−ycosβ​)(R−xcosβ−ysinβ)2R(xsinβ−ycosβ)2+(R−xcosβ−ysinβ)2 ​​dβ.​(7)(8)​
这两个几何图形可以通过变换坐标而相互联系起来。设(θ,s)表示平行梁几何中直线的坐标。同样地，设(β, u)表示上述扇形梁几何中直线的坐标。现在我们可以用平行梁坐标来表示扇梁坐标：
$$\beta =\theta +\frac{\pi }{2}-\arctan \left(\frac{s}{\sqrt{R^2-s^2}}\right)\mathrm{and}u=\frac{sR}{\sqrt{R^2-s^2}}.\text{\hspace{1em}(9)}$$

2.5 平面环境下的正弦图外推

目的是在考虑表征射线变换范围的矩条件的同时推断正弦图。该方法是[23,12]对扇形光束几何形状的一种适应(无需重新组合)。

平面射线变换的矩条件通常用平行光束坐标表示[23]：函数$(\theta ,s)\to g(\theta ,s)$在$\mathcal{P}$范围内当且仅当
$${\int }_{-\infty }^{\infty }g(\theta ,s)s^{n}ds\in \mathrm{span}\{e^{ik\theta }:|k|\le n,k+n\mathrm{even}\}\mathrm{for}n\in \mathbb{N}.\text{\hspace{1em}(10)}$$
为了得到其他采集几何形状的相应条件，比如扇形束，我们首先注意到，可以用任意一组多项式$U_n(s)$代替上面的多项式$T_n(s)=s^n$，其中$U_n$为n次，n为偶数时$U_n$为偶数，n为奇数时$U_n$为奇数。当选择$U_n$作为关于某个权重函数$W(s)$的多项式的标准正交族时，这就对g在完全基$U_n(s)W(s)e^{ik\theta }$中的级数展开施加了约束。因此，我们得到了范围的一个新的表征，详细参见[28,37]：$g$位于$\mathcal{P}$的范围内当且应当：
$$g(\theta ,s)=\sum _n\sum _k{c}_{n,k}{e}^{ik\theta }{U}_n(s)W(s)\text{\hspace{1em}(11)}$$
其中$c_{n,k}$仅当$|k|\le n$且$k+n$为偶数时才为非零。通过忽略项，其中$n>N$，我们得到基于有限系数集$c=(c_{0,0},c_{1,0},\dots ,c_{N，n}{)}^T$的$g$的近似。在这种形式下，力矩条件很容易转移到其他几何形状，例如扇形束。

我们可以利用这一点来推断一个次采样的正弦图。如果$g$在其定义域的一个子集上已知，我们估计系数$c_{n,k}$以最小化公式(11)中的误差，然后使用相同的系数在其定义域的其余部分估计$g$。我们可以在离散情况下将其写成矩阵形式，其中$g$是表示正弦图值的数组：
$$g=\mathsf{B}\cdot c.\text{\hspace{1em}(12)}$$
这里，$\mathsf{B}$是一个矩阵，其列由基函数$e^{ik\theta }{U}_n(s)W(s)$组成，这些基函数在扇形波束坐标下的格子的上采样(这些值可以使用第2.4节的坐标映射来计算)。设$g_k$为已知区域的正弦图样本，$g_u$为未知区域的正弦图样本。
g = [ g k g u ] and B = [ B k B u ] . (13) g=

$$\begin{bmatrix}g_k\\g_u\end{bmatrix}$$ \quad\text{and}\quad\mathsf B= $$\begin{bmatrix}\mathsf B_k\\\mathsf B_u\end{bmatrix}$$ . \tag{13} g=[gk​gu​​]andB=[Bk​Bu​​].(13)
我们求出最小二乘解
$${\mathsf{B}}_k^{\ast }\cdot {\mathsf{B}}_k\cdot \hat{c}={\mathsf{B}}_k^{\ast }\cdot g_k\text{\hspace{1em}(14)}$$
然后
$$g_u={\mathsf{B}}_u\cdot \hat{c}.\text{\hspace{1em}(15)}$$
这里${\mathsf{B}}^{\ast }$是$\mathsf{B}$的伴随函数，它的作用是将一个正弦图$g$变换成一个向量，这个向量的分量是$g$与不同基函数$e^{ik\theta }{U}_n(s)W(s)$的内积。

由于有限数据层析的严重病态性质，除非使用某种正则化，否则公式(14)的解是无用的。我们使用了Tikhonov正则化，这意味着执行正则化范围条件问题相当于使用公式(13)中的基函数的岭回归。作为多项式的标准正交族，我们尝试了Legendre多项式和Chebyshev多项式，发现Chebyshev多项式产生的结果明显更稳健。

注意，就像反向投影算子一样，我们不需要存储显式矩阵$\mathsf{B}$，我们只需要考虑它对系数向量$c$的作用以及它的伴随函数对正弦图$g$的作用。我们确实需要矩阵${\mathsf{B}}_k^{\ast }\cdot {\mathsf{B}}_k$，它的大小取决于我们选择的$N$。我们发现$N=50$产生了一个足够的近似值，结果是一个可管理的$650\times 650$矩阵。此外，当N = 50时，运算符$\mathsf{B}$和${\mathsf{B}}^{\ast }$的计算复杂度低于反向投影运算符，导致重建过程的计算负担与正常FBP重建的计算负担相同。

${\mathsf{B}}_k^{\ast }\cdot {\mathsf{B}}_k$的计算是一个相对昂贵的计算操作，然而，这个矩阵只需要在训练期间计算一次，然后可以存储在推理期间使用。因此，与训练网络的计算负荷相比，计算是便宜的。

2.6 网络结构

在数据分布没有先验条件的情况下，从某种意义上说，执行范围条件是我们能做的最好的事情。然而，由于数据来自现实世界的对象，因此数据中隐藏着我们可以利用的额外信息。这个想法是首先做出一个保留范围条件的初始外推，然后用一个从数据中学习的神经网络来弥补外推的不足。FBP重构算子的一个直接推广是用卷积网络代替数据空间中的卷积，该卷积网络可以通过针对大型训练集进行训练来适应。

由于FBP中使用的滤波器通常更容易在傅里叶域中表示，因此使用通过在频域中滤波来操作的网络似乎是合理的。这种架构在其他任务中取得成功的一个例子是傅里叶神经算子(FNO)[18]。在算子学习的背景下引入了FNO体系结构。更具体地说，这个想法是学习一个神经网络，该神经网络以网格独立的方式从其初始/边界条件输出偏微分方程的解(可以在一个网格上训练并在另一个网格上评估)。此外，FNO结构很适合学习傅里叶积分算子。射线变换的逆和FBP都是傅里叶积分算子的例子[16]。

FNO的关键组成部分是频谱卷积，它本质上是一个普通的卷积层，但每次卷积都在频域进行。给定一个有$C_{in}$通道的输入，作为函数$f:[0,1]\to {\mathbb{R}}^{C_{\mathrm{in}}}\in L_1([0,1],{\mathbb{R}}^{C_{\mathrm{in}}})$，频谱卷积层产生任意多个通道$C_{out}$的输出$h:[0,1]\to {\mathbb{R}}_{\mathrm{out}}^C\in L_1([0,1],{\mathbb{R}}^{C_{\mathrm{out}}})$。对于每个通道i，输出由
$$\hat{h}(\omega {)}_i=\sum _{j=1}^{C_{\mathrm{in}}}{\hat{k}}_{ij}(\omega )\hat{f}(\omega {)}_j\text{\hspace{1em}(16)}$$
这相当于将输入$f$与过滤器$k_{ij}={\mathcal{F}}^{-1}({\hat{k}}_{ij})$进行卷积。由于函数域是紧致的，我们考虑离散频率$\omega$。

在我们的FNO-BP网络中，FNO将正弦图的每一列，即从相同旋转角度收集的每一组测量值，视为单独的通道。在这种情况下，将从L个方向采样的正弦图视为函数$\mathbb{R}\to {\mathbb{R}}^L$的离散样本，该函数将空间坐标$u$映射到相应的样本集$g(:,u)$：
u ↦ [ g ( β 1 , u ) ⋮ g ( β L , u ) ] . (17) u\mapsto

$$\begin{bmatrix}g(\beta_1,u)\\\vdots\\g(\beta_L,u)\end{bmatrix}$$ . \tag{17} u↦ ​g(β1​,u)⋮g(βL​,u)​ ​.(17)
总之，该模型可以写成
$$g\mapsto \text{ReLU}({\mathcal{P}}^{\ast }(k_{\text{RL}}\ast g+\text{FNO}(g)\left)\right)\text{\hspace{1em}(18)}$$
其中$k_{RL}$是常用的Ram-Lak滤波器，相应的网络架构参见图1。FNO层本身由一组光谱卷积层组成，这些层之间有高斯误差线性单元(GELU)激活函数以及线性跳过连接层，参见[18]。在探索了一些超参数的配置之后，我们决定使用具有三个隐藏层的FNO，每个隐藏层由C = 60个通道组成，每个核- $k_{ij}$由280个值参数化，这些值对应于滤波器$k_{ij}$的前280个傅立叶系数。

FNO以1D和2D的形式存在，尽管我们的数据本质上是2D的(sinogram是关于线和角度的函数)，我们的FNO-BP采用1D的形式有两个原因。首先，低维傅里叶变换更有效，其次，FNO执行平移不变变换，这意味着它会错误地平等对待正弦图的所有区域。但事实并非如此，因为在需要特别注意的集中区域，有限角度断层扫描会导致数据缺失或质量低下。

2.7 实现

该模型是使用PyTorch和ODL在Python中实现的[2]，代码可以在GitHub上找到[[32]](https: //github. com /EricOldgren /KEX—CT reconstruction)。模型在2500个phantom的合成数据集上训练了30个epoch，使用Adam优化器，学习率为3·10−5，否则使用默认PyTorch参数。合成的数据由圆形的phantom组成，其中有不同类型的简单孔。请注意，phantom中心各不相同，孔的大小和几何形状各不相同，分为矩形和椭圆形。这与[HTC-2022](https://github. com/99991/HTC2022-TUD-HHU-version-1)的测试数据类似，尽管测试数据包含一些内部结构复杂得多的phantom。

3 Results

FNO-BP模型是根据模拟HTC-2022挑战[20]中可用的教学数据而构建的合成数据[32]进行训练的。然后在HTC-2022测试数据[20]上评估FNO-BP模型，以及使用第2.5节中描述的正弦图外推法的标准FBP和FBP。所得到的重建将所有负像素设置为零，并使用Otsu的方法进行阈值处理[26]。计分方法如表1所示，采用HTC-2022[17]比赛时采用的马修相关系数。

在表1中，可以看到FNO-BP模型的分数与常规FBP相比，有和没有正弦图外推，以及向HTC-2022提交的获奖作品[10]。对不同数量的可用角度(从90◦到30◦)进行了测试，得分作为每个难度级别给出的三个测试图像的平均值。应该注意的是，这两个表，已经90°的角跨度(所需的180°的一半)是相当有限的。

在表3中，我们看到了测试数据集中不同数据的结果重建。根据重构方法和可用的角跨度对图像进行排序。从这些图像中，我们可以看到，通过正弦图外推利用距离条件(FNO-BP和FBP+range cond)的两种方法都保持了椎间盘的一般形状。FNO-BP专门针对圆形phantom进行训练，这解释了这一事实，但FBP+range cond除外。不是，这意味着部分形状保持能力可能是由于正弦图外推。请注意，不同图像的phantom中心是不同的，这是在重建中准确捕获的。两种FBP方法都对重建进行 smear，但是当使用范围条件时， smear被包含在disc中。虽然FBP方法似乎夸大了较低角跨度的孔，但FNO-BP似乎使它们更小。

从表1中可以看出，FNO-BP方法的重建效果比FBP方法中的任何一种都好，但比中标的HTC-2022方法差[10]。与FBP方法相比，FNO-BP的能力在表2和表3中得到了证实，对于所选的地面事实，FNO-BP实现了与HTC-2022的获胜提交相似的精确重建，其角度跨度为50◦，而标准FBP和sinogram 外推法之后的FBP在角度跨度为80◦或90◦时开始落后。对于最小角度，与获胜的HTC-2022提交的方案相比，上述三种方法的表现似乎都不佳。

在表2和表3中，我们研究了不同方法重建特定对象的能力，但不同的角跨度。这里我们也放弃了任何阈值。从这些图像中，我们证实，即使没有阈值化，利用距离条件的方法在重建后也在很大程度上保持了物体的形状，特别是FNO-BP，即使在低角跨度的情况下，也非常接近地保持了物体的形状，而不会在物体半径外产生太多的smear。所有方法的孔的可识别性都是相似的，尽管有范围条件的FBP似乎引入了更多的伪影，在物体中心周围形成环的形状。

4 Conclusion

FNO-BP是一种用于层析重建的新型神经网络架构，其设计在推理过程中具有与FBP重建相同量级的计算复杂度。因此，在计算上，FNO-BP比基于模型的方法(迭代和变分)和当代领域信息学习方法(如LPD[3])要高效得多，这些方法需要多个正向和反向投影。