CN101718528A - 基于数字图像的快速求解圆参数方法 - Google Patents

Abstract

本发明是一种基于数字图像的圆参数的快速求解方法，包含步骤：1)利用数字图像传感器采集被测物体图像，建立平面坐标系；2)对图像进行直方图分析，确定适合分离待测圆边界的阈值。对采集得到的图像进行区域分割，设定一个包含待测圆的正方形区域；3)依据步骤2)中确定的阈值选定一个圆内点，并由该点搜索得到图像中待测圆形物体边缘上一个点；4)设定动态步长S，按照步骤3)中的边缘点搜索规则，得到待测圆形物体边缘的一个点组序列；5)求解圆心坐标及直径参数。本发明可应用于工业检测、医疗生物等领域中圆形零件或器具的几何参数测量，很好解决了数字图像中圆参数求解中速度与效果之间的矛盾，具有较高的实际应用价值。

Description

基于数字图像的快速求解圆参数方法

技术领域

本发明涉及一种对数字图像中边缘特征为圆的参数的快速求解方法。

背景技术

圆检测技术是一种对数字图像中物体边缘特征为圆的参数的求解技术。它在工业检测、医疗生物等领域中圆形零件或器具检测等方面有着广泛的应用。比较经典的圆检测方法有Hough变换圆检测法、最小二乘圆检测法。

圆检测问题可以定义为：对采集得到一幅灰度图像中圆形物体，边缘是最基本的特征之一。图像目标边缘是图像目标连通域与其它区域分界，在灰度图像中物体边缘的灰度值常常会发生急剧变化；因此，根据灰度值的大小变化关系可以搜索得到其边缘点；最后，根据这些点的位置关系求出圆的几何参数。

Hough变换法。Hough变换检测圆是在三维参数空间中的数据运算过程。

圆的Hough变换算法步骤如下：

(1)初始化一个变换域(a，b，r)空间的三维数组A(a，b，r)。a方向上的量化数目为被检测图像宽度上的像素数，b方向上的量化数目为被检测图像高度上的像素数，r方向上的量化数目为待测圆半径估计最大值rmax与最小值rmin之差(步长为1个像素)；

(2)顺序搜索图像中的所有的边缘轮廓点(x，y)(设图像中灰度值为255的是白色点)。对每一目标点，计算a＝|x-rcosθ|，b＝|y-rcosθ|，式中r＝r_min，...r_max，θ＝0，1，...359。在变换域(a，b，r)空间中的对应各点加1：A(a，b，r)＝A(a，b，r)+1；

(3)对A(a，b，r)进行搜索，找到一个三元组(a1，b1，r1)，使所有的A(a，b，r)中A(a1，b1，r1)最大，该三元组即对应待测圆；

(4)如果还存在其它的圆，则重复(1)～(3)步，直到所有的圆都检测完；

(5)画出所有的待测圆(a1，b1，r1)，(a2，b2，r2)，...，(an，bn，rn)。

Hough变换检测圆的最大优点是，对噪声点和粗大边缘像素点不敏感，可以完成对较高精度的圆的参数检测；并且可以同时完成对整个图像中的所有圆的检测。Hough变换检测圆的缺点是，对于Hough变换在三维数据空间的数组累加运算，需要同时判断圆心坐标和半径三个参数的变化，运算量和存储量都比较大，运算所耗费的时间较多，致使整体的检测周期延长，不适用于现实测量中。

最小二乘圆拟合法。

圆的最小二乘拟合方法步骤：

(1)初始化一个二维数组E(x_i，y_i)，搜所得到输入的二值边缘图像的圆的n个边缘像素点集E。(x_i，y_i)表示边缘像素的图像坐标。

(2)设圆的方程为：

(x-x₀)²+(y-y₀)²＝R²

取其残差为：

ε_i＝(x_i-x₀)²+(y_i-y₀)²-R²

其残差平方和为：

$\mathrm{\σ}(x_0,y_0,R)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[{(x_i-x_0)}^2+{(y_i-y_0)}^2-R^2]}^2$

(3)根据最小二乘原理，对上式求偏导方程组(偏导等于0)。最后化简整理得到圆得几何参数：

