论文笔记 圆轨道MDCT/CBCT锥角伪影及如何抑制[1]

标题

On the data acquisition, image reconstruction, cone beam artifacts, and their suppression in axial MDCT and CBCT – A review

摘要

近年来，多排CT（MDCT）和锥束CT（CBCT）在临床诊断中扮演了越来越重要的角色，但这些CT成像技术本身会导致一些伪影称为CB伪影。本文探究了圆轨道CT中CB伪影形成的根本原因，并针对这些伪影提出了一些可能的解决方案。

Introduction

介绍了CT的发展史，从单排到多排再到锥束。圆轨道MDCT和CBCT的一个不发避免的问题就是他们不满足数据充分条件（Tuy条件），这个问题也就导致了CB伪影的出现。

数据不充分——CB伪影的根本原因

Tuy-Smith数据充分性条件要求与待成像物体相交的每一个平面都与射线源的轨迹至少存在一个角点，满足这个条件的射线源轨迹才能重建出不含伪影的图像。
作者将圆轨迹CBCT不满足Tuy-Smith DSC解释为如果要从从圆轨迹上的的锥形束中找出两条共轭光线，除非这两条光线位于圆轨道所在平面上，否则这两条光线不会共线（也就是数据不一致data inconsistency）。这种解释在我看来比较奇怪，我的解释是与圆轨道平面平行的那些与物体相交的平面显然与射线源轨迹不存在交点，因此，圆轨道不满足Tuy条件。这里作者的解释或许更加接近伪影产生的本质。作者在这里提到的data insufficiency和data inconsistency其实是同一问题的两个方面，这个问题在Radon域中表现为整个Radon空间不能被投影数据完全覆盖（存在null space）。这里需要注意的是Tuy条件是在探测器为2-D理想探测器的前提下提出的，也就是大小要能包含整个物体的投影，采样率满足奈奎斯特采样定理。如果这两个条件任意一个不满足都会产生额外的伪影。数据在水平面上截断以及抑制方法在之前的文献中已经有所介绍（另外这种情况在扇形束中也会出现不能算是锥束伪影），在本文中主要探讨longitudinal方向的截断。

Radon transform and Radon shells on Radon space

这部分主要讲了一些三维Radon变换与反变换的公式，公式都是在三维球面坐标系下定义的。
三维X-ray变换就是沿直线L的线积分$$Xf(\overrightarrow{r})={\int }_{-\infty }^{\infty }f(\overrightarrow{r_0}+\overrightarrow{k}\cdot l)dl$$其中$\overrightarrow{r_0}$是直线$L$上的点，$\overrightarrow{k}$是单位方向向量。
三维Radon变换公式如下$$Rf(\overrightarrow{r})={\int }_0^{2\pi }{\int }_{-\pi /2}^{\pi /2}{\int }_0^{\infty }f(\overrightarrow{r})\delta (\overrightarrow{r}\cdot \overrightarrow{n}-\rho )\rho d\rho d\theta d\phi$$其中$$\overrightarrow{n}=(sin\theta cos\phi ,sin\theta sin\phi ,cos\theta )$$其实就是垂直于$\overrightarrow{n}$的与原点距离为$\rho$的平面上的面积分。
反变换的公式为$$f(\overrightarrow{r})=\frac{-1}{8\pi }{\int }_{-\pi /2}^{\pi /2}{\int }_0^{2\pi }{\int }_0^{\infty }\frac{{\partial }^2}{\partial {\rho }^2}Rf(\overrightarrow{r})rdrd\beta d\alpha$$其中$\overrightarrow{r}$定义为$$\overrightarrow{r}=r\cdot (sin\beta cos\alpha ,sin\beta sin\alpha ,cos\beta )$$即重建一个三维点需要所有包含该点的Radon平面，这些平面在Radon域构成的shell称为Radon shell。

3-D傅里叶切片定理

三维下的傅里叶切片定理可以理解为Radon域过原点直线上的值（相当于数据域垂直于该直线的所有平面面积分）的一维傅里叶变换与原数据三维傅里叶变换后对应直线上的值相等。即$$F_{\rho }(Rf(\rho ,\overrightarrow{n}))=F(u,v,w){|}_{(\rho ,\overrightarrow{n})}$$其中$(u,v,w)$为$(\rho ,\overrightarrow{n})$对应的直角坐标。

