CN103854271A - 一种平面摄像机标定方法 - Google Patents

Abstract

本发明公开了一种平面摄像机标定方法，具体为：估计变形前特征点在图像上的位置，并据此计算变形前图像上像素点与其所在剖分三角形顶点的插值关系；根据该插值关系利用特征点变形后的像素坐标计算变形前图像每个像素点在变形后图像上对应的位置，并利用该位置邻近的变形后像素点的像素值插值计算得到其像素值，从而完成图像校正。本发明利用局部线性内插值代替了本来需要通过高次矩阵求逆才能完成的多项式拟合；相比于现有方法，本发明在保证不降低图像校正效果与标定计算精度的前提下，计算复杂度更低、稳定性更高，同时有着很高的畸变抗性。

Description

一种平面摄像机标定方法

技术领域

本发明属于机器视觉领域，涉及一种用于平面特征测量、定位的摄像机标定方法。

背景技术

视觉测量技术因其非接触、高精度、无损伤等优点，在RFID标签封装机、装片机、划片机、引线键合机等IC封装设备上得到了广泛的应用。由于此类设备上的测量对象多为平面距离与角度且不需要得到摄像机内、外参数，故仅采用一副标定靶标图像进行标定计算校正各种线性、非线性图像变形并得到被拍摄平面图像坐标与世界坐标之间的转换关系。

目前主流的摄像机标定方法都是基于三维定位、测量的，其计算一般需要5~20幅从不同空间方位拍摄的标定图像，这类标定方法操作起来过于复杂，无法满足IC封装设备不加人工干预且仅使用一幅标定图像完成计算的在线运行需求。因此，开发精度、稳定性适应设备运行要求的平面标定算法具有重要的理论意义和实用价值。

针对上述要求，国内研究者提出了一些平面标定算法：文献“一种足球机器人视觉***非线性畸变的数字校正方法”（机器人技术与应用,2002,4:37-40）以二元三次多项式作为目标式，最小二乘拟合得到变形前后图像像素点的对应关系，并通过双线性插值得到变形前每个像素点的像素值。文献“一种光学镜头摄像机图像几何畸变的修正方法”（中国图象图形学报,2000,5(7):593-596）与文献“一种简单而精确的径向畸变标定方法”（中国图象图形学报,2004,9(4):429-434）以不同阶次的二元多项式作为目标式分析与实验，得出使用五次多项式能够得到最佳的效果。这类基于多项式的平面标定算法精度较高，但是高次多项式最小二乘拟合的过程涉及到高阶矩阵的求逆运算（以二元五次多项式为例，对应的求逆矩阵为21阶），也导致此类方法的计算量相对较大,这种算法实现起来复杂且稳定性相对不高。

发明内容

本发明提出一种平面摄像机标定方法，在保证标定精度的条件下减小了计算复杂度。

一种平面摄像机标定方法，具体为：创建与待校正图像相同大小的空图像，依据标定文件中存储的坐标映射关系找到空图像中每个像素点Q_j’（x_j’y_j’）对应于待校正图像的像素点Q_j（x_j，y_j），利用待校正图像像素点Q_j（x_j，y_j）的邻近像素点的像素值进行插值运算得到像素点Q_j’（x_j’y_j’）的像素值pix(Q_j’)，赋于像素值后的空图像即为待校正图像通过校正后的理想图像；校正后的理想图像的像素坐标值乘以标定文件中存储的坐标单位转换系数z得到对应的世界坐标值；

所述标定文件的建立包括以下步骤：

标定图像成像的步骤：

调整标定板成像位置，保证其所成的标定图像I中特征点网格没有残缺，利用标定板成像得到标定图像I；

在标定图像中提取特征点信息的步骤：

提取标定图像I中每个特征点P_i的图像坐标(X_i，y_i)，进一步确定特征点P_i(X_i，y_i)在特征点网格上对应的网格位置坐标(r_i，c_i)；遍历各特征点的网格位置坐标，找到其中行坐标的最大值与列坐标的最大值，从而得到特征点网格的行数R与列数C；估算变形前相邻特征点之间的像素间距d＝max((H+2)/(R-1),(W+2)/(C-1))，H和W分别为标定图像I的高和宽；记录相邻特征点的间距的世界坐标值D；

依据特征点的图像坐标对特征点进行三角剖分，并筛选有效剖分三角形集合的步骤：

对所有特征点按照其图像坐标进行得罗里（Delaunay）三角剖分，得到剖分三角形集合τ，计算各剖分三角形的形状偏差角η，从三角形集合τ中剔除η＞阈值η₀的剖分三角形得到有效剖分三角形集合T_e；

构建特征点变形前的图像坐标计算公式的步骤：

若d＝(H+2)/(R-1)，则按照式I估算特征点网格最左上角点P_lu在变形前对应点Q_1u的图像坐标（XX_1u，YY_1u）；若d＝(W+2)/(C-1)，则按照式II估算特征点网格最左上角点P_1u在变形前对应点Q_1u的图像坐标（XX_1u，YY_1u）

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=-1\\ {\mathrm{YY}}_{\mathrm{lu}}=(W-d\×(C-1))\÷2\end{array}\right.---I$

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=(H-d\×(R-1))\÷2\\ {\mathrm{YY}}_{\mathrm{lu}}=-1\end{array}\right.---\mathrm{II}$

构建每个特征点Q_i变形前的图像坐标（XX_i，yy_i）：