$\left\{\begin{array}{c}{x}_0=\frac{(\stackrel{\‾}{x^2}\stackrel{\‾}{x}+\stackrel{\‾}{x}\stackrel{\‾}{y^2}-\stackrel{\‾}{x^3}-\stackrel{\‾}{x{y}^2})({\stackrel{\‾}{y}}^2-\stackrel{\‾}{y^2})-(\stackrel{\‾}{x^2}\stackrel{\‾}{y}+\stackrel{\‾}{y}\stackrel{\‾}{y^2}-\stackrel{\‾}{x^{2}y}-\stackrel{\‾}{y^3})(\stackrel{\‾}{x}\stackrel{\‾}{y}-\stackrel{\‾}{\mathrm{xy}})}{2({\stackrel{\‾}{x}}^2-\stackrel{\‾}{x^2})({\stackrel{\‾}{y}}^2-\stackrel{\‾}{y^2})-2{(\stackrel{\‾}{x}\stackrel{\‾}{y}-\stackrel{\‾}{\mathrm{xy}})}^2}\\ y_0=\frac{(\stackrel{\‾}{x^2}\stackrel{\‾}{y}+\stackrel{\‾}{y}\stackrel{\‾}{y^2}-\stackrel{\‾}{y^3}-\stackrel{\‾}{{\mathrm{yx}}^2})({\stackrel{\‾}{y}}^2-\stackrel{\‾}{x^2})-(\stackrel{\‾}{x^2}\stackrel{\‾}{x}+\stackrel{\‾}{x}\stackrel{\‾}{y^2}-\stackrel{\‾}{y^{2}x}-\stackrel{\‾}{x^3})(\stackrel{\‾}{x}\stackrel{\‾}{y}-\stackrel{\‾}{\mathrm{xy}})}{2({\stackrel{\‾}{x}}^2-\stackrel{\‾}{x^2})({\stackrel{\‾}{y}}^2-\stackrel{\‾}{y^2})-2{(\stackrel{\‾}{x}\stackrel{\‾}{y}-\stackrel{\‾}{\mathrm{xy}})}^2}\\ R=\sqrt{{x_0}^2-2\stackrel{\‾}{x}{x}_0+{y_0}^2-2\stackrel{\‾}{y}{y}_0+\stackrel{\‾}{x^2}+\stackrel{\‾}{y^2}}\end{array}\right.$

其中：

$\left\{\begin{array}{ccc}\stackrel{\‾}{x}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i}{n}& \stackrel{\‾}{y}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i}{n}& \stackrel{\‾}{x^2}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^2}{n}\\ \stackrel{\‾}{y^2}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i^2}{n}& \stackrel{\‾}{x^3}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^3}{n}& \stackrel{\‾}{y^3}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i^3}{n}\\ \stackrel{\‾}{\mathrm{xy}}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i}{n}& \stackrel{\‾}{x{y}^2}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i^2}{n}& \stackrel{\‾}{x^{2}y}=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^2{y}_i}{n}\end{array}\right.$

最小二乘圆拟合法只需对边缘点循环一次即可完成运算，复杂的开根号运算只是在最后计算圆的半径时计算一次，因而该算法的计算效率很高，特别适合在嵌入式系统中应用。与Hough变换进行圆拟合比较起来，该法的缺点是抗噪声性能不够理想，小的噪声就会对算法的运算结果造成影响，噪声较大时甚至出现运算结果错误，因此在运用最小二乘法前需要对边缘二值化图像用开、闭运算进行降噪处理和边缘细化处理，在得到的理想的边缘图像的基础上，运用最小二乘拟合法检测圆，可以实现高精度的圆参数检测。

发明内容

本发明的目的是提供一种基于数字灰度图像快速求解其中圆形物体边缘的参数的方法。

基于数字图像的圆参数的快速求解方法，属于数字图像处理的范畴。圆参数指的是：圆直径和圆心坐标。本发明包含下列步骤：首先，利用数字图像采集传感器采集被测图像；对图像进行直方图分析，确定适合分离圆边界的阈值；然后，对采集得到的图像依据确定的阈值寻找一个圆内点，确定步长搜索图像中待测圆形物体边缘点；最后，求出圆心坐标及直径等参数。本发明可以直接对采集得到的原图像直接进行扫描求取圆直径和圆心坐标。包括下列步骤特点：

1)利用数字图像传感器采集被测物体图像；采用工业数字摄像机作为数字图像传感器，根据摄像机分辨率为：行列模式(如1280×1024)，摄像机拍摄的图像存储也为行列格式；故可以设定存储图像的第一个像素点为原点o，沿o点的行方向为x轴方向，沿o点的列方向为y轴方向，这样就可以得到一个平面坐标系xoy，单位为像素。

2)根据建立的平面坐标系xoy，对图像进行直方图分析，对采集得到的图像进行区域分割，设定一个包含待测圆的正方形区域，确定适合分离待测圆边界的阈值；

