CN109712212B - 一种工业ct伪影校正方法 - Google Patents

Abstract

本发明提供一种工业CT伪影校正方法，包括步骤：基于正投影法获得伪影场，对其进行阈值分割获得包含被测对象及受伪影强干扰的背景区域S(x,y)；实际CT图像Q(x,y)除去S(x,y)得到伪影弱干扰背景图像G(x,y)，获取其灰度直方统计图中最大峰值的灰度值；Q(x,y)中包括S(x,y)得到前景及伪影强干扰背景图像H(x,y)，获取[0,T]的最佳拟合高斯曲线；计算其概率密度函数并归一化处理；统计分析Q(x,y)中每个像素点的局部灰度，获取占比灰度值最大的灰度，计算其对应的归一化概率密度函数；对受伪影强干扰的背景区域灰度值校正及对Q(x,y)各个灰度值校正，与现有技术相比，本校正方法在针对单一材料复杂内部结构工件的CT扫描时，能有效降低伪影引起的灰度差异，实现精确、可靠、无损测量其内部结构尺寸。

Description

一种工业CT伪影校正方法

技术领域

本发明涉及工业CT成像技术领域，具体涉及一种工业CT伪影校正方法。

背景技术

工业CT(industrial computerized tomography)，全称是工业用计算机断层成像技术，是指应用于工业中的核成像技术，其基本原理是依据辐射在被检测物体中的减弱和吸收特性，同物质对辐射的吸收本领与物质性质有关。所以，利用放射性核素或其他辐射源发射出的、具有一定能量和强度的X射线或γ射线，在被检测物体中的衰减规律及分布情况，就有可能由探测器陈列获得物体内部的详细信息，最后用计算机信息处理和图像重建技术，以图像形式显示出来。工业CT能在对检测物体无损伤条件下，以二维断层图像或三维立体图像的形式，清晰、准确、直观地展示被检测物体的内部结构、组成、材质及缺损状况，是目前国内外对零部件内部尺寸测量所普遍采用的方法。

但是工业CT设备本身存在很多不可避免的缺陷，与常规的X射线摄像相比，更容易产生伪影，伪影是指实际物体被扫描时，重建数值与物体真实衰减系数之间的差异。从严格意义上讲，伪影是由非理想成像系统造成的受检测物体内部结构图像的伪迹特征。对于CT成像装置而言，各种原因导致的伪影类型众多，如：混叠伪影、部分容积效应伪影、散射伪影、射束硬化伪影、噪声伪影、运动伪影以及金属伪影等。CT图像中的伪影使重建图像的质量严重下降，易引起尺寸测量的误差。综上，伪影的校正研究对于工业CT成像具有极其重要的意义。

发明内容

本发明所要解决的技术问题是针对上述现有技术现状提供一种工业CT伪影校正方法。

本发明解决上述技术问题所采用的技术方案为：一种工业CT伪影校正方法，包括以下步骤：

(1)基于正投影法获得被测对象的伪影场F(x,y)并对其进行阈值分割，设所述伪影场F(x,y)中灰度值在[a,b]范围内的图像区域为S(x,y)，该区域包含被测对象及受伪影强干扰的背景区域；

(2)实际CT图像为Q(x,y)，提取所述实际CT图像Q(x,y)中除去S(x,y)的背景区域，得到伪影弱干扰背景图像G(x,y)，对G(x,y)进行灰度直方图h₂(Z)统计，并计算h₂(Z)的最大峰值的像素灰度值g；

(3)提取实际CT图像Q(x,y)中包括S(x,y)的背景区域，得到前景及伪影强干扰背景图像H(x,y)，对H(x,y)进行灰度直方图h₁(Z)统计，并通过最大类间方差法获取将H(x,y)中前景与伪影强干扰背景图像进行分割的阈值T，对h₁(Z)在[0,T]范围内的灰度直方图h₁(Z_i)进行最小残差法高斯拟合，获取最佳拟合高斯曲线f₁(i)；

(4)计算步骤(3)中f₁(i)的概率密度函数PDF(x)，归一化处理后得到PDF_nor(x)；

(5)对实际CT图像Q(x,y)的每个像素点(a,b)进行局部灰度统计分析，获取该像素点中占比像素灰度值最大的灰度V(a,b)，并计算该灰度V(a,b)对应的归一化的概率密度函数PDF_nor(V(a,b))；

(6)计算h₁(Z)在[0,T]范围内的最大峰值像素灰度值w，即为受伪影强干扰的背景区域的像素最大峰值灰度值；计算h₁(Z)在[T+1,L-1]范围内的最大峰值像素灰度值q，即为被测对象区域的像素最大峰值灰度值；其中，L-1为最大像素灰度值；

