赵 霞，赵金龙，赵荣格，陈 燕，桂志国，刘 祎+

(1.中北大学 山西省生物医学成像与影像大数据重点实验室，山西 太原 030051； 2.中北大学 信息与通信工程学院，山西 太原 030051)

0 引 言

计算机断层扫描(CT)已被广泛应用于各种临床应用。虽然CT[1]检查得到了广泛的认可，但是当过多的X射线照射剂量停留在人体内时，会给人体健康带来极大的隐患。因此CT检查应遵循辐射防护最优化的ALARA(as low as reasonably achievable)原则，尽可能地保护受检者。然而，X射线剂量的降低会使投影数据中产生大量的噪声，不能得到较为满意的重建图像[2]。因此，用最小的代价获得最佳CT重建图像，有着重要的实际应用价值。

图像域统计迭代重建算法是低剂量CT成像的一种方法，该算法是将所获得的投影数据进行统计建模，然后用某种迭代方法对目标函数进行求解，其最优解为待重建图像[3-12]。惩罚加权最小二乘算法[13-17]是最常见的一种统计迭代重建算法，在各种惩罚模型中，MRF是一种广泛使用的模型，但MRF没有明确地考虑非均匀性，不是一个令人满意的模型，因此，Zhang等对整个图像使用高斯混合MRF来考虑分布复杂度[18]。TV先验模型也是一种常用的模型，尤其是在投影数据缺失的情况下，He等基于加权方差的加权因子与TV模型相结合提出自适应加权全变分(AWTV)模型[19]，Zhang等先把梯度保真约束项和能够区分图像平滑区和细节区的边缘指示函数应用到TV模型中得到基于梯度保真项的自适应全变分模型[20]。Zhang等提出了基于绝对差值排序(ROAD)和中值先验(MP)模型的TV低剂量CT重建算法[21]。

TV正则化方法在图像去噪的同时能较好地保持图像的边缘，但是容易使图像的平坦区域产生阶梯效应。高斯MRF正则化方法使用固定的平滑系数导致图像边缘较为模糊。根据TV正则化和高斯MRF正则化的优缺点，本文提出了一种基于分区域处理的联合先验低剂量CT统计迭代重建算法。该算法首先对低剂量CT图像进行中值滤波，再利用滤波后图像的梯度值划分出图像的边缘区域和平坦区域，然后将平坦区域中的图像块提取出来用高斯MRF正则项处理，边缘区域中的图像块提取出来用TV正则项处理，最后将这两种正则项作为联合先验应用到惩罚加权最小二乘法中，使用SOR迭代重建算法进行求解。

1 噪声模型

在低剂量CT重建中，投影数据和噪声模型是至关重要的一个环节。以往的研究揭示了CT投影数据噪声产生的两个主要来源，即X射线量子噪声和系统电子噪声。探测器检测获得的对数变化前的CT投影数据可以用统计独立的泊松分布加上统计独立的高斯或正态分布来描述

(1)

文献[17]中指出，当投影数据中有电子噪声时，噪声方差可以定义为

(2)

2 基于分区域处理的低剂量CT统计迭代重建算法

2.1 惩罚加权最小二乘算法

参见文献[19]中的加权最小二乘(weighted least square，WLS)算法，其在图像域的代价函数为

Φ(u)=(y-Gu)′∑-1(y-Gu)

(3)

Φ(u)=(y-Gu)′∑-1(y-Gu)+βR(u)

(4)

式中：第一项为式(3)中的WLS，第二项R(u)表示惩罚项，β是平滑参数。一般的惩罚项有高斯MRF模型、TV模型等。

2.2 改进算法

在基于高斯MRF的惩罚加权最小二乘法中，由于目标函数中正则项的权重是固定的，没有考虑平坦区域和边缘区域具有不同的特点，使得图像边缘较为模糊。TV正则化方法虽然在去噪的同时能较好地保持图像的边缘，但是容易使图像的平坦区域产生阶梯效应。因此，本文提出了一种基于分区域处理的联合先验低剂量CT统计迭代重建算法，该算法首先对低剂量CT图像进行区域划分，再对不同区域的图像块进行不同的处理，然后将改进的正则项应用到惩罚加权最小二乘法中，最后使用SOR迭代重建算法进行求解。

2.2.1 区域划分并分区域处理