待测圆直径像素值R预估计是依据数字灰度图像中待测物体边缘像素点灰度值的急剧变化规律。编程搜索图像中的待测圆，在图像中设定一个包含待测圆的正方形区域L×L，目的是为了提高搜索速度。不同摄像机的像元值大小不同，不同物距时的像元值大小也不同，这里设像元值为k(mm/pixel)，设正方形区域边长为L(mm)，可以大致估算出待测圆直径的像素值(pixel)：R＝L/k。

3)依据确定的阈值寻找一个圆内点，搜索图像中待测圆形物体边缘点。具体搜索方法为：确定圆内一点p(x₀，y₀)，搜索得到点组(A₀，B₀，C₀，D₀)。当待测圆大致位置已知时，可以人工选定一个合适圆内点p(x₀，y₀)；当位置未知时，可以根据确定边界分离阈值进行判断，选定一个合适的圆内点p(x₀，y₀)。根据采集到的灰度图像中待测物体圆端面灰度值分布情况，决定由P点在步骤2)中设定的L×L区域中，沿其所在行方向向左、向右搜索边缘内侧点(或外侧点)A₀(x₁，y₁)、B₀(x₂，y₂)；接着由B₀点沿其所在列方向向下(或向上)搜索圆边缘的一个内侧点(或外侧点)C₀(x₃，y₃)；然后从点C₀(x₃，y₃)沿其所在行方向向左搜索圆边缘的一个内侧点(或外侧点)D(x₄，y₄)。假设所测圆图像为标准圆，根据摄像机坐标系xoy则有：A₀B₀和C₀D₀均为圆上平行于x轴的弦，B₀C₀和A₀D₀均为平行于y轴的弦；且A₀B₀⊥B₀C₀，B₀C₀⊥C₀D₀，几何关系如图1。可以把点A₀、B₀、C₀、D₀这样的四个点看成一个点组(A₀，B₀，C₀，D₀)。求出弦A₀B₀和C₀D₀中点坐标X₀，X₁，求出弦B₀C₀和A₀D₀的中点坐标Y₀，Y₁。

4)确定动态步长S，按照步骤3)中的搜索方法，得到待测圆形物体边缘点的一个点组序列，具体方法是：

动态步长S的选择。设一个步长变量S′，保持圆内点p(x₀，y₀)的x轴坐标不变，使点p(x₀，y₀)沿直线y＝y₀以S′为步长按递减(或递加)规则变化，并从该点按照步骤3)中方法搜索整个圆可以得到一个这样的点组序列如：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_i，B_i，C_i，D_i)}。

依次改变步长S′按照上面的规则搜索，就可以得到很多这样的点组序列如：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_j，B_j，C_j，D_j)}，{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_k，B_k，C_k，D_k)}等等。通过计算可以得到同样多组的中点坐标序列{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}，{(X₀，Y₀)，(X₁，Y₁)......(X_j，Y_j)}，{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}，{(X₀，Y₀)，(X₁，Y₁)......(X_k，Y_k)}等等，这里i，j，k＝1，2，3......n。根据数理统计学中分布的χ²-检验法判断各组弦中点坐标序列的分布情况。在保证测量精确性和较快运行时间的条件下得到一个步长范围[a，b]。在该区间[a，b]内选定一个合适的步长S’。计算由步骤2)中估算出来的直径与步长S’比值，可以得到一个比例因子，设为K；则有：K＝R/S’。应用到待测圆直径大小不同的情况下，就可以得到一个动态步长S＝R/K，这里，R为变量。

5)求解圆心坐标及直径参数，具体方法是：

圆参数求解。根据步骤4)中选定的步长S，按照步骤3)中搜索点的规则得到一个点组序列：{(A0，B0，C0，D0)，(A1，B1，C1，D1)......(Ai，Bi，Ci，Di)}，这里AiBi和CiDi均为圆上平行于x轴的弦，BiCi和AiDi均为平行于y轴的弦，且AiB⊥BiCi，BiCi⊥CiDi。由此可以得到一个中点序列：{(X0，Y0)，(X1，Y1)......(Xi，Yi)}。假设待测圆为标准圆，根据步骤1)中的摄像机坐标系xoy可知：平行于x轴的弦中点坐标即为待测圆圆心在xoy坐标系中的x轴向坐标，平行于y轴的弦中点坐标即为待测圆圆心在xoy坐标系中的y轴向坐标。又由AiBi⊥BiCi，BiCi⊥CiDi可知：弦AiCi和BiDi为圆的两条直径，几何关系如图2；因此，经过计算可以得到一个直径序列{R0，R1......Ri}。剔除其中的粗大误差点，对剩余点代入公式：