(7)对受伪影强干扰的背景区域的像素灰度值k进行校正，校正后的像素灰度值

(8)对Q(x,y)各个像素灰度值进行校正处理，校正后的CT图像

进一步地，步骤(1)具体包括：设被测对象为T(x,y)，其中0≤x≤m，0≤y≤n，

其中，

为被测对象内部及表面区域，图像中任意一点(a,b)在任意角度θ的线积分K_θ(a,b)为：

设伪影正投影数组为f(x,y)，利用K_θ(a,b)进行周向累积计算，

将f(x,y)进行灰度级映射成像F(x,y)，

伪影场

其中，Z表示所成图像的像素灰度级，S(x,y)＝{F(x_i,y_i)∈[a,b]}。

进一步地，步骤(2)中，G(x,y)＝Q(x,y)-Q(x,y)∩S(x,y)，对G(x,y)进行灰度直方图h₂(Z)统计，

其中，

Z＝0,1,...,L-1，Z表示所成图像的像素灰度值。

进一步地，步骤(3)中，H(x,y)＝Q(x,y)∩S(x,y)，对H(x,y)进行灰度直方图h₁(Z)统计，

通过最大类间方差法获取将H(x,y)中前景与伪影强干扰背景图像进行分割的阈值T包括：设任意分割值T_i，图像中像素灰度值小于阈值T_i的像素个数

其平均像素灰度值

图像中像素灰度值大于阈值T_i的像素个数

其平均像素灰度值

类间方差

采用遍历的方法获得类间方差g最大时的阈值T，T＝max(g)；

对h₁(Z)在[0,T]范围内的灰度直方图h₁(Z_i)进行最小残差法高斯拟合，拟合后的高斯曲线函数

其中，a是幅值，

是灰度值均值，σ是灰度值标准差，高斯拟合曲线迭代的初始灰度值均值

初始灰度值标准差

将

和σ_init代入公式

中，得到

其中，幅值a_init的初始值为

则得到幅值的初始值

其中，灰度值均值

幅值

灰度值标准差

选取不同的拟合后高斯曲线的灰度值均值

幅值a以及灰度值标准差σ，建立函数类Ψ＝{f₀(i),f₁(i),...,f_m(i)}(m<L)，求出使

时，最佳拟合高斯曲线灰度值均值

以及最佳灰度值标准差σ_opt。

进一步地，步骤(4)中，拟合后的高斯曲线函数f₁(i)的概率密度函数

其归一化概率密度函数

进一步地，步骤(5)中，对实际CT图像Q(x,y)的每个像素点(a,b)进行局部灰度统计分析的像素点范围为大小为3×3的像素矩阵。

与现有技术相比，本发明的优点在于：针对单一材料复杂内部结构工件的CT扫描时，能有效降低伪影引起的灰度差异，实现精确、可靠、无损测量其内部结构尺寸。

附图说明

图1为本发明优选实施例的校正方法的流程图；

图2为球阵列的实物图；

图3为球阵列的CT扫描图；

图4为球阵列位置1的实际CT图像；

图5为球阵列位置2的实际CT图像；

图6为球阵列位置1的正投影仿真CT图；

图7为球阵列位置2的正投影仿真CT图；

图8为去除受伪影强干扰的背景区域灰度直方图；

图9为被测对象及受伪影强干扰的背景区域图；

图10为被测对象及受伪影强干扰的背景区域灰度直方图；

图11为受伪影强干扰的背景区域灰度统计分布及高斯拟合结果图；

图12为球阵列校正后的CT图像；

图13为实际CT图像阈值分割效果图；

图14为校正后CT图像阈值分割效果图；

图15为实际CT局部的图像阈值分割效果图；

图16为校正后CT局部的图像阈值分割效果图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

参阅图1所示，为本优选实施例提供的一种工业CT伪影校正方法的流程图，包括以下步骤：

(1)基于正投影法获得被测对象的伪影场F(x,y)并对其进行阈值分割，设所述伪影场F(x,y)中灰度值在[a,b]范围内的图像区域为S(x,y)，该区域包含被测对象及受伪影强干扰的背景区域；

(2)实际CT图像为Q(x,y)，提取所述实际CT图像Q(x,y)中除去S(x,y)的背景区域，得到伪影弱干扰背景图像G(x,y)，对G(x,y)进行灰度直方图h₂(Z)统计，并计算h₂(Z)的最大峰值的像素灰度值g；

