浅谈CT图像重建的数理基础

原创 刘宜 XI区 收录于合集#原理 20 个 #科普 66 个

长按识别下方二维码关注我

大家好，欢迎大家访问XI区！

01

引言

想要知道一个西瓜的内部长什么样，皮厚不厚，有没有籽，最好的办法就是把西瓜切开来，对于临床来说，为了获得患者的体内信息而把病人切开是不现实的。传统的X线设备能够获得一系列密度叠加的像，但不能得到截面影像，对于空间位置的识别不够精细，如何不用动刀就获得准确的截面诊断信息成为了亟待解决的问题。CT的出现给我们提供了一种全新的成像手段，它能获取物体的截面信息，辅助临床诊疗。

CT通过扫描来获取投影值，根据投影值来求解成像剖面上衰减系数的分布，投影数据是用探测器采集的，探测器上所成的像是物体叠加后的图像。CT通常所说的正弦图是通过拉东变换得到的，它引入了一个狄拉克函数来当作点源，投影数据是对物体的线积分。相当于把函数f(x, y)通过线积分表示成了另外的一种直线参数的形式，把二维平面(x, y)坐标系映射到了直线参数（θ, P）坐标系。从物体的投影数据来得到物体的内部断层成像的过程就称之为图像重建。图像重建是一个数学工序，它的目的是把投影数据中的叠加效应除去，并营造出一个无重叠的原本物体的截面图像，这个数学工序叫做算法。

02

反投影与滤波反投影

2.1 解方程组

早期的算法研究以解方程组和反投影方法为主。

扫描得到的投影值

假设我们通过扫描得到了上面4个投影值，那么任务就是要还原出真实的衰减系数分布，对于上图的矩阵，我们可以列出一系列方程组：

尝试解这个方程组会发现不只一个解，从线性代数的角度上，这个方程组属于欠定方程组，4个方程不独立，把上面两个方程相加，把下面两个方程相加，这两个方程就变成一样的了，方程数小于未知数个数，它会有很多个解，但二维图像只能有一种可能，那么要能够得出唯一的解，就需要获得更多的投影值，也就是添加更多的方程组，对于CT来说就是变换不同的角度扫描，采集更多的投影数据。

而采集过程中不可避免的会引入噪声，导致方程组不相容，对于超定方程组来说是无解的，我们很多情况下只能得到它的近似解。而且CT的矩阵十分庞大，运算量非常大，这也导致方程的求解变得很不现实。

2.2 反投影

反投影提供一种求解衰减系数的新思路，直接反投影将获得投影值反投回去，反投其实就相当于把投影值分配给各个像素，由于事先并不知道真实的衰减系数分布，只能均匀的“涂抹”，包括原先像素值为零的点，这样就会导致星状伪迹，造成图像模糊，边界不清。

计算机模拟直接反投影导致的图像模糊、边缘失锐

2.3 中心切片定理

为了解决这个问题，Bracewell利用傅里叶空间得出了中心切片定理和傅里叶重构算法，从数学上证明了，根据投影值是完全可以重建出原始图像的。二维图像的中心切片定理指出：密度函数f(x, y)的投影函数��(��)的傅里叶变换��(��)等于函数f(x, y)的傅里叶变换F(��x, ��y)沿探测器平行方向过原点的片段。

图解中心切片定理

如果探测器围绕物体旋转至少180°，那么物体的二维傅里叶变换F(��x, ��y)就可以通过探测器探测得到的��(��)覆盖整个空间。一旦函数F(��x, ��y)为已知函数，原本的图像函数f(x, y)就可以通过二维傅里叶反变换获得。

探测器在每个探测方向上向傅里叶空间添入一条线。当这些线覆盖了整个傅里叶空间后，原本图像可由二维傅里叶反变换获得重建。

熟悉磁共振成像的学者看到这个图可能会很眼熟，联想到磁共振也有傅里叶空间的概念——“Κ空间”，“刀锋”技术和“螺旋桨”技术正是利用了这样的放射状Κ空间填充模式，Κ空间的中心决定了图像的对比，Κ空间的边缘决定了图像的细节。这种技术下Κ空间的中心区域有大量的信息重叠，图像有较高的信噪比，降低了运动伪影。这种技术也叫做“极射线扫描方法”，极是极坐标的意思。在这个方法中，x梯度和y梯度总是同时开启，同时关闭，没有相位编码，只有读出编码，读出编码的方向取决于x梯度和y梯度。理论上来讲可以用FBP算法来重建这个MR图像。

