CN102750693A

CN102750693A - Zernike矩对曲线边缘高精度定位的修正方法

Info

Publication number: CN102750693A
Application number: CN2012101768480A
Authority: CN
Inventors: 王珏; 王慧倩; 蔡玉芳; 刘明进; 邹永宁; 陶李
Original assignee: Chongqing University
Current assignee: Chongqing University
Priority date: 2012-05-31
Filing date: 2012-05-31
Publication date: 2012-10-24
Anticipated expiration: 2032-05-31
Also published as: CN102750693B

Abstract

基于Zernike矩对存在一定曲率的连续光滑的曲线边缘高精度定位的修正方法，涉及一种工业CT图像处理方法，首先采用一组标准工件工业CT图像的实际曲率与测量曲率半径误差值数据拟合得到曲率半径误差补偿因子函数，然后通过Zernike矩方法对工业CT图像进行亚像素边缘提取，并进行细化处理，再获取每一个边缘点的曲率半径和曲率圆心，最后应用曲率半径误差补偿因子函数，得到修正后的曲率半径，并根据曲线该点切线法线方向修正边缘点坐标，提高边缘定位精度。

Description

Zernike矩对曲线边缘高精度定位的修正方法

技术领域

本发明涉及一种工业CT定量化图像处理方法，尤其是涉及一种基于Zernike矩的曲线边缘高精度定位方法。

背景技术

工业CT(Industry Computerized Tomography)广泛应用于航空、航天、汽车等工业精密工件内部结构的测量及其缺陷检测领域。工业CT图像尺寸测量方法是一种以工业CT图像为信息载体并从中提取定量数据的测量方法，它以工业CT断层图像为研究对象，计算目标的长度、宽度、高度、壁厚等各种几何参数。与其它测量方法相比，其既能测得实体复杂的内部几何形状(尤其是封闭内腔)的参数，又不破坏实体。

目前，图像几何参数测量发展迅速，其测量的高精度一直是追求的目标。现有技术中先采用传统算子提取像素级边缘，初步确定椭圆参数再采用三级灰度矩和曲率不变原理提取亚像素最后通过最小二乘拟合确定椭圆参数(参见周拥军，寇新建.基于矩和曲率保持的圆形标志精确量测.武汉大学学报：信息科学版，2009，34(2)：203-206.)，这种方法考虑了边缘领域的灰度分布和几何特性，较像素级算子提高了测量精度，但亦有缺陷，其测量精度受到灰度矩边缘检测的局限且需要求矩和超越方程，测量时间较长。

现有技术中基于空间矩的圆标记中心亚像素定位算法(参见郭玉波，姚郁，遆晓光.基于空间矩的圆标记中心亚像素定位算法.吉林大学学报：工学版，2009，39(1)：160-163.)通过改进的空间矩圆形边缘模板进行边缘检测，再通过拟合方法确定圆参数，该测量方法重新推导了空间距的边缘模型，使其适用于圆形边缘的边缘检测，但亦存在不足，其测量精度受到空间矩边缘检测的局限且需要计算六个矩模板和迭代求解非线性方程组，耗时较长。

Ghosal和Mehrotal在空间矩方法基础上首次提出了利用Zernike正交矩来检测亚像素边缘(参见Ghosal S,Mehrotra R.Orthogonal moment operators for sub-pixel edgedetection.Pattern Recognition,1993,26(2):295-306.)，在他们的算法中建立了理想的阶跃灰度模型，通过图像的三个不同阶次Zernike矩计算模型的4个参数，以这4个参数为判断边缘的依据来确定图像中物体边缘的位置。Ghosal方法利用Zernike矩旋转不变性将空间矩需要计算的6个矩模版减小至3个，减少了计算量且精度比灰度矩和空间矩高。

但是，分析Zernike矩计算模型和在工程应用的结果，不难发现Ghosal提出的Zernike矩方法是基于直线边缘的，当检测目标为曲线时，算法将每个模板内的曲线边缘作为直线处理，且当待检测的曲线边缘曲率增大时，其边缘检测误差也随之增大。