(3)提取实际CT图像Q(x,y)中包括S(x,y)的背景区域，得到前景及伪影强干扰背景图像H(x,y)，对H(x,y)进行灰度直方图h₁(Z)统计，并通过最大类间方差法获取将H(x,y)中前景与伪影强干扰背景图像进行分割的阈值T，对h₁(Z)在[0,T]范围内的灰度直方图h₁(Z_i)进行最小残差法高斯拟合，获取最佳拟合高斯曲线f₁(i)；

(4)计算步骤(3)中f₁(i)的概率密度函数PDF(x)，归一化处理后得到PDF_nor(x)；

(5)对实际CT图像Q(x,y)的每个像素点(a,b)进行局部灰度统计分析，获取该像素点中占比像素灰度值最大的灰度V(a,b)，并计算该灰度V(a,b)对应的归一化的概率密度函数PDF_nor(V(a,b))；

(6)计算h₁(Z)在[0,T]范围内的最大峰值像素灰度值w，即为受伪影强干扰的背景区域的像素最大峰值灰度值；计算h₁(Z)在[T+1,L-1]范围内的最大峰值像素灰度值q，即为被测对象区域的像素最大峰值灰度值；其中，L-1为最大像素灰度值；

(7)对受伪影强干扰的背景区域的像素灰度值k进行校正，校正后的像素灰度值

(8)对Q(x,y)各个像素灰度值进行校正处理，校正后的CT图像

进一步地，步骤(1)具体包括：设被测对象为T(x,y)，其中0≤x≤m，0≤y≤n，

其中，

为被测对象内部及表面区域，图像中任意一点(a,b)在任意角度θ的线积分K_θ(a,b)为：

设伪影正投影数组为f(x,y)，利用K_θ(a,b)进行周向累积计算，

将f(x,y)进行灰度级映射成像F(x,y)，

伪影场

其中，Z表示所成图像的像素灰度级，S(x,y)＝{F(x_i,y_i)∈[a,b]}。

进一步地，步骤(2)中，G(x,y)＝Q(x,y)-Q(x,y)∩S(x,y)，对G(x,y)进行灰度直方图h₂(Z)统计，

其中，

Z＝0,1,...,L-1，Z表示所成图像的像素灰度值。

进一步地，步骤(3)中，H(x,y)＝Q(x,y)∩S(x,y)，对H(x,y)进行灰度直方图h₁(Z)统计，

通过最大类间方差法获取将H(x,y)中前景与伪影强干扰背景图像进行分割的阈值T包括：设任意分割值T_i，图像中像素灰度值小于阈值T_i的像素个数

其平均像素灰度值

图像中像素灰度值大于阈值T_i的像素个数

其平均像素灰度值

类间方差

采用遍历的方法获得类间方差g最大时的阈值T，T＝max(g)；

对h₁(Z)在[0,T]范围内的灰度直方图h₁(Z_i)进行最小残差法高斯拟合，拟合后的高斯曲线函数

其中，a是幅值，

是灰度值均值，σ是灰度值标准差，高斯拟合曲线迭代的初始灰度值均值

初始灰度值标准差

将

和σ_init代入公式

中，得到

其中，幅值a_init的初始值为

则得到幅值的初始值

其中，灰度值均值

幅值

灰度值标准差

选取不同的拟合后高斯曲线的灰度值均值

幅值a以及灰度值标准差σ，建立函数类Ψ＝{f₀(i),f₁(i),...,f_m(i)}(m<L)，求出使

时，最佳拟合高斯曲线灰度值均值

以及最佳灰度值标准差σ_opt。

进一步地，步骤(4)中，拟合后的高斯曲线函数f₁(i)的概率密度函数

其归一化概率密度函数

进一步地，步骤(5)中，对实际CT图像Q(x,y)的每个像素点(a,b)进行局部灰度统计分析的像素点范围为大小为3×3的像素矩阵。

校正试验：

为了验证本校正方法的可靠性，提高工业CT尺寸测量的精度，通过本校正方法将校正后的CT图像与未校正的CT图像采用行业惯用的VG(VGSTUDIO MAX)进行复杂点阵镂空结构的尺寸测量实验。

本优选实施例采用球阵列的复杂结构模拟件作为实验被测对象，高精度不锈钢球排列成空间阵列，外部采用PMMA(polymethyl methacrylate，聚甲基丙烯酸甲酯)固定件密封，该模拟件具有规律的内部结构，且空间范围内各不锈钢球的直径基本一致，可较为全面的反应空间各位置尺寸测量的精度，图2为球阵列的实物图，图3为球阵列的CT扫描图。

本优选实施例使用的CT系统为YXLon FXE 225.48射线源和Varian 4030E探测器组成，实验使用的扫描参数如表1所示。

