一种基于Delaunay三角剖分的特征点坐标自动对应方法

摘要

本发明属于机器视觉领域，具体为一种基于Delaunay三角剖分的特征点坐标自动对应方法。该方法包括以下步骤：（1）二值分割，筛选有效特征点图像区域；（2）最小二乘椭圆拟合计算特征点的像素坐标；（3）对特征点集进行Delaunay三角剖分；（4）筛选有效剖分三角形；（5）遍历有效剖分三角形的非最长边，建立特征点集合的四邻域信息表；（6）查询四邻域信息表，四邻域递归搜索遍历特征点集合，在遍历过程中计算每个特征点的世界坐标，完成其像素坐标与世界坐标的对应。当图像存在拍摄倾角、畸变以及特征点缺失时，本发明提出的方法仍可有效运行，特别适合IC封装视觉定位系统的在线标定。

说明书

技术领域

本发明属于机器视觉领域，涉及一种用于摄像机在线标定过程中圆形 阵列靶标特征点坐标的自动对应方法。特征点坐标对应指的是特征点像素 坐标与世界坐标的对应。

背景技术

视觉定位因其非接触、高精度、无损伤等优点，在IC封装设备上得到 广泛的应用。摄像机标定是视觉定位的重要环节之一，用于建立图像坐标 系与世界坐标系之间的联系，校正各种线性、非线性图像变形。目前，常 用的标定靶标有棋盘格型与圆点阵列型。其中圆点阵列型靶标采用圆形图 像区域中心点作为特征点，具有更高的标定精度，在IC封装设备上被广泛 应用。靶标特征点阵列分布且间距（世界坐标）已知，其世界坐标容易得 到，然而特征点像素提取顺序往往与世界坐标的顺序不一致，为了实现标 定算法，必须保证标定特征点的图像坐标与世界坐标正确对应。在拍摄标 定图像时，现有方法大多采用人工调整调整标定板位置的方式使之放正以， 并人工统计特征点网格的行列数。这违背了IC封装设备在线、全自动的运 行要求，影响了设备的工作效率。另外，在特征点网格分布密集的情况下， 人工统计特征点网格的行列数也极易出现错误。因此，开发无需人工干预 与人工统计、适应IC封装设备在线运行要求的特征点坐标对应方法具有重 要的理论意义和实用价值。

针对摄像机在线标定减少人工干预的需求，国内外研究者提出了多种 特征点坐标的自动对应方法:文献“Fully automatic algorithm for region of  interest location in camera calibration”(Optical Engineering,2002, 41(6):1220-1226)提出一种基于RADON变换的智能感兴趣图像区域方法实 现阵列圆点靶标上特征点坐标的自动对应，但它只能在拍摄倾角小、畸变 （图像非线性变形）小的前提下有效；文献“Robust recognition or  checkerboard pattern for camera calibration”(Optical Engineering, 2006,45(9):1-9)与文章“基于圆点阵列靶标的特征点自动提取方法”(中国机 械工程,2010，21(16):1906-1910)提出的方法能够适应标定板存在旋转的情 况，却需要在靶标图像区域中加上特殊的三角标记；Delaunay三角剖分是 计算几何理论中一种“使剖分三角形的最小角最大”的角度最优三角剖分。设 ρ为平面上的任一点集，则在ρ的每个Delaunay剖分三角形的外接圆的内 部，都不包含ρ中的任何一点。它常常被用于平面点集合的划分，通过遍 历剖分三角形的端点，可以得到点集合中任一点的“附近点”信息。文献 “Automatic Grid Finding in Calibration Patterns Using Delaunay  triangulation”(Technical Report,NRC-46487/ERB-1104,National Research  Council,Canada,August,2003)首次引进了Delaunay三角剖分方法实现特征点 坐标的自动对应，但该方法只对棋盘格型标定靶标有效且依赖于三个标记 圆提供的方位信息。文献“摄像机标定中特征点的一种自动对应方法”(光电 子·激光,2011,22(5):736-739)提出的对应方法基于Delaunay三角剖分，它允 许特征点网格存在一定程度的残缺，但仍需保证特征点网格外围至少有一 条边是完整的。

发明内容

本发明所要解决的技术问题是：提供一种基于Delaunay三角剖分的特 征点坐标自动对应方法，该方法提出的特征点坐标对应方法使得摄像机标 定可以在无需人工干预的情况下自动完成，简化了摄像机标定过程，满足 了IC封装设备对于摄像机在线标定的需求。

本发明提供的一种基于Delaunay三角剖分的特征点坐标自动对应方 法，其特征在于，该方法包括下述步骤：