$X=\underset{0}{\overset{i}{\mathrm{\Σ}}}{X}_i/(i+1),Y=\underset{0}{\overset{i}{\mathrm{\Σ}}}{Y}_i/(i+1),R=\underset{0}{\overset{i}{\mathrm{\Σ}}}{R}_i/(i+1)$

求取平均值得到(X，Y)和R即为待测圆的圆心坐标和和直径。

本发明具有速度快，精度高的优点。经典的Hough变换圆检测法，是在三维空间进行计算。最小二乘圆检测法计算公式很复杂，且搜索得到的数据多；因而具有很高的精度，但同时也大大增加了计算量，运行时间较长。

本发明只在二维空间进行计算，设计了一种原理简单，计算量小，求解精度高的快速圆参数求解技术来解决实际的需要。它的圆参数求解原理非常简单，巧妙运用摄像机坐标，进行直径像素值预估计，提出动态步长S，并给出了相应的计算方法。

表格1列出了不同圆检测方法对同一幅直径为30.02mm的标准圆图像在相同条件下求得的圆参数数据。与国际上经典的Hough变换圆检测法和最小二乘圆检测法相比，在相同条件下本发明求解的圆参数具有很高的精度。

表格1不同方法在相同条件下求解得到的圆参数数据

	本发明	最小二乘法	Hough变换法
				机器视觉测量直径(pixel)	561.23	560.19	564
机器视觉测量直径(mm)	30.044	29.988	30.191
				人工测量直径值(mm)	30.02	30.02	30.02
绝对误差(mm)	0.024	-0.032	0.171
				相对误差(％)	0.080	-0.106	0.570

表格2列出了不同圆检测方法对同一幅直径为30.02mm的标准圆图像在不同机器上所需要的处理时间。与国际上经典的Hough变换圆检测法和最小二乘圆检测法相比，从表格2中在四台不同配置的计算机上三种方法所需的处理时间的实验数据可以看出：本发明的速度有了显著的提高。

表格2不同方法在不同环境下的运行时间(表中ms表示毫秒)

	CPU：Pentium(R)43.20GHzRAM：512MB	CPU：Core4200+2.20GHzRAM：1.75GB	CPU：Pentium(R)43.00GHzRAM：512MB	CPU：Pentium(R)41.80GHzRAM：768MB
					本发明	3.10(ms)	4.77(ms)	3.25	5.70(ms)
最小二乘法	7.98(ms)	14.37(ms)	8.30	15.85(ms)
					Hough变换法	92.30(ms)	91.42(ms)	91.36	192.15(ms)

应用本发明可以对数字灰度图像中圆形物体边缘的参数进行快速有效的求解。大量的实验证明，本发明很好解决了数字图像中圆参数求解中速度与效果之间的矛盾，具有较高的实际应用价值。

附图说明

图1是点组(A₀，B₀，C₀，D₀)在摄像机坐标系中的几何关系图；

图2是点组(A_i，B_i，C_i，D_i)在摄像机坐标系中的几何关系图；

图3是本发明流程示意图。

具体实施方式

基于数字图像的快速求解圆参数方法的原理是：建立摄像机平面坐标系xoy，选定圆内一点；并由上述的搜索规则，在采集的灰度图像中待测圆的边缘上搜索得到一个点组(A₀，B₀，C₀，D₀)。A₀B₀和C₀D₀均为圆上平行于x轴的弦，B₀C₀和A₀D₀均为平行于y轴的弦；且B₀C₀⊥A₀B₀，B₀C₀⊥C₀D₀。选定一个合适的步长S，按照相同的方法搜索得到一个点组序列：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_i，B_i，C_i，D_i)}，并计算得到一个相应的弦中点坐标序列：{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}。同理可知：B_iC_i⊥A_iB_i，B_iC_i⊥C_iD_i。根据几何学知识可知：弦A_iC_i和B_iD_i为圆的两条直径，因此可以得到一个直径序列{R₀，R₁......R_i}。剔除其中的粗大误差点，对剩余点求取平均值得到(X，Y)和R即为待测圆的圆心坐标和和直径。本发明无需建立复杂的空间模型，没有复杂的计算公式，原理非常简单；因此，本发明是一种快速的圆参数求解方法。