为了提高存在一定曲率的连续光滑的曲线边缘检测定位精度，满足工业CT图像的高精度测量要求，针对以上发现的问题，提出了对基于Zernike矩的曲线边缘定位之曲率半径补偿因子修正方法。

发明内容

有鉴于此，为了解决上述问题，本发明的提供了一种Zernike矩对工业CT图像曲线边缘定位的修正方法。本发明实现容易、效果好、精度高,是一种能有效、快速、高精度测量具有曲线轮廓工件的方法。

本发明的目的是这样实现的：

本发明提供的Zernike矩对曲线边缘高精度定位的修正方法,包括以下步骤：

1）计算曲率半径误差补偿因子函数τ。

2）获取图像并提取图像的边缘；

3）对图像中的所有边缘点进行细化处理；

4）获取曲线每个边缘点的曲率半径和曲率圆心；

5）应用曲率半径误差补偿因子函数τ和定位点切线的法线方向，分别得到修正后的曲线曲率半径值和修正方向，从而修正边缘点坐标。

进一步，所述曲率半径误差补偿因子函数τ的计算步骤如下：

11）获取一组已知实际曲率的质量均匀的标准钢制工件的工业CT图像；

12）对工业CT图像采用Zernike矩方法进行边缘检测；

13）获取每个点的曲率半径，并计算曲率半径绝对误差；

14）对曲率半径绝对误差与实际曲率数据按折线段分段拟合，可表达如下：