Radon和傅里叶空间下的数据不充分的表现

这部分首先要理解一次锥束投影对应的Radon shell，这个可以与上图对应着理解。$S$为射线源，$O$为坐标原点，这一次投影可以填充满Radon域以$SO$为直径的球形范围。射线源旋转一圈可以得到如下图所示的Radon shell.

可以看到不管旋转多少圈，在Radon域都会存在一个null space。这种null space就是数据不充分的根本原因。在傅里叶域中，这种null space表现为一个锥角为$si{n}^{-1}(R/D)$的圆锥，其中$R$为重建的半径，$D$为$DSO$。重建的半径越大，傅里叶域中的null space越大，CB伪影也就越严重。

射线源轨迹与数据充分性条件

对于MDCT/CBCT，满足条件且最有效的方式是螺旋扫描。此外，还有circle-plus-line, circle-plus-arc, circle-plus-two-arcs, dual orthogonal circles, dual orthogonal circles plus line, dual elliptical等。螺旋CT的问题在于不适合处理存在不自觉运动的投影数据（例如心脏CT），对于这种情况，医生一般会选择使用大锥角的圆轨迹MDCT/CBCT。

CB伪影的表现形式和形态学

数据不充分导致的CB伪影

这种伪影主要表现为几何扭曲和亮度降低，并且这种伪影会随着与射线源轨迹平面的距离增大而变得越来越严重。具体情况如下图。

数据不一致导致的CB伪影

如之前所说，数据不一致指的是重建平面上共轭光线其实并不共线，如下图所示。

这种现象的后果如下图所示，主要是靠近圆轨道的骨骼附近存在一些doubled structure。（之前不是说不充分跟不一致是同一问题的不同表现吗？这里怎么分开了）

longitudinal方向数据截断导致的CB伪影

前两种伪影即便是在理想的探测器下依旧会出现，而数据截断的伪影是由探测器不够长造成的。目前SOTA的MDCT大概能重建160mm的物体，即足够重建大多数病人的心脏，而CBCT可以重建的范围更大一些，大概能包含肺部。数据截断造成伪影的原因如下图所示，主要是物体中部分点与射线源的连线延伸到了探测器之外。

longitudinal方向采样率太低导致的CB伪影

对于32排以下的MDCT，longitudinal aliasing伪影（例如windmill，pinwheel和bear-claw）仅在螺旋扫描中出现（这个是为什么？？）。但随着排数的增加（超过64排），axial CT的反投影必须要根据3D几何来完成（之前都是当做多排2D的来做的吗？），这也使longitudinal插值变得必要起来。下图是windmill伪影的例子。

这种aliasing伪影在形态学上通常并不与任何病例现象相似，但在诊断中也很让人头疼。

抑制伪影的方法

基于锥角的3D加权

这种方法针对数据不一致导致的伪影。我们认为锥角小的光线比锥角大的光线更加可靠，所以在重建中，我们给小锥角光线一个较大的权值，大锥角光线小的权值，这样就可以一定程度上抑制伪影。这种方法在目标物体较小（整体锥角比较小）的时候可以很好的抑制伪影，但当锥角本身很大的时候效果不是很明显。

启发式的探测器扩展

主要针对数据截断造成的伪影。可能的方法有复制最外侧一排的数据。但这种方法并不能保证扩展出来的数据是正确的（不知道有没有人用深度学习做的）。

处理longitudinal采样不足的方法

最有效的方法显然是硬件上增加采样率，但这种方法受限于制造水平和成本，不太现实。其他的基于硬件的方法主要有飞焦点，可以近似采样率翻倍。效果如下图。

此外，还有基于软件的方法如低通滤波，但这种方法不像飞焦点那样既保留特征又去除伪影。最近，有人提出了一种探测器交错扫描的算法可以有效的抑制aliasing伪影[1]。这种方法使用基于晶格的采样框架，抑制伪影和图像分辨率之间的trade-off可以分析优化。效果如下图所示。

用插值来填满null space

2D的Radon变换可以将3D空间域的剧烈变化转换成3DRadon域的平缓变化，这使用插值来填充null space变得可能。但目前为止，这方面并没有取得太大的成功。