极射线扫描的K空间采样示意图

从中心切片定理的角度上，反投影可以理解为：在一个方向上的反投影等价于向即傅里叶平面添加一个中心片段，做180°的反投影就可以完整得到二维傅里叶变换函数。接下来只需进行二维傅里叶逆变换即可将F(��x, ��y)还原到f(x, y)。

当我们往傅里叶平面里一条一条地添加“中心片段”时，中心片段在傅里叶平面原点的密度高于在远离原点的区域的密度。在傅里叶空间原点附近的区域是低频区域，图像低频信息表示图像中灰度值变化缓慢的区域，它形成了图像基本的灰度等级，是图像的大致概貌和近似信息。感性的来讲，对低频分量的过分加权使得会使得反映图像边缘、噪声和细节的高频分量减少，导致图像变得边界不清，过分柔和，也就是看起来模糊。这也解释了前面所说为什么光靠反投影是不能得到原本图像的。

2.4滤波反投影

为了消除模糊的效果，需要对傅里叶空间要进行加权校正，使其高低频分量密度均匀。一个典型的做法是把投影数据p(s, θ)的一维傅里叶变换��(��, θ)乘以|��|。然后，再对乘积|��|��(��, θ)做一维傅里叶反变换。对投影数据做了处理(即滤波)之后，把处理后(即滤波后)的数据做反投影。这样就得到了重建的图像f(x, y)。这个方案通常被称作先滤波后反投影(Filtered back projection，FBP)算法，函数|ω|被称作斜坡滤波器。斜坡滤波器其实是个高通滤波器，它抑制低频成分，增强高频成分。

斜坡滤波有不只一种实现方法，实际上并不需要用傅里叶变换来做斜坡滤波。根据傅里叶变换理论，在��域中做乘法等价于在s域中做卷积。也就是说斜坡滤波既可以在频域中用乘法实现，也可以在空间域中用卷积实现

其中h(s)是卷积积分中的卷积核，它是H(��)=|��|的一维傅里叶反变换。我们在重建参数里见到的Hr40，Bv36，Br60等一系列卷积核都与之相关。

计算机模拟不同角度数对重建结果的影响

并且我们可以自由地更换斜坡滤波和反卷积的次序，斜坡滤波又可以进一步分解成求导运算和希尔伯特变换，这样，我们利用不同的排列可以得到更多的重建算法。所有算法都包含了反投影的步骤，那么不进行反投影可以得到图像吗，答案是肯定的。利用中心切片定理可以让反投影在傅里叶域中实现，但是这种方法并不常用，因为这个方法是从极坐标系(��, θ)到直角坐标系(��x, ��y)的插值运算，而插值运算往往会引入较大的误差，导致结果欠佳。

平行束图像重建解析算法一览表

直接反投影会带来星状伪迹，导致图像边缘失锐，通过高通滤波后又会导致高频噪声的放大，实际应用中，我们通过乘以一个窗函数来抑制噪声。窗函数是通过控制带宽来把高频信息和和高频噪声一起舍弃，引入窗函数伴随着高频信息的丢失，这又有可能带来振铃效应，解决办法是选用合适的窗函数（韩窗、汉宁窗等）。

引入窗函数导致的振铃效应

添加了不同窗函数的滤波反投影重建效果

2.5扇形束重建与三维平行光束重建

需要指出的是，上面我们所说的都是针对平行束重建而言，平行束图像重建的理论基础中心切片定理在扇形束图像重建中并没有能与之对应的。解决办法利用了点扩散函数的移动不变性，通过引入雅可比因子变量代换来得到扇形束的图像重建算法。

如果扇形束的焦点轨道是整个圆周，这就叫全扫描。如果扇形束的焦点轨道是部分圆周，这就叫短扫描。扫描角度β可以小于360°，这种扫描方式叫做短扫描。角度β的最小取值范围取决于数据采集时物体与探测器之间的几何关系，可能大于180°，可能等于180°，也可能小于180°，通常所说的半扫描重建就是扫描覆盖180°来采集数据进行重建。

三维平行光束重建和锥形束重建算法较为复杂，限于篇幅这里不再作解释，值得注意的是对于大的锥形束张角，锥形束重建算法的伪影很严重，在远离中心平面（也就是轨道平面）的区域尤为明显。使用宽体探测器需要用到更大张角，不可避免的会带来这种伪影，导致图像的前后一部分不可用，这也是宽体探测器CT的一个限制，为了消除这种影响，锥形束伪影的校正技术一度成为研究热点。

