一种基于一致性采样算法的标靶单应性矩阵计算方法

1.本发明涉及相机标定技术领域，具体地说，涉及一种基于一致性采样算法的标靶单应性矩阵计算方法。


背景技术：

2.单应性矩阵在图像校正、图像拼接、相机标定和视觉slam(simultaneous localization and mapping，即时定位与地图创建)等计算机视觉领域有着非常重要的作用，单应性矩阵是二维世界坐标系到三维图像坐标系的映射转换关系矩阵，它的好坏决定了对应应用领域的计算结果。
3.相机标定在三维世界和相机成像平面之间建立了对应关系，此映射关系通过单应性矩阵体现。单应性矩阵的结果精度会直接影响相机内参、外参和图像的精确性。
4.张正友平面标定法是以棋盘格点为标定模板的。但该平面标定法所使用的标定点过多，故标定点检测与匹配过程中的误差源较多。
5.经检索，中国专利号zl2016109183420，公开了一种高精度的摄像机动态标定方法；该申请案利用模板图像匹配的方法确定标定图像中的棋盘格靶标区域，通过棋盘格x型标定点检测方法剔除边缘畸变较大的棋盘格外标定点与噪声点；基于高斯曲面优化算法提高棋盘格靶标x型标定点的定位精度；基于非量测畸变校正的摄像机标定方法实现摄像机标定，引入畸变模型与畸变系数优化模型提高摄像机标定精度；最后通过采集到的平面运动靶标序列图像完成摄像机的动态标定，确定平面运动图像的图像坐标与世界坐标对应的最佳单应性矩阵h；该申请案整个标定过程只需求解与优化单个畸变系数，标定精度得到了提高；但该申请案仍然存在标定点多，标定点检测与匹配过程中的误差源较多的问题。


技术实现要素：

