混合高斯/泊松最大似然函数下的CBCT图像重建

郑蓉珍，赵 芳，李 波，田 昕*

(1.武汉大学 电子信息学院，湖北 武汉 430072;2.武汉大学 中南医院心血管内科，湖北 武汉 430071;3.武汉大学 口腔医院放射科，湖北 武汉 430079)

1 引 言

锥形束计算机断层成像技术(Cone Beam Computed Tomography，CBCT)因其X射线辐射剂量低、分辨率高、数据采集时间短等优点，近年来逐渐成为CT领域的研究热点。在低辐射剂量下，CBCT图像重建质量容易受到噪声等因素的干扰。加大辐射剂量可以提升图像质量，但同时也会增加患者患癌和遗传病变的风险[1]。所以，如何在低辐射剂量下实现高质量CBCT图像重建是一个非常有意义的研究问题。

CBCT三维重建算法包括解析算法和迭代算法。以Feldkamp-Davis-Kress (FDK)算法三维图像重建算法，一种易于实现的基于圆轨迹扫描的近似三维图像重建算法为代表的解析算法，重建速度快，在小锥角时具有良好的重建效果，但如果投影数据采样不充分或是测量噪声较大时，解析法重建出的图像效果不甚理想，且对噪声敏感。而各种迭代算法，例如同时型迭代重建算法(Simultaneous Algebraic Reconstruction Technique，SART)，可以恢复出高质量的图像，但其缺点同样是噪声会在迭代过程中被逐渐放大，严重影响重建图像质量。目前对CBCT图像去噪的主要思路是首先对投影数据去噪，然后通过去噪后的投影数据实现CBCT图像的三维重建。在去噪方法中，非局部均值算法具有良好的性能，但算法计算量大，去噪时间比较长[2]，因此，在CBCT去噪中应用较为受限。2010年，Yong Yin等对小波算法提出改进，提出了适合CBCT图像特征的去噪算法[3]，但该方法不具有自适应性。近年来，三维块匹配滤波(Block-Matching and 3D filtering, BM3D)算法具有较好的去噪性能，被广泛应用于CBCT的图像去噪中[4]。

综上所述，传统方法中对于低辐射剂量CBCT图像重建中的噪声去除是先去噪，再重建。但是，这种处理方式存在如下问题：(1)虽然可以去除投影数据中的噪声，但是同时也降低了投影数据中图像边缘的质量，从而使得由投影数据重建的三维图像往往会出现过平滑的现象；(2)在去噪过程中，未充分考虑低辐射剂量下CBCT投影数据中真实数据和噪声的分布特性。一方面，探测器在工作过程中会不可避免的产生服从高斯分布的加性系统噪声；另一方面由于辐射剂量的减少，所以在X射线穿过被测人体时会产生更严重的光电效应、光子散射等现象，在光电转换过程中光子数的减少会产生严重的泊松噪声[5]，因此，在低辐射剂量情况下CBCT系统投影数据噪声分布近似为复杂的高斯/泊松混合分布模型，而不是传统的高斯噪声模型。

本文将噪声去除与图像重建过程建立在统一的重建模型中，它包含一个基于混合高斯/泊松最大似然函数的保真项和一个基于三维全变分正则化方法的约束项。保真项用于描述混合高斯/泊松噪声环境下重建值与观测值尽可能的相似，约束项要求在重建过程中去除噪声的同时尽可能地保证重建图像的边缘信息。我们通过可分离近似法和扩展拉格朗日法对上述模型进行了求解，并通过仿真数据和真实数据验证了算法的有效性。

2 图像重建算法

2.1 CBCT三维图像重建原理及数学描述

典型的基于等距算法原理的CBCT三维图像重建过程描述如下：

如图1(a)所表示的笛卡尔坐标系xyz旋转β角度可得图1(b)tsz表示的坐标系图，其中p′和ξ分别表示经过待重建点P(t,s,z)映射到虚拟探测器上的横纵坐标，D表示射线源S到到虚拟探测器坐标原点O的距离。

图1 FDK扫描投影结构示意图Fig.1 FDK scan projection structure diagram

FDK算法的重建过程可以通过式(1)进行表示：

f(t,s,ξ)=

(1)

其中：U表示加权因子，Rβ(p,ξ)表示在角度β下经过被测物体的投影，具体过程可以参见文献[6]。

2.2 基于混合高斯泊松最大似然函数的CBCT三维图像重建算法

2.2.1 混合高斯/泊松噪声环境下的CBCT图像重建模型

基于以上推导本文所提出的CBCT重建模型可以表示为：

(2)

