一种工业CT散射校正方法与流程

一种工业ct散射校正方法
技术领域
1.本发明涉及ct图像校正技术领域，特别涉及一种工业ct散射校正方法。


背景技术：

2.工业ct检测技术是一种在x射线透射成像技术上发展起来的实用性无损检测手段，适用于复杂结构产品及零部件内部结构及缺陷的检测，具有成像直观，定量、定位、定性准确等优势，广泛应用于工业无损探伤、医疗卫生等领域。
3.x射线与物质发生关系主要是：透射和散射两种。其中透射是指x射线穿透物质时其能力、强度发生了变化，产生了射线衰减；射线的衰减规律基本上符合lamber-beer定律，即射线在很小的厚度范围内，其强度的衰减量正比于入射射线强度和穿透物体的厚度，工业ct既是基于该定律的基础上结合卷积反投影等算法实现内部结构及缺陷的可视化；散射是指x射线在射入物体后，与物质原子的相互作用，包括：光电效应、康普顿效应、电子对效应和瑞利散射，其主要特征就是x射线偏离了入射方向。
4.在工业ct扫描过程中，由于x射线接收用探测器具有较高的射线敏感性，可以同时接收到透射和散射x射线。但是在ct重建算法中忽略了散射线的存在，这就造成了ct图像中存在散射伪影，在ct图像上表现为阴影(或条纹)，这对内部缺陷的检测造成了严重的干扰。因此，研究面阵工业ct散射校正方法非常必要的。
5.现阶段，工业ct的散射校正大致可分为硬件校正和软件校正两大类。硬件校正是在x射线成像系统的各个仪器组件上添加一些校正工具，通过物理上阻挡散射线被探测器接收到，以达到散射校正的目的。比如：在探测器前面加装准直器，阻挡通过物体产生的散射x射线进入探测器内以及经典的空气隙(air gap)方法，即增加被照物体与探测器的距离,也能够减小散射成分；软件校正方法是指对已经得到的x射线投影图像，采用数字图像处理的方法，在计算机中通过对图像本身的分析和对被照物体的性质的估计，得出一个散射分布图，来进行散射校正。其中包括卷积法、反卷积法、蒙特卡罗模拟法、模型估计法等等。但是软件校正方法算法复杂、计算量大、效率较低。近年来，国内外相关研究人员采用散射校正板法取得了较好的效果，通过二维铅块点阵列(二维孔阵列铅板)置于被检工件与探测器之前，通过扫描2次，获得物体+散射校正板图像a、单独物体图像b，利用图像a二维插值运算得到散射强度分布图，再用图像b减去散射强度分布图得到散射校正后图像，该方法(二维孔阵列)要求每个孔指向焦点，这就使得每个孔的角度方向均不一样，造成加工难度及成本极高。难以实现广泛的应用普及。


技术实现要素：