第1步：读取标定图像I，再对其进行二值分割，得到二值图像I’；

第2步：对二值图像I’进行连通域分析，得到图像区域集合N，计算每 个图像区域的面积值S；

第3步：提取图像区域集合N中每个图像区域的轮廓点坐标，并据此 计算得到每个图像区域的圆度值ΔR，及其半径值R；将ΔR与R带入式2 计算得到N中每个图像区域的相对圆度值Δr；

第4步：在图像区域集合N中筛选满足“S在[S₁,S₂]范围内且Δr＜Δr₀” 的图像区域，将得到的图像区域组成有效特征点图像区域集合N_e，其中， S₁、S₂分别为预先设定的图像区域面积高、低阈值，Δr₀为预先设定的图像 区域相对圆度阈值；

第5步：对有效特征点图像区域集合N_e中每个有效特征点图像区域n_ei， i表示特征点序号，根据第3步提取的轮廓点坐标通过拟合计算得到一个椭 圆方程，由该方程参数计算得到的椭圆中心点即为图像区域n_ei对应的特征 点P_i，其坐标为(x_i，y_i)；

第6步：设由所有特征点P_i组成的点集合为P，对P进行Delaunay三 角剖分，得到剖分三角形集合τ；

第7步：对于剖分三角形集合τ中每个剖分三角形，按照公式 η＝min(|ω₁-90°|,|ω₂-90°|,|ω₃-90°|)计算剖分三角形的形状偏差角η，其中，ω₁、 ω₂、ω₃为剖分三角形的角度值；预先设定的形状偏差角阈值η₀，在剖分三 角形集合τ中筛选满足“η≤η₀”的三角形，得到有效剖分三角形集合τ_e；

第8步：遍历有效剖分三角形集合τ_e中每个三角形的非最长边，将其 两个端点的序号以及方向角分别加入到对方的四邻域信息结构体中，到遍 历结束时，可以建立特征点集合的四邻域信息表；

第9步：查询特征点集合的四邻域信息表，任选一个四邻域点完整的 特征点P₀作为种子点，四邻域递归搜索遍历特征点集合；在遍历过程中， 确定每一层递归当前点P_k相对于种子点P₀的特征点网格坐标(r_k,c_k)，k表示 递归层数，到遍历结束时得到特征点集合中每个特征点P_i的相对网格坐标 (r_i,c_i)；

第10步：找到特征点集合P中所有特征点相对网格坐标的行坐标最小 值r_min=min(r_i)与列坐标最小值c_min=min(c_i)，计算每个特征点P_i的世界坐标 (X_i,Y_i)，

$\left(\begin{array}{c}{X}_i=L\times (c_i-c_{\min })\\ Y_i=L\times (r_i-r_{\min })\end{array}\right)$

其中，L为特征点间距的世界坐标值。将每个特征点P_i的像素坐标(x_i，y_i) 与世界坐标(X_i,Y_i)一一对应，完成特征点坐标的对应。

本发明可应用于IC封装设备在线标定过程中圆形阵列靶标特征点坐标 自动对应。在IC封装设备的在线标定过程中，标定板是由机器自动放置的， 由于特征点的间距往往很小（一般不超过3mm），这样一来很容易导致部分 特征点图像区域被遮挡形成残缺的情况发生，当遮挡发展到特征点网格外 围没有一条边是完整的时候，背景技术提到的方法就不再适用。另外，背 景技术提到的方法还需要人工统计特征点网格的行列数，这不仅违背了设 备自动化运行的要求而且也极易出错。本发明提供的ACDT(Automatic  Correspondence based on Delaunay Triangulation)方法在继承以前方法优点的 前提下着力解决了上述的问题，无需任何人工干预即可实现特征点坐标的 对应。当因为标定板定位误差出现特征点网格的旋转、残缺问题时，ACDT 方法仍可正常发挥功能，另外，该方法对于不超过一般测量用镜头畸变最 高要求（不超过1%）的畸变不敏感。

本发明提出的特征点坐标对应方法使得摄像机标定可以在无需人工干 预与人工统计的情况下自动完成，简化了摄像机标定过程，满足了IC封装 设备对于摄像机在线标定的需求；由于利用了特征点集的阵列分布特点与 局部抗畸变性质，该方法对于畸变、拍摄倾角不敏感，在特征点网格存在 由靶标旋转、平移造成的遮挡缺陷时仍然可以正常发挥功能；相比于背景 技术里提到的两种同样基于Delaunay三角剖分的特征点对应方法，本发明 提出ACDT方法不依赖于标记圆提供的方位信息，无需人工统计特征点网 格个数，并且在遮挡很严重（特征点网格外围没有一条边是完整的）的情 况下仍然可以正常工作，适应性、鲁棒性更好。