其中：x代表需要重建的三维CBCT图像；f(x)是最大似然函数的保真项，用于约束重建值与观测值的相似程度；g(x)是三维全变分正则化项，用于去除重建过程中的噪声；λ是正则化参数。

对于高斯噪声而言，假设其分布服从标准正态高斯分布而言，最大化似然函数f(x)可以得到：

(3)

其中：Ai代表服从式(1)的第i个角度下的投影矩阵，yi代表第i个角度的投影数据(观测值)，M代表总的投影数目。

(a)混合高斯/泊松最大似然函数的保真项

(4)

其中：[bi]k和[yi]k分别代表向量bi和yi中第k个数据。式(4)较为复杂，不宜直接求解其最大似然函数。采用GAT算法的类似思路[7]，对输入投影数据进行Anscombe变换，其过程如下：

(5)

在此情况下yi的概率密度函数可以表示为[8]：

p(yi|bi)=

(6)

已知bi=Aix，那么关于x的最大似然函数可以表示为[9]：

(7)

(b)三维全变分正则化项

g(x)是正则化项，考虑到重建过程中噪声的影响，这里可以引入三维全变分正则化项(3D Total Variation，3DTV)[10]来对重建图像进行约束，使得重建过程在去除噪声的同时能够较好地保留图像的边缘和细节信息。这里采用各项同性3D全变分正则化，此时g(x)可以表示为：

g(x)=‖x‖TV=

(8)

其中：βx,βy和βz是沿着不同方向的权重系数,Dx,Dy和Dz分别代表x,y和z方向的差分计算符号。假设x是三维重建图像x(x,y,z)对应的一维矢量(将三维图像按照首尾链接的方式变成一维矢量)，该转化过程用vec(·)表示，则Dxx，Dyx和Dzx的计算过程如下：

(9)

边界处的值可以采用补零的方式进行填充。

2.2.2 模型优化求解

为了描述方便，在下文中变量均转化为矢量进行描述。基于可分离近似的方法[11]，式(2)可以转化为如下两个迭代步骤：

(10)

其中：

(11)

f(x)=

(12)

(13)

其中：ρr是正则化参数，u是中间变量，y是拉格朗日乘子。与D类似，u和y可以表示为：

分别对x,u和y进行求解，从而可以将式(13)转化为求解如式(14)所示的三个子问题：

(14)

可以得到：

这里拟采用Barzilar-Borwein方法来更新公式(11)中的迭代步长a(k)[14]。

令g(k)为f(x)在第k次迭代时的梯度值，则：

(16)

其中：s(k-1)=x(k)-x(k-1),y(k-1)=g(k)-g(k-1)。

2.2.3 算法描述

总结前文步骤，所提出算法的流程如下：

算法:混合高斯/泊松最大似然函数下的CBCT图像重建输入:投影数据yi,投影矩阵Ai,循环次数TO和TI,正则化参数λ,其他参数ρr,βx,βy,βz,σ2变量初始赋值:x(1),a(1),y(1),u(1) 1.根据公式(5),计算yi;2.For k=1: TO3.根据公式(12),计算f(x);4.根据公式(11),计算z(k); For t=1: TI5.根据公式(15),计算x(t);6.根据公式(16),计算u(t);7.根据公式(17),计算y(t); End8.令x(k)=x(TI);9.根据公式(18),计算a(k); End输出:x

3 实验结果与分析

3.1 仿真数据实验

仿真实验数据为head phantom模体，模体大小为128×128×128，每个体素大小为4×4×4 mm3,投影规格为256×200×211，投影的角度为0°～210°，射线源到探测器的距离为1 500 mm，射线源到模体的距离为1 100 mm。通过调整仿真数据中平均入射光子数I0仿真不同剂量的X射线下得到的不含噪声的投影数据，对投影数据添加不同程度的混合高斯/泊松噪声进行后续实验，对其添加泊松噪声，泊松噪声的大小与I0大小成负相关关系，添加的高斯噪声通过调整期均方差σ调整其大小，σ越大添加的高斯噪声越大。

(a)主观结果

对模体在I0=50 000的辐射剂量情况下得到的投影数据添加混合高斯/泊松噪声后各算法重建head phantom模体取第45个切片，其中泊松噪声的大小与I0大小成负相关关系，高斯噪声：σ=50。重建图像中红色方框区域为目标放大区域，简称ROI区域。图2第1行分别表示FDK，SART，预先用BM3D对投影数据去噪然后用FDK算法重建和本文算法的重建结果，图2第2行分别表示各算法重建图像相同ROI区域的放大效果图。之所择图2(a)中标红的区域为目标放大区域，因为此区域包含了丰富的图像细节和边缘信息，放大之后便于区别各算法在保持重建图像边缘和细节信息方面的能力(彩图见期刊电子版)。对比其他算法的重建图像和本文算法的重建图像可以看本文算法在去除噪声的同时很好地保留了图像的细节和边缘信息。

