近年来中国经济和贸易的快速发展，随着货运量和旅检量的飞速增长，各类违禁品层出不穷，检出难度越来越大。传统的DR扫描设备已经无法完全满足海关、民航等领域内的安全检查需求，因此在民航、旅检的场景下，CT检查方案^[1]及其衍生的技术手段（如基于X射线衍射分析的二级检查方案^[2]）已经成为主流的发展方向。但在海关的集装箱检查领域，由于扫描客体尺寸过大，CT检查设备的应用案例很少，目前仅有中国海关在深圳莲塘口岸部署了一台可以实现集装箱车辆整体CT扫描的安检设备^[3]。但对于大型货物/车辆CT扫描系统，其成像视野直径达6 m，探测器覆盖长度约7.5 m，庞大系统带来的复杂性使其环形伪影的情况十分严重。不仅存在一般CT系统中的探测器通道缺陷、能谱响应非线性、探测暗场和亮场特性不稳定等因素所带来的环形伪影，还存在加速器系统不稳定、成像系统振动、机械变形等问题引入的伪影，而且大型货物/车辆CT扫描系统的重量、尺寸都使得机械部件的振动、变形无法通过调试解决，从而导致大尺寸CT系统中的环形伪影不可避免地要采用算法进行校正。

对于CT图像中的环形伪影，常见的校正方法有3类：数据采集校正法、投影正弦图校正法和重建图像校正法^[4-5]。其中，数据采集校正法主要依靠硬件校正，常用的有暗场校正和增益不一致性校正^[6-8]；投影正弦图校正法一般采用傅里叶变换、小波变换转换到频域进行高频伪影滤波，或Hough变换滤除直线^[9-11]；传统重建图像校正法则是将笛卡尔坐标系下的图像变换至极坐标，再进行直线伪影滤除，此外最近还发展了滑动窗滤波法、TV正则化法、迭代法、深度学习法等方法^[12-15]。

在大型货物/车辆CT安全检查系统中，被检物种类繁多，结构差异性大，图像中心区域的环形伪影可能受到物体结构影响，难以从投影正弦图中识别，且终端用户一般无法获得投影正弦图进行前处理校正，因此，不适合采用投影正弦图校正法。另一方面，实际安检系统获得理想无伪影的图像难度大，迭代校正法处理时间长，因而基于深度学习^[16-19]和迭代校正^[20-22]的方法均无法较好的解决实际应用中的问题，需要寻找更适合的后处理校正算法。Hough变换法^[23]是一种应用较为广泛的直线检测方法，可为环形伪影校正提供较好的检测效果。但是，目前其应用方向主要在投影正弦图校正，对于其在后处理中应用的相关描述较少。此外，实际重建的大型货物/车辆图像一般含有较大噪声，单纯用Hough变换检测效率低，且容易误检，校正后会对物体自有结构带来较大影响，难以达到理想的效果。

对此，本文创新性地提出一种基于图像梯度特征聚类分析的Hough变换的后处理校正方法，对重建图像进行极坐标变换，利用变换后图像的梯度特征进行聚类分析与数据初筛，再对重组数据采用Hough变换直线检测，映射找出原图像中的伪影数据，进而插值消除极坐标图像的直线伪影，达到重建图像环形伪影消除的目的。

1.   算法介绍

利用重建图像进行后处理的校正算法，主要基于坐标变换法。通过极坐标变换将笛卡尔坐标系中的环状伪影变换为线状伪影，找出并去除线状伪影，再变换回笛卡尔坐标系（图1）。

图  1  重建图像从（a）笛卡尔坐标系转化为（b）极坐标

Figure  1.  Reconstructed images in (a) Cartesian coordinate and (b) polar coordinate

下载: 全尺寸图片 幻灯片

1.1   坐标变换

为了尽量不丢失图像信息，不直接按照原图像像素点$ \left(x,y\right) $进行变换，而是先确定极坐标图像各个点的像素值，按照公式：

从极坐标$ \left(r,\theta \right) $映射到笛卡尔坐标，再通过插值，确定极坐标图像各个像素点的值。由此，原图像中的环形伪影在极坐标中就呈现为与$ \theta $轴相互平行的直线。问题转换为如何寻找极坐标图像上的直线伪影。

1.2   Hough变换直线检测

Hough变换是直线检测应用较广泛的方法，基本过程：将直线上每个数据点变换为参数平面中的一条直线或曲线，利用共线数据点对应的参数曲线相交于参数空间中的一点的关系，将直线的提取问题转换为计数问题^[24]。设已知图像空间$ x\text{-}y $上的一条直线$ l $，在极坐标系下可表示为：

其中，$\rho \ge 0 $，$ 0\le \varphi \le \text{π} $，$ \rho $为直线$ l $到原点的距离，$ \varphi $为直线与$ y $轴正方向的夹角。直线$ l $上不同的点在参数平面$ \rho \text{-} \varphi $中转化为一簇相交的正弦曲线。通过统计特性找到峰值点，根据Hough变换的对偶性，便可确定参数平面上的峰值点对应的直线。

