基于Mean Shift算法的CT图像金属伪影消除方法及消除系统

摘要

本发明提供了一种基于Mean Shift算法的CT图像金属伪影消除方法和消除系统。通过预处理CT图像去除噪声并平滑轻度伪影，分割出CT图像中的金属区域，对金属区域的投影数据作线性插值生成插值图像，并修复组织信息后得到先验图像，再将先验图像的投影数据替代原投影数据，得到校正后的CT图像。与相关技术相比，本发明提供的基于Mean Shift算法的CT图像金属伪影消除方法及消除系统可以去除CT图像中的轻度金属伪影，平滑效果好，不丢失细节信息；金属区域分割效果好，计算效率更高；去除伪影效果明显，而且还保护了原CT图像中人体组织结构信息。

说明书

技术领域

本发明涉及用于计算机断层摄影(CT)中减少金属伪影的技术领域，尤其涉及一种基于Mean Shift算法的CT图像金属伪影消除方法及消除系统。

背景技术

随着医学技术的发展，CT图像已经在临床医学这一领域得到了广泛的应用，极大地提高了医学诊断水平。然而，CT技术利用X射线束扫描人体层面后接收再成像，对于一些体内含有金属移植物的特殊病人，如金属补牙、人工关节、心脏起搏器等等，由于射线硬化、散射、噪声以及部分容积效应等原因，很有可能得到含有金属伪影的CT图像。金属伪影使得CT图像质量变差，组织结构难以判断，可能导致误诊的结果。因此，CT图像的金属伪影消除(MAR,Metal Artifact Reduction)对提高临床诊断的准确性有很重要的意义。MAR技术是医学图像处理领域的一个热门研究方向，许多学者为此做出贡献，所涉及到的图像处理技术可大致分为三大类，分别是投影修补法、迭代法以及以上两种方法的混合法。投影修补法主要通过插值、先验图像等方法生成新的投影数据，替代被腐蚀的金属区域的投影数据，从而消除CT图像的金属伪影，一般而言，算法简单复杂度较低，但容易丢失组织信息。迭代法利用代数或者统计学的方法，根据设定的优化准则来不断迭代校正图像，在噪声较大时也能得到较好的结果，但是算法复杂度相对较高。混合法是以上两种方法的结合，互相取长补短，发现显著提高图像质量，且计算效率较高的算法。其中投影修补法是MAR技术的有别于其他图像处理的特色处理方法，很多研究者都提出了自己的投影修补算法，但仍未得到理想的结果，这个方向仍然有探索的空间。

发明内容

本发明需要解决的技术问题是提供一种基于Mean Shift算法的CT图像金属伪影消除系统，其去除伪影效果明显，能够保护人体组织结构的特点。

本发明提供了一种基于Mean Shift算法的CT图像金属伪影消除方法，该方法包括如下步骤：

S1、通过自适应Mean Shift平滑算法预处理CT图像，对含有金属伪影的CT图像进行平滑滤波，消除CT图像中的噪声并平滑轻度伪影；

S2、采用简化的Mean Shift分割算法，通过分析CT图像的特征空间和聚类方法得到CT图像中的金属区域；

S3、将金属区域和预处理CT图像前投影得到投影数据，对于坐标x处的投影数据，采用线性插值算法得到校正后的投影数据其中，β为投影角度，{p^β,q^β}为插值空间；

S4、根据线性插值后的投影数据反投影重建，得到校正伪影后不含金属区域的CT图像v_li＝{v_li(i)|i∈I_li}，其中，i为CT图像中的像素点，I_li代表线性插值处理后的图像；

S5、用线性插值得到的图像替换原CT图像中的金属区域，得到不含金属区域、组织信息完好的CT图像v_nom＝{v_nom(i)|i∈I_nometal}，进而得到差异图像v_diff＝{v_nom(i)-v_li(i)|i∈I_diff}，则最终的先验图像v_prior可表示为：v_prior(i)＝v_li(i)+v_diff(i)*w(i)，其中，I_nometal代表不包含金属的CT图像，I_diff代表原CT图像和线性插值图像的差值图像，权重系数对式中的v_diff(i)归一化处理，σ取值区间为(0,1)；

