光栅相衬CT中吸收信号环形伪影的去除方法

杜天宇, 高昆, 吴朝

(中国科学技术大学国家同步辐射实验室, 安徽 合肥 230029)

0 引 言

在传统的X 射线吸收成像中,X 射线穿过不同吸收系数的物质,探测器上会记录不同强度的信号。当物质具有相近的吸收系数时,传统的吸收成像不能得到令人满意的结果。光栅X 射线相衬成像是一类新兴的多模态成像技术,在数据采集过程中可以同时获取吸收、折射率和散射系数三种信号。样品折射率的差异可以显示某些吸收衬度无法看到的结构,有望在临床医学、生物学、材料学等领域得到应用[1−3]。然而,研究人员对于相衬计算机断层扫描(CT)设备的伪影去除算法却鲜有讨论。

在过去几十年里,研究人员提出了很多针对传统吸收成像的伪影去除算法。1998 年,Raven[4]利用正弦图灰度值在垂直于角度方向的剧烈变化,使用低通滤波器去除伪影。2004 年,Sijbers 和Postnov[5]提出了一种两步算法,通过形态学操作提取感兴趣区域,在极坐标中通过滑动窗口生成伪影模板,并用原图像减去伪影模板获得校正后的图像。2009 年,Prell 等[6]提出了两种基于滤波的伪影去除算法,并且证明了在极坐标系下的伪影去除算法比在笛卡尔坐标系下的更加高效。同年,M¨unch 等[7]提出了一种基于小波分解和傅里叶滤波的去噪滤波器,可以在保留原始图像信息的同时去除伪影。2010 年,Abu Anas 等[8]阐述了环形伪影在正弦图域中表现为条纹伪影,去除条纹伪影比环形伪影更加简单。2011 年,Wang 等[9]提出一种改进Canny 算法用于去除环形伪影,该算法对正弦图中的条纹伪影进行检测,并采用分段B 样条拟合伪影数据。2017 年,Liang 等[10]采用了一种基于迭代的伪影去除算法,通过迭代残差图像提取细节。2018 年,Vo 等[11]分析了不同形貌伪影的产生原因,并针对不同伪影提出了相应的解决方案。

在基于光栅的相衬成像装置中,光栅结构缺陷会影响某些探测器通道中的信号,使正弦图中引入明显的亮条纹或暗条纹,在重建图中会出现严重的环形伪影[12]。根据计算机断层扫描的特点,靠近探测器中心通道的伪影具有更高的强度[13]。因此,对于中心通道周围的样本区域,受光栅缺陷影响的伪影较强,通过现有算法难以去除。本文提出一种针对光栅相衬CT 设备中吸收信号伪影的去除算法,并通过评价指标证明了该算法在强伪影条件下有很好的表现。

1 实验装置

图1 为典型的三光栅相衬X 射线成像装置的原理图,该装置由三个光栅组成:一个源光栅G0,一个分束光栅G1 和一个分析光栅G2。入射光束穿过G0,在G1 处分裂成小束,在G2 处产生条纹,由G2 解析为莫尔条纹图样,被探测器记录[14−16]。通过横向步进G1(相位步进),可以使后方摩尔条纹产生变化,并通过重复这个过程提取出每个像素对应的吸收、折射率和散射系数信息[17−20]。光栅相衬CT 成像设备结合了计算机层析技术,在旋转过程中的每个投影角度重复相位步进过程,然后分别对三种数据进行重构[21]。

为了验证所提出的伪影去除算法,根据图1 构建了一个光栅相衬成像装置,使用60 kV 管电压、200 mA 管电流的常规X 射线管作为光源。系统中所用三块金材料吸收光栅G0、G1、G2 的周期分别为30、24、120µm,在光栅均匀区域的位移曲线可见度平均可达22%。G0 到G1 的距离为300 mm,G1 到G2 的距离为1200 mm,CT 转轴位于G1 下游140 mm 处。成像数据的采集基于“周进步采集”方案,即光学系统连续旋转一周采集投影数据后再步进一步,然后进行第二步的旋转采集,如此重复直至所有步进数据采集完毕,这样可以大大提升采集速度。CT 成像实验过程中每周拍摄800 个投影角并均匀分布于360◦上,单周曝光时间为6 s。实验中,使用样品为聚苯乙烯(PS)圆柱体,浸泡在质量百分比为1%的氯化钠水溶液中。图2 为样品的成像原始数据,图2(b)所示的重建图说明了样品的层析截面,其中较亮的月牙形区域是氯化钠水溶液,较暗的圆盘区域是PS 圆柱体。

