一种去除CT图像中同心椭圆伪影的变分模型.docx

一种去除 CT 图像中同心椭圆伪影的变分模型

摘要：在计算机断层成像（CT）中，伪影会降低重建图像的质量。针对该问题，本文提出一 种去除 CT 图像中同心椭圆伪影的方法。该方法基于方向全变分（DTV）的思想，将椭圆伪影 的去除问题建模为能量最小化问题，并通过对椭圆伪影边缘特征的分析，建立适应于椭圆伪 影的变分模型。由于提出的模型为具有可分裂结构的非光滑凸优化问题，因此本文使用交替 方向乘子法（ADMM）。最后，通过仿真实验验证该模型对于去除椭圆伪影的有效性。
关键词：CT 影像；交替方向乘子法；椭圆伪影；非光滑问题

自从 1971 年 Hounsfield 发明了世界上首台 X 射线计算机断层成像扫描（computed tomography，CT）机，CT 技术在扫描方式和重建方法上经历了 5 代发展[1] ，在医疗和工业 检测等方面产生了深刻影响。由于系统不完善，CT 成像所得到的图像中可能存在多种伪影，
对临床诊断与故障诊断存在负面影响[2-3]。
随着 CT 技术的深入应用，以及快速成像、定量成像、低剂量成像[4]的新需求，非标准 轨迹扫描及分布式多光源成像在解决伪影问题上体现出了巨大的潜力[5] ，但也伴随着新型 伪影等其他问题的出现。伪影严重影响临床诊断的效果，尤其是伪影可能被认为是医学图像 （例如肿瘤）的组成部分，这可能会导致错误的诊断[6] 。例如在放射治疗中，严重的伪影 可能被识别为骨骼，导致定位错误。正是由于这些伪影使得后处理和定量评估变得复杂， 因此开发出高效且可靠的伪影去除方法至关重要。
由于探测器损坏或探测器响应不一致，标准圆轨迹扫描 CT 图像中环形伪影是最常见的 伪影之一[7] 。近年来，许多研究人员一直在开发用于去除 CT 中环形伪影的方法[8] ，这些方 法大致可以分为投影域的前处理[9]和图像域的后处理[10]两类。投影域的前处理需要访问原 始投影数据，这限制了算法的适用性。而针对图像域的后处理技术，有研究者提出了一种 基于形态算子的环形伪影校正方法[11] ，该方法中参数的选取对伪影校正的精度至关重要。 随后，研究人员提出了一种使用独立分量分析（ICA）的有效方法[12] ，但该方法并不能有效 地保持图像细节。
在假设环形伪影中的灰度值特征是局部极值的前提下，研究人员使用中值滤波来去除

极坐标下的条纹伪影[13]；最近，研究人员提出了一种稀疏约束下基于变分的环状伪影去除 方法[14] ，该方法具有有效地校正效果，但坐标变换的过程降低了图像的空间分辨率；很多 研究者利用环状伪影在极坐标系中表现为条状伪影[15-17] ，在极坐标系中进行图像处理后， 图像被变换回笛卡尔坐标系，进而提出新型的环状伪影去除方法。
在成像过程中，对于非标准扫描轨迹，如变螺距、变半径螺旋轨迹，则可能会出现椭 圆伪影，而上述去除伪影的方法并不能有效地应用到去除该类伪影，并且针对去该类伪影 的研究工作相对缺乏。针对上述问题，在结合椭圆伪影几何特性和利用加权型全变分能有 效保持图像边缘特性的基础上，本文提出一种去除椭圆伪影的新方法。由于提出的模型是 含有线性算子的非光滑凸优化问题，因此可以利用经典的交替方法乘子法（alternating direction mutipler method，ADMM）将初始问题转化为多个易求解的子问题，进而可以高 效求解，并能保证算法的收敛性。最终数值实验表明，本文提出的模型和算法在不影响图 像细节的情况下，能有效地去除 CT 图像中的椭圆伪影。
1 模型建立与求解
假定 g 为观测到的具有椭圆伪影的图像，为了去除图像中的椭圆伪影，本文提出下述

各向异性全变分模型：





（1）

其中拟合项 Ⅱf - g Ⅱ2(2)是为了保证复原图像f 不偏离观测图像f太远，正则项 ⅡT▽fⅡ2; 1 （注：
Ⅱ ·Ⅱ2; 1 表示内范数为 l2 范数，外范数为 l1 范数）是为了有效刻画图像的细节特征，参数 λ 用 来平衡这两者关系的权重。另外，这里 ▽f = (▽xf ; ▽yf )T 表示梯度向量，T = diag(t1; t2 ) = ΛτR-θ 用来刻画椭圆伪影和图像边缘等结构特征，其中

τ 为扩散参数。
模型（1）是非光滑凸优化问题，显然目标函数也是强制的，因此存在唯一的最优解。 虽然针对非光滑问题可以有各种光滑化方法求解，但仅能得到近似解。由于其非光滑项具 有可分性结构，因此本文利用交替方向法求解，即：将初始问题转化为多个易求解的子问 题。为此，通过引入两个辅助变量 w 和 v ，可将问题（1）转化为约束优化问题：
（3）
为了有效求解问题（3），在增广拉格朗日方法框架下，可将其转化为鞍点问题：
。 （4）