首先，在本发明的步骤中，建立的摄像机平面坐标系xoy简单实用。待测圆直径像素值R预估计的过程一方面提高搜索速度；另一方面对后面的动态步长的确定提过了依据。如果待测图像简单且对比非常明显，待测圆直径值大致范围已知时，可以不进行直径像素值的预估计这一过程。

其次，在圆内点的选择方面：当待测物大致位置已知时，可以人工选定一个合适圆内点p(x₀，y₀)；当位置未知时，可以根据确定边界分离阈值进行判断，选定一个合适的圆内点p(x₀，y₀)。根据采集得到的灰度图像中待测物体圆端面灰度值分布情况，决定由P点沿其所在行方向向左、向右搜索边缘内侧点(或外侧点)A₀(x₁，y₁)、B₀(x₂，y₂)；接着由B₀点沿其所在列方向向下(或向上)搜索圆边缘的一个内侧点(或外侧点)C₀(x₃，y₃)；最后从点C₀(x₃，y₃)沿其所在行方向向左搜索圆边缘的一个内侧点(或外侧点)D₀(x₄，y₄)。假设所测圆图像为标准圆，根据摄像机坐标系xoy则有：A₀B₀和C₀D₀均为圆上平行于x轴的弦，B₀C₀和A₀D₀均为平行于y轴的弦；且A₀B₀⊥B₀C₀，B₀C₀⊥C₀D₀，几何关系如图1。求弦中点坐标(X，Y)：

X₀＝(x₁+x₂)/2，X₁＝(x₃+x₄)/2

Y₀＝(y₁+y₂)/2，Y₁＝(y₃+y₄)/2

动态步长的确定。动态步长的提出是对不同直径大小的圆来说的。根据待测圆边缘占有的全部像素点数选择一个合适的搜索步长，具有较普遍的意义。动态步长S的确定。首先，可以把步骤3)中的点A₀、B₀、C₀、D₀看成一个组(A₀，B₀，C₀，D₀)；然后，设定一个步长S，由圆内点p(x₀，y₀)，按照步骤3)中的点搜索方法搜索得到一个点组序列{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_i，B_i，C_i，D_i)}。相应的有A_iB_i和C_iD_i均为圆上平行于x轴的弦，B_iC_i和A_iD_i均为圆上平行于y轴的弦。经计算可求出弦中点序列{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}。把步长S改作变量，并按照一定的差值改变步长S，按照同样的方法搜索可以得到很多这样的点组序列如：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_i，B_i，C_i，D_i)}，{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_j，B_j，C_j，D_j)}，{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_k，B_k，C_k，D_k)}等等。同理，可以得到相应的弦中点坐标序列：{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}，{(X₀，Y₀)，(X₁，Y₁)......(X_j，Y_j)}，{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}，{(X₀，Y₀)，(X₁，Y₁)......(X_k，Y_k)}等等，这里i，j，k＝1，2，3......n。根据数理统计学中分布的χ²-检验法判断各组弦中点坐标序列的分布情况。根据测量精确性及运行时间来选择合适的步长S。

本发明对同一幅直径为30.02mm的标准圆图像按不同步长进行实验取得多组实验数据，直径像素预估计值R＝502。首先，对这些数据利用拉伊达法则进行误差点剔除；然后，对剩余点按照数理统计学中非参数估计中的样本点分布的χ²-检验法，在显著水平a＝0.05时判断各组中点坐标序列的分布情况如表格3：

表格3是样本数据分析结果及运行时间(表中μs表示微秒)

步长(pixel)	点数(pixel)	正态分布	运行时间(μs)
				20	26	无法判断	3420.27
15	34	无法判断	3540.28
				12	40	无法判断	3585.95
10	48	服从	3826.07

步长(pixel)	点数(pixel)	正态分布	运行时间(μs)
				8	58	服从	4227.80
6	76	服从	4431.81

由表格3可以得到一个步长区间[10，6]。根据测量精确性及运行时间来选择合适的步长S′＝8。比较步长与直径的关系可以得到一个比例因子，设为K；则有：K＝R/S’＝502/8＝62。应用到待测圆直径不同的情况下，就可以得到一个动态步长S＝R/K＝R/62(R为变量)。