计算机模拟下不同的锥形束张角的图像重建结果

2.6小结

至此，我们可以总结出滤波反投影算法的一些特点：

1) FBP是一种解析算法，重建速度快，对计算设备要求低，几乎扫描结束的同时就出来图像了，我们看到的CT上扫描实时预览图多为这种方法重建得到；

2) FBP算法依赖于投影数据的完整性，因此必须有足够的剂量才能重建出质量较好的图像；

3) FBP没有考虑噪声模型对图像的影响，对噪声的控制只能通过在频域上加窗函数来实现，也无法去伪影。

03

迭代重建

3.1 商用迭代重建算法及其实现

正是由于诸多局限，且随着计算机的发展和人们对低剂量扫描需要的日益增长，迭代重建重新兴起，各个厂家都开发出了各自的商用迭代重建算法。

商用迭代重建算法及其出现的时间

迭代重建算法的通用流程

迭代重建的算法繁多，但万变不离其宗，图像重建的解析算法与图像重建的迭代算法之间的主要区别在于对图像的模型化。在解析算法中，我们假设图像是连续的，每个像素只是一个点，这些离散的点是以图像显示为目的的。这些点的选择可以是随意的，与图像重建无关。但是，在迭代算法中，每个像素是个小面积。这个小面积在计算当前图像的投影数据时要用到，像素模型对重建图像的质量好坏影响很大。

图像重建的解析算法与迭代算法之间的另一个区别是解析算法着力对一个积分方程求解，而迭代算法着力对一个线性方程组求解，我们可以用一个线性方程组把真实的成像几何结构与成像物理效应模型化，解一个线性方程组比解一个积分方程要容易，意味着迭代算法比解析算法得到的重建图像更为准确。

迭代算法是用来对目标函数做极小化的，目标函数可以有效地考虑测量值的噪声影响。比起解析算法只能通过加窗控制噪声，迭代算法控制噪声的办法更为丰富，诸如：

1) 算法收敛以前就提前停止迭代运算；

2) 在投影/反投影对中对成像的几何结构和物理效应做更精确的建模；

3) 用平滑而互相重叠的像素来代替平坦而不重叠的像素；

4) 应用正确的噪声模型来建立目标函数；

5) 使用先验知识（贝叶斯法则）。

3.2 各类商用迭代重建算法特点

各种迭代重建算法及其性能对比

上表中按时间顺序列举了目前14种商用CT的重建算法，分类并简要介绍了它们的特点。首个商业迭代算法直到2009年才出现，尽管长期以来人们都知道它的理论优势，但由于运算能力的限制，第一代商用迭代算法，包括自适应统计迭代（ASiR；GE Healthcare，Waukesha，Wis）、图像空间迭代（IRIS；Siemens Healthineers，Forchheim，Germany）、自适应剂量减少迭代（AIDR；Canon Medical Systems，Otawara，Japan）和iDose（Philips Healthcare，Bothell，Wash），都不是真正意义上的迭代算法，仅是对投影数据或图像数据进行处理的迭代降噪方法。

滤波反投影、加权滤波反投影、迭代重建算法的效果对比

由于图像数据和投影数据的相互性，出现了混合迭代算法，这种算法在图像空间上迭代以消除二次噪声，在投影空间上迭代减少伪影（例如，Sinogram Affirmed Iterative Reconstruction[SAFIRE；Siemens Healthineers]）。随着计算机的发展，运算能力提高，许多CT制造商已经过渡到基于模型的IR，如ASiR-V（GE Healthcare）、Iterative Model Reconstruction（IMR；Philips Healthcare）或Forward Projected Model-Based Iterative Reconstruction SoluTion（FIRST；Canon Medical Systems）。

对于混合迭代算法，这类IR算法在原始数据域去伪影，在图像数据域降噪声，结合了光子统计、CT系统的详细模型和扫描对象的模型，重建速度较快。

使用SAFIRE迭代重建的图像相比FBP具有更低的噪声，允许更低的辐射剂量

对于基于模型的迭代算法，这种IR算法考虑更多CT系统的细节，如探测器和焦点的形状，X射线束的能量谱等，而付出的代价是重建时间相对较长。