S6、将先验图像前投影得到校正后的投影数据，替换原CT图像中金属区域的投影数据，完成投影修补；

S7、将得到的投影数据反投影重建，并添加步骤S2中分割得到的金属区域，生成最终的校正图像。

优选的，所述自适应Mean Shift平滑算法进行平滑的步骤包括：

(1)将CT图像中的每个像素点的位置和像素值构成特征空间向量x；

(2)从CT图像左上方的第一个像素点i开始，利用公式计算偏移均值向量M(x)，其中w(x)为样本点的权重，G(x)为高斯核函数，h为带宽；

(3)设置一个停止条件ε，如果M(x)>ε，将M(x)+x的值赋给x，迭代执行步骤(2)，直至M(x)≤ε，或者达到迭代次数上限而结束迭代，开始下一个像素的计算，直至遍历整个CT图像。

优选的，所述带宽h的确定步骤包括：

(1)给定初始带宽h₀，并设置步长Δh、阈值δ；

(2)若带宽内的点与中心点的CT值之差|p-p₀|≤δ，则判定两点是相似的，统计相似点n₀和总采样点n；

(3)若n₀≥75％*n，则h＝h+Δh，迭代执行步骤(2)，直至满足n₀<65％*n，或者达到迭代次数上限而结束迭代。

优选的，所述迭代次数上限为五次。

优选的，在步骤S2中，对于d维空间中n个样本点x_i,i＝1,...,n，在x点的偏移均值向量的基本形式为：其中S_h是半径为h的高维球，S_h(x)＝{y:(y-x)^T(y-x)≤h²}，k是n个样本点中落入高维球的点的个数，则M_h(x)是落入高维球区域的样本点相对于x点的偏移量再取平均值。

优选的，所述Mean Shift分割算法的实现步骤包括：

(1)将每个像素点的位置和像素值构成特征空间向量x，并设置一个带宽h；

(2)随机选取初始的聚类中心x，计算每个像素点到聚类中心的距离，若(x_i-x)^T(x_i-x)≤h²，标记x_i并给该聚类投票；

(3)利用上述公式计算M_h(x)，得到新的聚类中心x'，若x'与x的距离小于h/2，将两个聚类合并，否则生成一个新的聚类，迭代执行步骤(2)，直至所有的点都被标记；

(4)每个像素点属于票数最多的聚类，从而得到每个聚类的所有像素点，像素值最高的聚类为金属区域。

本发明还提供了一种基于Mean Shift算法的CT图像金属伪影消除系统，包括：

预处理模块，用于通过自适应Mean Shift平滑算法预处理CT图像，对含有金属伪影的CT图像进行平滑滤波，消除CT图像中的噪声并平滑轻度伪影；

金属分割模块，用于采用简化的Mean Shift分割算法，通过分析CT图像的特征空间和聚类方法得到CT图像中的金属区域；以及

投影修补模块，用于对金属区域的投影数据作线性插值生成插值图像，并修复组织信息后得到先验图像，用先验图像的投影数据替代原投影数据，得到校正后的图像。

与相关技术相比，本发明提供的基于Mean Shift算法的CT图像金属伪影消除方法及消除系统通过对CT图像进行预处理，可以去除CT图像中的轻度金属伪影，另外，通过自适应选取平滑窗口大小，平滑效果好，不丢失细节信息；通过对CT图像采用简化的Mean Shift分割算法，不但金属区域分割效果好，而且计算效率更高；通过对插值图像的组织信息的修复，去除伪影效果明显，而且还保护了原CT图像中人体组织结构信息。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图，其中：

图1为本发明基于Mean Shift算法的CT图像金属伪影消除方法的流程图；

图2(a)为利用本发明基于Mean Shift算法的CT图像金属伪影消除方法对CT图像进行平滑处理前的效果图；