s = \{\begin{matrix} a \cdot Ω + b & 0 < Ω < m 1 \\ c \cdot Ω + d & m 1 \leq Ω < m 2 \end{matrix}

其中，s为曲率半径测量值绝对误差，Ω为曲率，a,b,c,d,m1,m2为常值，a,c分别为两条直线的斜率，m1为折线段转折点对应的曲率，m2为可测的最大曲率，考虑了模板大小的影响。由此，根据曲率与曲率半径的关系，曲率半径误差补偿因子函数τ可表达为：

τ = \{\begin{matrix} - \frac{c}{R} - d & 1 / m 2 < R \leq 1 / m 1 \\ - \frac{a}{R} - b & R > 1 / m 1 \end{matrix}

其中，τ表示曲率半径误差补偿因子函数，R表示曲率半径值，a,b,c,d,m1,m2为常值。

进一步，所述边缘点细化步骤如下：

31）Zernike矩算法亚像素边缘检测中保留每个点的边缘方向φ和边缘距离l；

32）对图像中的所有边缘点在边缘方向φ上做边缘距离l的非极小值抑制。所谓的l的非极小值抑制，即将边缘点的边缘方向φ分为四个范围0~π/4，π/4~π/2，π/2~3π/4，π/2~π，再比较每个边缘点的边缘距离l₀在边缘方向上相邻的两个像素点的边缘距离值l₁,l₂，判断l₀是否比l₁,l₂小，若不是，则认为此点不是边缘点。

33）对图像边缘进行8方向链码边缘跟踪，去除伪边缘。

进一步，所述获取曲线每个边缘点的曲率半径和曲率圆心步骤如下：

41）根据已知三个相邻边缘点P_i-1，P_i，P_i+1的坐标，可以计算出拟合圆的圆心坐标P_y(x_c,y_c)，

其中圆心坐标位置计算公式如下：

x_{c} = \frac{(y_{i - 1} - y_{i + 1}) (y_{i - 1}^{2} - y_{i}^{2} + x_{i - 1}^{2} - x_{i}^{2}) - (y_{i - 1} - y_{i}) (y_{i - 1}^{2} - y_{i + 1}^{2} + x_{i - 1}^{2} - x_{i + 1}^{2})}{2 [(y_{i - 1} - y_{i + 1}) (x_{i - 1} - x_{i}) - (y_{i - 1} - y_{i}) (x_{i - 1} - x_{i + 1})]}

y_{c} = \frac{(x_{i - 1} - x_{i + 1}) (x_{i - 1}^{2} - x_{i}^{2} + y_{i - 1}^{2} - y_{i}^{2}) - (x_{i - 1} - x_{i}) (x_{i - 1}^{2} - x_{i + 1}^{2} + y_{i - 1}^{2} - y_{i + 1}^{2})}{2 [(x_{i - 1} - x_{i + 1}) (y_{i - 1} - y_{i}) - (x_{i - 1} - x_{i}) (y_{i - 1} - y_{i + 1})]}

42）求拟合圆及其半径R

R = \sqrt{{(x_{i} - x_{c})}^{2} + {(y_{i} - y_{c})}^{2}}

(x-x_c)²+(y-y_c)²＝R²

43）求过P_i点的法线L的方向。

y = \frac{y_{i} - y_{c}}{x_{i} - x_{c}} \cdot x + y_{i} - \frac{y_{i} - y_{c}}{x_{i} - x_{c}} \cdot x_{i}

进一步，所述曲率修正是沿曲线该点的切线法线方向通过所得曲线曲率半径测量值和曲率半径误差补偿因子函数τ按以下方式形成：

R′＝R+τ，

其中，R′表示补偿后的曲率半径测量值。

本发明的优点在于：本发明针对Zernike矩方法理想直线边缘模型对曲线边缘定位误差较大的不足，采用一组已知曲率的标准工件工业CT图像进行Zernike矩方法亚像素边缘提取并细化处理，然后通过获取每个边缘点的曲率半径和其切线法线方向，通过数据拟合得到实际曲率与曲率半径测量值的绝对误差的关系，再由此获得曲率半径误差补偿因子函数，利用曲率半径误差补偿因子函数对Zernike矩方法检测的曲线曲率进行了修正；本发明的方法可以应用于任意具有曲线轮廓工件ICT图像边缘定位和目标测量中，实现了具有曲线轮廓工件工业CT图像的高精度定位和测量。

本发明的其它优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到指导。本发明的目标和其它优点可以通过下面的说明书，权利要求书，以及附图中所特别指出的结构来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述，其中：

图1为本发明实施例提供的Zernike矩对圆形边缘高精度定位的修正方法流程图；

图2为本发明实施例的一组已知半径的高精度标准圆形工件工业CT图片；

图3为本发明实施例的圆测量中圆半径测量值绝对误差与实际曲率关系图；

图4为本发明实施例的圆形工件工业CT图片。

具体实施方式

以下将结合附图，对本发明的优选实施例进行详细的描述；应当理解，优选实施例仅为了说明本发明，而不是为了限制本发明的保护范围。

图1为本发明实施例提供的Zernike矩对圆边缘高精度定位的修正方法流程图，如图所示,包括以下步骤：

1）计算圆半径误差补偿因子函数τ。

2）获取图像并提取图像的边缘；

3）对图像中的所有边缘点进行细化处理；

4）获取半径测量值和圆心；

5）计算圆半径误差补偿因子函数τ；

6）计算修正后的半径值，并与精确值比较。

步骤一、圆半径误差补偿因子函数的计算步骤如下：

1、获取一组已知半径的标准圆形钢制工件的工业CT图片；图2为本发明实施例的圆形工件工业CT图片，其中20个标准圆形工件用序号标出，如图2中的[1]-[20]。

建立标准实验模型：选取20个已知实际半径的质量均匀的标准圆形钢制工件，其实际半径是通过游标卡尺直径测量值计算所得，采用6MeV高能工业CT系统扫描：扫描分度数为2048，采样时间为0.021s，准直器宽度为0.4mm。扫描图像中1个像素对应实际工件0.2612mm；

2、对圆形工件ICT图像采用Zernike矩方法进行边缘检测，并进行细化处理；

3、通过最小二乘拟合法得到圆半径，并计算半径测量绝对误差；

4、通过实际曲率与半径测量值绝对误差数据拟合得到半径误差补偿因子函数：

通过实际曲率与半径测量值绝对误差数据分段拟合得到折线段，其公式可表达为：

s = \{\begin{matrix} a \cdot Ω + b & 0 < Ω < m 1 \\ c \cdot Ω + d & m 1 \leq Ω < m 2 \end{matrix}

其中，s为曲率半径测量值绝对误差，Ω为曲率，a,b,c,d,m1,m2为常值，a,c分别为两条直线的斜率，m1为折线段转折点对应的曲率，m2为可测的最大曲率，考虑了模板大小的影响。本实施例中a约为12.6153，b约为-0.2981，c约为1.6889，d约为-0.0421，m1约为0.0234，m2约为0.2857。其圆半径测量结果如图3和表1所示，图3为本发明实施例的圆测量中圆半径测量值绝对误差与曲率关系图。图3中横坐标为实际曲率，纵坐标为圆半径测量值绝对误差。表1为图2中20个圆形工件ICT图片（序号标出）半径测量结果。

表1

圆序号	实际半径/像素	Zernike矩测量值/像素	绝对误差/像素
				1	1.91	2.76	0.85
2	3.79	4.17	0.38
				3	5.74	6.01	0.27
4	7.66	7.84	0.18
				5	9.57	9.71	0.14
6	11.50	11.61	0.11
				7	13.28	13.36	0.08

8	15.33	15.42	0.09
				9	17.13	17.17	0.04
10	19.16	19.22	0.06
				11	22.99	23.02	0.03
12	30.65	30.65	0.00
				13	42.11	42.10	-0.01
14	47.86	47.81	-0.05
				15	57.43	57.35	-0.08
16	64.30	64.24	-0.06
				17	67.00	66.92	-0.08
18	80.40	80.29	-0.11
				19	94.68	94.48	-0.20
20	114.38	114.17	-0.21

5、通过实际曲率和半径测量绝对误差的关系，所述曲率半径误差补偿因子函数通过以下公式形成为：

τ = \{\begin{matrix} - \frac{c}{R} - d & 1 / m 2 < R \leq 1 / m 1 \\ - \frac{a}{R} - b & R > 1 / m 1 \end{matrix}

其中，τ表示曲率半径误差补偿因子函数，R表示曲率半径值，a,b,c,d,m1,m2为常值，单位为像素。本实施例中a约为12.6153，b约为-0.2981，c约为1.6889，d约为-0.0421，m1约为0.0234，m2约为0.2857。

步骤二、采用Zernike矩方法提取边缘轮廓，Zernike矩亚像素边缘检测的具体实现步骤如下：

1、计算模板Re[M₁₁]，Im[M₁₁]和M₂₀；

其中M_pq表示p阶q次Zernike矩模板，M₁₁表示一阶一次Zernike矩模板，M₂₀表示二阶零次Zernike矩模板，Re[M₁₁]表示M₁₁模板实部，Im[M₁₁]表示M₁₁模板虚部。

2、利用模板Re[M₁₁]，Im[M₁₁]分别与图像卷积运算得到对应的Zernike矩实部Re[A₁₁]，和虚部Im[A₁₁]，M₂₀与图像卷积运算得到二阶零次Zernike矩A₂₀。p阶q次Zernike矩模板M_pq与图像卷积运算后得到p阶q次Zernike矩A_pq。

3、根据以下公式计算边缘模型参数φ，l，k，h：

φ = \tan^{- 1} (\frac{Im [A_{11}]}{Re [A_{11}]});

l = \frac{A_{20}}{A_{11}^{'}};

k = \frac{3 A_{11}^{'}}{2 {(1 - l^{2})}^{3 / 2}};

h = \frac{A_{00} - \frac{kπ}{2} + k \sin^{- 1} (l) + kl \sqrt{(1 - l^{2})}}{π};

其中，φ表示边缘中心与对应像素中心的连线和水平方向的夹角，l表示边缘中心到像素中心距离，k表示边缘相对于背景的灰度阶跃值，h表示背景灰度值，A′_pq为A_pq旋转φ角后的Zernike矩。

4、选取阈值δ和τ，对每个像素点进行边缘点判定，若该点的参数满足条件l≤δ，k≥τ，则认为此点边缘点，利用式子

[\begin{matrix} x_{s} \\ y_{s} \end{matrix}] = [\begin{matrix} x \\ y \end{matrix}] + \frac{Nl}{2} [\begin{matrix} \cos (φ) \\ \sin (φ) \end{matrix}],

计算亚像素边缘点坐标，N为模版的大小，(x_s,y_s)表示亚像素位置，(x,y)表示中心像素位置。

步骤三、边缘细化

由于l≤δ,k≥τ的判定具有全局性，会产生提取局部边缘较粗的现象，所以进行了边缘细化处理，其步骤如下：

1、Zernike矩算法亚像素边缘检测中保留每个点的边缘方向φ和边缘距离l；

2、对图像中的所有边缘点在边缘方向φ上做边缘距离l的非极小值抑制。所谓的l的非极小值抑制，即将边缘点的边缘方向φ分为四个范围0~π/4，π/4~π/2，π/2~3π/4，π/2~π，再比较每个边缘点的边缘距离l₀在边缘方向上相邻的两个像素点的边缘距离值l₁,l₂，判断l₀是否比l₁,l₂小，若不是，则认为此点不是边缘点。

3、对图像边缘进行8方向链码边缘跟踪，去除伪边缘。

步骤四、获取圆半径和圆心；

由于本实施例是针对圆形工件，因此直接采用最小二乘拟合法求半径和圆心。圆的亚像素边缘数据计算出来以后，将这些数据对圆方程进行最小二乘拟合，由方程系数可直接求出圆半径。圆的一般方程为：

F(X,P)＝x²+y²+Ax+By+C＝0

P＝[A,B,C]表示曲线参数，X＝[x²,y²,x，y,1]^T表示坐标向量。对于圆上任意边缘点

则F(X_i,P)表示该点到F(X,P)＝0的代数距离，根据最小二乘原理，目标函数为：

\min Σ_{i = 1}^{n} F^{2} (X_{i}, P)

求解出圆参数P，则圆半径值为：

R = \sqrt{- C + \frac{A^{2}}{4} + \frac{B^{2}}{4}},

步骤五、应用曲率半径误差补偿因子函数，得到修正后的半径值：

根据圆半径测量值R，计算圆形工件工业CT图像的Zernike矩曲率半径误差补偿因子函数：

τ = \{\begin{matrix} - \frac{c}{R} - d & 1 / m 2 < R \leq 1 / m 1 \\ - \frac{a}{R} - b & R > 1 / m 1 \end{matrix}

其中，τ表示曲率半径误差补偿因子函数，R表示曲率半径值，a,b,c,d,m1,m2为常值。本实施例中a约为12.6153，b约为-0.2981，c约为1.6889，d约为-0.0421，m1约为0.0234，m2约为0.2857。单位为像素。

步骤六、计算修正后的半径值，并与精确值比较：

根据半径误差补偿因子函数和半径测量值，修正后的半径值应为：

R′＝R+τ

其中，R′表示补偿后的半径测量值。

图4中圆形工件工业CT图片的修正后半径测量结果如表2所示。

表2圆形工件ICT图片半径测量值

图序号	a	b	c	d
					实际半径/像素	5.72	9.51	20.96	58.02
修正后半径测量值/像素	5.70	9.50	20.95	58.02

修正后测量值绝对误差/像素

-0.02

-0.01

0.00

最后应说明的是，以上实施例仅用以说明本发明技术方案而非限制，尽管参照较佳实施例对本发明进行了说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者同等替换，而不脱离本发明方案的精神和范围。本发明就每个步骤所涉及的技术来说，都是可以用本技术领域的常识做到的，但这些步骤组合起来形成一种基于Zernike矩的曲线边缘高精度定位的修正方法。

Claims

1.Zernike矩对曲线边缘高精度定位的修正方法,其特征在于，包括以下步骤：

1）计算曲率半径误差补偿因子函数τ；

2）获取图像并提取图像的边缘；

3）对图像中的所有边缘点进行细化处理；

4）获取曲线每个边缘点的曲率半径和曲率圆心；

5）应用曲率半径误差补偿因子函数τ，得到修正后的曲线曲率半径，并根据曲线切线法线方向，修正边缘点坐标值。

2.根据权利要求1所述的Zernike矩对曲线边缘高精度定位的修正方法,其特征在于，所述曲率误差补偿因子函数τ的计算步骤如下：

12）对工业CT图像采用Zernike矩方法进行边缘检测；

13）获取每个点的曲率半径测量值，并计算曲率半径绝对误差；

s = \{\begin{matrix} a \cdot Ω + b & 0 < Ω < m 1 \\ c \cdot Ω + d & m 1 \leq Ω < m 2 \end{matrix}

其中，s为曲率半径测量值绝对误差，Ω为曲率，a,b,c,d,m1,m2为常值，a,c分别为两条直线的斜率，m1为折线段转折点对应的曲率，m2为可测的最大曲率，由此，根据曲率与曲率半径的关系，曲率半径误差补偿因子函数τ可表达为：

τ = \{\begin{matrix} - \frac{c}{R} - d & 1 / m 2 < R \leq 1 / m 1 \\ - \frac{a}{R} - b & R > 1 / m 1 \end{matrix}

3.根据权利要求2所述的Zernike矩对曲线边缘高精度定位的修正方法，其特点在于，所述边缘点细化是对图像中的所有边缘点做边缘梯度方向上边缘距离参数l的非极小值抑制来细化边缘。

4.根据权利要求3所述的Zernike矩对曲线边缘高精度定位的修正方法，其特点在于，所述的获取曲线每个边缘点的曲率半径和曲率圆心，其步骤如下：

41）根据已知三个边缘点P_i-1，P_i，P_i+1的坐标，计算出拟合圆的圆心坐标P_y(x_c,y_c)，其中圆心坐标位置计算公式如下：

x_{c} = \frac{(y_{i - 1} - y_{i + 1}) (y_{i - 1}^{2} - y_{i}^{2} + x_{i - 1}^{2} - x_{i}^{2}) - (y_{i - 1} - y_{i}) (y_{i - 1}^{2} - y_{i + 1}^{2} + x_{i - 1}^{2} - x_{i + 1}^{2})}{2 [(y_{i - 1} - y_{i + 1}) (x_{i - 1} - x_{i}) - (y_{i - 1} - y_{i}) (x_{i - 1} - x_{i + 1})]}

y_{c} = \frac{(x_{i - 1} - x_{i + 1}) (x_{i - 1}^{2} - x_{i}^{2} + y_{i - 1}^{2} - y_{i}^{2}) - (x_{i - 1} - x_{i}) (x_{i - 1}^{2} - x_{i + 1}^{2} + y_{i - 1}^{2} - y_{i + 1}^{2})}{2 [(x_{i - 1} - x_{i + 1}) (y_{i - 1} - y_{i}) - (x_{i - 1} - x_{i}) (y_{i - 1} - y_{i + 1})]}

42）通过以下公式求拟合圆及其半径R：

R = \sqrt{{(x_{i} - x_{c})}^{2} + {(y_{i} - y_{c})}^{2}}

(x-x_c)²+(y-y_c)²＝R²

43）通过以下公式求过P_i点的法线L的方向：

y = \frac{y_{i} - y_{c}}{x_{i} - x_{c}} \cdot x + y_{i} - \frac{y_{i} - y_{c}}{x_{i} - x_{c}} \cdot x_{i} .

5.根据权利要求4所述的Zernike矩对曲线边缘高精度定位的修正方法,其特征在于，所述曲率半径修正是沿曲线该点的切线法线方向通过所得曲线曲率半径和曲率半径误差补偿因子函数τ按以下方式形成：

R′＝R+τ，

其中，R表示曲率半径测量值，R′表示补偿后的曲率半径测量值。