表1

参数名称	数值
		扫描电压	160kV
扫描电流	300uA
		探测器帧速率	2fps
探测器像素个数	3200×2304

选取实际CT图像和正投影仿真图像进行对比，图4和图5为球阵列位置1和位置2的实际CT图像，图6和图7为球阵列位置1和位置2的正投影仿真CT图，基于对尺寸测量精度的验证，位置1和位置2的选取主要关注伪影的分布情况及对被测对象边缘区域的干扰情况，从图中可以看出，实际CT图像背景中存在的条状伪影与正投影仿真CT图中存在的伪影基本保持一致，因此可证明正投影伪影仿真法的有效性和可靠性。

对正投影仿真CT图像F(x,y)进行阈值分割，取灰度范围在[0,3]内的图像区域S(x,y)；

提取实际CT图像Q(x,y)中除去被测对象及受伪影强干扰的背景区域S(x,y)的背景区域G(x,y)，对其对应的灰度直方图统计如图8所示，该灰度直方图中最大峰值像素灰度值g＝95；

提取实际CT图像Q(x,y)中S(x,y)对应的像素灰度值H(x,y)，该区域包含被测对象及受伪影强干扰的背景区域，如图9所示，其对应的灰度直方图h₁(Z)如图10所示；

对该灰度直方图h₁(Z)采用最大类间方差法进行自动阈值分割，确定分割值T₁＝136将前景与受伪影强干扰的背景进行分割。以T₁＝136为分界线，计算灰度直方图在[0,136]范围内的最大峰值像素灰度值w＝114，即为受伪影强干扰的背景区域的像素最大峰值灰度值；计算h₁(Z)在[137,255]范围内的最大峰值像素灰度值q＝165，即为被测对象区域的像素最大峰值灰度值；

针对H(x,y)的灰度直方图h₁(Z)在[0,136]范围内的灰度直方图进行最小残差法高斯拟合，如图11所示，离散点为h₁(Z)在[0,136]范围内真实的统计分布，曲线为高斯拟合结果；

计算获得最佳拟合高斯曲线灰度值均值

以及最佳灰度值标准差σ_opt＝9.599；

计算拟合后的高斯曲线函数f₁(i)的概率密度函数

其归一化概率密度函数

对受伪影强干扰的背景区域的像素灰度值进行校正，设受伪影强干扰的任意一点背景像素灰度值为k，校正后像素灰度值的转换公式为

对实际CT图像Q(x,y)每个像素点进行局部灰度统计分析，设任意像素点为(a,b)，分析范围为以该像素点为中心的大小为3×3的像素矩阵，取该区域占比像素灰度值最大的灰度V(a,b)，计算V(a,b)对应的归一化概率密度函数PDF_nor(V(a,b))；

对实际CT图像Q(x,y)进行各个像素灰度值校正处理，校正后图像的处理公式为：

实际CT图像Q(x,y)的校正后图像Q_cal(x,y)如图12所示。

参阅图13至图16，分别为实际CT图像阈值分割效果图、校正后CT图像阈值分割效果图、实际CT局部的图像阈值分割效果图、校正后CT局部的图像阈值分割效果图，采用相同的阈值对实际CT图像和灰度校正后CT图像进行边缘识别，从局部放大效果图来看，减少了伪影对边缘分割的干扰。

以直径为2mm的不锈钢球阵列为例，对其原始CT图像集和经过正投影伪影影响区统计的灰度校正后的CT图像集均采用固定阈值分割法，查看每个球体的直径分布情况，如表2所示，校正前最大直径为2.1mm，这是因为一部分伪影被计入不锈钢球体直径内造成最大直径被放大。

比较校正前后球体直径的平均值，可以发现校正后更接近于2mm直径，且校正前最大球体直径测量值为2.1mm，校正后该最大直径为2.0mm更接近实际尺寸，从所有球体直径偏差也可以看出从0.26mm降低至0.14mm。

本优选实施例提供的工业CT伪影校正方法，针对单一材料复杂内部结构工件的CT扫描时，能有效降低伪影引起的灰度差异，实现精确、可靠、无损测量其内部结构尺寸。

需要说明的是，本技术领域的普通技术人员应当认识到，以上的实施方式仅是用来说明本发明，而并非用作为对本发明的限定，只要在本发明的实质精神范围之内，对以上实施方式所作的适当改变和变化都落在本发明要求保护的范围之内。

Claims (6)

1.一种工业CT伪影校正方法，其特征在于，包括以下步骤：