图2(b)为利用本发明基于Mean Shift算法的CT图像金属伪影消除方法对CT图像进行平滑处理后的效果图；

图3(a)为未进行金属分割的原CT图像图；

图3(b)为对图3(a)所示原CT图像采用阈值分割出的金属区域图；

图3(c)为对图3(a)所示原CT图像采用区域生长法分割出的金属区域图；

图3(d)为对图3(a)所示原CT图像采用Mean Shift算法分割出的金属区域图；

图4(a)为含金属伪影的原CT图像图；

图4(b)为对图4(a)所示含金属伪影的原CT图像图进行预处理的CT图像效果图；

图4(c)为对图4(a)所示含金属伪影的原CT图像图进行金属分割的CT图像效果图；

图4(d)为对图4(a)所示含金属伪影的原CT图像图进行线性插值的CT图像效果图；

图4(e)为对图4(a)所示含金属伪影的原CT图像图进行处理后的先验图像效果图；

图4(f)为对图4(a)所示含金属伪影的原CT图像图进行处理的校正图像效果图；

图5(a)为本发明未进行基于Mean Shift算法的CT图像金属伪影消除方法处理的原CT图像图；

图5(b)为对图5(a)所示原CT图像进行传统线性插值算法处理的效果图；

图5(c)为本发明对图5(a)所示原CT图像进行基于Mean Shift算法的CT图像金属伪影消除方法处理的效果图；

图6为本发明基于Mean Shift算法的CT图像金属伪影消除系统的框架图。

具体实施方式

下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

请参阅图1，本发明提供了一种基于Mean Shift算法的CT图像金属伪影消除方法，包括如下步骤：

S1、通过自适应Mean Shift平滑算法预处理CT图像，对含有金属伪影的CT图像进行平滑滤波，消除CT图像中的噪声并平滑轻度伪影；

S2、采用简化的Mean Shift分割算法，通过分析CT图像的特征空间和聚类方法得到CT图像中的金属区域；

S3、将金属区域和预处理CT图像前投影得到投影数据，对于坐标x处的投影数据，采用线性插值算法得到校正后的投影数据其中，β为投影角度，{p^β,q^β}为插值空间；

S4、根据校正后的投影数据反投影重建，得到校正后不含金属区域的CT图像v_li＝{v_li(i)|i∈I_li}，其中，i为CT图像中的像素点，I_li是线性插值处理后的图像；

S5、用线性插值得到的图像替换原CT图像中的金属区域，得到不含金属区域、组织信息完好的CT图像v_nom＝{v_nom(i)|i∈I_nometal}，进而得到差异图像v_diff＝{v_nom(i)-v_li(i)|i∈I_diff}，则最终的先验图像v_prior可表示为：v_prior(i)＝v_li(i)+v_diff(i)*w(i)，其中，I_nometal代表不包含金属的CT图像，I_diff代表线性插值图像和原CT图像的差值图像，权重系数对式中的v_diff(i)归一化处理，σ取值区间为(0,1)；

S6、将先验图像前投影得到校正后的投影数据，替换原CT图像中金属区域的投影数据，完成投影修补；

S7、将得到的投影数据反投影重建，并添加步骤S2中分割得到的金属区域，生成最终的校正图像。

步骤S1用于对CT图像进行图像预处理，应用了自适应Mean Shift平滑算法，不但可以去除程度较轻的伪影，避免出现过度校正的情况，而且不会破坏原有CT图像的组织结构。

所述自适应Mean Shift平滑算法的基本原理是根据初始点计算偏移均值向量，移动初始点到其偏移均值，重新计算偏移均值向量，不断迭代此过程，直至满足结束条件，是一种基于特征空间的概率密度梯度函数的估计算法。

基于以上原理，引入核函数和样本点权重的概念，核函数使得样本与中心点的距离不同，其偏移量对偏移均值向量的贡献也不同，而权重则体现了每个样本点重要性的不同。综上，偏移均值向量可阐述为：