圆参数求解。根据选定的步长S按照同样的方法搜索得到一个点组序列{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_i，B_i，C_i，D_i)}，并求其相应的弦中点坐标序列{(X₀，Y₀)，(X₁，Y₁)......(X_j，Y_j)}，这里i，j＝1，2......n。由A₀B₀和C₀D₀均为圆上平行于x轴的弦，B₀C₀和A₀D₀均为平行于y轴的弦；同理可知：A_iB_i和C_iD_i均为圆上平行于x轴的弦，B_iC_i和A_iD_i均为平行于y轴的弦。依据几何学中弦与圆心的关系及摄像机坐标系xoy可知：中点(X₀，Y₀)，(X₁，Y₁)即为圆心点坐标，即可得到圆心点序列{(X₀，Y₀)，(X₁，Y₁)......(X_j，Y_j)}。又由B₀C₀⊥A₀B₀，B₀C₀⊥C₀D₀，有几何知识及勾股定理可知：弦A₀C₀和B₀D₀均为圆的直径。

根据勾股定理即可求出圆的直径：

$R_0=\sqrt{A_0{B_0}^2+B_0{C_0}^2},R_1=\sqrt{B_0{C_0}^2+C_0{D_0}^2}$

$R_0=\sqrt{{(x_3-x_1)}^2+{(y_3-y_1)}^2},R_1=\sqrt{{(x_4-x_2)}^2+{(y_4-y_2)}^2}$

同理可知：A_iB_i⊥B_iC_i，B_iC_i⊥C_iD_i。由此可以得到一个相应的直径序列{R₀，R₁......R_j}。

在实际测量中，由于偶然误差得客观存在，所得得数据总存在着一定得离散性。但也可能由于粗大误差出现个别离散较远的数据，这通常称为坏值或可疑值。下面采用拉伊达准则来剔除坏点：

①在比较好的测量条件下，能获得较清晰的图像。待测物边缘均匀容易区分，获得的数据波动范围小时，可以采用均值区间法剔除误差较大点。首先，把所求得圆心坐标数据代入公式：

$X=\underset{0}{\overset{i}{\mathrm{\Σ}}}{X}_i/(i+1),Y=\underset{0}{\overset{i}{\mathrm{\Σ}}}{Y}_i/(i+1),$ $R=\underset{0}{\overset{i}{\mathrm{\Σ}}}{R}_i/(i+1)$ 计算求得圆心坐标、直径的平均值。接着设定一个合适的条件：-δ≤X-X_i≤δ，-δ≤Y-Y_i≤δ(δ可以根据测量要求设定)。根据该条件从所求得的圆心参数数据值中剔除误差较大点。分别求出剩余数据的圆心坐标、直径的平均值，并以此时的X、Y、R为圆参数的精确值。也可根据情况，依次减小δ，(X，Y为每次计算时，所有满足要求数据的平均值)进行迭代计算，直至达到测量精度要求；并以此时的圆心坐标平均值X、Y为待测圆圆心坐标的精确值，以此时的直径平均值为待测圆直径的精确值。

②在测量条件有变化时，采集得到的图像中物体边缘会有少量凸点、凹点或缺陷，因而可能导致搜索得到的数据中含有坏点。我们知道坏点的存在对测量结果的平均值的影响很大。这里采用拉伊达法则进行坏点剔除。该方法按正态分布理论，以最大误差范围3s为依据进行判别，偏差Δx_i＝x_i-x， $s=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2}=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\mathrm{\Δ}{x}_i\right)}^2}$ 指的是样本偏差。如果偏差|Δx_i|＞3s时，则认为x_i是含有粗大误差的坏值。根据拉伊达法则对所得到的数据进行分析处理剔除误差较大点，按照公式：

$X=\underset{0}{\overset{i}{\mathrm{\Σ}}}{X}_i/(i+1),Y=\underset{0}{\overset{i}{\mathrm{\Σ}}}{Y}_i/(i+1),R=\underset{0}{\overset{i}{\mathrm{\Σ}}}{R}_i/(i+1)$

对剔除坏点后所有符合要求的剩余数据分别求出X、Y、R平均值并以此作为圆参数的精确数值。

图1是本发明中第一个点组(A₀，B₀，C₀，D₀)的搜索方法：弦A₀B₀和C₀D₀分别在直线x＝x₀和x＝x₃上，弦B₀C₀在直线y＝y₂上。如果待测圆为标准圆，则有：弦A₀B₀和C₀D₀关于圆心对称；所以有弦A₀B₀和C₀D₀平行于x轴，弦B₀C₀和A₀D₀平行于y轴，且B₀C₀⊥A₀B₀，B₀C₀⊥C₀D₀。

图2是本发明中搜索得到的一个点组序列中的任意一个点组(A_i，B_i，C_i，D_i)。同理，有弦A_iB_i和C_iD_i平行于x轴，弦B_iC_i和A_iD_i平行于y轴，且B_iC_i⊥A_iB_i，B_iC_i⊥C_iD_i。根据几何学知识易知：平行于x轴的弦中点坐标即为待测圆圆心在xoy坐标系中的x轴向坐标，平行于y轴的弦中点坐标即为待测圆圆心在xoy坐标系中的y轴向坐标，弦A_iC_i和B_iD_i为圆的两条直径。

