工业CT图像中弱边缘尺寸测量方法与流程
本发明涉及一种工业CT图像中弱边缘尺寸测量方法。
背景技术：
工业计算机断层成像(Computed Tomography，CT)是通过采集不同角度下的X射线投影重建物体CT图像，其优势体现在可以无损地测量物体的外部和内部结构及缺陷的尺寸，有效弥补传统测量方法只能测量物体外表面结构、解剖观测内部缺陷尺寸等不足。但是，为了保证CT系统可以穿透一定厚度的物体，射线源的功率需要足够大，随之而来地焦斑尺寸也需要足够大，这就会导致重建后的CT图像物体边界模糊退化，清晰度降低；在工程应用中，缺陷种类十分复杂，其缺陷类型(如铸件气孔、疏松等)就导致边界不清晰；或者因为缺陷表面与CT扫描平面不垂直等因素，加上CT图像的容积效应，造成CT图像中缺陷边缘有较宽的灰度过渡区，形成了渐变边缘或称弱边缘；这就影响到工件内外部尺寸测量和内部缺陷测量的精度。现阶段，基于工业CT图像的测量方法基本上以检测人员手工测量为主，大量采用传统的“半高宽”法，这样获得的缺陷尺寸与解剖验证实验所获得的结果往往是脱节的，这种方式费时费力且造成缺陷尺寸测定存在明显的“放大效应”。既大大降低了工业CT缺陷测量的有效性，又对工业CT拓展其应用领域设置了严重的阻碍。
针对上述问题，经过对现有技术的文献检索发现，重庆大学曾理教授对工业CT缺陷自动识别及检测开展过研究工作，其运用Facet模型与分形维数相结合的方法对目标进行分割，利用分形维数自动定位出缺陷大致位置，再用Facet模型进行边缘检测，以降低该模型时间复杂度；该方法应用到铁路货车摇枕和侧架的工业CT检测中，实验结果表明能够获得较为准确的缺陷位置和形状，自动化程度较高；但是该方法对模拟图像的分割效果不理想，有待于进一步研究。
清华大学的张丽研究员等人对基于工业CT的工件缺陷快速定位和识别方法进行了研究，其对工业CT图像处理思路和方法做了阐述，采用统计滤波对CT图像进行噪声处理，利用level set方法(用于缺陷跟踪)和Gabor小波方法(用于缺陷特征提取)来实现对缺陷的识别，指出两种方法在CT缺陷识别领域中比较有潜力，但并未就上述两种方法的实现过程进行简要阐释，也未对方法所能实现效果加以描述。
在缺陷快速定位方面，有人提出用覆盖(Blanket)算法计算分形面面积，并用最小二乘法进行直线拟合得到局部盒维数的方法实现定位，该方法简单省时，尤其对于结构复杂缺陷具有明显优势，缺点对于内部灰度均匀的缺陷应用上受到限制，需要与传统分割方法结合使用，在过渡区域缺陷存在算法选择误区。
南华大学的李林升对工业CT铸件缺陷开展了自动识别技术的研究工作，其提出在工业CT扫描图像上采用划分区域和重复运用C-V模型的方法实现对多灰度级缺陷进行目标识别，并在实践中尝试；该方法不足之处在于当灰度级较多时，重复进行C-V模型的计算，非常耗时。
山东大学的孙同景教授对工业CT图像自适应缺陷识别开展了研究工作，其针对X射线投影图像特点，提出在Kalman滤波算法新息正交性的基础上，进行图像数据自适应补偿，定义了灰度动态评价函数和相关阈值，并通过实验验证了方法的有效性；然而该方法只是在工业CT投影(DR)图像上实验验证，并未在工业CT图像上进行有效性研究。
重庆大学王珏教授所带领的团队就此进行了多方面的探索性研究，其提出采用形态学开-闭重建滤波，抹去图像上小于结构元素的细节，抑制图像噪声保证原有的特征，采用C-V方法或分水岭方法进行图像分割，再进行最小二乘法拟合得到了较好的结果；不足之处在于单纯使用形态学开-闭重建运算对噪声进行滤波，没有从根本上解决“过分割”问题，局限性较大。
大量实验研究均是基于规则形状的工件边缘，对于不规则缺陷引起的弱边缘只是从面积上进行了计算，对边缘定位精度没有准确性分析，仅处于理论研究阶段。
中北大学的潘晋孝教授采用边缘退化模型来对CT图像进行边缘恢复，实验结果显示该方法较常用图像尺寸测量方法相对误差减少了1.5％；但是研究人员重点研究有规则边缘的CT图像，对不规则边缘或缺陷边缘未作进一步的研究。
福州大学的王美清教授针对图像中弱边缘分割技术开展了研究，提出改进LI的保持距离水平集方法，提出自适应分割弱边缘的活动轮廓模型，在模型中将图像灰度均值加入到自适应力的系数中，并采用相似性系数和错误分割率评价对算法进行有效性评估；实验结果表明该方法能够分割多目标图像，对常规噪声(高斯噪声、椒盐噪声)具有较强的抗噪性，但未就工业CT图像特有噪声开展相关研究。
暨南大学的蔡利栋教授对图像中弱边缘检测技术开展研究工作，提出一种基于非线性灰度变换的弱边缘检测方法，其原理首先对图像灰度作非线性变换，平滑后计算梯度值并进行非最大化抑制，最后利用梯度直方图选取合适的梯度阈值来标识边缘点；通过对医用CT图像进行实验分析，在弱边缘检出能力上较常规算法有明显提高；由于医用CT与工业CT从射线源能量强度等系统性能上就有显著不同，导致所形成的CT图像质量也有较大差异，所以该算法在工业CT图像上使用效果有所抑制。
技术实现要素：
本发明所要解决的技术问题是针对由于工件或缺陷表面与CT扫描平面不垂直，加上CT图像的容积效应导致图像中边界模糊退化，清晰度降低使得工业CT图像测量误差大等问题，提供一种速度快、精度高、可用于人工测量的工业CT图像中弱边缘尺寸测量方法。
本发明解决上述技术问题所采用的技术方案为：一种工业CT图像中弱边缘尺寸测量方法，其特征在于：包含如下步骤：
步骤一、分别采集被测物体和标准试块的CT图像，其中标准试块的密度和厚度与被测物体一致；
步骤二、获取标准试块的CT图像垂直界面一维点扩散函数：
所述标准试块采用圆形标准试块或方形标准试块；
当标准试块为圆形标准试块时，采用GB/T 29069-2012《无损检测工业计算机层析成像(CT)系统性能测试方法》中5.3圆盘卡法获得圆形标准试块CT图像垂直界面一维点扩散函数PSF[k]，求取PSF[k]的累计幅值A；
当标准试块为方形标准试块时，在标准方形试块的CT图像中框选一个长方形区域，使得标准方形试块的一条直线边缘图像包含在该长方形区域中，对该区域进行二值化分割，获取该直线边缘点，再根据直线边缘点进行拟合获得该直线边缘的斜率k，计算与该直线边缘垂直直线的斜率k′＝1/k，根据该垂直直线的斜率在标准试块边缘上均匀提取N条相隔距离为h的与该直线边缘相垂直的剖线：yn＝k′xn+Bn，Bn＝d+nh,n∈[1,N]，d为起始斜截距；根据每条剖线yn在标准方形试块的CT图像上经过的位置，提取与之对应的直线边缘点的对应灰度值ln，获取所有直线边缘上所有像素的灰度加以合并平均，得到边缘响应函数ERF，然后对边缘响应函数ERF进行拟合求导获取方形标准试块CT图像垂直界面一维点扩散函数PSF[k]，具体步骤如下：
边缘响应函数ERF是长度为N的一维数组，设为ERFN，设E[m]为ERFN的子数组，其中a为数组起点，b为数组终点。0≤a<b≤N，其中数组长度为m＝b-a，m∈[15～25]；对E[m]进行3次方最小二乘法拟合，获得拟合后的E[m]′，取中点值E[a+m/2]′作为该段数组拟合后新的Enew[m]值，再对Enew[m]进行3次方最小二乘法拟合，获得拟合后的Enew[m]′，对Enew[m]′进行求导dEnew[m]′，取中点值dEnew[a+m/2]′作为该段数组的PSF值，其中W是归一化系数，W＝max(dEnew[a+m/2]′)；
该标准试块的CT图像垂直界面一维点扩散函数为PSF[s]，对PSF[s]进行累加获得累计幅值
步骤三、获取被测物体待测尺寸的灰度曲线：在被测物体的CT图像中画一条通过待测尺寸的直线，使得该直线包含待测尺寸的两端边界，获取该直线上对应的灰度值l[n]；
步骤四、对被测物体待测尺寸的两端边界分别进行复原重构：
对步骤三获得的灰度直线上对应的灰度值l[n]以两端边界中心位置为界进行分割分段处理，分别获得2组CT图像垂直界面一维点扩散函数，它分别对应被测物体待测尺寸的两端边界CT图像垂直界面一维点扩散函数，假设左边界对应的CT图像垂直界面一维点扩散函数为PSF[L]，右边界对应的CT图像垂直界面一维点扩散函数为PSF[R]；设h[L]、h[R]分别为PSF[L]和PSF[R]对标准试块的CT图像垂直界面一维点扩散函数PSF[s]作反卷积后得到的复原边界系数，为左侧复原边界系数，为右侧复原边界系数；
在这里，对被测物体待测尺寸的左右两端边界类型进行分类：
1、垂直于CT扫描平面的垂直平面边界；
2、倾斜于CT扫描平面的斜平面边界；
3、不规则面的边界；
然后对获取的h数据进行分类，这里的h数据为h[L]或h[R]：
1、如果h数据中有且仅有一个最大值max[h]，并且max[h]附近数据急剧下降，h数据曲线呈山峰状，即num{max[h]}＝1，则判定为第1种边界；
2、如果h数据显示为方波，即num{max[h]}为多个且连续；则判定为第2种边界；
3、如果h数据显示为曲线，则判定为第3种边界；
对步骤三获得的灰度值l[n]求导，有两个极值分别对应待测尺寸的左右边界的位置，设为j、k，j表示左边界位置，k表示右边界位置；取j、k中间点，即将步骤三获得的灰度值l[n]在这点上一分为二，分开进行边界复原处理：
设左边段的灰度值为L(w)，假定一个长度估计值x0，x0初值小于测量长度，其退化前的幅值为其中A为步骤二中求取的标准试块的CT图像垂直界面一维点扩散函数为PSF[s]的累计幅值；然后根据判断出的左侧边界类型分别进入如下处理：
1、当左侧边界类型判断为1，即垂直平面边界，初始值x0不变，然后计算Xold(L)：
2、当边界类型判断为2，即斜平面边界，计算左侧复原边界系数h[L]中大于的点数R，即计算方波的宽度，然后计算Xold(L)：
3、当边界类型判断为3，即不规则面的边界，设左侧复原边界系数h[L]中max{h[L]}对应的位置为y，取其后的所有数据设为h[y]，然后计算Xold(L)：
将Xold(L)与步骤二获得的标准试块的CT图像垂直界面一维点扩散函数为PSF[s]进行卷积，计算与左边段的灰度值为L(w)的标准偏差S(w)；然后将x0＝x0+1，循环迭代计算，直至最终回计算时的x0值大于测量长度；获得S(w)的曲线，取数组S(w)中的最小值对应的x0，此时的x0即为最佳长度值，然后将最佳长度值x0代入Xold(L)中，获得左侧最佳重构曲线Xnew(L)；
右侧最佳重构曲线Xnew(R)与左侧最佳重构曲线Xnew(L)的获取方式相同，将右侧最佳重构曲线Xnew(R)与左侧最佳重构曲线Xnew(L)合并，得到被测物体待测尺寸的最佳重构曲线Xnew(w)；
步骤五、根据像素对应实际距离对被测物体待测尺寸长度进行测量：根据CT扫描设备的实际成像范围a×a，单位mm，被测物体的CT图像尺寸n×n，单位为像素数，计算出每个像素之间对应的实际距离为所述步骤三拖动的直线起点与终点已知，分别为(x1,y1)和(x2,y2)，根据步骤四得到的被测物体待测尺寸的最佳重构曲线Xnew(w)，计算N值：被测物体待测尺寸应实际长度为
与现有技术相比，本发明的优点在于：采用本发明提供的方法获取的被测物体待测尺寸速度快、精度高。
附图说明
图1为本发明实施例中标准方形试块的CT图像中框选的长方形区域示意图；
图2为本发明实施例中标准方形试块的边缘响应函数ERF；
图3为本发明实施例中标准方形试块的CT图像垂直界面一维点扩散函数；
图4为本发明实施例中圆形被测物体的CT图像；
图5为在图4中的CT图像中画一条通过圆形被测物体直径的直线后的示意图。
具体实施方式
以下结合附图实施例对本发明作进一步详细描述。
如图1所示的工业CT图像中弱边缘尺寸测量方法，其包含如下步骤：
步骤一、分别采集被测物体和标准试块的CT图像，其中标准试块的密度和厚度与被测物体一致；
步骤二、获取标准试块的CT图像垂直界面一维点扩散函数：
所述标准试块采用圆形标准试块或方形标准试块；
采用GB_T 29069-2012《无损检测工业计算机层析成像(CT)系统性能测试方法》中5.3圆盘卡法获得圆形标准试块CT图像垂直界面一维点扩散函数PSF[k]，求取PSF[k]的累计幅值A；
当标准试块为方形标准试块时，在标准方形试块的CT图像中框选一个长方形区域，使得标准方形试块的一条直线边缘图像包含在该长方形区域中，参见图1所示，对该区域进行二值化分割，获取该直线边缘点，再根据直线边缘点进行拟合获得该直线边缘的斜率k，计算与该直线边缘垂直直线的斜率k′＝1/k，根据该垂直直线的斜率在标准试块边缘上均匀提取N条相隔距离为h的与该直线边缘相垂直的剖线：yn＝k′xn+Bn，Bn＝d+nh,n∈[1,N]，d为起始斜截距；根据每条剖线yn在标准方形试块的CT图像上经过的位置，提取与之对应的直线边缘点的对应灰度值ln，获取所有直线边缘上所有像素的灰度加以合并平均，得到边缘响应函数ERF，参见图2所示；然后对边缘响应函数ERF进行拟合求导获取方形标准试块CT图像垂直界面一维点扩散函数PSF[k]，参见图3所示，具体步骤如下：
边缘响应函数ERF是长度为N的一维数组，设为ERFN，设E[m]为ERFN的子数组，其中a为数组起点，b为数组终点。0≤a<b≤N，其中数组长度为m＝b-a，m∈[15～25]；对E[m]进行3次方最小二乘法拟合，获得拟合后的E[m]′，取中点值E[a+m/2]′作为该段数组拟合后新的Enew[m]值，再对Enew[m]进行3次方最小二乘法拟合，获得拟合后的Enew[m]′，对Enew[m]′进行求导dEnew[m]′，取中点值dEnew[a+m/2]′作为该段数组的PSF值，其中W是归一化系数，W＝max(dEnew[a+m/2]′)；
该标准试块的CT图像垂直界面一维点扩散函数为PSF[s]，对PSF[s]进行累加获得累计幅值
步骤三、获取被测物体待测尺寸的灰度曲线：在被测物体的CT图像中画一条通过待测尺寸的直线，使得该直线包含待测尺寸的两端边界，获取该直线上对应的灰度值l[n]；以圆形被测物体为例，需要测量该圆形被测物体的直径，其CT图像参见图4所示，在圆形被测物体CT图像中画一条通过该圆形被测物体的直径的直线，参见图5所示；
步骤四、对被测物体待测尺寸的两端边界分别进行复原重构：
对步骤三获得的灰度直线上对应的灰度值l[n]以两端边界中心位置为界进行分割分段处理，分别获得2组CT图像垂直界面一维点扩散函数，它分别对应被测物体待测尺寸的两端边界CT图像垂直界面一维点扩散函数，假设左边界对应的CT图像垂直界面一维点扩散函数为PSF[L]，右边界对应的CT图像垂直界面一维点扩散函数为PSF[R]；设h[L]、h[R]分别为PSF[L]和PSF[R]对标准试块的CT图像垂直界面一维点扩散函数PSF[s]作反卷积后得到的复原边界系数，为左侧复原边界系数，为右侧复原边界系数；
在这里，对被测物体待测尺寸的左右两端边界类型进行分类：
1、垂直于CT扫描平面的垂直平面边界；
2、倾斜于CT扫描平面的斜平面边界；
3、不规则面的边界；
然后对获取的h数据进行分类，这里的h数据为h[L]或h[R]：
1、如果h数据中有且仅有一个最大值max[h]，并且max[h]附近数据急剧下降，h数据曲线呈山峰状，即num{max[h]}＝1，则判定为第1种边界；
2、如果h数据显示为方波，即num{max[h]}为多个且连续；则判定为第2种边界；
3、如果h数据显示为曲线，则判定为第3种边界；
对步骤三获得的灰度值l[n]求导，有两个极值分别对应待测尺寸的左右边界的位置，设为j、k，j表示左边界位置，k表示右边界位置；取j、k中间点，即将步骤三获得的灰度值l[n]在这点上一分为二，分开进行边界复原处理：
设左边段的灰度值为L(w)，假定一个长度估计值x0，x0初值小于测量长度，其退化前的幅值为其中A为步骤二中求取的标准试块的CT图像垂直界面一维点扩散函数为PSF[s]的累计幅值；然后根据判断出的左侧边界类型分别进入如下处理：
1、当左侧边界类型判断为1，即垂直平面边界，初始值x0不变，然后计算Xold(L)：
2、当边界类型判断为2，即斜平面边界，计算左侧复原边界系数h[L]中大于的点数R，即计算方波的宽度，然后计算Xold(L)：
3、当边界类型判断为3，即不规则面的边界，设左侧复原边界系数h[L]中max{h[L]}对应的位置为y，取其后的所有数据设为h[y]，然后计算Xold(L)：
将Xold(L)与步骤二获得的标准试块的CT图像垂直界面一维点扩散函数为PSF[s]进行卷积，计算与左边段的灰度值为L(w)的标准偏差S(w)；然后将x0＝x0+1，循环迭代计算，直至最终回计算时的x0值大于测量长度；获得S(w)的曲线，取数组S(w)中的最小值对应的x0，此时的x0即为最佳长度值，然后将最佳长度值x0代入Xold(L)中，获得左侧最佳重构曲线Xnew(L)；
右侧最佳重构曲线Xnew(R)与左侧最佳重构曲线Xnew(L)的获取方式相同，将右侧最佳重构曲线Xnew(R)与左侧最佳重构曲线Xnew(L)合并，得到被测物体待测尺寸的最佳重构曲线Xnew(w)；
步骤五、根据像素对应实际距离对被测物体待测尺寸长度进行测量：根据CT扫描设备的实际成像范围a×a，单位mm，被测物体的CT图像尺寸n×n，单位为像素数，计算出每个像素之间对应的实际距离为所述步骤三拖动的直线起点与终点已知，分别为(x1,y1)和(x2,y2)，根据步骤四得到的被测物体待测尺寸的最佳重构曲线Xnew(w)，计算N值：被测物体待测尺寸应实际长度为
采用上述原理，还可以对面积型尺寸周向边界复原测量：首先，对测量图像进行传统的二值化分割，获取感兴趣对象的初始轮廓线B，B是所有初始轮廓线构成的集合；然后对二值化初始轮廓线采用图像处理中的形态学膨胀运算，进行多次运算，膨胀次数一般大于20，形成一个新的外轮廓线E，E围绕B；对二值化初始轮廓线采用图像处理中的形态学腐蚀运算，进行多次运算，膨胀次数一般大于20，形成一个新的内轮廓线F。为内轮廓线F的所有点在外轮廓线E中寻找最近似的点，搜索方法：内轮廓线F中的任意一点f1计算它与外轮廓线E中所有点的距离，当外轮廓线E中存在唯一的一点e1，它与f1距离最小，则确定e1为f1的近似点；当外轮廓线E中两点e2和e3与f1的距离最小，在外轮廓线E弧线e2e3上的中点e1与f1的距离为弧线e2e3范围内的局部极大值，则确定e1为f1的近似点；连接e1和f1，必然与初始轮廓线B相交于一点b1；设直线e1f1是点b1的剖面直线，获取该直线上所有像素的灰度，对初始轮廓线B上所有点计算剖面直线EF，对每条剖面直线进行所述步骤四中的单边复原方法，复原出感兴趣对象实际边界C；在感兴趣对象复原后轮廓线位置已知的条件下，可通过所述步骤五获得对象任意两边界之间的实际距离。对于面积测量，统计边界C包含的像素点数M＝num(C)，根据CT实际成像范围a×a，单位(mm)，图像尺寸n×n，单位(像素数)。我们可以计算出每个像素之间对应的实际距离为面积计算方法近似采用统计边界C包含的像素点数与相邻像素之间实际距离的平方相乘，即面积
- 还没有人留言评论。精彩留言会获得点赞！