附图说明

图1表示ACDT方法的流程图。

图2(a)表示标定图像，图2(b)表示经过筛选得到的有效特征点图像区 域。

图3表示本发明实例提出的有效特征点图像区域筛选算法的流程图。

图4(a)表示特征点集的Delaunay三角剖分结果，图4(b)表示有效剖分 三角形，图4(c)表示有效剖分三角形的非最长边。

图5表示本发明实例提出的特征点网格四邻域信息表建立算法的流程 图。

图6表示选取10号点作为种子点，四邻域递归遍历特征点集合的路径。

具体实施方式

本发明利用特征点集合Delaunay剖分三角形的性质得到各点的邻域点 信息，并利用特征点网格的局部抗畸变性质四邻域递归搜索遍历特征点集 合，确定各点在特征点网格上所处的位置。

本发明联合图像区域面积与图像区域相对圆度两个条件进行有效特征 点图像区域的筛选。其中，图像区域面积可近似等于组成该图像区域的像 素点个数，图像区域相对圆度可以根据GB1598-80《形状和位置公差检 测规定》中关于圆度评定的最小图像区域圆模型计算得到。现有的类似方 法大多只通过图像区域面积进行筛选，另外一些加入了简单的形状判断， 相同条件下的筛选效果不如本方法精确。

本发明利用特征点与图像像素点共同的阵列分布特点，将常用于图像 连通域分析的四邻域递归搜索遍历算法应用于特征点集合的遍历，并在遍 历过程中根据搜索路径确定每个特征点相对于种子点的位置。

本发明利用特征点网格的局部抗畸变性，确定当前点局部行列方向角 与当前点四邻域点位置。

特征点网格的局部抗畸变性：设特征点网格ξ中一点P_i（i表示特征点 序号）与其四邻域点P_i1、P_i2、P_i3、P_i4构成的向量中 相邻向量间的夹角分别为α_i12、α_i23、α_i34、α_i41，按照式1定义特征点网格ξ 的最大畸变角β_max。仿真实验与应用实践证明，在图像畸变不超过测量用图 像采集系统最高要求（不大于1%）的情况下，β_max可以保持在一个较小的 范围（不超过10°）内，这种性质称作特征点网格的局部抗畸变性。

β_max=max(|α_i12-90°|,|α_i23-90°|,|α_i34-90°|,|α_i41-90°|)    式1

如图1所示，本发明方法的具体步骤包括：

第1步：读取标定图像I，用二值分割法（如OSTU自适应阈值分割方 法等）分割I，得到二值图像I’；

第2步：对图像I’进行连通域（Blob）分析得到图像区域集合N，计算 每个图像区域的面积值S。

第3步：使用轮廓提取算法（如八邻域跟踪算法等）提取图像区域集 合N中每个图像区域的轮廓点，并据此计算得到每个图像区域的圆度值Δ R，及其半径值R。将ΔR与R带入式2计算得到N中每个图像区域的相对 圆度值Δr。

$\mathrm{\Delta r}=\frac{\mathrm{\Delta R}}{R}$式2

计算圆度值ΔR以及半径值R可以根据GB1598-80《形状和位置公差 -检测规定》中关于圆度评定的最小图像区域圆法，也可以采用其它圆度 建模计算方法，如最小二乘圆法、最小外接圆法、最小内切圆法等。

第4步：，在图像区域集合N中筛选满足“S在[S₁,S₂]范围内且Δr＜Δr₀” 的图像区域，将得到的图像区域组成有效特征点图像区域集合N_e。图像区 域面积高、低阈值S₁、S₂以及图像区域相对圆度阈值Δr₀可以根据标定板 的成像实际情况设定，S₁、S₂一般分别取预估计特征点图像区域面积值得 1.2倍与0.8倍，Δr₀的取值范围一般在0.05到0.3之间。

第5步：对N_e中每个有效特征点图像区域n_ei（i表示特征点序号），根 据其轮廓点坐标（第3步中提取得到）通过拟合计算（如最小二乘法）得 到一个椭圆方程，由该方程参数计算得到的椭圆中心点即为图像区域n_ei对 应的特征点P_i，其图像坐标为(x_i，y_i)。

第6步：设由所有特征点P_i组成的点集合为P，对P进行Delaunay三 角剖分，得到剖分三角形集合τ。

第7步：设定形状偏差角阈值η₀（取值范围一般在15°~25°），按照式 3计算τ中每个剖分三角形的形状偏差角η，在剖分三角形集合τ中筛选满 足“η≤η₀”的三角形，得到有效剖分三角形集合τ_e，其中，ω₁、ω₂、ω₃为 剖分三角形的角度值；