$\left\{\begin{array}{c}{\mathrm{xx}}_i={\mathrm{XX}}_{\mathrm{lu}}+c_i\×d\\ {\mathrm{yy}}_i={\mathrm{YY}}_{\mathrm{lu}}+r_i\×d\end{array}\right.;$

利用特征点变形前的图像坐标计算公式计算有效剖分三角形集合变形前的三角形集合的步骤：

对于有效剖分三角形集合T_e中每个三角形△_k，该三角形三顶点P_k1、P_k2、P_k3对应的特征点网格位置坐标带入特征点变形前的图像坐标计算公式，计算得到对应的变形前点Q_k1、Q_k2、Q_k3的图像坐标，从而得到△_k对应的变形前三角形△_k’，由△_k’组成的三角形集合即为变形前的有效剖分三角形集合T_e’；

利用有效剖分三角形集合变形前的三角形集合求解三角形内插值系数的步骤：

定义标定图像I在变形前的图像为I’，遍历I’的每一个像素点B_j，在有效剖分三角形集合变形前的三角形集合τ_e’中找到包含B_j的三角形△_k’；设△_k’的三顶点Q_k1、Q_k2、Q_k3的图像坐标分别为（xx_k1,yy_k1）、（xx_k2,yy_k2）、（xx_k3,yy_k3），像素点B_j在变形前的图像坐标为（xx_j，yy_j），求解二元一次方程组III，得到三角形内插值系数m_j、n_j：

$\left\{\begin{array}{c}{\mathrm{xx}}_{k1}+m_j\×({\mathrm{xx}}_{k2}-{\mathrm{xx}}_{k1})+n_j\×({\mathrm{xx}}_{k3}-{\mathrm{xx}}_{k1})={\mathrm{xx}}_j\\ {\mathrm{yy}}_{k1}+m_j\×({\mathrm{yy}}_{k2}-{\mathrm{yy}}_{k1})+n_j\×({\mathrm{yy}}_{k3}-{\mathrm{yy}}_{k1})={\mathrm{yy}}_j\end{array}\right.---\mathrm{III}$

依据三角形内插值系数计算标定图像的像素点图像坐标的步骤：

遍历标定图像I在变形前的图像I’的每一个像素点B_j，根据B_j对应的三角形△_k’在有效三角形集合τ_e中找到对应的三角形△_k，△_k的三顶点P_k1、P_k2、P_k3的图像坐标分别为（x_k1,y_k1）、（x_k2,y_k2）、（x_k3,x_k3）；计算B_j在变形后的标定图像I中对应的像素点A_j的图像坐标（x_j，y_j）：

$\left\{\begin{array}{c}{x}_j=x_{k1}+m_j\×(x_{k2}-x_{k1})+n_j\×(x_{k3}-x_{k1})\\ y_j=y_{k1}+m_j\×(y_{k2}-y_{k1})+n_j\×(y_{k3}-y_{k1})\end{array}\right.;$

存储坐标单位转换系数和标定图像的像素点变形前后的图像坐标映射关系的步骤：

计算坐标单位转换系数z=D/d，建立标定图像I在变形前的图像I’中每个像素点B_j对应在I中像素点A_j的图像坐标映射关系（xx_j，yy_j）→（x_j，y_j）。

进一步地，所述三角剖分采用得罗里（Delaunay）三角剖分方法。

进一步地，所述插值运算采用双线性插值法、最近邻域插值、双三次插值或分形插值。

本发明的技术效果体现在：

本发明利用了特征点集的阵列分布特点与局部抗畸变性质，该方法用局部线性内插值代替了本来需要通过高次矩阵求逆才能完成的多项式拟合；相比于背景技术里提到的其它平面测量标定方法，本发明提出的方法在保证不降低图像校正效果与标定计算精度的前提下，计算复杂度更低、稳定性更高，同时有着很高的畸变抗性。由于其计算过程中只用到线性插值（三角形内插值），故需要求逆的矩阵仅为2阶，实验证明，对于有着明显非线性畸变的图像，该方法仍然能够达到理想的校正效果与测量精度。

附图说明

图1为标定文件制作方法的流程图。

图2为应用标定文件进行图像校正与坐标转换方法的流程图。

图3为用于制作标定文件的圆点型标定图像。

图4为对附图3中的标定特征点按照其图像坐标进行Delaunay三角剖分，并以形状偏差角阈值为15°筛选得到有效剖分三角形集合示意图。

图5-1为用于测试标定效果的电路板图。其与附图3所示的标定图像是由同一相机在相同条件(物距、焦距、光圈系数等)下拍摄得到。

图5-2为图5-1电路板图像中的待测宽度特征尺寸E的示意图。

图6-1为使用本发明所示标定算法通过附图3所示的标定图像制作标定文件，并对图5-1所示的电路板图像进行校正处理后的图像。

图6-2为使用MIL8.0（加拿大Matrix公司开发的商用图像处理软件包）的标定模块包通过图3所示的标定图像制作标定文件，并对图5-1所示的电路板图像进行校正处理后的图像。

具体实施方式

本发明提出的基于剖分三角形内插值的平面标定算法的核心思想在于利用标定特征点集有效剖分三角形分割标定图像所形成的局部区域内进行线性插值计算完成摄像机标定，从而取代了传统方法需要通过二元高次多项式的拟合计算才能完成的任务。

特征点集的有效Delaunay剖分三角形充满并均匀分布在标定图像的整个幅面上，且由于图像整体畸变对于每个剖分三角形局部区域的影响可以忽略不计（局部畸变抗性）。这两个性质为完成基于平面特征尺寸测量的摄像机标定提供了一种方法：估计变形前所有标定特征点在图像上的位置，并据此计算（三角形内插值）得到变形前图像上每个像素点与其所在剖分三角形顶点的插值关系；根据该插值关系利用提取得到变形后标定特征点图像坐标计算变形前图像每个像素点在变形后图像上对应的位置，并利用该位置邻近的变形后像素点的像素值插值计算得到其像素值，从而完成图像校正；变形前图像中图像坐标到世界坐标的转换，可以乘以转换系数（特征点间距世界坐标值/变形前图像上特征点间距像素坐标值）的方式进行。

利用剖分三角形将变形前的图像进行分成很多个面积不大的（因此要求特征点个数不应太少）三角形区域，通过三角形内插值计算得到其中每个像素点图像坐标与其所落在的剖分三角形顶点图像坐标的内插值关系；并认为此插值关系在图像变形后保持不变（因为每个三角形区域的面积足够小），据此得到变形前图像中每个像素点对应在变形后图像上的位置。变形前的图像可认为其图像坐标与世界坐标呈完全的线性比例关系，故只需在上述基础上再存储该比例关系（后面称为坐标单位转换系数）即可完成标定文件的制作，具体步骤如图1所示，具体为：

第1步：通过平移调整标定板成像位置，保证其所成的标定图像I中，标定特征点网格没有残缺，即每一行的点个数相等。标定板中标定特征点的间距的世界坐标值D（单位一般为mm）。标定特征点是标定板上以一定的标准距离阵列分布的点，在棋盘格型标定板中特征点为每个棋盘格四周的角点，在圆点型标定板中为每个阵列圆形区域的几何中心点。

第2步：由摄像机采集标定图像I，采用适当的标定特征点提取方法（圆点型标定板可采用“基于椭圆拟合”的提取方法，棋盘型标定板可采用“基于Harris角点”的提取方法）计算得到I中每个标定特征点P_i的图像坐标(x_i，y_i)（i为特征点的提取顺序序号），并将P_i的图像行坐标x_i和图像列坐标y_i分别存入N×1（N为从I提取得到的特征点个数）矩阵X与N×1矩阵Y的第i行。

第3步：对第1步得到的特征点进行排序定位（例如，基于得罗里（Delaunay）三角剖分的特征点坐标自动对应方法、基于圆点阵列靶标的特征点坐标自动提取方法、棋盘格模板角点的自动识别与定位方法等等），找到点(x_i，y_i)在特征点网格上对应的位置坐标(r_i,c_i)，并将其中的网格行坐标r_i和网格列坐标c_i分别存入N×1矩阵R与N×1矩阵C的第i行。

第4步：遍历各特征点的网格位置坐标，找到其中行坐标的最大值与列坐标的最大值，从而得到特征点网格的行数R与列数C，如式（1）所示。

$\begin{array}{cc}R=\max \underset{i}{\mathrm{\Σ}}{r}_i,& C=\max \underset{i}{\mathrm{\Σ}}{c}_i\end{array}---\left(1\right)$

第5步：对第2步得到的特征点P_i按照其像素坐标(x_i，y_i)进行得罗里（Delaunay）三角剖分，得到剖分三角形集合τ。平面点集的三角剖分将由由f个点组成的平面点集的凸包平面划分为2f-2-g个三角形，其中，g为点集中落在凸包边界上的点的个数，且这些三角形的顶点均为点集中的点。Delaunay三角剖分是各类三角剖分中一种满足三角形最小内角值最大的最优化剖分方法。

第6步：设定形状偏差角阈值η₀（取值范围一般在5°~25°），按照式(2)计算各剖分三角形的形状偏差角η，从三角形集合τ中剔除其中η>η₀的三角形得到有效剖分三角形集合τ_e；

ω₁、ω₂、ω₃分别为三角形的三个顶角的角度值。

第7步：设标定图像的高、宽分别为H、W,按照式（3）估算未变形情况下特征点之间的像素间距d。

d＝max((H+2)/(R-1),(W+2)/(C-1))        (3)

第8步：根据第7步计算得到的像素间距d，估算特征点网格最左上角点P_lu在变形前对应点（记作Q_lu）的图像坐标（XX_lu,YY_lu）。

如果d=(H+2)/(R-1),则按照式（4）进行计算；

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=-1\\ {\mathrm{YY}}_{\mathrm{lu}}=(W-d\×(C-1))\÷2\end{array}\right.---\left(4\right)$

如果d=(W+2)/(C-1),则按照式（5）进行计算。

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=(H-d\×(R-1))\÷2\\ {\mathrm{YY}}_{\mathrm{lu}}=-1\end{array}\right.---\left(5\right)$

第9步：根据点Q_lu的图像坐标（XX_lu,YY_lu）与每个特征点P_i的网格坐标(r_i,c_i)，按照式（6）计算变形前每个特征点Q_i（此处i的含义与第1步中的含义一致）的图像坐标（xx_i，yy_i）。

$\left\{\begin{array}{c}{\mathrm{xx}}_i={\mathrm{XX}}_{\mathrm{lu}}+c_i\×d\\ {\mathrm{yy}}_i={\mathrm{YY}}_{\mathrm{lu}}+r_i\×d\end{array}\right.---\left(6\right)$

第10步：将第6步得到的有效剖分三角形集合τ_e中每个三角形△_k（k为有效剖分三角形序号）的顶点P_k1、P_k2、P_k3根据其对应的序号i找到对应的网格坐标，带入式（6）计算得到其对应的变形前点Q_k1、Q_k2、Q_k3的图像坐标。从而得到△_k对应的变形前三角形△_k’。由△_k’组成的三角形集合称为变形前的有效剖分三角形集合，记作τ_e’。

第11步：设未变形条件下的标定图像为I’,遍历I’的每一个像素点B_j(j为图像像素点的顺序序号），在τ_e’中找到并记录包含B_j的三角形△_k’（τ_e’中有且仅有一个三角形包含B_j）。设△_k’顶点Q_k1、Q_k2、Q_k3的图像坐标分别为（xx_k1,yy_k1）、（xx_k2,yy_k2）、（xx_k3,yy_k3）,解式（7）所示的二元一次方程组，得到三角形内插值系数m_j、n_j并分别存入WH×1矩阵M、N的第j行中，其中（xx_j，yy_j）为像素点B_j的像素坐标。

$\left\{\begin{array}{c}{\mathrm{xx}}_{k1}+m_j\×({\mathrm{xx}}_{k2}-{\mathrm{xx}}_{k1})+n_j\×({\mathrm{xx}}_{k3}-{\mathrm{xx}}_{k1})={\mathrm{xx}}_j\\ {\mathrm{yy}}_{k1}+m_j\×({\mathrm{yy}}_{k2}-{\mathrm{yy}}_{k1})+n_j\×({\mathrm{yy}}_{k3}-{\mathrm{yy}}_{k1})={\mathrm{yy}}_j\end{array}\right.---\left(7\right)$

第12步：遍历I’的每一个像素点B_j，根据B_j对应的三角形△_k’的序号k在三角形集合τ_e中找到三角形△_k，设其顶点P_k1、P_k2、P_k3的图像坐标分别为（x_k1,y_k1）、（x_k2,y_k2）、（x_k3,x_k3）；根据像素点序号j在矩阵M、N中三角形内插值系数m_j、n_j。按照式（8）计算B_i在变形后图像I中对应的像素点A_j的图像坐标（x_j，y_j）。

$\left\{\begin{array}{c}{x}_j=x_{k1}+m_j\×(x_{k2}-x_{k1})+n_j\×(x_{k3}-x_{k1})\\ y_j=y_{k1}+m_j\×(y_{k2}-y_{k1})+n_j\×(y_{k3}-y_{k1})\end{array}\right.---\left(8\right)$

第13步：依据第1步中输入的特征点间距世界坐标值D和第7步计算得到的特征点间距像素坐标值d，计算坐标单位转换系数z=D/d；构建I’中每个像素点B_j对应在I中像素点A_j的图像坐标映射关系“（xx_j，yy_j）→（x_j，y_j）”，将坐标单位转换系数和坐标映射关系一起存储到硬盘上，便完成了标定文件F_CALIBRATION的制作。

总结上述步骤，第1和2步为标定图像成像的步骤，第3和4步为在标定图像中提取特征点信息的步骤，第5和6步为依据图像坐标对特征点进行三角剖分，并筛选有效剖分三角形的步骤，第7、8和9步为构建特征点变形前的图像坐标计算公式的步骤，第10步为利用特征点变形前的图像坐标计算公式计算有效剖分三角形集合变形前的三角形集合的步骤，第11步为利用有效剖分三角形集合变形前的三角形集合求解三角形内插值系数的步骤，第12步为依据三角形内插值系数计算标定图像的像素点图像坐标的步骤，第13步为存储坐标单位转换系数和标定图像的像素点变形前后的图像坐标映射关系的步骤。

上述步骤完成后，对于任意一副与标定图像在相同条件下拍摄的图像，都可以载入标定文件校正拍摄过程中引入的各类线性、非线性变形，并建立图像坐标系与世界坐标系之间的对应关系，完成各种图像测量应用。具体步骤如附图2所示，具体为：

第1步：载入待处理图像II（II与标定图像I由相同的摄像机在相同条件下拍摄得到，故其高度、宽度也分别为H、W）以及标定文件F_CALIBRATION。

第2步：创建一副高度H、宽度W的空图像II’，用于存放校正后的图像。

第3步：遍历II’中每个像素点F_j’，按照标定文件中存储的对应关系找到II’中每个像素点Q_j’（x_j’y_j’）变形后对应于图像II的位置图像坐标（x_j,y_j）。利用（x_j，y_j）对应的四个临近像素点Q_jlu（int(x_j,),int(y_j)）、Q_jld（int(x_j,)+1,int(y_j)）、Q_jru（int(x_j,),int(y_j)+1）、Q_jrd（int(x_j,)+1,int(y_j)+1）的像素值pix(Q_jlu)、pix(Q_jld)、pix(Q_jru)、pix(Q_jrd)（函数int()表示向下取整操作），按照式（9）所示的双线性插值法得到像素点Q_j’的像素值pix(Q_j’)。遍历结束后，II’中储存的就是图像II通过校正后的理想图像。本发明不局限于上述插值方式，还可以采用“最近邻域插值”、“双三次插值”、“分形插值”等插值算法。

pix(Q_j')＝(1-u)(1-v)pix(Q_jlu)+(1-u)vpix(Q_jld)+u(1-v)pix(Q_jru)+uvpix(Q_jrd)(9)

式（9）中u＝x_j-int(x_j)、v=y_j-int(y_j)。

第4步：读取标定文件F_CALIBRATION中储存的特征点间距世界坐标值D和特征点间距像素坐标值d，计算坐标单位转换系数z=D/d。对于在II’计算得到的特征尺寸像素值，都可以通过乘以该系数的方式转换为世界坐标值。

下面结合一实例作进一步的描述，应用本发明提出的基于方法根据标定图像制作标定文件，并载入标定文件校正一副同样条件拍摄的变形的电路板图像，具体步骤如下：

第1步：读取标附图3所示的标定图像I（图像规格为636×480），采用基于区域边缘椭圆拟合的算法，提取I中各标定特征点P_i的像素坐标(x_i，y_i)。并由用户输入标定特征点间距的世界坐标值D=1mm。

第2步：依据第1步提取得到的特征点坐标，找到各标定特征点P_i(x_i，y_i)在特征点网格上对应的位置坐标(r_i,c_i)，结果见表1（点个数太多，在此仅显示一部分用以说明结果的数据格式）。

表1变形后标定特征点图像坐标及其网格位置坐标

特征点序号	点Pi的图像坐标(xi,yi)	点Pi网格位置坐标(ri,ci)
			1	(21.0043,33,1918)	（1,1）
2	(43.4125,31.1803)	（1,2）
			3	(64.4419,29.1543)	（1,3）
4	(87.0882,27.6020)	（1,4）

5	(110.4163,25.8313)	（1,5）
			6	(134.2161,24.5569)	（1,6）
7	(158.6630,22.9776)	（1,7）
			8	(183.3105,21.7021)	（1.8）
9	(208.5767,20.7308)	（1,9）
			10	(234.3952,19.8137)	(1,10)
...	...	...
			441	(386.2931,464.4960)	(18,16)
442	(411.7172,464.1817)	(18,17)
			443	(437.1926,463.6335)	(18,18)
444	(462.0515,463.0496)	(18,19)
			445	(486.5858,462.1726)	(18,20)
446	(510.8604,461.0551)	(18,21)
			447	(534.7117,459.8381)	(18,20)
448	(558.0476,458.3972)	(18,23)
			449	(580.7406,456.8767)	(18,24)
450	(602.9831,454.9909)	(18,25)

第3步：遍历各点的网格位置坐标，找到其中行坐标的最大值与列坐标的最大值，得到特征点网格的行数R与列数C分别为R=18、C=25。

第4步：对第2步得到的特征点P_i按照其图像坐标(x_i,y_i)进行Delaunay三角剖分，得到剖分三角形集合τ。设定形状偏差角阈值为15°筛选得到有效剖分三角形集合τ_e，结果如图4所示。

第5步：标定图像的高、宽分别为H＝436、W＝480，按照式（3）计算变形前特征点之间的像素间距d＝28.35294118。按照式（4）计算得到特征点网格最左上角点P_lu在变形前对应点Q_lu的图像坐标（XX_lu,YY_lu）=（-1,-22.2353）。

第7步：根据Q_lu的图像坐标以及变形前特征点像素间距d，按照式（6）计算得到变形前每个特征点Q_i的图像坐标（xx_i，yy_i），结果见表2（点个数太多，在此仅显示一部分）。

表2变形前标定特征点的网格位置坐标及其图像坐标

特征点序号	点Qi的网格位置坐标(ri,ci)	点Qi的图像坐标(xi,yi)
			1	(1,1)	(-22.2353,-1)
2	(1,2)	(6.1176,-1)
			3	(1,3)	(34.4706,-1)
4	(1,4)	(62.8235,-1)
			5	(1,5)	(91.1765,-1)
6	(1,6)	(119.5294,-1)
			7	(1,7)	(147.8824,-1)
8	(1,8)	(176.2353,-1)
			9	(1,9)	(204.5882,-1)
10	(1,10)	(232.9412,-1)
			...	...	...
441	(18,16)	(403.0588,481)
			442	(18,17)	(431.4118,481)
443	（18,18）	(459.7647,481)
			444	（18,19）	(488.1176,481)
445	(18,20)	(516.4706,481)
			446	(18,21)	(544.8235,481)
447	(18,22)	(573.1765,481)
			448	(18,23)	(601.5294,481)
449	(18,24)	(629.8824,481)
			450	(18,25)	(658.2353,481)

第8步：将τ_e中每个三角形△_k（k为有效剖分三角形序号）的顶点P_k1、P_k2、P_k3根据其对应的序号i在步骤2的计算结果中找到对应的网格坐标，带入式（6）计算得到其对应的变形前点Q_k1、Q_k2、Q_k3的图像坐标，从而得到△_k对应的变形前三角形△_k’。由△_k’组成的三角形集合称为变形前的有效剖分三角形集合记作τ_e’。

第9步：设未变形条件下的标定图像为I’,遍历I’中每个像素点B_j，在τ_e’中找到并记录包含B_j的三角形△_k’。设△_k’顶点Q_k1、Q_k2、Q_k3的图像坐标分别为（xx_k1,yy_k1）、（xx_k2,yy_k2）、（xx_k3,yy_k3）,按照式（7）列出二元一次方程组求解得到三角形内插值系数m_j、n_j。根据三角形内插值系数m_j、n_j，按照式（8）计算得到B_j在变形后图像I中对应的像素点A_j的图像坐标（x_j，y_j），结果见表3（像素点个数太多，仅显示一部分变形）。

表3变性前图像像素点对应在变形后图像上的图像坐标

像素点序号	变形前像素点图像坐标	变形后对应位置图像坐标
			1	(1,1)	(38.4227,33.2211)
2	(2,1)	(39.1779,33.1501)
			3	(3,1)	(39.9330,33.0792)
4	(4,1)	(40.6882,33.0083)
			5	(5,1)	(41.4434,33.9373)
6	(6,1)	(42.1985,33.8664)
			7	(7,1)	(42.9832,33.8063)
8	(8,1)	(43.7717,32.7476)
			9	(9,1)	(44.5500,32.6776)
10	(10,1)	(45.3268,32.6062)
			11	(11,1)	(46.1037,32.5347)
12	(12,1)	(46.8806,32.4633)
			13	(13,1)	(47.6575,32.3918)
14	(14,1)	(48.4343,32.3203)
			15	(15,1)	(49.2112,32.2489)
...	...	...

305265	(622,480)	(574.5041,456.4620)
			305266	(623,480)	(575.3045,456.4084)
305267	(624,480)	(576.1049,456.3548)
			305268	(625,480)	(576.9052,456.3011)
305269	(626,480)	(577.7056,456.2475)
			305270	(627,480)	(578.5060,456.1939)
305271	(628,480)	(579.3064,456.1402)
			305272	(629,480)	(580.1067,456.0876)
305273	(630,480)	(580.9047,456.0398)
			305274	(631,480)	(581.6892,455.9733)
305275	(632,480)	(582.4737,455.9068)
			305276	(633,480)	(583.2582,455.8403)
305278	(634,480)	(584.0427,455.7738)
			305279	(635,480)	(584.8217,455.7073)
305280	(636,480)	(585.6116,455.6408)

第11步：将第10步计算得到的标定特征点变形前后图像坐标对应关系“（xx_j，yy_j）→（x_j，y_j）”，第1步由用户输入的标定特征点间距物理值D以及第5步计算得到的变性前特征点像素间距d，存储起来得到标定文件“F_CALIBRATION_1”。

第12步：载入电路板图像II（如附图5-1所示，其拍摄条件与附图3所示的标定图像相同）与标定文件“F_CALIBRATION_1”。并创建一副高度H＝480、宽度W＝636的空图像II’，用于存放校正后的图像。

第13步：遍历II’中每个像素点F_j’，按照标定文件中存储的对应关系找到II’中每个像素点Q_j’（x_j’y_j’）变形后对应于图像II的位置像素坐标（x_j，y_j）。利用（x_j,y_j）对应的四个临近像素点Q_jlu（int(x_j,),int(y_j)）、Q_jld（int(x_j,)+1,int(y_j)）、Q_jru（int(x_j,),int(y_j)+1）、Q_jrd（int(x_j,)+1,int(y_j)+1）的像素值pix(Q_jlu)、pix(Q_jld)、pix(Q_jru)、pix(Q_jrd)，按照式（9）所示的双线性插值法得到像素点Q_j’的像素值pix(Q_j’)。遍历结束后，II’中储存的就是图像II通过校正后的理想图像，如附图6-1所示。

第14步：对于附图5-2所示的电路板宽度特征尺寸E，可在校正后的理想图像中测得E的像素尺寸E_PIX_MINE=463.3528373。按照特征点间距世界坐标值D和特征点间距像素坐标值d，计算坐标单位转换系数z=D/d=0.03526971mm/pix，并通过该系数转换得到E的世界坐标尺寸E_WORLD_MINE=z×E_PIX_MINE=16.342(mm)。

第15步：使用MIL8.0的标定模块通过标定附图1所示的标定图像制作标定文件，校正附图3-1所示的电路板图像，结果如附图6-2所示。同时，使用MIL8.0的检测模块测量特征尺寸E得到其世界坐标值E_WORLD_MIL=16.376(mm)。

由此可见，使用本发明提出的方法与使用成熟商用软件包MIL8.0对图像变形（尤其是非线性变形）进行校正的效果接近，且基于二者进行的实际图像特征尺寸测量结果亦可相差在0.2%左右。通过这一对比可以说明，本发明提出的方法在降低计算复杂度、提高算法稳定性的同时保持了与现有标定算法相同的图像校正效果以及标定计算精度。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种平面摄像机标定方法，具体为：创建与待校正图像相同大小的空图像，依据标定文件中存储的坐标映射关系找到空图像中每个像素点Q_j’（x_j’y_j’）对应于待校正图像的像素点Q_j（x_j,y_j），利用待校正图像像素点Q_j（x_j,y_j）的邻近像素点的像素值进行插值运算得到像素点Q_j’（x_j’,y_j’）的像素值pix(Q_j’)，赋于像素值后的空图像即为待校正图像通过校正后的理想图像；校正后的理想图像的像素坐标值乘以标定文件中存储的坐标单位转换系数z得到对应的世界坐标值；

所述标定文件的建立包括以下步骤：

标定图像成像的步骤：

调整标定板成像位置，保证其所成的标定图像I中特征点网格没有残缺，利用标定板成像得到标定图像I；

在标定图像中提取特征点信息的步骤：

提取标定图像I中每个特征点P_i的图像坐标(x_i，y_i)，进一步确定特征点P_i(x_i，y_i)在特征点网格上对应的网格位置坐标(r_i，c_i);遍历各特征点的网格位置坐标，找到其中行坐标的最大值与列坐标的最大值，从而得到特征点网格的行数R与列数C；估算变形前相邻特征点之间的像素间距d＝max((H+2)/(R-1),(W+2)/(C-1))，H和W分别为标定图像I的高和宽；记录相邻特征点的间距的世界坐标值D；

依据特征点的图像坐标对特征点进行三角剖分，并筛选有效剖分三角形集合的步骤：

对所有特征点按照其图像坐标进行得罗里（Delaunay）三角剖分，得到剖分三角形集合τ，计算各剖分三角形的形状偏差角η，从三角形集合τ中剔除η＞阈值η₀的剖分三角形得到有效剖分三角形集合τ_e；

构建特征点变形前的图像坐标计算公式的步骤：

若d＝(H+2)/(R-1)，则按照式I估算特征点网格最左上角点P_lu在变形前对应点Q_lu的图像坐标（XX_lu，YY_lu）；若d＝(W+2)/(C-1)，则按照式II估算特征点网格最左上角点P_lu在变形前对应点Q_lu的图像坐标（XX_lu，YY_lu）

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=-1\\ {\mathrm{YY}}_{\mathrm{lu}}=(W-d\×(C-1))\÷2\end{array}\right.---I$

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=(H-d\×(R-1))\÷2\\ {\mathrm{YY}}_{\mathrm{lu}}=-1\end{array}\right.---\mathrm{II}$

构建每个特征点Q_i变形前的图像坐标（xx_i，yy_i）：

$\left\{\begin{array}{c}{\mathrm{xx}}_i={\mathrm{XX}}_{\mathrm{lu}}+c_i\×d\\ y{y}_i={\mathrm{YY}}_{\mathrm{lu}}+r_i\×d\end{array}\right.;$

利用特征点变形前的图像坐标计算公式计算有效剖分三角形集合变形前的三角形集合的步骤：

对于有效剖分三角形集合τ_e中每个三角形Δ_k，该三角形三顶点P_k1、P_k2、P_k3对应的特征点网格位置坐标带入特征点变形前的图像坐标计算公式，计算得到对应的变形前点Q_k1、Q_k2、Q_k3的图像坐标，从而得到Δ_k对应的变形前三角形Δ_k’，由Δ_k’组成的三角形集合即为变形前的有效剖分三角形集合τ_e’；

利用有效剖分三角形集合变形前的三角形集合求解三角形内插值系数的步骤：

定义标定图像I在变形前的图像为I’，遍历I’的每一个像素点B_j，在有效剖分三角形集合变形前的三角形集合τ_e’中找到包含B_j的三角形Δ_k’；设Δ_k’的三顶点Q_k1、Q_k2、Q_k3的图像坐标分别为（xx_k1,yy_k1）、（xx_k2,yy_k2）、（xx_k3,yy_k3），像素点B_j在变形前的图像坐标为（xx_j，yy_j），求解二元一次方程组III，得到三角形内插值系数m_j、n_j：

$\left\{\begin{array}{c}{\mathrm{xx}}_{k1}+m_j\×({\mathrm{xx}}_{k2}-{\mathrm{xx}}_{k1})+n_j\×({\mathrm{xx}}_{k3}-{\mathrm{xx}}_{k1})={\mathrm{xx}}_j\\ {\mathrm{yy}}_{k1}+m_j\×({\mathrm{yy}}_{k2}-{\mathrm{yy}}_{k1})+n_j\×({\mathrm{yy}}_{k3}-{\mathrm{yy}}_{k1})={\mathrm{yy}}_j\end{array}\right.---\mathrm{III}$

依据三角形内插值系数计算标定图像的像素点图像坐标的步骤：

遍历标定图像I在变形前的图像I’的每一个像素点B_j，根据B_j对应的三角形Δ_k’在有效三角形集合τ_e中找到对应的三角形Δ_k，Δ_k的三顶点P_k1、P_k2、P_k3的图像坐标分别为（x_k1,y_k1）、（x_k2,y_k2）、（x_k3,x_k3）；计算B_j在变形后的标定图像I中对应的像素点A_j的图像坐标（x_j，y_j）：

$\left\{\begin{array}{c}{x}_j=x_{k1}+m_j\×(x_{k2}-x_{k1})+n_j\×(x_{k3}-x_{k1})\\ y_j=y_{k1}+m_j\×(y_{k2}-y_{k1})+n_j\×(y_{k3}-y_{k1})\end{array}\right.;$

存储坐标单位转换系数和标定图像的像素点变形前后的图像坐标映射关系的步骤：

计算坐标单位转换系数z=D/d，建立标定图像I在变形前的图像I’中每个像素点B_j对应在I中像素点A_j的图像坐标映射关系（xx_j，yy_j）→（x_j，y_j）。

2.根据权利要求1所述的平面摄像机标定方法，其特征在于，所述三角剖分采用得罗里（Delaunay）三角剖分方法。

3.根据权利要求1或2所述的平面摄像机标定方法，其特征在于，所述插值运算采用双线性插值法、最近邻域插值、双三次插值或分形插值。

4.一种平面摄像机标定***，包括：

校正模块，用于创建与待校正图像相同大小的空图像，依据标定文件中存储的坐标映射关系找到空图像中每个像素点Q_j’（x_j’y_j’）对应于待校正图像的像素点Q_j（x_j，y_j），利用待校正图像像素点Q_j（x_j，y_j）的邻近像素点的像素值进行插值运算得到像素点Q_j’（x_j’,y_j’）的像素值pix(Q_j’)，赋于像素值后的空图像即为待校正图像通过校正后的理想图像；校正后的理想图像的像素坐标值乘以标定文件中存储的坐标单位转换系数z得到对应的世界坐标值；

所述标定文件的构建***包括以下模块：

标定图像成像的第一模块，用于调整标定板成像位置，保证其所成的标定图像I中特征点网格没有残缺，利用标定板成像得到标定图像I;

在标定图像中提取特征点信息的第二模块，用于提取标定图像I中每个特征点P_i的图像坐标(x_i，y_i)，进一步确定特征点P_i(x_i，y_i)在特征点网格上对应的网格位置坐标(r_i，c_i);遍历各特征点的网格位置坐标，找到其中行坐标的最大值与列坐标的最大值，从而得到特征点网格的行数R与列数C;估算变形前相邻特征点之间的像素间距d＝max((H+2)/(R-1),(W+2)/(C-1))，H和W分别为标定图像I的高和宽；记录相邻特征点的间距的世界坐标值D；

依据特征点的图像坐标对特征点进行三角剖分，并筛选有效剖分三角形集合的第三模块，用于对所有特征点按照其图像坐标进行得罗里（Delaunay）三角剖分，得到剖分三角形集合τ，计算各剖分三角形的形状偏差角η，从三角形集合τ中剔除η＞阈值η₀的剖分三角形得到有效剖分三角形集合τ_e；

构建特征点变形前的图像坐标计算公式的第四模块，用于若d＝(H+2)/(R-1)，则按照式I估算特征点网格最左上角点P_lu在变形前对应点Q_lu的图像坐标（XX_lu,YY_lu）；若d＝(W+2)/(C-1)，则按照式II估算特征点网格最左上角点P_lu在变形前对应点Q_lu的图像坐标（XX_lu，YY_lu）

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=-1\\ {\mathrm{YY}}_{\mathrm{lu}}=(W-d\×(C-1))\÷2\end{array}\right.---I$

$\left\{\begin{array}{c}{\mathrm{XX}}_{\mathrm{lu}}=(H-d\×(R-1))\÷2\\ {\mathrm{YY}}_{\mathrm{lu}}=-1\end{array}\right.---\mathrm{II}$

构建每个特征点Q_i变形前的图像坐标（xx_i，yy_i）：

$\left\{\begin{array}{c}{\mathrm{xx}}_i={\mathrm{XX}}_{\mathrm{lu}}+c_i\×d\\ y{y}_i={\mathrm{YY}}_{\mathrm{lu}}+r_i\×d\end{array}\right.;$

利用特征点变形前的图像坐标计算公式计算有效剖分三角形集合变形前的三角形集合的第五模块：对于有效剖分三角形集合τ_e中每个三角形△_k，该三角形三顶点P_k1、P_k2、P_k3对应的特征点网格位置坐标带入特征点变形前的图像坐标计算公式，计算得到对应的变形前点Q_k1、Q_k2、Q_k3的图像坐标，从而得到Δ_k对应的变形前三角形Δ_k’，由Δ_k’组成的三角形集合即为变形前的有效剖分三角形集合τ_e’；

利用有效剖分三角形集合变形前的三角形集合求解三角形内插值系数的第六模块，用于定义标定图像I在变形前的图像为I’，遍历I’的每一个像素点B_j，在有效剖分三角形集合变形前的三角形集合τ_e’中找到包含B_j的三角形Δ_k’；设Δ_k’的三顶点Q_k1、Q_k2、Q_k3的图像坐标分别为（xx_k1,yy_k1）、（xx_k2,yy_k2）、（xx_k3,yy_k3），像素点B_j在变形前的图像坐标为（xx_j,yy_j)，求解二元一次方程组III，得到三角形内插值系数m_j、n_j：

$\left\{\begin{array}{c}{\mathrm{xx}}_{k1}+m_j\×({\mathrm{xx}}_{k2}-{\mathrm{xx}}_{k1})+n_j\×({\mathrm{xx}}_{k3}-{\mathrm{xx}}_{k1})={\mathrm{xx}}_j\\ {\mathrm{yy}}_{k1}+m_j\×({\mathrm{yy}}_{k2}-{\mathrm{yy}}_{k1})+n_j\×({\mathrm{yy}}_{k3}-{\mathrm{yy}}_{k1})={\mathrm{yy}}_j\end{array}\right.---\mathrm{III}$

依据三角形内插值系数计算标定图像的像素点图像坐标的第七模块，用于遍历标定图像I在变形前的图像I’的每一个像素点B_j，根据B_j对应的三角形Δ_k’在有效三角形集合τ_e中找到对应的三角形Δ_k，Δ_k的三顶点P_k1、P_k2、P_k3的图像坐标分别为（x_k1,y_k1)、（x_k2,y_k2）、（x_k3,x_k3）；计算B_j在变形后的标定图像I中对应的像素点A_j的图像坐标（x_j,y_j）：

$\left\{\begin{array}{c}{x}_j=x_{k1}+m_j\×(x_{k2}-x_{k1})+n_j\×(x_{k3}-x_{k1})\\ y_j=y_{k1}+m_j\×(y_{k2}-y_{k1})+n_j\×(y_{k3}-y_{k1})\end{array}\right.;$

存储坐标单位转换系数和标定图像的像素点变形前后的图像坐标映射关系的第八模块，用于计算坐标单位转换系数z=D/d，建立标定图像I在变形前的图像I’中每个像素点B_j对应在I中像素点A_j的图像坐标映射关系（xx_j,yy_j）→（x_j,y_j）。

Images