基于模型的迭代考虑了更多的CT系统细节，图像质量较FBP更好，但重建时间也增加了

百花齐放终将殊途同归，任何一种算法都着力于提高图像质量，降低噪声、去除伪影、优化剂量，目前没有证据表明混合迭代算法的性能低于基于模型的迭代算法。

各类商用迭代重建算法简要解释：

IRIS——Iterative Reconstruction in Image Space（图像空间迭代）

SAFIRE——Sinogram Affirmed Iterative Reconstruction（正弦波图形法迭代）

ASIR——Adaptive Statistical Iterative Reconstruction（自适应统计迭代）

MBIR——Model Based Iterative Reconstruction（基于模型的迭代，商品名Veo）

iDose——Double Model Based Iterative Reconstruction（基于双模式迭代）

AIDR——Adaptive Iterative Dose Reduction（自适应剂量减少迭代）

ADMIRE——Advanced Model-Based IterativeReconstruction（高级模型迭代）

IMR——Iterative Model Reconstruction（模型迭代)

FIRST——Forward Projected Model-Based Iterative Reconstruction SoluTion（前向投影模型迭代）

3.3 迭代重建算法的优势

迭代算法有很多优点：

1) 允许对X射线源和探测器进行建模，从而提高空间分辨率和重建准确性；

2) 充分考虑光子的统计特性，增加低噪声投影权重并降低高噪声投影权重，从而减少伪影，提高剂量利用率；

3) 考虑真实物理模型，可以在保持解剖边缘锐利的同时减低图像噪音；

4) 迭代算法允许使用不完整的数据进行重建。

3.4 迭代重建算法的噪声

迭代重建算法通常通过对图像特定区域内信号变异程度采取非线性方程变换进行降噪。这也是为什么IR图像特定区域内的空间分辨率依赖于周围结构的对比和噪声水平，IR算法使用的参数决定了这种依赖性。具体来说，算法会设定一个阈值区分高对比度区域与低对比度区域，低于这个阈值被认为是噪声，重建时被平滑，因而降低空间分辨率；高于这个阈值被认为是结构边界，重建时最小程度平滑，因而保持空间分辨率。

这解释了为何IR算法的图像质量和对应可接受的辐射剂量与诊断任务的类型高度相关，对于低对比度的诊断任务，IR算法很难鉴别低对比度物体的真实边界和噪声，因此很难保持空间分辨率。事实上，很多研究已经证明：尽管IR算法可以使得图像质量看起来更好，但过度降低辐射剂量仍会影响低对比度病灶的检出能力。

对比噪声比经常被用作衡量图像质量的指标，因为它易于测量，而且放射科医生认识到，许多疾病的诊断依赖于肉眼观察组织结构之间的CT值差异。然而，由于迭代重建图像中的非线性噪声特性，对比噪声比用来衡量迭代重建图像质量便不再合适了，因为其并不能很好地反映迭代重建图像在低对比度病灶检出能力上的缺陷。

迭代算法的非线性会改变噪声的纹理，造成所谓的“蜡像样伪影”，放射科医生并不喜欢这种现象。噪声功率谱（NPS）可以描述噪声纹理的特点。这种“蜡像样伪影”实质上是噪声功率谱的频峰左移（与FBP相比，IR算法中低频噪声要多得多）带来的视觉结果。

不同层厚和管电流下改变卷积核重建图像数据ＮＰＳ处理结果

图2 Ａ．H30s重建；Ｂ．H60s重建；Ｃ．U90s重建

图3 Ａ．H30s重建；Ｂ．H60s重建；Ｃ．U90s重建

上图卷积核从H30到U90，图像视觉效果“平滑”到“锐利”，噪声功率谱表现为峰值右移，也就是说噪声从低频分量占据主导到高频分量占据主导，图像边界更清晰，量子噪声也更大，空间分辨率提升而密度分辨率下降。对于迭代重建来说，情况则是反过来的，随着峰值的左移，噪声由高频分量多转为低频分量多，部分学者把这称为“噪声纹理”的改变，形成视觉上的“蜡像感”、“卡通感”。非线性运算会影响图像内噪声的空间分布，导致噪声纹理和空间分辨率改变，后者依赖于目标物体和周围背景之间的对比。