η＝min(|ω₁-90°|,|ω₂-90°|,|ω₃-90°|)    式3

第8步：遍历有效剖分三角形集合τ_e中每个三角形的非最长边，对每 条边的两个端点Q₁、Q₂进行如下操作：将Q₁的点序号以及方向角加入到 Q₂的四邻域信息结构体，并将Q₂的点序号以及方向角加入到Q₁的四邻域 信息结构体。到遍历结束时，可以建立特征点网格（即特征点集合P）的四 邻域信息表。

四邻域信息表是一个自定义类型的结构体数组，结构体在数组中的索 引序号表示了其所代表的特征点序号，其中四个long型变量记录了当前索 引序号对应的点的四邻域点序号，另外四个double型变量记录了各四邻域 点对应的方向角。

第9步：查询特征点网格的四邻域信息表，任选一个四邻域点完整的 特征点P₀作为种子点，四邻域递归搜索遍历特征点集。在遍历过程中，确 定每一层递归当前点P_k（k表示递归层数）相对于种子点P₀的特征点网格 坐标(r_k,c_k)（后面简称“相对网格坐标”）。到遍历结束时，可以得到特征点 集合中每个特征点P_i的相对网格坐标(r_i,c_i)。

四邻域递归搜索遍历的具体方法如下：

（9.1）计算第k层递归当前点P_k的局部行列方向角：当k=0时，规定 种子点的局部行正、列正、行负、列负方向角δ_0R+、δ_0C+、δ_0R-、δ_0C-分别为 0°、90°、180°、270°。当k＞0时，设从第k-1层递归当前点P_k到P_k的向量 与水平向右方向的夹角为γ，根据P_k相对P_k-1所处的位置，按照式4 计算P_k的局部行列方向角δ_kR+、δ_kC+、δ_kR-、δ_kC-，并将它们归一化到0~360° 的范围内。

式4