本文算法首先对低剂量CT图像进行平坦区域和边缘区域的划分，由于低剂量CT图像信噪比低，进行准确划分比较困难，所以先对其进行中值滤波，将滤波后的图像进行梯度变换，然后根据图像的像素值确定一个阈值，将大于阈值的像素都等于255，小于等于阈值的像素都等于0，即可划分出图像的平坦区域和边缘区域。

对于二维CT图像ui,j来说，梯度变换公式为

(5)

图1 Shepp-Logan模型区域划分

图1(a)是Shepp-Logan模型，图1(b)是对Shepp-Logan模型进行FBP算法重建后的图像，图1(c)是划分区域后的图像。

根据划分区域后的图像，将低剂量CT图像中所对应的边缘区域的图像块提取出来进行TV正则化处理，将平坦区域的图像块提取出来进行高斯MRF正则化处理。高斯MRF正则项的表达式如下

(6)

对于二维CT图像ui,j来说，它的全变分表达式为

(7)

将图像值非负作为约束条件，且可通过如下最优化问题来得到图像u

(8)

采用梯度下降法来求解该最优化问题，其中全变分的计算公式为

(9)

式中:ε是一个非常小的正参数，防止分母为0。

2.2.2 基于分区域处理的惩罚加权最小二乘算法

本文算法在图像域的目标函数为

Φ(u)=(y-Gu)′∑-1(y-Gu)+β1R(u1)+β2TV(u2)

(10)

在重建实验中，两个惩罚系数β1和β2是可以独立设置的，首先固定β1，并搜索β2使得RMSE最小，然后固定β2再搜索最优的β1来确定β1和β2的取值，u1和u2分别是从图像平坦区域和边缘区域提取的图像块，R(u1)表示根据式(6)对平坦区域的图像块进行高斯MRF正则化处理，TV(u2)表示根据式(9)对边缘区域的图像块进行TV正则化处理。

(11)

(12)

其具体实现步骤如下：

步骤1 使用FBP算法对低剂量CT投影数据进行重建后的图像作为初始化图像，记为u0。

步骤2 将u0先进行中值滤波，然后利用滤波后图像的梯度值将其分为平坦区域和边缘区域。

步骤3 分别将u0所对应的平坦区域和边缘区域的图像块提取出来。

步骤4 根据式(6)对平坦区域的图像块进行高斯MRF正则化处理，根据式(9)对边缘区域的图像块进行TV正则化处理。

步骤5 利用式(12)计算出迭代重建图像u。

重复步骤2～步骤5，经过一定次数目标函数的解达到收敛后的重建结果作为最终的重建图像。

3 实验结果与分析

为了验证本文所提算法的有效性和可行性，本文用模拟数据和临床数据来进行实验，然后将本文算法与FBP算法、基于TV先验的惩罚加权最小二乘法(PWLS-TV)以及基于高斯MRF先验的惩罚加权最小二乘法(PWLS-MRF)的重建图像进行了比较。为了进一步评价重建图像的质量，本文采用相关系数、信噪比、归一化均方误差对重建图像进行定量描述分析，定义如下：

(1)相关系数(correlation coefficient，CORR)

(13)

(2)信噪比(signal to noise ratio，SNR)

(14)

信噪比越高表明重建图像质量越好。

(3)归一化均方误差(root mean squared error，RMSE)

(15)

RMSE越小表明重建结果越接近原始图像。

3.1 模拟数据实验

图2 Shepp-Logan模型在不同投影角度下 不同重建算法结果

图2中第1列为FBP算法结果，第2列为PWLS-MRF算法结果，第3列为PWLS-TV算法结果，第4列为本文算法结果。为了更加清楚看清本文算法重建的图像比其它重建算法重建的结果更好，选取图2的两个感兴趣区(region of interest，ROI)，即ROI1和ROI2进行放大，结果如图3所示。

由图2、图3可知，FBP重建结果中出现条形伪影，且角度越少越严重；采用PWLS-MRF算法进行重建，虽然重建图像中的条形伪影得到了有效的抑制，但没有很好保护图像的边缘信息，其边缘模糊(如图3箭头所指边缘)；PWLS-TV算法虽然保持了图像的边缘信息，但同时在平坦区域带来了阶梯效应(如图3框中区域)，而本文算法在去除阶梯伪影的同时保护了图像的边缘信息。