图3是本发明的详细流程图。首先，由图像传感器获取待测物体图像，对图像进行直方图分析取得合适的分离所求边缘的阈值。其次，建立摄像机平面坐标系，设定待测圆所在的大致区域。选定圆内一点，设定步长S，按照一定规律搜索得到待测圆边缘上的一个点组序列{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)......(A_i，B_i，C_i，D_i)}。计算求得相应的圆参数序列{(X₀，Y₀)，(X₁，Y₁)......(X_i，Y_i)}和{R₀，R₁......R_i}，对这些圆参数序列进行坏点剔除；最后，将符合要求的剩余数据代入公式计算得到圆心坐标(X，Y)和直径R。

基于数字图像的快速圆参数求解方法的要点是：

1、摄像机平面坐标系xoy的建立，要保证x轴和y轴分别与摄像机扫描图像时的行列方式的行和列的平行关系。

2、设定待测圆在图像中的大致区域。要求采集得到的图像中待测边缘的内侧或外侧像素点要均匀，且其灰度值能明显区分出来，以保证更好地识别待测圆目标，更好地设定搜索区域。为了做到这一点，要保证待测零件尺寸具有较高的精度，利用照明系统或人工方法尽量改善测量条件。选择合适的阈值，能更好突出分离待测边缘。

3、圆内点的选择。要求：保证该点在待测圆内，位于边缘内侧点外附近，最好能使其y轴坐标大致在圆心y轴坐标附近(或x轴坐标大致在圆心x轴坐标附近)。确定圆内点的方法：①当待测物圆端面在屏幕中的大致位置已知时，可以人工设定一个圆内点。②当待测圆端面位置不好判断时，根据选定的阈值搜索得到待测圆边缘上一个内侧点或外测点，再通过几何关系确定一个圆内点。

4、动态步长选择。根据步长计算公式，要保证选择的步长不能过小或过大；否则，要根据直径预估计值适当减小或增大比例因子。

5、把剔除误差点后剩余数据代入公式计算出圆心坐标(X，Y)和直径R。

本发明中计算用到的点序列仅为待测圆边缘上的一部分像素点，这样使搜索的样本点个数大大减小，提高了计算速度。

本发明中圆参数求解方案是整个方法的重点。它的核心就是合理建立摄像机平面坐标系，确定合适的动态步长，按照摄像机扫描原理搜索点序列，计算得到相应的弦中点序列，剔除其中误差较大点。根据摄像机平面坐标系，建立几何关系，运用简单的几何原理求解得到精确的圆心坐标(X，Y)和直径R。

Claims (8)

1.一种基于数字图像的圆参数的快速求解方法，所述圆参数为直径和圆心坐标，其特征在于：

1)利用数字图像传感器采集被测物体图像，建立基于摄像机行列模式的平面坐标系；

2)对图像进行直方图分析，确定适合分离待测圆边界的阈值。对采集得到的图像进行区域分割，设定一个包含待测圆的正方形区域；

3)依据步骤2)中确定的阈值选定一个圆内点，并由该点搜索得到图像中待测圆形物体边缘上一个点；

4)设定动态步长S，按照步骤3)中的搜索规则，得到待测圆形物体边缘的一个点组序列；

5)求解圆心坐标及直径参数。

2.根据权利要求1所述的圆参数的快速求解方法，其特征在于，步骤1)所述的数字图像传感器为工业数字摄像机，根据摄像机分辨率为行列模式，将摄像机拍摄的图像也按行列格式存储。

3.根据权利要求2所述的圆参数的快速求解方法，其特征在于，步骤1)所述的建立平面坐标系是设定存储图像的第一个像素点为原点o，沿o点的行方向为x轴方向，沿o点的列方向为y轴方向，得到的平面坐标系xoy，单位为像素。

4.根据权利要求1所述的圆参数的快速求解方法，其特征在于，步骤2)所述确定待测圆边界的阈值是依据数字灰度图像中待测物体边缘像素点灰度值的急剧变化规律。编程搜索图像中的待测圆，在图像中设定一个包含待测圆的正方形区域L×L，估算出待测圆直径R的像素值(pixel)：R＝L/k

其中，k(mm/pixel)为像元值，L(mm)为正方形区域边长。