（9.2）确定第k层递归当前点P_k的四邻域点位置：按照式5确定P_k的 四邻域点P_kj(j表示四邻域点序号）相对于P_k的位置，其中ε₀为预先设定 的方向偏差阈值（一般在20°内）。

$\left(\begin{array}{c}|{\theta }_{\mathrm{kj}}-{\delta }_{\mathrm{kR}+}|\le {ϵ}_0\Rightarrow P_{\mathrm{kj}}=P_{\mathrm{kR}+}\\ |{\theta }_{\mathrm{kj}}-{\delta }_{\mathrm{kC}-}|\le {ϵ}_0\Rightarrow P_{\mathrm{kj}}=P_{\mathrm{kC}-}\\ |{\theta }_{\mathrm{kj}}-{\delta }_{\mathrm{kR}-}|\le {ϵ}_0\Rightarrow P_{\mathrm{kj}}=P_{\mathrm{kR}-}\\ |{\theta }_{\mathrm{kj}}-{\delta }_{\mathrm{kC}+}|\le {ϵ}_0\Rightarrow P_{\mathrm{kj}}=P_{\mathrm{kC}+}\end{array}\right)$式5

θ_kj表示向量与水平向右方向的夹角，P_kR+、P_kR+、P_kR+、P_kR+分别表示P_k在行正、列正、行负、列负方向上的四邻域点。

（9.3）计算第k层递归当前点P_k的相对网格坐标：设第k-1层递归当前 点P_k-1相对于种子点P₀的网格坐标为(r_k-1,c_k-1)，根据P_k相对于P_k-1所处的位 置，按照式6计算P_k相对于种子点P₀的网格坐标(r_k,c_k)。

$\left(\begin{array}{c}{P}_k=P_{k-1R+}\Rightarrow [r_k,c_k]=[r_{k-1},c_{k-1}]+\left[\mathrm{0,1}\right]\\ P_k=P_{k-1C+}\Rightarrow [r_k,c_k]=[r_{k-1},c_{k-1}]+\left[\mathrm{1,0}\right]\\ P_k=P_{k-1R-}\Rightarrow [r_k,c_k]=[r_{k-1},c_{k-1}]+[0,-1]\\ P_k=P_{k-1C-}\Rightarrow [r_k,c_k]=[r_{k-1},c_{k-1}]+[-\mathrm{1,0}]\end{array}\right)$式6

第10步：找到特征点集合P中所有特征点相对网格坐标的行坐标最小 值r_min=min(r_i)与列坐标最小值c_min=min(c_i)，按照式7计算每个特征点P_i的 世界坐标(X_i,Y_i)，其中，L为特征点间距的世界坐标值。将每个特征点P_i的 像素坐标(x_i，y_i)与世界坐标(X_i,Y_i)一一对应，完成特征点坐标的对应。

$\left(\begin{array}{c}{X}_i=L\times (c_i-c_{\min })\\ Y_i=L\times (r_i-r_{\min })\end{array}\right)$式7

实例：

本发明实例提出的ACDT方法的核心思想在于将图像连通域分析中常 用的四邻域递归搜索遍历方法类比地用来遍历同样具有阵列分布特点的特 征点集合，并在遍历过程中确定所有特征点在特征点网格上的位置，从而 完成特征点坐标的对应。

由特征点集Delaunay三角剖分结果中的有效剖分三角形非最长边组成 的线段集合与特征点网格四邻域点对构成的线段集合相等，这一性质为遍 历特征点提供了一种方法：通过遍历有效剖分三角形非最长边，建立特征 点网格的四邻域信息表，它是一个元素个数等于特征点个数的自定义结构 体数组，结构体的C语言定义如下：

在特征点集合遍历过程中，可以利用特征点网格的局部抗畸变性，查 询四邻域信息表更新当前特征点对应的局部行、列方向，并基于此确定当 前特征点的四邻域点相对于它的位置，得到特征点相对于种子点的网格坐 标。

下面结合附图对本发明的具体实施方式作进一步的描述，应用本发明 提出的ACDT方法对图2(a)所示的标定图像（特征点间距0.8mm）进行特 征点像素坐标与世界坐标的对应，具体步骤如下：

第1步：使用OSTU自适应阈值方法对图2(a)进行二值分割，在分割 得到的二值图像中进行有效特征点筛选。观察发现特征点成像的直径在85 个像素左右，设定图像区域面积高、低阈值分别为7000、5000，图像区域 相对圆度阈值为0.25，进行有效特征点图像区域筛选，结果如图2(b)所示， 筛选过程的流程图如图3所示。

第2步：根据有效特征点图像区域的轮廓点坐标，最小二乘椭圆拟合 计算该图像区域对应的特征点像素坐标，具体过程如下：

假设有效特征点图像区域集合N_e中某个图像区域n_ei（每个图像区域的 中心点对应一个特征点）的轮廓点坐标为(x_il，y_il)(i表示特征点序号，l表示 轮廓点序号），将其带入椭圆方程的一般式中，可以组成以椭圆方程一般式 参数A_i～F_i为未知数的方程组，该方程组残差的平方和Z_i如下式所示；

$Z_i=\underset{l}{\Sigma }{({x_{\mathrm{il}}}^2{A}_i+{2x}_{\mathrm{il}}{y}_{\mathrm{il}}{B}_i+{y_{\mathrm{il}}}^2{C}_i+x_{\mathrm{il}}{D}_i+y_{\mathrm{il}}{E}_i+F_i)}^2$

计算使得Z_i取最小时值时，方程组参数的最小二乘解a_i～f_i，带入下式 计算可以得到n_ei对应的区域中心像素坐标(x_i，y_i)，该坐标即为有效特征点图 像区域n_ei对应的特征点坐标；

$x_i=\frac{b_i{e}_i-c_i{d}_i}{a_i{c}_i-{b_i}^2},$$y_i=\frac{b_i{d}_i-a_i{e}_i}{a_i{c}_i-{b_i}^2}$

对N_e中的每个区域都执行上述操作，可以得到所有特征点的像素坐标， 结果见表2。

第3步：对特征点集合进行Delaunay三角剖分，结果如图4(a)所示。

第4步：设定形状偏差角阈值为5°筛选有效剖分三角形，结果如图4(b) 所示。

第5步：去掉有效剖分三角形的最长边，结果如图4(c)所示，遍历有效 剖分三角形的非最长边，将其两个端点的序号以及方向角分别加入到对方 的四邻域信息结构体中，建立特征点网格的四邻域信息表，结果见表1。该 步骤的算法流程图如图5所示。

表1特征点集四邻域信息表

第6步：选取10号点作为种子点四邻域递归搜索遍历特征点集合，遍 历过程如图6所示。

第7步：在特征点集合的遍历过程中，根据搜索路径确定每个特征点 相对于种子点的网格坐标，结果见表2。

第8步：根据特征点间距，对特征点相对于种子点的网格坐标进行转 换，得到特征点的世界坐标，完成特征点坐标对应，结果见表2

表2特征点对应的中间结果与最终结果

以上所述为本发明的较佳实施例而已，但本发明不应该局限于该实施例和附图所公开的内容。 所以凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。