下面在交替方向乘子法框架下，即固定其中 4 个变量求解剩余的 1 个变量，求解问题

（4）即有下述高斯赛德尔迭代格式：

下面考虑如何解决（5）中的子问题。
问题（5）的第 1 个子问题是经典的 l2 -l1 问题，整理后可得：

通过使用软阈值算子，其闭合解可以写成：

问题（5）的第 2 个子问题是一个光滑优化问题，整理后可得：

其最优解 vk+1 满足以下线性方程组：

通过简单的计算，可得：
（10）
问题（5）的第 3 个子问题是一个最小二乘优化问题，整理后为：

该问题是凸优化问题，其对应的欧拉-拉格朗日方程为：
(λI -Y1Δ)fk+1 = λg -div (αk )-Y1 div (vk+1) ， （12） 其中 I 为单位算子，对于线性方程（12），不同的边界条件对应着不同的数值方法。当使用 诺依曼边界条件或狄利克雷边界条件时，拉普拉斯算子-Δ 是半正定的。在这种情况下，可
以使用预条件共轭梯度（PCG）法来求解它，因为 λ 和 Y1 都是正标量，上面方程组左侧的 矩阵是对称且正定的。本文假设边界条件是周期性的，故右侧算子是块循环矩阵，方程

（12）可以利用快速傅里叶变换和逆变换快速求解：

其中 F (·)和 F−1 (·) 分别是快速傅立叶变换（FFT）及其逆变换。 综上，求解模型（1）的算法框架如表 1：
表 1 算法 1
Table 1 Algorithm 1

算法 1：使用ADMM 解决问题（4）
1. 初始化： λ > 0; Y1 > 0; Y2 > 0 ，选择 f 0 ; v0 ; α0 和 β0 的初始值
2. 对于 k = 1; 2; · · · 通过（7），（10）和（13）依次获得 wk+1 ; vk+1 ; fk+1
3. 直到满足 ≤ 10−5 或 k ≥ 500 时终止
4. 将 f := fk+1 作为输出图像
针对算法 1，有如下收敛性结果，其证明类似于文献 [18]。
定理： 由算法 1 产生的序列{(wk ; vk ; fk ; αk ; βk )}收敛到问题（4）的鞍点，且序列{fk } 收敛到问题（1）的解。
2 实验与讨论
实验使用 Shepp-Logan 模体数据对所提出的方法进行评估，展示所提出的去除椭圆伪 影模型及算法的有效性。本文有两类参数需要选取，即模型参数和算法参数，根据多次实 验可知，正则化参数 λ 为 0.01，权重参数 τ 取 2，惩罚参数 Y1 和 Y2 均取 3，旋转矩阵参数 θ 选取如下：
根据 Bayram 和等的观点[19] ，旋转矩阵中参数 θ 表示沿该倾斜角度的像素变化具有更
高的权重。如图 1，根据椭圆几何关系可以 导出：
椭圆长轴在y轴上(椭圆长轴在x轴上) ，（14）
（a） （b）
图 1 椭圆几何关系示意图，其中子图（a）中
椭圆长轴在 x 轴上，子图（b）中椭圆长 轴在y 轴上
Fig.1 Schematic diagram of the geometric
relationship of the ellipse, the long axis of the ellipse in subfigure (a) is on the x-axis and the long axis of the ellipse in subfigure (b) is on the y-axis
其中k 是同心椭圆伪影的长半轴与短半轴之 比， (x;y)是 CT 图像中像素点的坐标。在本 文实验中，图 2 中椭圆伪影 k 取2，图 4 中 椭圆伪影 k 取5/3。
图 2 显示了包含椭圆伪影的原始和去除 伪影后的 Shepp-Logan 模体。第 1 行中的图 像是原始图像，第 2 行中的图像是带有椭圆 伪影的图像，第 3 行中的图像是使用本文模

型去除了椭圆伪影的图像。（a）、（g）、（m）是整个图像，（b）、（h）、（n）中红框圈出的 是所选的第 1 组 ROI，蓝框圈出的是所选的第 2 组 ROI，（c）和（d）、（i）和（j）、（o） 和（p）分别是对应红色实线框中的放大视图，（e）和（f）、（k）和（l）、（q）和（r）分 别是对应蓝色实线框中的放大视图 。图 2（ g）中的椭圆伪影很严重，而校正后的图像 （图 2（m）） 几乎与参考图像（图 2（a）） 相同。显然，在经过本文模型处理之后（图 2 第 3 行中），图像中几乎看不到伪影。

（a）

（g）

（m）


（b）

（h）

（n）


（c）

（i）

（o）


（d）

（j）

（p）


（e）

（k）

（q）


（f）

（l）

（r）