5.根据权利要求1所述的圆参数的快速求解方法，其特征在于，步骤3)和4)中所述圆内点是在所述平面坐标系xoy内确定的点p(x₀，y₀)。当待测圆大致位置已知时，可以人工选定该点p；当位置未知时，可以根据确定边界分离阈值进行判断，选定该点p。

6.根据权利要求1所述的圆参数的快速求解方法，其特征在于，在步骤3)中求待测圆边缘点的搜索方法是依据采集到的灰度图像中待测圆端面灰度值分布情况，决定由P点在设定的所述L×L区域中，沿其所在行方向向左、向右搜索边缘内侧点或外侧点A₀(x₁，y₁)、B₀(x₂，y₂)；接着由B₀点沿其所在列方向向下或向上搜索圆边缘的一个内侧点或外侧点C₀(x₃，y₃)；然后从点C₀(x₃，y₃)沿其所在行方向向左搜索圆边缘的一个内侧点或外侧点D₀(x₄，y₄)；最后得到待测圆边缘的点组(A₀，B₀，C₀，D₀)。其中，A₀B₀和C₀D₀均为圆上平行于x轴的弦，B₀C₀和A₀D₀均为平行于y轴的弦，且A₀B₀⊥B₀C₀，B₀C₀⊥C₀D₀；求出弦A₀B₀和C₀D₀中点坐标X₀，X₁，求出弦B₀C₀和A₀D₀的中点坐标Y₀，Y₁。

7.根据权利要求6所述的圆参数的快速求解方法，其特征在于，在步骤4)中所述动态步长S的选择包括步骤：

设一个步长变量S′，保持圆内点p(x₀，y₀)的x轴坐标不变，使所述点p(x₀，y₀)沿直线y＝y₀以S′为步长按递减或递加规则变化，并从该点按照所述测圆边缘点的搜索方法搜索整个圆得到点组序列：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)……(A_i，B_i，C_i，D_i)}；

依次改变步长S′，按照上述的规则搜索，得到多个点组序列：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)……(A_j，B_j，C_j，D_j)}，{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)……(A_k，B_k，C_k，D_k)}。通过计算得到同样多的中点坐标序列：{(X₀，Y₀)，(X₁，Y₁)……(X_i，Y_i)}，{(X₀，Y₀)，(X₁，Y₁)……(X_j，Y_j)}，{(X₀，Y₀)，(X₁，Y₁)……(X_i，Y_i)}，{(X₀，Y₀)，(X₁，Y₁)……(X_k，Y_k)}，这里i，j，k＝1，2，3……n；

根据数理统计学中分布的χ²-检验法判断各组弦中点坐标序列的分布情况。在精度要求许可的情况下得到一个步长范围[a，b]，在该区间[a，b]内选定一个合适的步长S’；

计算由所述待测圆直径R与步长S’比值，得到一个比例因子K，则有：K＝R/S’。应用到不同直径的待测圆，得到一个动态步长S＝R/K。这里，R为变量。

8.根据权利要求7所述的圆参数的快速求解方法，其特征在于，在步骤5)中所述求解圆心坐标及直径参数包括步骤：

根据选定的动态步长S，按照步骤3)中搜索点的规则得到一个点组序列：{(A₀，B₀，C₀，D₀)，(A₁，B₁，C₁，D₁)……(A_i，B_i，C_i，D_i)}和中点序列：{(X₀，Y₀)，(X₁，Y₁)……(X_i，Y_i)}；

在所述摄像机坐标系xoy中，平行于x轴的弦中点坐标即为待测圆圆心在xoy坐标系中的x轴向坐标，平行于y轴的弦中点坐标即为待测圆圆心在xoy坐标系中的y轴向坐标；又由A_iB_i⊥B_iC_i，B_iC_i⊥C_iD_i可知：弦A_iC_i和B_iD_i为圆的两条直径，由此经过计算可以得到一个直径序列{R₀，R₁……R_i}；

根据所述的圆参数序列：{(X₀，Y₀)，(X₁，Y₁)……(X_i，Y_i)}和{R₀，R₁……R_i}，采用拉伊达准则来剔除其中的粗大误差点。对剩余点按下式求取平均值，得到圆心坐标(X，Y)和直径R：

$X=\underset{0}{\overset{i}{\mathrm{\Σ}}}{X}_i/(i+1),$ $Y=\underset{0}{\overset{i}{\mathrm{\Σ}}}{Y}_i/(i+1),$ $R=\underset{0}{\overset{i}{\mathrm{\Σ}}}{R}_i/(i+1).$

Images