图1 三光栅相衬成像装置示意图Fig.1 Illustration of a three-grating phase contrast imaging device

图2 原始样品数据。(a)原始正弦图;(b)原始重建图Fig.2 Original sample data. (a)Original sinogram;(b)Original reconstruction image

2 方法与结果

根据图像重建算法,吸收信号是通过斜坡滤波核重建的,它会放大高频噪声。在强伪影的情况下,重建之前应该先进行正弦图域的处理。因此将算法分成两部分,分别是正弦图域的处理和重建图域的处理。由于样本的吸收系数在空间上连续变化,期望相邻像素服从相似的分布,使用如下正弦图域的算法[11]:1)将正弦图的灰度值沿着角度方向由小到大进行排序,并在排序前记录每个像素的位置;2)对排序后的正弦图沿着探测器通道方向进行平滑滤波;3)根据步骤1)中记录的像素位置,将平滑滤波后的结果重新排列到原始正弦图上。

将图2(a)中的正弦图沿着角度方向排序,并绘制投影角度为300◦时的轮廓图,如图3(a)、(b)所示。对排序后的图像沿着探测器通道方向进行平滑滤波,选择相同的投影角度绘制轮廓图,如图3(c)、(d)所示,其中平滑滤波选择大小为5 的均匀滤波器,如果尺寸较大会导致正弦图中的信息模糊,尺寸较小会导致平滑效果受到影响。最后将滤波后的图像重新排列到原始正弦图上并进行重建,可得到如图4(a)、(b)所示的正弦图和重建图。

图3 正弦图域的算法步骤。(a)排序后的正弦图;(b)投影角度为300◦时图3(a)的部分轮廓图;(c)滤波后的正弦图;(d)投影角度为300◦时图3(c)的部分轮廓图Fig.3 Algorithm steps of sinogram-domain.(a)Sorted sinogram;(b)The partial profile of Fig.3(a)at 300◦projection angle;(c)Filtered sinogram;(d)The partial profile of Fig.3(c)at 300◦projection angle

对比图2 和图4 可以看出,正弦图域算法去除了原始重建图中距离中心通道较远的部分弱伪影。对于距离中心通道较近的强伪影,采用如下重建图域的算法进行去除:1)将重建图从笛卡尔坐标系转换到极坐标系,并沿径向使用平滑滤波;2)计算滤波前后的差值,并将结果转换回笛卡尔坐标系;3)利用图像分割算法对正弦图域处理后的重建图进行分割,得到每一类样品的分布;4)通过形态学腐蚀操作获得每一类样品的内部区域,保护样品的边界像素信息;5)根据残差图像中样品区域的灰度分布特点定位伪影像素,并用临近的非伪影像素均值替代。

图4 正弦图域算法的处理结果。(a)处理后的正弦图;(b)处理后的重建图Fig.4 Results of sinogram-domain algorithm. (a)Sinogram after processing;(b)Reconstruction image after processing

环形伪影是以探测器中心通道为圆心的同心圆环,将重建图从笛卡尔坐标系转换到极坐标系,环形伪影表现为垂直于极径方向的条纹伪影,如图5(a)所示。因为样品的边界和伪影在极径方向都具有较强的梯度,沿径向平滑滤波,并计算滤波前后的差值,得到如图5(b)所示的残差图像,其中平滑滤波和正弦图域算法中的平滑滤波采用相同的策略。将残差图像转换到笛卡尔坐标系中,可以清晰地看到样品边界和伪影的轮廓,如图6(a)所示。