(1)基于正投影法获得被测对象的伪影场F(x,y)并对其进行阈值分割，设所述伪影场F(x,y)中灰度值在[a,b]范围内的图像区域为S(x,y)，该区域包含被测对象及受伪影强干扰的背景区域；

(2)实际CT图像为Q(x,y)，提取所述实际CT图像Q(x,y)中除去S(x,y)的背景区域，得到伪影弱干扰背景图像G(x,y)，对G(x,y)进行灰度直方图h₂(Z)统计，并计算h₂(Z)的最大峰值的像素灰度值g；

(3)提取实际CT图像Q(x,y)中包括S(x,y)的背景区域，得到前景及伪影强干扰背景图像H(x,y)，对H(x,y)进行灰度直方图h₁(Z)统计，并通过最大类间方差法获取将H(x,y)中前景与伪影强干扰背景图像进行分割的阈值T，对h₁(Z)在[0,T]范围内的灰度直方图h₁(Z_i)进行最小残差法高斯拟合，获取最佳拟合高斯曲线f₁(i)；

(4)计算步骤(3)中f₁(i)的概率密度函数PDF(x)，归一化处理后得到PDF_nor(x)；

(5)对实际CT图像Q(x,y)的每个像素点(a,b)进行局部灰度统计分析，获取该像素点中占比像素灰度值最大的灰度V(a,b)，并计算该灰度V(a,b)对应的归一化的概率密度函数PDF_nor(V(a,b))；

(6)计算h₁(Z)在[0,T]范围内的最大峰值像素灰度值w，即为受伪影强干扰的背景区域的像素最大峰值灰度值；计算h₁(Z)在[T+1,L-1]范围内的最大峰值像素灰度值q，即为被测对象区域的像素最大峰值灰度值；其中，L-1为最大像素灰度值；

(7)对受伪影强干扰的背景区域的像素灰度值k进行校正，校正后的像素灰度值

(8)对Q(x,y)各个像素灰度值进行校正处理，校正后的CT图像

2.根据权利要求1所述的工业CT伪影校正方法，其特征在于，步骤(1)具体包括：设被测对象为T(x,y)，其中0≤x≤m，0≤y≤n，

其中，

为被测对象内部及表面区域，图像中任意一点(a,b)在任意角度θ的线积分K_θ(a,b)为：

设伪影正投影数组为f(x,y)，利用K_θ(a,b)进行周向累积计算，

将f(x,y)进行灰度级映射成像F(x,y)，

伪影场

其中，Z表示所成图像的像素灰度级，S(x,y)＝{F(x_i,y_i)∈[a,b]}。

3.根据权利要求1所述的工业CT伪影校正方法，其特征在于，步骤(2)中，G(x,y)＝Q(x,y)-Q(x,y)∩S(x,y)，对G(x,y)进行灰度直方图h₂(Z)统计，

其中，

Z表示所成图像的像素灰度值。

4.根据权利要求1所述的工业CT伪影校正方法，其特征在于，步骤(3)中，H(x,y)＝Q(x,y)∩S(x,y)，对H(x,y)进行灰度直方图h₁(Z)统计，

通过最大类间方差法获取将H(x,y)中前景与伪影强干扰背景图像进行分割的阈值T包括：设任意分割值T_i，图像中像素灰度值小于阈值T_i的像素个数

其平均像素灰度值

图像中像素灰度值大于阈值T_i的像素个数

其平均像素灰度值

类间方差

采用遍历的方法获得类间方差g最大时的阈值T，T＝max(g)；

对h₁(Z)在[0,T]范围内的灰度直方图h₁(Z_i)进行最小残差法高斯拟合，拟合后的高斯曲线函数

其中，a是幅值，

是灰度值均值，σ是灰度值标准差，高斯拟合曲线迭代的初始灰度值均值

初始灰度值标准差

将

和σ_init代入公式

中，得到

其中，幅值a_init的初始值为

则得到幅值的初始值

其中，灰度值均值

幅值

灰度值标准差

选取不同的拟合后高斯曲线的灰度值均值

幅值a以及灰度值标准差σ，建立函数类Ψ＝{f₀(i),f₁(i),...,f_m(i)}(m<L)，求出使

时，最佳拟合高斯曲线灰度值均值

以及最佳灰度值标准差σ_opt。

5.根据权利要求1所述的工业CT伪影校正方法，其特征在于，步骤(4)中，拟合后的高斯曲线函数f₁(i)的概率密度函数

其归一化概率密度函数

6.根据权利要求1所述的工业CT伪影校正方法，其特征在于，步骤(5)中，对实际CT图像Q(x,y)的每个像素点(a,b)进行局部灰度统计分析的像素点范围为大小为3×3的像素矩阵。

Images