1.3   图像梯度特征提取

虽然Hough变换直线检测受直线中残缺部分影响小，但容易将被检物上原本的线状结构或不连续噪点误判为伪影直线，给图像造成信息丢失。因此本文采用一种修正的Hough变换法，在Hough变换直线检测之前，通过图像列梯度剔除不存在伪影的区域，减少误判。对于正常无环形伪影的图像，极坐标图像中列灰度值是渐变的；而环形伪影区域，图像列的灰度值存在突变现象。故而可以根据图像列灰度的一阶和二阶梯度数值来估计是否存在伪影。

为了降低图像噪声的影响，提高算法的鲁棒性，对极坐标图像沿着$ \theta $方向计算均值：

其中，$ {N}_{\theta } $为角度总数。

图2分别是无环形伪影和含环形伪影图像数据列均值的一阶和二阶梯度分布，容易看出伪影部分的梯度值存在明显的波动。

图  2  （a）无环形伪影和（b）含环形伪影图像的梯度分布

Figure  2.  Gradient distribution of images (a) without ring artifacts and (b) with ring artifacts

下载: 全尺寸图片 幻灯片

1.4   阈值选择

由于大型货物/车辆CT安全检查系统所面对的被检物种类繁多，图像列的一阶和二阶梯度的变化范围大，无法采用固定阈值的方式进行伪影提取。因此我们采用k-means算法对图像列梯度数值进行聚类分析，将可能存在伪影的列和正常图像进行区分。

k-means算法的处理过程：首先，随机地选择k个样本，每个样本初始地代表一个簇的平均值或中心；对剩余的每个样本，根据其与各簇中心的距离，赋给最近的簇；然后重新计算每个簇的平均值。不断重复此过程，直到簇内具有较高的相似度，而簇间的相似度较低时，准则函数收敛（图3）。

图  3  图像梯度值聚类分布

Figure  3.  Clustering distribution of image gradient values

下载: 全尺寸图片 幻灯片

将梯度数据分为两簇，基于k-means聚类定位簇中心后，以两个簇中心的中点$ {G}_{\mathrm{t}\mathrm{h}}=\left({C}_{1}+{C}_{2}\right)/2 $为阈值。

图2（b）虚线部分为$ \left[{-G}_{\mathrm{t}\mathrm{h}},{G}_{\mathrm{t}\mathrm{h}}\right] $范围，此范围外便认为存在粗大误差，可能受到伪影影响。将这些数据筛选出，重新组成一幅图像（图4），对其进行去噪增强等操作后，再利用1.2节的Hough变换法进行检测直线。

图  4  重组前后的极坐标系图像

Figure  4.  The images in polar coordinate before and after reorganization

下载: 全尺寸图片 幻灯片

1.5   处理流程

以列方向伪影为例，具体实施过程：首先计算极坐标图像列均值的梯度，然后根据k-means聚类确定阈值，选出误差较大的数据列，组成一幅新图像；将重组图像Hough变换至参数平面，进行直线检测，再筛选出角度接近垂直的直线；最后，插值修正伪影数据，并反变换回原图像空间（图5）。

图  5  算法处理流程

Figure  5.  The flow of algorithm processing

下载: 全尺寸图片 幻灯片

2.   测试结果

为验证上述方法，在Shepp-Logan模型上添加环形伪影和噪声，进行伪影消除效果对比：将无伪影和噪声的Shepp-Logan重建图像进行极坐标变换（图6（a）），在其上随机添加垂直线段和噪声（图6（b）），其中垂直线段在笛卡尔坐标系则对应环形伪影（图6（f））。在传统Hough变换法中，直接通过Hough变换进行直线检测，并筛选出角度接近垂直的线段（图6（c））。随后对找出的线段插值处理，修正其所在位置的灰度值，再反变换回笛卡尔坐标系，得到校正图像（图6（g））。在基于图像梯度特征提取的Hough变换法中，对待处理图像计算列均值的一阶梯度和二阶梯度，根据k-means聚类确定阈值，选出误差较大的数据列，组成一幅重组图像；将重组图像Hough变换至参数平面，进行直线检测，然后筛选出角度接近垂直的线段，并映射回原尺寸图像（图6（d））；再插值修正直线段，并反变换回原图像空间，得到最终图像图6（h）。

图  6  环形伪影消除结果对比

注：（a）（e）无伪影和噪声的图像，（b）（f）添加环形伪影和噪声的图像，（c）（g）基于传统Hough变换法的环形伪影消除图像，（d）（h）基于图像梯度聚类分析的环形伪影消除图像。两行分别对应极坐标系（a）～（d）和笛卡尔坐标系（e）～（h）的图像。

Figure  6.  Comparison of ring artifact elimination results

下载: 全尺寸图片 幻灯片