图 2 去伪影前后 Shepp-Logan 模体及其 ROI，第 1 行为参考图像，第 2 行为被伪影污染的图像，
第 3 行为去伪影后的图像
Fig.2 Shepp-Logan phantom and its ROI before and after artifact removal, the first
row is the reference image, the second row is the image contaminated with artifacts, and the third row is the artifact-removed image
本文实验通过计算 PSNR 和 SSIM 评估去伪影后图像的质量，在表 2 中显示这些数据。 从表格中可以清楚地看到，通过本文模型和算法去除椭圆伪影后，SSIM 和 PSNR 在整体图 像和ROI 区域都有所增加，图像质量在去噪后明显提升。
表 2 去伪影前后 Shepp-Logan 模体的 SSIM 和 PSNR 对比
Table 2 Comparison of SSIM and PSNR of Shepp-Logan phantom before and
after artifact removal

对应图像
Shepp-Logan 模体
SSIM（含伪影）
SSIM（不含伪影）
PSNR（含伪影）
PSNR（不含伪影）
（g）、
（m）
0.990
0.997
49.052
51.442
（i）、
（o）
0.961
0.984
44.178
48.122
（j）、
（p）
0.981
0.992
46.369
49.611
（k）、
（q）
0.989
0.997
46.150
47.742
（l）、
（r）
0.963
0.988
43.389
47.108

0.220

灰度值
0.215
0.21
0.20
0.20
0.19
0.19
0.18
原始图像
含伪影图像
去伪影后图像

350 360 370 380 390 400 410 420 像素值
图 3 像素对比图，其中红色折线表示原始图像灰度值，蓝色折线表示含伪影图像灰度值，
绿色折线表示去伪影后图像灰度值
Fig.3 Pixel comparison chart, the red line indicates the gray value of the
reference image, blue line indicates the gray value of the image with artifacts, and green line indicates the gray value of the image after removing artifacts

（a）

（g）

（m）

（b）

（h）

（n）

（c）

（i）

（o）

（d）

（j）

（p）

（e）

（k）

（q）

（f）

（l）

（r）

图 4 去伪影前后胸腔 CT 及其 ROI，第 1 行为原始图像，第 2 行为被伪影污染的图像，第 3 行
为去伪影后的图像
Fig.4 Chest CT and its ROI before and after artifact removal. The first row is the
reference image, the second row is the image contaminated with artifacts, and the third row is the artifact-free image

实验中还取出上述图像的一行显著变化的像素点（第 204 行中第 350 到第 420 个像素 点）来直观地展示伪影去除效果（图 3）。图中通过 3 种不同颜色的实线来显示原始图像 的灰度值变化，其中红色实线表示原始图像该行的灰度值，蓝色虚线表示含伪影图像该行 的灰度值，绿色虚线表示去除伪影后图像该行的灰度值。从图 3 中可以明显观察到，椭 圆伪影导致含伪影图像灰度值显著偏离原始图像，而本文提出的模型和算法有效修复了这 种偏差。
为进一步测试模型与算法的鲁棒性，本文使用图像结构更复杂的胸腔 CT 进行仿真实验。 图 4 显示胸腔 CT 图像去除椭圆伪影的结果（ROI相应的图像位置与图 2 相同）。从图 4 中 可见，即使在图像结构更复杂的 CT 图像中，大部分伪影都能得到去除，且图像结构仍能被 较好保存。
表 3 展示该胸腔 CT 图像去伪影前后 SSIM 和 PSNR 的变化，可见经过本文模型处理后， CT 图像结构相似度和峰值信噪比普遍提高，同心椭圆伪影得到有效去除。
表 3 去伪影前后胸腔 CT 的 SSIM 和PSNR 对比
Table 3 Comparison of SSIM and PSNR of chest CT before
and after artifact removal

对应图像
Shepp-Logan 模体
SSIM（含伪影）
SSIM（不含伪影）
PSNR（含伪影）
PSNR（不含伪影）
（g）、
（m）
0.619
0.655
32.472
34.238
（i）、
（o）
0.917
0.940
32.984
33.371
（j）、
（p）
0.833
0.828
34.254
35.229
（k）、
（q）
0.933
0.954
33.328
34.801
（l）、
（r）
0.859
0.874
34.766
34.156

3 结论
本文利用椭圆伪影的结构特征，结合方向全变分可以有效描述该结构特征的性质，建 立一个新型的去椭圆伪影的模型，并利用算子分裂型方法求解，数值实验验证了模型和方 法的有效性。
然而，CT 图像中通常可能存在不止一种形态的伪影，本文只针对同心椭圆伪影进行 去除，在接下来地研究中，将建立起适用于多种伪影的模型，以求在临床诊断上取得更 高的性能。
参考文献
[1] WANG G, YU H, de MAN B. An outlook on X-ray CT research and development[J]. Medical Physics, 2008, 35(3): 1051-1064.
[2] 汤少杰, 张砚博, 汪友明, 等. 可同时抑制多种图像伪影的最优骨校正[J]. CT 理论与应用研究, 2018，27(3): 301-314. DOI:10.15953/j.1004-4140.2018.27.03.03.
TANG S J, ZHANG Y B, WANG Y M, et al. An optimal bone correction capable of simultaneously suppressing various types of image artifacts[J]. CT Theory and Applications, 2018, 27(3): 301-314. DOI:10.15953/j.1004-4140.2018.27.03.03. (in Chinese).