6.本发明所要解决的技术问题是针对上述现有技术，而提供一种成本低的工业ct散射校正方法。
7.本发明解决上述技术问题所采用的技术方案为：一种工业ct散射校正方法，其特征在于包括如下步骤：
8.步骤1、加工点阵板，并在点阵板上加工垂直于板面的通孔阵列；
9.其中通孔阵列包括m行n列等间距分布的通孔，m和n均为正整数；
10.步骤2、通过工业ct系统对被检测样品进行周向扫描，获得s幅被检测样品的第一周向dr图像，将第s幅第一周向dr图像中像素点为(x,y)的像素值记为as(x,y)，s∈1、2
…
s；上述的工业ct系统包括依次间隔布置的x射线机、用于供被检测样品放置的置物台和探测器；
11.步骤3、在工业ct系统的探测器前面竖向放置步骤1中的点阵板，采用与步骤2中相同的扫描工艺对被检测样品进行周向扫描，获得s幅被检测样品的第二周向dr图像，将第s幅第二周向dr图像中像素点为(x,y)的像素值记为bs(x,y)，bs(x,y)和as(x,y)的扫描角度相同；
12.步骤4、采用与步骤2中相同的扫描工艺对点阵板进行dr扫描，获得一幅点阵板投影图像；
13.步骤5、计算散射线分布图像；
14.散射线分布图像计算的具体步骤为：
15.步骤5-1、将点阵板投影图像进行二值化处理，得到点阵板投影二值化图像；
16.步骤5-2、对点阵板上的每个通孔进行编号，其中i和j分别对应为点阵板上行和列所对应的序号，i∈{1、2、
…
m}，j∈{1、2、
…
n}；
17.步骤5-3、搜索点阵板上每一个通孔在点阵板投影二值化图像中所占的像素位置，并将每个通孔所占的像素位置归为第一子类ck
i,j
，计算每个通孔在点阵板投影二值化图像上成像后的质心位置；
18.步骤5-4、计算每个通孔的质心位置对应在第二周向dr图像和第一周向dr图像中的图像灰度值，分别记为bs(x,y)
i,j
和as(x,y)
i,j
；
19.步骤5-5、根据点阵板投影二值化图像的像素尺寸v
cal
以及点阵板的通孔孔径d计算出通孔孔径占据的像素数量
20.步骤5-6、获取以每个通孔的质心位置为圆心，为半径的圆，并将在该圆范围内的像素位置归为第二子类cy
i,j
，找出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点，在第二周向dr图像中提取出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点所对应位置的灰度值，计算所有点的灰度值均值，设为k
si,j
；
21.步骤5-7、计算每一幅第二周向dr图像中各个通孔质心位置的散射值l
si,j
；l
si,j
的计算公式为：
22.l
si,j
＝as(x,y)
i,j-bs(x,y)
i,j-k
si,j
；
23.步骤5-8、根据每一幅第二周向dr图像中各个通孔质心位置的散射值，采用二维插值法计算出每幅第一周向dr图像的散射线分布图像，将第s幅第一周向dr图像的散射线分布图像中像素点为(x,y)的像素值记为es(x,y)；
24.步骤6、对每一幅第一周向dr图像进行散射线校正，即获得s幅散射校正后的图像；其中第s幅散射校正后的图像中像素点为(x,y)的像素值为as(x,y)-es(x,y)；
25.步骤7、采用s幅散射线校正的图像进行ct重建，获得校正后的ct图像。
26.本方案中，所述步骤1中点阵板厚度是根据工业ct系统中x射线机最高能量而定。
27.优选地，所述点阵板材料为钢；
28.当x射线机最高能量为50kev-150kev时，点阵板厚度为5mm；
29.当x射线机最高能量为150kev-250kev时，点阵板厚度为10mm；
30.当x射线机最高能量为250kev-350kev时，点阵板厚度为15mm；
31.当x射线机最高能量为350kev-450kev时，点阵板厚度为20mm。
32.进一步地，所述点阵板材料为除钢以外的其他材料时，则点阵板厚度d
x
计算公式为：
[0033][0034]
其中，d1表示点阵板材料为钢时的板厚；ρ1表示钢的密度；ρ
x
表示其他材料的密度。
[0035]
进一步的，所述步骤3中的扫描工艺包括管电压、管电流、放大比、焦点尺寸和积分时间。
[0036]
本方案中，所述步骤5-3中每个通孔在点阵板投影二值化图像上成像后的质心位置计算方式为：计算通孔在点阵板投影二值化图像中所占像素位置的均值，则该均值即为通孔的质心位置。
[0037]
与现有技术相比，本发明的优点在于：传统方法必须采用放射形状的孔阵列板，加工难度大、成本高、时间长，而本发明中在点阵板上加工垂直于板面的通孔阵列，解决了上述问题，并通过图像处理的方式，在散射线分布图像计算时计算各个通孔质心位置的散射值，以实现散射校正的目的。因此本方法成本低，容易实现且由于散射校正精度高。
附图说明
[0038]
图1为本发明实施例中点阵板的示意图；
[0039]
图2为本发明实施例中获得某一角度的被检测样品的第一周向dr图像；
[0040]
图3为与图2中同一角度的被检测样本的第二周向dr图像；
[0041]
图4为本发明实施例中点阵板投影图像；
[0042]
图5为本发明实施例中点阵板投影二值化图像；
[0043]
图6为本发明实施例中图4的局部放大图；
[0044]
图7为本发明实施例中属于第二子类但不属于第一子类的点示意图；
[0045]
图8为本发明实施例中散射校正前后的图像；其中(a)为散射校正前的图像，(b)为散射校正后的图像。
具体实施方式
[0046]
以下结合附图实施例对本发明作进一步详细描述。
[0047]
本实施例中的工业ct散射校正方法包括如下步骤：
[0048]
步骤1、加工点阵板，并在点阵板上加工垂直于板面的通孔阵列；
[0049]
其中通孔阵列包括m行n列等间距分布的通孔，m和n均为正整数；点阵板的通孔孔径d越小越好，一般取1-3mm(直径)；通孔孔距(两通孔的圆心距)范围为通孔孔径的1.5-3倍之间，纵向及横向孔距保持一致，等间距分布；采用机械加工等方式制造垂直于板面的孔阵列；
[0050]
本实施例中点阵板厚度根据工业ct系统中x射线机最高能量(管电压)而定，当点阵板材料选取为钢；
[0051]
当x射线机最高能量为50kev-150kev时，点阵板厚度为5mm；
[0052]
当x射线机最高能量为150kev-250kev时，点阵板厚度为10mm；
[0053]
当x射线机最高能量为250kev-350kev时，点阵板厚度为15mm；
[0054]
当x射线机最高能量为350kev-450kev时，点阵板厚度为20mm；
[0055]
本实施例中，实验用的工业ct系统中x射线机最高能量(管电压)为90kv，采用不锈钢制作厚度为5mm的点阵板，通孔的直径为3mm，通孔中心间距为9mm，点阵板的具体实物图如图1所示；
[0056]
点阵板材料为除钢以外的其他材料时，则点阵板厚度d
x
计算公式为：
[0057][0058]
其中，d1表示点阵板材料为钢时的板厚；ρ1表示钢的密度；ρ
x
表示其他材料的密度；点阵板的具体材料根据实际需要更换；
[0059]
步骤2、通过工业ct系统对被检测样品进行周向扫描，获得s幅被检测样品的第一周向dr图像，将第s幅第一周向dr图像中像素点为(x,y)的像素值记为as(x,y)，s∈1、2
…
s；上述的工业ct系统包括依次布置的x射线机、用于供被检测样品放置的置物台和探测器；工业ct系统的具体成像方式为本领域技术人员所公知的现有技术，在此不再展开赘述；
[0060]
其中在周向扫描前根据被检测对象材料、密度、尺寸等特征调整检测工艺参数；获得某一角度的被检测样品的第一周向dr图像如图2所示；
[0061]
步骤3、在工业ct系统的探测器前面竖向放置步骤1中的点阵板，采用与步骤2中相同的扫描工艺对被检测样品进行周向扫描，获得s幅被检测样品的第二周向dr图像，将第s幅第二周向dr图像中像素点为(x,y)的像素值记为bs(x,y)，bs(x,y)和as(x,y)的扫描角度相同；
[0062]
点阵板与探测器之间保持平行，点阵板尽量靠近探测器放置，但对于其与探测器之间的距离并没有特别要求；本实施例中扫描工艺包括管电压、管电流、放大比、焦点尺寸和积分时间；获取与步骤2中相同角度的第二周向dr图像，如图3所示；
[0063]
步骤4、采用与步骤2中相同的扫描工艺对点阵板进行dr扫描，获得一幅点阵板投影图像，将点阵板投影图像中像素点为(x,y)的像素值记为d(x,y)；点阵板投影图像如图4所示；
[0064]
步骤5、计算散射线分布图像；
[0065]
散射线分布图像计算的具体步骤为：
[0066]
步骤5-1、将点阵板投影图像进行二值化处理，得到点阵板投影二值化图像；
[0067]
图像二值化的阈值选择，可采用自动或手动设置方法，其中自动方法优选最大类间方差法，获得的点阵板投影二值化图像如图5所示；
[0068]
步骤5-2、对点阵板上的每个通孔进行编号，其中i和j分别对应为点阵板上行和列所对应的序号，i∈{1、2、
…
m}，j∈{1、2、
…
n}；
[0069]
步骤5-3、搜索点阵板上每一个通孔在点阵板投影二值化图像中所占的像素位置，并将每个通孔所占的像素位置归为第一子类ck
i,j
，计算每个通孔在点阵板投影二值化图像
上成像后的质心位置；
[0070]
其中每个通孔在点阵板投影二值化图像上成像后的质心位置计算方式为：计算通孔在点阵板投影二值化图像中所占像素位置的均值，则该均值即为通孔的质心位置；
[0071]
步骤5-4、计算每个通孔的质心位置对应在第二周向dr图像和第一周向dr图像中的图像灰度值，分别记为bs(x,y)
i,j
和as(x,y)
i,j
；
[0072]
步骤5-5、根据点阵板投影二值化图像的像素尺寸v
cal
以及点阵板的通孔孔径d计算出通孔孔径占据的像素数量
[0073]
本实施例中，像素尺寸v
cal
＝0.15mm，点阵板孔径d＝3mm，计算出孔径占据的像素数量为20个；
[0074]
步骤5-6、获取以每个通孔的质心位置为圆心，为半径的圆，并将在该圆范围内的像素位置归为第二子类cy
i,j
，找出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点，在第二周向dr图像中提取出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点所对应位置的灰度值，计算所有点的灰度值均值，设为k
si,j
；
[0075]
如图6所示为点阵板投影图像的局部放大图，通过该图6可以看出：当通孔处于边缘位置时，通孔会发生变形，从而形成如图6中所示的斜孔；为此，通过步骤5-6中的计算方式找出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点(如图7中所示的网格区域)，计算网格区域中所有点在第二周向dr图像中对应位置的灰度值均值，这样能有效提高散射校正的准确性；
[0076]
步骤5-7、计算每一幅第二周向dr图像中各个通孔质心位置的散射值l
si,j
；l
si,j
的计算公式为：
[0077]
l
si,j
＝as(x,y)
i,j-bs(x,y)
i,j-k
si,j
；
[0078]
步骤5-8、根据每一幅第二周向dr图像中各个通孔质心位置的散射值，采用二维插值法计算出每幅第一周向dr图像的散射线分布图像，将第s幅第一周向dr图像的散射线分布图像中像素点为(x,y)的像素值记为es(x,y)；
[0079]
上述的二维插值法为现有技术，在此不再展开；
[0080]
步骤6、对每一幅第一周向dr图像进行散射线校正，即获得s幅散射校正后的图像；其中第s幅散射校正后的图像中像素点为(x,y)的像素值为as(x,y)-es(x,y)；
[0081]
步骤7、采用s幅散射线校正的图像进行ct重建，获得校正后的ct图像。
[0082]
步骤7中的ct重建方式为现有技术，在此不再展开赘述。另外如图8(a)为散射校正前的效果图；如图8(b)为散射校正后的效果图。
[0083]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