表1为各算法在不同角度下重建出的Shepp-Logan图像相较于原始图像的质量评价参数对比。从表1可知，本文算法和其它两种重建算法相比，CORR、SNR值更大，RMSE更小。表明本文算法重建出的图像质量更好，因此，

表1 Shepp-Logan模型在不同角度下 各算法的客观评价参数

本文算法的改进是可行的。

3.2 临床数据实验

由于无法获得实际的投影数据，所以本文实验二使用了从The Cancer Image Archive(TCIA)中获得的一张高剂量的实际肺部图像切片(如图4所示)进行实验的仿真模拟来验证本文所提算法的有效性。

图4 高剂量肺部切片

实验二采用仿真平行束扫描高剂量肺部切片来获取投影数据，稀疏角度的个数分别为80、120、150，每个角度上有888个探测器，实验二获取原始投影数据和给投影数据中加噪声来模拟低剂量CT稀疏投影数据的实验参数设置同实验一，然后对仿真的低剂量CT稀疏投影数据进行实验，实验结果如图5和图6所示。

图5 肺部切片在不同投影角度下不同重建算法结果

图6 肺部切片在不同角度下各重建图像局部放大图

图6是将图5中的ROI1和ROI2进行放大。从图5、图6可以看出，稀疏角度数越多重建结果越好，由于PWLS-MRF算法中权重是固定的，使得图像边缘较为模糊(如图6箭头所指边缘)。PWLS-TV算法重建效果明显改善，但是在欠采样的情况下重建图像产生了阶梯伪影(如图6框中区域)。而本文算法在平坦区域平滑的同时保持了图像的边缘信息，投影角度数越多越接近原始图像。

表2为肺部切片采用PWLS-MRF、PWLS-TV和本文算法重建的图像相对于模板图像的RMSE、CORR、SNR等质量评价参数的对比。由表2可知，对比其它两种重建算法，本文算法的CORR、SNR值更大，RMSE更小。表明本文算法重建出的图像质量更好，而且去噪能力更强。

表2 肺部切片在不同角度下各算法的客观评价参数

本文实验一和实验二使用向模板和高剂量CT图像的投影数据加噪声的方式来模拟低剂量CT投影数据，本文实验三直接使用从The Mayo Clinic(USA)获取的一张低剂量的实际腹部图像切片(如图7所示)，采用仿真平行束对其扫描来获取投影数据，稀疏角度的个数分别为80、120、150，每个角度上有888个探测器，然后对仿真的低剂量CT稀疏投影数据进行实验，实验结果如图8和图9所示。

图7 低剂量腹部切片

图8 腹部切片在不同投影角度下不同重建算法结果

图9 腹部切片在不同角度下各重建图像局部放大图

图9是选取图8的3个ROI，即ROI1、ROI2、ROI3进行放大。由图8、图9可以看出，第一列的FBP算法投影角度数越少，重建结果中噪声越多，细节信息失去的越多。第二列的PWLS-MRF算法重建的图像边缘信息被过度平滑(如图9箭头所指边缘)。第三列的PWLS-TV算法虽然较好地保持了图像的边缘，但是在噪声多的地方产生阶梯伪影(如图9框中区域)。第四列的本文算法消除了PWLS-TV引入的阶梯伪影，同时保持了图像的边缘细节信息。从视觉效果上来看，本文算法明显优于其它3种重建算法。

表3为不同算法对投影角度数为80、120、150的噪声投影数据进行重建结果的客观评价参数，在条件相同的情况下，与其它两种算法相比，本文算法重建图像的SNR值和CORR值更大，RMSE值更小，表明本文算法的去噪能力更强且与原始图像更相似。

4 结束语

本文提出了一种基于分区域处理的低剂量CT统计迭代重建算法，该算法能够充分地利用高斯MRF和TV各自的优点对平坦区域和边缘区域进行处理，使得重建出的图像在平滑平坦区域的同时保留图像的边缘信息。将本文算法与其它重建算法比较可以发现，本文算法不仅去噪能力强，而且还能有效地保护重建图像的边缘细节信息。但是，本文所提的算法中有一些参数是凭经验确定的，可能会影响重建图像的质量。因此，可以将自适应参数选择作为未来的研究方向。

表3 腹部切片在不同角度下各算法的客观评价参数