图2 仿真数据主观实验结果.Fig.2 Subjective experimental result of simulation data

(b) 客观结果

表1，表2利用各种客观图像质量评价指标RMSE，CC，SSIM，UQI，PSNR及各算法的收敛次数ITE来对重建图像性能进行辅助评价。

表1 仿真数据客观实验结果(I0=100 000)

Tab.1 Objective experimental results of simulation data.(I0=100 000)

算法RMSECCUQIPSNRITEFDK0.067 80.9560.95471.501 41SART0.076 10.9470.93870.501 734BM3D+FDK0.068 00.9560.95371.476 71文中算法0.052 90.9730.97273.664 37

表2 仿真数据客观实验结果(I0=50 000)

Tab.2 Objective experimental result of simulation data.(I0=50 000)

算法RMSECCUQIPSNRITEFDK0.076 80.9430.94270.426 41SART0.084 30.9360.92169.612 425BM3D+FDK0.071 10.9520.94971.095 91文中算法0.061 70.9650.96172.322 16

表1，表2分别为在I0=100 000和I0=50 000的辐射剂量情况下的投影数据添加混合高斯/泊松噪声后各算法重建head phantom模体中取第45个切片的各项客观评价指标。其中两组数据所添加混合噪声中泊松噪声的大小与I0大小成负相关关系，高斯噪声相同：σ=50。对比各算法重建图像的各客观评价指标可以看出，不同辐射剂量不同强度混合高斯/泊松噪声情况下本文算法重建图像的同一切片误差最小，例如，相对于其他方法而言，PSNR最高可以提升2.1 dB；同时和原图的相关系数最高，结构相似性最高，UQI指标和信噪比也都是最高的。且图像质量随着噪声强度的增强退化越来越严重的情况下本文算法仍然能保持较高的重建图像质量。从计算复杂度而言，由于BM3D需要进行图像块匹配，因此，BM3D+FDK具有最高的计算复杂度；本文算法的计算复杂度次之，FDK具有最低的计算复杂度。从算法收敛性而言，对比迭代算法SART达到收敛时需要几十次的迭代次数，本文算法收敛的迭代次数均在10次以内，有很高的重建效率。虽然本文算法相对于经典的FDK算法而言，算法的计算复杂度有所提升，迭代次数有所增加，会对快速重建造成一定的影响；但是，快速重建问题可以通过基于GPU的并行加速算法来解决，因此并不会对算法的实际使用造成太大的影响。

3.2 真实数据实验

实验数据采用美国瓦利安(Varian)公司新一代250×200 mm2非晶硅平板探测器PaxScan 2520DX所获取的投影数据，其分辨率为768×960，共有360张投影数据。重建三维图像后取其不同的切片如图3所示。对比重建方法包括FDK，SART，BM3D+FDK及本文算法。对比可以看出，传统的FDK及SART算法在重建过程中容易受到噪声的干扰，且随着噪声的增强重建图像质量退化严重，因此，重建图像质量较差。BM3D+FDK算法对投影数据进行去噪，可以有效地提升重建图像的质量，但是同时也会使得重建图像边缘模糊。而本文算法在重建过程中有效地去除了噪声，并在很大程度上保留了图像的细节信息，具有较高重建图像质量。

图3 真实数据实验结果Fig.3 Reconstructed data of real projection data

4 结 论

本文针对低剂量下CBCT图像重建过程中的噪声有效去除问题，提出了一种基于高斯/泊松最大似然函数的CBCT三维图像重建算法。首先提出了一种结合图像去噪与三维重建的统一模型，包括一个保真项和一个约束项。在该模型中，考虑到实际存在的高斯/泊松混合噪声，保真项通过最大化高斯/泊松似然函数，来约束重建图像与观测数据尽可能的相似。由于重建图像是三维图像，因此，可以通过三维全变分作为约束项，来减少噪声对重建图像质量的影响。本文通过head phantom仿真数据和真实投影数据进行三维重建，并和其他的比较经典的重建算法FDK，SART和BM3D+FDK算法进行了主客观对比分析，结果证明本文算法具有较好的去噪效果，例如，相对于其他方法而言，PSNR最高可以提升2.1 dB，同时本文方法能在很大程度上保留重建图像的细节信息，收敛速度快，重建效率高。本文算法可以在低剂量下重建出较好质量的CBCT图像，可以被广泛应用到各种医学研究领域，例如牙科诊断。下一步的主要工作方向是研究基于GPU的并行加速算法，从而更进一步提升本文算法的运行速度。