技术特征：
1.一种工业ct散射校正方法，其特征在于包括如下步骤：步骤1、加工点阵板，并在点阵板上加工垂直于板面的通孔阵列；其中通孔阵列包括m行n列等间距分布的通孔，m和n均为正整数；步骤2、通过工业ct系统对被检测样品进行周向扫描，获得s幅被检测样品的第一周向dr图像，将第s幅第一周向dr图像中像素点为(x,y)的像素值记为a
s
(x,y)，s∈1、2
…
s；上述的工业ct系统包括依次间隔布置的x射线机、用于供被检测样品放置的置物台和探测器；步骤3、在工业ct系统的探测器前面竖向放置步骤1中的点阵板，采用与步骤2中相同的扫描工艺对被检测样品进行周向扫描，获得s幅被检测样品的第二周向dr图像，将第s幅第二周向dr图像中像素点为(x,y)的像素值记为b
s
(x,y)，b
s
(x,y)和a
s
(x,y)的扫描角度相同；步骤4、采用与步骤2中相同的扫描工艺对点阵板进行dr扫描，获得一幅点阵板投影图像；步骤5、计算散射线分布图像；散射线分布图像计算的具体步骤为：步骤5-1、将点阵板投影图像进行二值化处理，得到点阵板投影二值化图像；步骤5-2、对点阵板上的每个通孔进行编号，其中i和j分别对应为点阵板上行和列所对应的序号，i∈{1、2、
…
m}，j∈{1、2、
…
n}；步骤5-3、搜索点阵板上每一个通孔在点阵板投影二值化图像中所占的像素位置，并将每个通孔所占的像素位置归为第一子类ck
i,j
，计算每个通孔在点阵板投影二值化图像上成像后的质心位置；步骤5-4、计算每个通孔的质心位置对应在第二周向dr图像和第一周向dr图像中的图像灰度值，分别记为b
s
(x,y)
i,j
和a
s
(x,y)
i,j
；步骤5-5、根据点阵板投影二值化图像的像素尺寸v
cal
以及点阵板的通孔孔径d计算出通孔孔径占据的像素数量步骤5-6、获取以每个通孔的质心位置为圆心，为半径的圆，并将在该圆范围内的像素位置归为第二子类cy
i,j
，找出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点，在第二周向dr图像中提取出属于第二子类cy
i,j
但不属于第一子类ck
i,j
的点所对应位置的灰度值，计算所有点的灰度值均值，设为k
si,j
；步骤5-7、计算每一幅第二周向dr图像中各个通孔质心位置的散射值l
si,j
；l
si,j
的计算公式为：l
si,j
＝a
s
(x,y)
i,j-b
s
(x,y)
i,j-k
si,j
；步骤5-8、根据每一幅第二周向dr图像中各个通孔质心位置的散射值，采用二维插值法计算出每幅第一周向dr图像的散射线分布图像，将第s幅第一周向dr图像的散射线分布图像中像素点为(x,y)的像素值记为e
s
(x,y)；步骤6、对每一幅第一周向dr图像进行散射线校正，即获得s幅散射校正后的图像；其中第s幅散射校正后的图像中像素点为(x,y)的像素值为a
s
(x,y)-e
s
(x,y)；步骤7、采用s幅散射线校正的图像进行ct重建，获得校正后的ct图像。2.根据权利要求1所述的工业ct散射校正方法，其特征在于：所述步骤1中点阵板厚度
是根据工业ct系统中x射线机最高能量而定。3.根据权利要求2所述的工业ct散射校正方法，其特征在于：所述点阵板材料为钢；当x射线机最高能量为50kev-150kev时，点阵板厚度为5mm；当x射线机最高能量为150kev-250kev时，点阵板厚度为10mm；当x射线机最高能量为250kev-350kev时，点阵板厚度为15mm；当x射线机最高能量为350kev-450kev时，点阵板厚度为20mm。4.根据权利要求3所述的工业ct散射校正方法，其特征在于：所述点阵板材料为除钢以外的其他材料时，则点阵板厚度d
x
计算公式为：其中，d1表示点阵板材料为钢时的板厚；ρ1表示钢的密度；ρ
x
表示其他材料的密度。5.根据权利要求1～4任一项所述的工业ct散射校正方法，其特征在于：所述步骤3中的扫描工艺包括管电压、管电流、放大比、焦点尺寸和积分时间。6.根据权利要求5所述的工业ct散射校正方法，其特征在于：所述步骤5-3中每个通孔在点阵板投影二值化图像上成像后的质心位置计算方式为：计算通孔在点阵板投影二值化图像中所占像素位置的均值，则该均值即为通孔的质心位置。

技术总结
本发明涉及一种工业CT散射校正方法，包括：加工点阵板，并在点阵板上加工垂直于板面的通孔阵列；通过工业CT系统对被检测样品进行周向扫描，获得S幅被检测样品的第一周向DR图像；在工业CT系统的探测器前面竖向放置点阵板，获得S幅被检测样品的第二周向DR图像；接着对点阵板进行DR扫描，获得点阵板投影图像；并根据点阵板投影图像计算散射线分布图像，并根据散射线分布图像对每一幅第一周向DR图像进行散射线校正，最后采用散射线校正的图像进行CT重建，获得校正后的CT图像。该方法成本低，容易实现。易实现。易实现。


技术研发人员：赵亮
受保护的技术使用者：奥影检测科技（上海）有限公司
技术研发日：2023.02.22
技术公布日：2023/7/13