6.1.发明要解决的技术问题
7.鉴于上述现有技术中存在的问题，本发明提供了一种基于一致性采样算法的标靶单应性矩阵计算方法；本发明是一种求相机内参矩阵和外参矩阵相乘得到的单应性矩阵的方法，可实现有效靶点的快速提取，去除了无效匹配点对求解单应性矩阵的影响。
8.2.技术方案
9.为达到上述目的，本发明提供的技术方案为：
10.本发明的一种基于一致性采样算法的标靶单应性矩阵计算方法，其步骤为：
11.(1)设计合理标靶；
12.(2)通过特征点图像像素坐标和世界坐标构建齐次矩阵方程组剔除误匹配点，进行标靶点非三点共线筛选，采用dlt直接线性变换法构建单应性矩阵方程和超定方程，由最小二乘法对h矩阵参数进行求解；
13.(3)通过平方代价函数均方误差计算预测图像像素坐标和真实值图像像素坐标的误差，并根据阈值k划分第i点属于内点集合还是外点集合；
14.(4)根据外点是服从高斯分布，设为标准正态分布n～(0,1)，则其均值为0，方差为σ；理论阈值和实际阈值的偏离程度为残差符合n维的χ2分布的特点，采用卡方单边独立性检验获取阈值k；
15.(5)平方代价函数均方误差小于阈值k，得出最优单应性矩阵。
16.更进一步地，步骤(1)中定性分析标靶取点个数n取值，由ransac算法计算：
[0017][0018]
式中，w为内点概率，wn为选择的n个点都是内点的概率，(1-wn)k表示重复k次不是全部的n个点都是内点的概率，p为置信度。
[0019]
更进一步地，步骤(1)中置信度取p＝0.95，n取满足计算单应性矩阵h的最小匹配数，取n＝9，即标靶取9个点。
[0020]
更进一步地，步骤(2)中，从步骤(1)设计的9个标靶点随机抽出5个样本并保证这5个样本之间不共线，具体为：
[0021]
从标靶中随机选取的5个特征点有10个约束单应性矩阵方程，矩阵方程形式：
[0022]
ah＝b
[0023]
a为图像坐标系坐标值和世界坐标系坐标值组成的系数矩阵，h是单应性矩阵；
[0024]
依据非线性方程组解的情况判断标靶点是否共线；判断标准为：
[0025]
(1)若r(a)＝r(a|b)《9，非线性方程组为欠定方程,h矩阵有n-r(a)个非零特解，该非零特解任意线性组合是h的解，此时标靶五点有三点共线点，为废靶；
[0026]
(2)若r(a)＝r(a|b)＝9，非线性方程组为超定方程，标靶为正常标靶，可使用5个不共线特征点组成的10个约束单应性矩阵方程求h矩阵。
[0027]
更进一步地，步骤(2)中，超定方程使用线性最小二乘法计算h，构造欧几里2-范数则最小二乘解为经推导得：
[0028]at
ah＝λh
[0029]
h是矩阵a
t
a特征值为λ的特征向量，此时欧几里2-范数当特征值λ取最小λ
min
时，最小，此时最小二乘解是a
t
a最小特征值λ
min
对应的特征向量；特征值λ
min
对应的特征向量组成h矩阵。
[0030]
更进一步地，步骤(3)中，将正常标靶中9个点剔除三共线点后，筛选出6个不共线点集，记为m0；从m0中取五不共线点，共个点集组成m'＝{m1,m2,m3,m4,m5}，用超定方程的最小二乘法，由点集m'＝{m1,m2,m3,m4,m5}计算出单应性矩阵对应的h'＝{h1,h2,h3,h4,h5}，h'称为临时模型参数；用该模型h'去测试图像像素坐标系的匹配点点集，计算平方误差代价函数：
[0031]
[0032][0033]
k是阈值，若则该i点可划分到内点集合{θ0}，若则该i点可划分到外点集合{θ1}。
[0034]
更进一步地，提出零假设p0：“阈值k与内点的数目无关”，在基于零假设p0成立的情况下，创建的分类器取得极大值时，由pearsonχ2计算公式知：此时卡方值也取极大值，表示理论阈值和实际阈值的偏离程度最大，可以拒绝零假设p0：“阈值k与内点的数目无关”；可接受备择假设p1：“阈值k与内点数目有关”，说明阈值k能最有效区分内点和外点。
[0035]
更进一步地，阈值k的计算过程如下：
[0036]
1)阈值k随机存在于[0,n]中，用二分法将阈值范围分成2个部分[0,k]和[k,n]，将二分阈值得到的内点和外点两种独立要素视为2个自由度，判断k值对内点和外点是否有显著区分，作为k为阈值的标准；确认联合自由度，选择显著性水平α，服从联合自由度为1的卡方分布，对应卡方分布表得抽样数据表；
[0037]
2)通过抽样数据表分别计算各行各列理论期望值α1,β1,α2,β2；
[0038]
3)由于自由度存在偏差，用yates修正计算公式求卡方值χ2；用各行各列理论期望值α1,β1,α2,β2计算ei(i＝1,2,3,4)，用内点和外点观察频数a1(k),a2(k),b1(k),b2(k)计算ai(i＝1,2,3,4)；
[0039]
4)在不同阈值中选出最大表示理论阈值和实际阈值的偏离程度最大，可以拒绝零假设p0：“阈值k与内点的数目无关”；可接受备择假设p1：“阈值k与内点数目有关”，说明阈值k能最有效区分内点和外点，此时为最优阈值k。
[0040]
3.有益效果
[0041]
采用本发明提供的技术方案，与已有的公知技术相比，具有如下显著效果：
[0042]
(1)本发明的一种基于一致性采样算法的标靶单应性矩阵计算方法，对圆形标靶取点进行优化，设计标靶取9个点，兼顾了求单应性矩阵h最小匹配数和运算复杂度；同时，进行标靶点非三点共线筛选，剔除三点共线的冗余点，可实现有效靶点的快速提取，去除了无效匹配点对求解单应性矩阵的影响。
[0043]
(2)本发明的一种基于一致性采样算法的标靶单应性矩阵计算方法，用平方代价
函数的误差计算图像像素坐标系预测值和真实值的误差，误差和阈值k比较分类成外点和内点；并采用卡方独立性检验获取，判断阈值k能否有效区分内点和外点，对内点和外点的区分更加精确，从而提高了标靶单应性矩阵计算精确性。
附图说明
[0044]
图1为本发明的一种基于一致性采样算法的标靶单应性矩阵计算方法流程图；
[0045]
图2为本发明中置信度为95％的正态分布内外点分布图；
[0046]
图3为设计的九点标靶图；
[0047]
图4为基于卡方χ2检验的平方误差代价函数的阈值k计算方法流程图。
具体实施方式
[0048]
为进一步了解本发明的内容，结合附图和实施例对本发明作详细描述。
[0049]
实施例1
[0050]
结合附图，本实施例的一种基于一致性采样算法的标靶单应性矩阵计算方法，对圆形标靶取点进行优化，剔除三点共线的冗余点，求出单应性矩阵。并且合理设计阈值判断平方误差代价函数。
[0051]
为实现上述目的，本实施例包括如下步骤：
[0052]
步骤一、设计合理标靶
[0053]
假设每个点是内点的概率为w＝内点数目/(内点数目 外点数目)。通常不知道w的具体值，wn是选择的n个点都是内点的概率，1-wn是选择的n个点至少有一个点不是内点的概率，(1-wn)k表示重复k次不是全部的n个点都是内点的概率，假设算法执行k次以后成功的概率是p，p也称为置信度，表示个体对命题真实性相信的程度。经验设置一般置信度设置在[0.95,0.99]的范围内，本实施例中置信度取p＝0.95，表示有多次抽样试验中有95％的可能性认为算法执行k次会成功，用matlab产生二元正态分布随机点数矩阵，绘制样本数据的95％置信椭圆区域，则有图2所示内外点分布图。
[0054]
定性分析n取值，由ransac算法计算公式：
[0055][0056]
本实施例中n取满足计算单应性矩阵h的最小匹配数，可以取n＝9，即标靶取9个点。考虑到h矩阵最多有9个独立变量，理论上最多需要9个点就能计算出单应性矩阵h。k表示计算单应性矩阵h所需实验次数，k与n成正相关，为简化运算复杂度，n取值并非越大越好，本实施例取n＝9。
[0057]
步骤二、通过圆形标靶特征点图像像素坐标和世界坐标构建齐次矩阵方程组剔除误匹配点，进行标靶点非三点共线筛选，为得到稳定的结果，从匹配数据集即步骤一设计的标靶9个点随机抽出5个样本并保证这5个样本之间不共线，具体为：
[0058]
令空间平面测量标靶的3d点世界坐标系对应的坐标为[xi,yi,zi]
t
,在图像像素坐标系对应的2d点坐标为[xi',yi',zi']
t
，两者映射关系表示为：
[0059][0060]
其中，h
ij
(1≤i＝j≤3,∈z)是单应性矩阵分量。
[0061][0062]
u＝[h
11 h
12 h
13
]
t
，v＝[h
21 h
22 h
23
]
t
，w＝[h
31 h
32 h
33
]
t
，由于是齐次坐标系，可令z
′i＝zi＝1。
[0063]
用dlt直接线性变换法同除第三行得:
[0064][0065]
变换写成矩阵形式：
[0066][0067]
知h
t
＝[u
t v
t w
t
]设h
t
＝[u v w]，上矩阵改写为关于未知数u,v,w构成的h的方程：
[0068][0069]
ab＝0
[0070]
即ah＝b
[0071]
一个特征点对应2个单应性矩阵方程即为2个约束条件，故从标靶中随机选取的5个特征点有10个约束单应性矩阵方程。可写为矩阵方程形式：ah＝b，a为图像坐标系坐标值和世界坐标系坐标值组成的系数矩阵，h是单应性矩阵。依据非线性方程组解的情况判断标
靶点是否共线。
[0072]
r(a)是系数矩阵a的秩，r(a|b)是增广矩阵的秩。
[0073]
(1)若r(a)＝r(a|b)《9，非线性方程组为欠定方程。h矩阵有n-r(a)个非零特解，该非零特解任意线性组合是h的解。此时标靶五点有三点共线点，为废靶。
[0074]
(2)若r(a)＝r(a|b)＝9，非线性方程组为超定方程。标靶为正常标靶，可使用5个不共线特征点组成的10个约束单应性矩阵方程求h矩阵。
[0075]
标靶五点坐标约束条件大于h矩阵未知数个数。可求h矩阵
[0076]
超定方程使用线性最小二乘法计算h，构造欧几里2-范数则最小二乘解为由于s尺度因子不变性，设约束条件范数||h||＝t，由拉格朗日乘数法：一阶偏导
[0077][0078]
l'(h)＝2a
t
ah-2λh＝0
[0079]
l'(λ)＝t-h
t
h＝0
[0080]
解得：a
t
ah＝λh，即h是矩阵a
t
a特征值为λ的特征向量，此时欧几里2-范数当特征值λ取最小λ
min
时，最小，此时最小二乘解是a
t
a最小特征值λ
min
对应的特征向量。特征值λ
min
对应的特征向量组成h矩阵。
[0081]
步骤三：通过平方代价函数均方误差计算图像像素坐标系预测值和真实值的误差，并根据阈值k划分第i点属于内点集合还是外点集合。
[0082]
将步骤二中正常标靶中9个点中剔除{(x0,y0),(x3,y3),(x7,y7)}三共线点后，筛选出6个不共线点集，记为m0＝{(x1,y1),(x2,y2),(x4,y4),(x5,y5),(x6,y6),(x8,y8)}，从m0中取五不共线点，共个点集组成的m'＝{m1,m2,m3,m4,m5}，用步骤二中(2)情况下超定方程的最小二乘法，由点集m'＝{m1,m2,m3,m4,m5}计算出单应性矩阵对应的h'＝{h1,h2,h3,h4,h5}，h'称为临时模型参数。用该模型h'去测试图像像素坐标系的匹配点点集，计算平方误差代价函数：
[0083][0084][0085]
k是阈值，若
则该i点可划分到内点集合{θ0}，若则该i点可划分到外点集合{θ1}。
[0086]
步骤四：求k阈值，根据外点是高斯白噪声，由中心极限定理知内外点均服从标准正态分布n～(0,1)，则其均值为0，方差为σ。理论阈值和实际阈值的偏离程度为残差符合n维的χ2分布的特点，其概率密度函数为阈值k采用卡方独立性卡方检验获取，判断阈值k能否有效区分内点和外点。创建一个分类器f:阈值k可划分{θ0}，{θ1}两个内点和外点集合，内点集合{θ0}表示分类器用来刻画理论阈值和实际阈值的偏离程度。可提出零假设p0：“阈值k与内点的数目无关”，在基于零假设p0成立的情况下，创建的分类器取得极大值时，由pearsonχ2计算公式知：此时卡方值也取极大值时，表示理论阈值和实际阈值的偏离程度最大，可以拒绝零假设p0：“阈值k与内点的数目无关”可接受备择假设p1：“阈值k与内点数目有关”。说明阈值k能最有效区分内点和外点。
[0087]
χ2表示理论阈值和实际阈值的偏离程度，a表示内点和外点观察频数，e表示基于p0的计算的期望频数，a与e的差值为残差。阈值k的计算过程如下：
[0088]
(1)零假设p0：“阈值k与内点的数目无关”，备择假设p1：“阈值k与内点数目有关”。
[0089]
(2)阈值k随机存在于[0,n]中，用二分法将阈值范围分成2个部分[0,k]和[k,n]，也视为2个自由度。同时将二分阈值得到的内点和外点两种独立要素也视为2个自由度，判断k值对内点和外点是否有显著区分作为k为阈值的标准。确认联合自由度(2-1)
×
(2-1)＝1。选择显著性水平α＝0.05，服从联合自由度为1的卡方分布，对应卡方分布表
[0090]
有下抽样数据表
[0091][0092]
(3)通过抽样数据表分别计算各行各列理论期望值α1,β1,α2,β2：
[0093]
阈值范围[0,k)对应的外点数目的理论期望值：同理β1,α2,β2也可通过横行总和和竖列总和计算得到。
[0094][0095]
(4)将二分阈值得到的内点和外点两种独立要素也视为2个自由度，联合自由度(2-1)
×
(2-1)＝1，由于自由度1存在偏差，用yates修正计算公式求卡方值χ2。用各行各列理论期望值α1,β1,α2,β2计算ei(i＝1,2,3,4)，用内点和外点观察频数a1(k),a2(k),b1(k),b2(k)计算ai(i＝1,2,3,4)。
[0096][0097]
卡方表的总行或总列累加得：
[0098]
(5)在不同阈值中选出最大表示理论阈值和实际阈值的偏离程度最大，可以拒绝零假设p0：“阈值k与内点的数目无关”；可接受备择假设p1：“阈值k与内点数目有关”。说明阈值k能最有效区分内点和外点，此时为最优阈值k。
[0099]
步骤5：取平方误差代价函数最小值和最优阈值k比较，若满足此时对应的
组成的是最优单应性矩阵。
[0100][0101]
以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的方式及实施例，均应属于本发明的保护范围。