其中，w(x)为样本点的权重，离x越近，采样点权重越大；G(x)为核函数，常用的有两种，分别是单位核函数和高斯核函数，本系统采用高斯核函数；h为带宽，h越大，核函数越平坦，图像越平滑。

对于d维空间的n个采样点，密度估计函数为：

其中，核函数K(x)＝c_k,dk(||x||²)，k(||x||²)称为核函数的剖面函数，h为带宽，w(x)为采样点的权重。

由上式可得，设k(x)的负导函数为g(x)，其对应的核函数为G(x)＝c_g,dg(||x||²)，则有：

综上可得，说明M_h,G(x)总是指向概率密度增加最快的方向，是一个归一化的概率密度梯度函数。Mean>

所述自适应Mean Shift平滑算法进行平滑的主要步骤包括：

(1)将CT图像中的每个像素点的位置和像素值构成特征空间向量x；

(2)从CT图像左上方的第一个像素点i开始，利用公式计算偏移均值向量M(x)，其中w(x)为样本点的权重，G(x)为高斯核函数，h为带宽；

(3)设置一个停止条件ε，如果M(x)>ε，将M(x)+x的值赋给x，迭代执行步骤(2)，直至M(x)≤ε，或者达到迭代次数上限而结束迭代，开始下一个像素的计算，直至遍历整个CT图像。

在用自适应Mean Shift平滑算法对CT图像平滑的过程中，带宽h是一个很重要的参数。如果h偏大，CT图像可能过于平滑而丢失细节，如果h偏小，可能达不到平滑的效果。人为选取固定带宽时，很容易出现上述情况，影响最终图像的处理结果。因此，本发明自适应选择带宽h，所述带宽h的确定步骤包括：

(1)给定初始带宽h₀，并设置步长Δh、阈值δ；

(2)若带宽内的点与中心点的CT值之差|p-p₀|≤δ，则判定两点是相似的，统计相似点n₀和总采样点n；

(3)若n₀≥75％*n，h＝h+Δh，迭代执行步骤(2)，直至满足n₀<65％*n，或者达到迭代次数上限而结束迭代。所述迭代次数上限为五次。

自适应带宽选择使得在细节信息不多的区域带宽较大，提高计算效率，而在细节和边缘信息较为密集的区域带宽较小，不丢失细节信息，以达到更好的平滑效果。

请一并参阅图2(a)和(b)，图2(a)为利用本发明基于Mean Shift算法的CT图像金属伪影消除方法对CT图像进行平滑处理前的效果图，图2(a)为利用本发明基于Mean Shift算法的CT图像金属伪影消除方法对CT图像进行平滑处理后的效果图。可见，对于1、2号位程度较轻的伪影，Mean Shift平滑算法可以达到比较理想的平滑效果，而且不会破坏骨骼等组织结构。对于3号位的严重伪影，也有一定的效果但不够明显，还需要进一步的处理。

所述步骤S2用于对所述CT图像进行金属区域分割，金属区域的分割是决定最终校正图像质量的关键因素，传统的阈值分割法虽然简单高效，但是容易将CT值高的骨骼组织误分类成金属，不能得到正确的分割结果。区域生长法一定程度上解决了误分类的问题，比阈值法的分割结果更加准确，但是基于由种子点根据相似性向外区域扩充的原理，分割出的金属边缘扩大不清晰，导致分割结果大于正确的金属区域。本系统采用一种简化的MeanShift算法进行金属分割，能够锐化金属边界，解决远离金属的骨骼误分类的问题，得到较满意的分割结果。

因此，本发明采用简化的Mean Shift分割算法，通过分析图像的特征空间和聚类方法来达到分割的目的。对于d维空间中n个样本点x_i,i＝1,...,n，在x点的偏移均值向量的基本形式为：其中S_h是半径为h的高维球，S_h(x)＝{y:(y-x)^T(y-x)≤h²}，k是n个样本点中落入高维球的点的个数，则M_h(x)是落入高维球区域的样本点相对于x的偏移量再取平均值。

基于上述原理，可以实现所述Mean Shift分割算法简化，实现步骤包括：