在上述对比结果中，两种方法所添加的环形伪影数量和位置相同，均为50条，其中传统Hough变换法校正31条，基于图像梯度特征提取的Hough变换法校正44条，多于传统Hough变换法。另一方面，图6（c）中传统Hough变换法在靠近旋转中心区域检测出多段归属于物体的自有结构，而基于图像梯度特征提取的Hough变换法由于对图像梯度进行先期筛选，大大减少了这些自有结构的误判。

更进一步地，将此方法应用于大型货物/车辆CT检查图像中。图7（a）为同方威视公司的FG9000 DT集装箱/车辆CT检查系统^[25]。图7（b）和图7（c）为用此系统扫描的断层图像，图7（b）为一幅风机叶片扫描图像，图7（c）为一个集装箱扫描图片，不难发现，重建图中存在大量环形或弧形伪影，对叶片结构和集装箱货物的分析判断造成严重干扰。

图  7  集装箱/车辆CT检查系统

Figure  7.  Cargo/Vehicle CT inspection system

下载: 全尺寸图片 幻灯片

根据算法过程，首先将重建图像进行极坐标变换。变换后如图8（a）和图8（d）所示。若采用传统Hough变换检测直线，图7（b）可以检测出39条斜率较低（$ \left|\theta \right| < $0.5° ）的线段（图8（b）），图7（c）检出了1578条线段（图8（e）），而通过图像列梯度做筛选再进行检测，风机叶片图像可以检测245条线段（图8（c）），集装箱货物可以检测249条线段（图8（f））。

图  8  风机叶片和集装箱图像的检测结果

注：（a）（d）极坐标变换图，（b）（e）Hough变换直线标记结果，（c）（f）基于图像梯度特征提取后Hough变换直线标记结果。

Figure  8.  Detection results of images of fan blade and container

下载: 全尺寸图片 幻灯片

得到伪影图像数据后，插值修正并反变换回原图像空间。图9展示了环形伪影的消除效果，其中图9（a）和图9（c）采用了直接Hough变换法，图9（b）和图9（d）采用了基于图像梯度特征聚类分析的Hough变换法。

图  9  风机叶片和集装箱货物在传统Hough变换法（a）（c）和基于图像梯度特征提取的Hough变换法（b）（d）下环形伪影消除效果

Figure  9.  The ring artifacts elimination effects of fan blade and contaner under traditional Hough transform method (a) (c) and Hough transform method based on image gradient feature extraction (b) (d)

下载: 全尺寸图片 幻灯片

3.   讨论

通过在Shepp-Logan模型、风机叶片断层图和集装箱断层图上的测试可以看出，基于图像梯度特征聚类分析的Hough变换校正方法对断层图像的环形伪影能够进行有效的修正。

从图6标记的检测结果看，在Shepp-Logan模型上新算法识别出了50条伪影中的44条，而且没有出现误判的情况，而传统Hough变换法只检测出了31条，且在旋转中心附近出现了许多误判的情况。在风机叶片断层图和集装箱断层图上，新算法分别检测出了245条和249条环形伪影，而传统Hough变换法却检测出了39条和1578条。由于实物图像没有无伪影的标准图，无法判断检测结果的准确性，但从图8（b）和图8（e）中标记的处理结果可以明显看出，传统Hough变换法在被检物较为简单时（风机叶片）有大量的漏检伪影，在被检物较为复杂时（集装箱货物）却又有大量的误检。

另一方面，新算法在图8（b）和图8（e）上标记的检测结果基本合理，而且两次的伪影检测数量相当，位置相近，说明新算法能够有效检测出由于系统固定偏差产生的较明显的环形伪影，但从图9（b）和图9（d）的处理效果看，在实际图像较为复杂时，极坐标图像列梯度受到物体本身结构的影响后，梯度的k-means聚类结果无法准确分类出较不明显的环形伪影，因此校正后的图像还存在一些环形伪影。

4.   总结

本文提出一种基于图像梯度特征聚类分析的环形伪影后处理校正方法。利用环状伪影位置灰度值突变的特点，通过极坐标系中图像列均值的一阶和二阶梯度进行k-means聚类分析，将可能存在环状伪影的图像列和无伪影的图像列进行区分，再对可能存在环状伪影的区域利用Hough变换进行直线检测，定位伪影数据后插值消除图像的伪影。由于通过图像列梯度对图像的环形伪影情况进行了预先分类，本文的方法可以提高环形伪影的检测准确性，减少对正常图像的误报。

在Shepp-Logan模型上的模拟测试和大型货物/车辆CT检查系统图像的实际测试都表明新算法消除环状伪影的效果能够明显提升，同时减少了对被测物结构的误处理，图像的细节得到有效保存，但实测的图像效果也表明由于图像列均值的一阶和二阶梯度的k-means聚类结果受到被测物结构形状的影响，部分环形伪影无法被检出，但通过优化极坐标图像的分类方法，有望进一步改进算法效果。