举个不恰当的例子，我们形容一个人长得不好看（图像质量不好），因为他脸上有麻子（量子噪声，高频），通过整形手术（迭代重建）以后他还是不好看，因为手术把麻子去掉了，也把面部细节改变（噪声纹理改变，低频），五官不再清晰分明（蜡像感）了。

3.5 小结

IR算法既可以降噪声，又可以去伪影，某些场合允许更低的辐射剂量来达到一致性的诊断结果。降噪效果会根据诊断任务的不同而发生变化，诊断任务可大致分为两类：

1) 高对比度诊断任务，即目标物体与周围背景的CT值之差较大（如胸部CT肺结节的检出，CTA钙化斑块的检出，CT结肠造影息肉的检出等），IR算法对低剂量图像的空间分辨率影响微乎其微，即使辐射剂量降低50%对于诊断效能的影响也很小。

2) 低对比度诊断任务，辐射剂量的轻微降低也会极大地影响细微结构的检出能力，空间分辨率的降低会对诊断产生负面影响（如腹部CT肝转移的检测和头颅CT脑梗死的检出）。

04

深度学习重建

深度学习（DL, Deep Learning）是机器学习（ML, Machine Learning）领域中一个新的研究方向，随着人工智能和计算机视觉的高速发展，基于深度学习的卷积神经网络在医学成像领域的大放异彩，解决了很多复杂的模式识别难题。这些卷积神经网络技术已应用于低剂量CT图像用来减少图像噪声，允许使用不完整或稀疏数据以实现无伪影重建，例如减轻金属植入物伪影。

最近，CT制造商开发了使用深度学习卷积神经网络来进行图像重建或降噪的算法，而非传统的迭代重建，包括TrueFidity（GE Healthcare）和Advanced intelligent Clear IQ Engine（AiCE；Canon Medical Systems）。其高度非线性噪声和空间分辨率特性需要进一步的技术验证，在临床应用推广之前需要进行严格评估。

右盆腔肿块、腹部肿块和肝转移瘤患者的增强CT图像（箭头），使用A 滤波反投影、B 迭代重建（ASiR-V 50%，GE Healthcare）和C 深度学习重建（TrueFidelity；GE Healthcare）。注意与细微肝转移瘤边缘相关的边界（容积CT剂量指数6.73mGy，层厚0.625mm）。

现有阶段的深度学习算法还处于探索期，代表了未来的一种可能的方向，有很多方面需要优化，结果可能并没有想象的那么完美，迭代重建算法依旧是现在的最可靠的CT重建算法，深度学习重建并不完全优于迭代重建算法，这种新的图像重建模式是否比当前的迭代算法具有更好的性能还有待观察。

更多干货，关注XI区！

作者简介

刘宜

刘宜，西门子医疗CT临床应用专家；2020年加入西门子，有多年临床工作经验，掌握CT成像原理和各部位扫描技术。目前负责南区西门子CT全产品及影像后处理系统培训工作。

参考文献：

[1] Mileto A, Guimaraes LS, McCollough CH, Fletcher JG, Yu L. State of the Art in Abdominal CT: The Limits of Iterative Reconstruction Algorithms. Radiology. 2019 Dec;293(3):491-503. doi: 10.1148/radiol.2019191422. Epub 2019 Oct 29. PMID: 31660806.

[2] Solomon J, Marin D, Roy Choudhury K, Patel B, Samei E. Effect of Radiation Dose Reduction and Reconstruction Algorithm on Image Noise, Contrast, Resolution, and Detectability of Subtle Hypoattenuating Liver Lesions at Multidetector CT: Filtered Back Projection versus a Commercial Model-based Iterative Reconstruction Algorithm. Radiology. 2017 Sep;284(3):777-787. doi: 10.1148/radiol.2017161736. Epub 2017 Feb 7. PMID: 28170300; PMCID: PMC5702911.

[3] 曾更生.医学图像重建入门[M].北京:高等教育出版社,2010

[4] 余晓锷,高海英,蔡凡伟,赵全.基于噪声功率谱的CT图像噪声评价[J].中国医学影像技术,2014,30(08):1243-1246.DOI:10.13929/j.1003-3289.2014.08.035.

[5] [6] [7] 谭健,刘占.IRIS,SAFIRE,ADMIRE,你都知道吗？[J/OL].爱影联盟,2017-06-29

仅供专业人士交流目的，不用于商业用途。

2022年10月10日

如果你觉得写得还不错，请分享、在看和打赏！

原标题：《浅谈CT图像重建的数理基础》

阅读原文