(1)将每个像素点的位置和像素值构成特征空间向量x，并设置一个带宽h；

(2)随机选取初始的聚类中心x，计算每个像素点到聚类中心的距离，若(x_i-x)^T(x_i-x)≤h²，标记x_i并给该聚类投票；

(3)利用公式计算M_h(x)，得到新的聚类中心x'，若x'与x的距离小于h/2，将两个聚类合并，否则生成一个新的聚类。迭代执行步骤(2)，直至所有的点都被标记；

(4)每个像素点属于票数最多的聚类，从而得到每个聚类的所有像素点，像素值最高的聚类为金属区域。

如图3(a)所示，为原CT图像，阈值法、区域生长法、Mean Shift算法对同一幅如图3(a)所示的CT图像处理得到的金属物体分割结果如图3(b)～(d)所示，图(b)所示的阈值分割和图(c)所示的区域生长法均出现了骨骼误分类的情况，而且边缘不够清晰，细节较差，而图(d)所示的Mean Shift算法能够更准确地把金属物体从背景中分割出来。

所述步骤S3～S7用于投影修补，以得到最终的校正图像。原始投影数据中金属投影值比其他人体组织的投影值金属高很多，所以在重建CT图像时容易生成伪影，影响图像质量，为了解决这一问题，需要用到插值法来修补投影数据。CT图像金属伪影校正中用到的插值方法有很多，如线性插值、非线性插值、三次样条插值等等。经过比较分析，本系统采用线性插值修补投影数据，线性插值能够有效消除大部分的金属伪影，而且简单高效，容易实现。

请参阅图4(a)～(f)所示的各阶段的效果图，可以看出，采用自适应Mean Shift平滑算法平滑图像，简化的Mean Shift分割算法分割金属区域，线性插值技术修补伪影，由不含金属物体的原CT图像的和插值图像的差异图像处理后生成先验图像，并由先验图像得到最终校正图像的去除CT图像中金属伪影的算法流程。

采用本系统提出的CT图像金属伪影消除算法处理一系列带金属伪影的CT图像，其中包括轻度伪影、重度伪影，图像中单个金属物体、多个金属物体等三组数据，实验结果及与传统线性插值算法的处理结果对比如图5(a)～(c)所示。

分析这三组数据，当伪影较轻时，传统线性插值可能造成过度校正的结果，而本系统提出的平滑算法处理得到的预处理图像，可以去除金属伪影，并保护图像中原有的组织信息。当图像仅含有一个金属物体时，两种图像处理算法均能有效去除伪影，但是相比较线性插值方法，本系统提出的算法能够避免次级伪影的产生，有更好的处理效果。当图像中含有多个金属，并且金属周围有骨骼组织时，线性插值很有可能破坏组织结构，但是本系统提出的算法在去除伪影的同时考虑了原有图像中人体组织信息的保护，因此有更好的处理效果。

请参阅图6，本发明还提供了一种基于Mean Shift算法的CT图像金属伪影消除系100，包括：

预处理模块1，用于通过自适应Mean Shift平滑算法预处理CT图像，对含有金属伪影的CT图像进行平滑滤波，消除CT图像中的噪声并平滑轻度伪影；

金属分割模块2，用于采用简化的Mean Shift分割算法，通过分析CT图像的特征空间和聚类方法得到CT图像中的金属区域；以及

投影修补模块3，用于对金属区域的投影数据作线性插值生成插值图像，并修复组织信息后得到先验图像，用先验图像的投影数据替代原投影数据，得到校正后的图像。

与相关技术相比，本发明提供的基于Mean Shift算法的CT图像金属伪影消除方法及消除系统通过对CT图像进行预处理，可以去除CT图像中的轻度金属伪影，另外，通过自适应选取平滑窗口大小，平滑效果好，不丢失细节信息；通过对CT图像采用简化的Mean Shift分割算法，不但金属区域分割效果好，而且计算效率更高；通过对插值图像的组织信息的修复，去除伪影效果明显，而且还保护了原CT图像中人体组织结构信息。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的专利保护范围内。