图5 极坐标系中的图像。(a)重建图像;(b)残差图像Fig.5 Image in polar coordinate system. (a)The reconstruction image;(b)The residual image

图6 重建图域算法的处理结果。(a)残差图像;(b)每一类样品的内部区域;(c)残差图像中的PS 区域;(d)PS 区域的灰度直方图;(e)PS 区域的伪影像素;(f)重建图域算法的处理结果Fig.6 The result of reconstruction image-domain algorithm. (a)The residual image;(b)The internal region of each kind of sample;(c)PS region of the residual image;(d)Grayscale histogram of PS region;(e)Artifact pixels of PS region;(f)The result of reconstruction image-domain algorithm

采用基于高斯混合模型(GMM)的机器学习算法,以灰度值和坐标作为特征,可以得到每一类样品的区域分布[22]。因为样品边界是图像中至关重要的视觉信息,为了保留样品的边界,通过形态学腐蚀操作得到样品的内部非边界像素区域,如图6(b)为PS 圆柱体和氯化钠水溶液的内部像素分布。对于每一类样品,可以获得残差图像中对应的区域。以PS 圆柱体为例,可以得到如图6(c)所示的残差图像区域,其中的像素值是在极坐标系下用原始图像减去滤波后的结果,期望它们服从高斯分布。图6(d)为该区域的灰度分布直方图,将满足

的像素定义为伪影像素,式中:x为该点像素值,µ和σ 为所有像素的均值和标准差,k为参数。当伪影严重时,k可以选择1 或更小;当伪影轻微时,可以选择2 或更大。在实验中选择k为1,可以得到图6(e)所示的伪影像素。对于图中的每个伪影像素,该值将被周围非伪影像素的平均值所替代。最终将各部分样品区域处理后的图像结合,得到的结果如图6(f)所示。

3 评价指标

为了验证所提出算法的有效性,将其与目前广泛使用的基于平滑滤波的方法[23]进行对比。在该方法中,平滑滤波器分别选择窗口大小为5 的标准高斯滤波器和窗口大小为5 的均值滤波器,并记作M1 和M2。使用Pfeiffer 等[24]的投影数据作为标准正弦图,将重建后的结果作为标准图像,并在正弦图中随机生成强条纹作为伪影数据,如图7 所示,最后通过算法将伪影去除并和标准图像进行比较。

图7 评价数据。(a)标准正弦图;(b)标准重建图;(c)具有伪影的正弦图;(d)具有伪影的重建图Fig.7 Evaluation data.(a)Standard sinogram;(b)Standard reconstruction image;(c)Sinogram with artifacts;(d)Reconstruction image with artifacts

图8 为使用不同算法对伪影数据校正得到的重建图像,使用峰值信噪比(PSNR)和结构相似度指数(SSIM)作为算法的评价标准,如表1 所示,其中PSNR 越大,SSIM 越接近1,说明图像越接近标准图像[25]。从图8 和表1 可以看出,所提出方法和基于平滑滤波的方法对图像中的环形伪影都有一定的抑制效果,但是基于平滑滤波的方法没有定位伪影,在滤波时会使用伪影数据,因此校正结果仍然存在伪影信息。所提出方法通过定位伪影,并针对伪影像素进行校正,因此具有更好的表现。

表1 不同算法的评价指标Table 1 Evaluation index of different algorithms

图8 不同算法的评价结果。(a)M1 算法;(b)M2 算法;(c)所提出算法Fig.8 Evaluation results of different algorithms. (a)M1 algorithm;(b)M2 algorithm;(c)The proposed algorithm

4 结 论

在基于光栅的X 射线相衬CT 中,采用正弦图域和重建图域相结合的图像处理算法,去除了吸收信号中的强伪影,并通过实验数据验证了所提出算法的有效性。该算法的主要优点是保留了样本的边界和细节,并且算法的框架可以随着图像分割算法的发展而不断改进。该算法有助于在未来的临床和安全相衬成像设备中实现快速和精确的检测。