two-pass重建

CB伪影在图像中表现为在高密度或低密度结构附近的阴影或者条状伪影。进一步的研究发现这种伪影在形态上与beam hardening伪影极为相似，而beam hardening伪影可以用iterative bone operation来抑制，因此，我们认为可以用类似的方法来抑制CB伪影。two-pass的重建可以分为两步，第一步是用传统的FDK或者CB-FBP算法来进行重建，第二步是将高密度与低密度结构分割出来以后进行正投影再用正投影的结果进行重建。如果将第一步重建的结果标为$f$，则第二步的重建可以表示为如下形式$$f_m=S(f)$$$$g_m=P(f_m)$$$$f_m^{\prime }=B(g_m)$$$$f_a=f_m^{\prime }-F_1(f_m)$$$$f_c=f-F_2(f_a)$$ 其中$S$为分割算子，$P$为投影算子，$B$为反投影算子。$F_1$用于处理由于正反投影造成的空间分辨率下降，$F_2$是一个低通滤波器用于去除高频的伪影。
这种two-pass的重建方法可以有效的去除极端密度物体附近的条状伪影，但对远离射线源平面处的形变和CT值降低并没有太大的作用，此外，分割过程中使用的阈值对算法的效果有较大的影响。具体效果如下图。

文章后面分析了一大堆two-pass方法有效的原因，这里就不写了。

short scan数据的锥束伪影及抑制

short scan指的是CT扫描不扫满360度的数据。主要是为了应对人体脏器或组织周期性的运动（心跳，呼吸）以及一些不自主的运动（肠收缩等）或者是跟踪一些运动（肺癌肿瘤等），也就是提高时间分辨率。

half/short scan

一般最常用的扫描方式是扫描平行光180度的范围，对应扇束中$180°+{\gamma }_m$的范围，并通过Parker权重来处理数据冗余。Parker权重后来也拓展到partial scan的情况下，也就是大于180度小于360度的扫描。
MDCT+Parker weighting的做法现在在心血管成像中非常常见，其中还可以分为单源跟双源两种形式，后者可以达到比前者更高的时间分辨率。圆轨道CT相比螺旋CT与心脏门控结合更加容易，可以根据门控来降低病人受到的辐射剂量。
但是，理论上Parker weighting仅仅在轨道平面上是准确的，因此，会导致由于处理数据冗余不当而造成的伪影。对此，学者们也提出了一些算法，这些算法使用平移不变的滤波取代了加权来应对数据冗余的问题[2-4]。
由于扫描不满360度不存在共轭光线，3D加权的方法在这里无法应用。因此，重建出来的图像的质量会降低。但是，在心脏和肺部CT中，CB伪影往往不是最主要的问题，因为骨骼结构比较少。

super-short scan

这种扫描方式主要针对只有一小部分图像ROI需要进行重建的情况。如下图所示，如果只重建物体的一个ROI，且该ROI偏离圆轨道中心，那么即便是小于$180°+{\gamma }_m$的扇形扫描范围也能满足重建的需求。

用half way平移的探测器来扩大FOV

由于360度的扫描会有一倍的数据冗余，我们可以将扫描分成两部分，在每一部分中探测器放在不同的位置。以此，在探测器大小有限的情况下扩大扫描的FOV。

[1] Xie H, Tang X. Optimization of data acquisition in axial CT under the framework of sampling on lattice for suppression of aliasing artifacts with algorithmic detector interlacing. Med Phys. 2017;44:101914.
[2] Zhuang T, Zambelli J, Nett B, Leng S, Chen G. Exact and approximate cone-beam reconstruction algorithms for C-arm based cone-beam CT using a two-concentric-arc source trajectory. SPIE Proc. Vol. 6913, 691321; 2008:12. https://doi.org/10.1117/12.772390
[3] Nett BE, Zhuang TL, Leng S, Chen GH. Arc-based cone-beam recon-struction algorithm using an equal weighting scheme. J X-ray Sci Tech. 2007;15:19–48.
[4] Zhu L, Yoon S, Fahrig R. A short-scan reconstruction for cone-beam CT using shift-invariant FBP and equal weighting. Med Phys. 2006;34:4422–4438.