计算机视觉教程1-2：单应性矩阵估计

举报
Mr.Winter 发表于 2022/03/21 21:45:41 2022/03/21
【摘要】 透视空间所有变换都是投影变换的特例，本节进一步研究投影变换矩阵(单应性矩阵)的估计

3a42366089df410590ace510e81229fd.png

1 导论

计算机视觉系列教程1-1：透视空间与透视变换中提到，透视空间所有变换都是投影变换的特例，本节进一步研究投影变换矩阵(单应性矩阵)的估计。

透视变换的核心是单应性矩阵 H H 或单应性向量 h h 。

H = [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] ⇔ h = [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] T H=\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] \Leftrightarrow h=\left[ \begin{matrix} h_{11}& h_{12}& h_{13}& h_{21}& h_{22}& h_{23}& h_{31}& h_{32}& h_{33}\\\end{matrix} \right] ^T

设 p s r c = [ x y 1 ] T p_{src}=\left[ \begin{matrix} x& y& 1\\\end{matrix} \right] ^T 与 p d s t = [ x ′ y ′ 1 ] T p_{dst}=\left[ \begin{matrix} x'& y'& 1\\\end{matrix} \right] ^T 是二维透视空间 P 2 \mathbb{P}^2 中，一次透视变换前后的对应点，因此其满足

p d s t = H p s r c ⟺ [ x ′ y ′ 1 ] = [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] [ x y 1 ] p_{dst}=Hp_{src}\Longleftrightarrow \left[ \begin{array}{c} x'\\ y'\\ 1\\\end{array} \right] =\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] \left[ \begin{array}{c} x\\ y\\ 1\\\end{array} \right]

若将单应性矩阵进行尺度缩放后作用于 p s r c p_{src} ，则

k H p s r c = k [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] [ x y 1 ] = k p d s t kHp_{src}=k\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] \left[ \begin{array}{c} x\\ y\\ 1\\\end{array} \right] =kp_{dst}

而透视空间中， k p d s t kp_{dst} 与 p d s t p_{dst} 实际上对应同一点，因此 k H kH 与 H H 相当于同一次透视变换，故单应性矩阵 H H 仅有8个自由度，通常通过设置 h 33 = 1 h_{33}=1 或 ∥ h ∥ 2 2 = 1 \lVert h \rVert _{2}^{2}=1 来约束冗余的参数。

下面详细阐述单应性矩阵估计方法。

2 基本直接线性变换(Basic DLT)

将上式改写为齐次形式

[ 0 0 0 − x − y − 1 y ′ x y ′ y y ′ x y 1 0 0 0 − x ′ x − x ′ y − x ′ − y ′ x − y ′ y − y ′ x ′ x x ′ y x ′ 0 0 0 ] [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] = [ 0 0 0 ] \left[ \begin{matrix} 0& 0& 0& -x& -y& -1& y'x& y'y& y'\\ x& y& 1& 0& 0& 0& -x'x& -x'y& -x'\\ -y'x& -y'y& -y'& x'x& x'y& x'& 0& 0& 0\\\end{matrix} \right] \left[ \begin{array}{c} h_{11}\\ h_{12}\\ h_{13}\\ h_{21}\\ h_{22}\\ h_{23}\\ h_{31}\\ h_{32}\\ h_{33}\\\end{array} \right] =\left[ \begin{array}{c} 0\\ 0\\ 0\\\end{array} \right]

其中系数矩阵的秩为2，因此一对变换点仅能确定2个自由度。因此需要无三点共线的四对变换点才能确定单应性矩阵 H H 。

[ 0 0 0 − x 1 − y 1 − 1 y 1 ′ x 1 y 1 ′ y 1 y 1 ′ x 1 y 1 1 0 0 0 − x 1 ′ x 1 − x 1 ′ y 1 − x 1 ′ 0 0 0 − x 2 − y 2 − 1 y 2 ′ x 2 y 2 ′ y 2 y 2 ′ x 2 y 2 1 0 0 0 − x 2 ′ x 2 − x 2 ′ y 2 − x 2 ′ ⋮ ] [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] = [ 0 0 0 ] ⇔ A h = 0 \left[ \begin{matrix} 0& 0& 0& -x_1& -y_1& -1& y_{1}^{'}x_1& y_{1}^{'}y_1& y_{1}^{'}\\ x_1& y_1& 1& 0& 0& 0& -x_{1}^{'}x_1& -x_{1}^{'}y_1& -x_{1}^{'}\\ 0& 0& 0& -x_2& -y_2& -1& y_{2}^{'}x_2& y_{2}^{'}y_2& y_{2}^{'}\\ x_2& y_2& 1& 0& 0& 0& -x_{2}^{'}x_2& -x_{2}^{'}y_2& -x_{2}^{'}\\ & & & & \vdots& & & & \\\end{matrix} \right] \left[ \begin{array}{c} h_{11}\\ h_{12}\\ h_{13}\\ h_{21}\\ h_{22}\\ h_{23}\\ h_{31}\\ h_{32}\\ h_{33}\\\end{array} \right] =\left[ \begin{array}{c} 0\\ 0\\ 0\\\end{array} \right] \Leftrightarrow {Ah=0}

对于形如 A x = b Ax=b 的非齐次线性方程可以通过伪逆形式计算 x x ；对于形如 A x = 0 Ax=0 的齐次线性方程的解则对应矩阵 A A 右奇异向量 v p v_p ，其中 v p v_p 对应的奇异值 σ p ≈ 0 \sigma _p\approx 0 或不对应奇异值。

3 归一化直接线性变换(Normalized DLT)

基本DLT估计算法的缺陷是：

(1) 单应性估计不具有相似不变性

假设在第一次估计下有 x d s t = H x s r c x_{dst}=Hx_{src} 。现对两张图像分别进行相似变换并重新进行单应性估计，得到 ( T d s t x d s t ) = H ′ ( T s r c x s r c ) \left( T_{dst}x_{dst} \right) =H'\left( T_{src}x_{src} \right) ，改写为 x d s t = ( T d s t − 1 H ′ T s r c ) x s r c x_{dst}=\left( T_{dst}^{-1}H'T_{src} \right) x_{src} ，大部分情况下 H ≠ T d s t − 1 H ′ T s r c H\ne T_{dst}^{-1}H'T_{src} ，这表明基本DLT算法无法抵抗相似变换的干扰。

(2) 估计的单应性矩阵容易产生病态条件，鲁棒性差

由于默认透视空间的尺度变换因子 w = 1 w=1 ，所以齐次坐标下很可能产生分量幅度差异大的情况，例如某特征点 X = [ 100 101 1 ] T X=\left[ \begin{matrix} 100& 101& 1\\\end{matrix} \right] ^T 。在这种情况下估计出的单应性矩阵，各个元素数值数量级可能会相差 1 0 4 10^4 以上，导致病态条件——特征点的轻微变化都会造成单应性矩阵的剧变。

基于(1)(2)两种缺陷，需要将基本DLT算法进行优化,优化的核心就是特征点坐标的归一化。设原图像特征点集合为 ，目标图像特征点集合为 ，则具体的算法为：

① 将特征点集合 X s r c X_{src} 、 X d s t X_{dst} 归一化

使用相似变换矩阵 T s r c T_{src} 、 T d s t T_{dst} 将特征点集合中心移至原点，且与原点平均距离为 2 \sqrt{2} 。由于默认尺度因子为 w = 1 w=1 ，所以归一化到 2 \sqrt{2} 可以保持齐次坐标的三个分量有相同的幅度，例如 X = [ 100 100 1 ] T ⇒ X n o r m a l = [ 1 1 1 ] T X=\left[ \begin{matrix} 100& 100& 1\\\end{matrix} \right] ^T\Rightarrow X^{normal}=\left[ \begin{matrix} 1& 1& 1\\\end{matrix} \right] ^T 。

这里给出一种相似变换矩阵的计算方式，设

X n o r m a l = T X ⇔ [ x ~ i y ~ i 1 ] = [ s t x s t y 1 ] [ x i y i 1 ] X^{normal}=TX\Leftrightarrow \left[ \begin{array}{c} \tilde{x}_i\\ \tilde{y}_i\\ 1\\\end{array} \right] =\left[ \begin{matrix} s& & t_x\\ & s& t_y\\ & & 1\\\end{matrix} \right] \left[ \begin{array}{c} x_i\\ y_i\\ 1\\\end{array} \right]

令

{ 1 N ∑ i = 1 N x ~ i = 1 N ∑ i = 1 N ( s x i + t x ) = 0 1 N ∑ i = 1 N y ~ i = 1 N ∑ i = 1 N ( s y i + t y ) = 0 1 N ∑ i = 1 N x ~ i 2 + y ~ i 2 = 2 \begin{cases} \frac{1}{N}\sum_{i=1}^N{\tilde{x}_i}=\frac{1}{N}\sum_{i=1}^N{\left( sx_i+t_x \right)}=0\\ \frac{1}{N}\sum_{i=1}^N{\tilde{y}_i}=\frac{1}{N}\sum_{i=1}^N{\left( sy_i+t_y \right)}=0\\ \frac{1}{N}\sum_{i=1}^N{\sqrt{\tilde{x}_{i}^{2}+\tilde{y}_{i}^{2}}}=\sqrt{2}\\\end{cases}

解得

{ t x = − s 1 N ∑ i = 1 N x i = − s x ˉ t y = − s 1 N ∑ i = 1 N y i = − s y ˉ s = 2 1 N ∑ i = 1 N x ~ i 2 + y ~ i 2 = 2 1 N ∑ i = 1 N ( x i − x ˉ ) 2 + ( y i − y ˉ ) 2 \begin{cases} t_x=-s\frac{1}{N}\sum_{i=1}^N{x_i}=-s\bar{x}\\ t_y=-s\frac{1}{N}\sum_{i=1}^N{y_i=-s\bar{y}}\\ s=\frac{\sqrt{2}}{\frac{1}{N}\sum_{i=1}^N{\sqrt{\tilde{x}_{i}^{2}+\tilde{y}_{i}^{2}}}}=\frac{\sqrt{2}}{\frac{1}{N}\sum_{i=1}^N{\sqrt{\left( x_i-\bar{x} \right) ^2+\left( y_i-\bar{y} \right) ^2}}}\\\end{cases}

② 运用基本DLT算法由 与 估计单应性矩阵

③ 解归一化，映射回实际图像

4 鲁棒单应性估计(Robust Homography Estimation)

结合基本DLT算法、归一化DLT算法及计算机视觉系列教程6-2：图像匹配算法概述.的RANSAC算法进行单应性矩阵估计，具体流程如下：

(1) 设置迭代次数 K = ∞ K=\infty ，内点集 S i n = ⊘ S_{in}=\oslash ，模型参数 H = H 0 H=H_0 ；
(2) 随机从样本数据集 S S 中选取4对特征点，并通过基本DLT算法确定测试模型 H t e s t H_{test} ；
(3) 用 H t e s t H_{test} 遍历样本数据集 S S ，估计误差 ε \varepsilon 在距离阈值 t t 内的点加入内点集 S i n S_{in} 。其中阈值 t = 5.99 σ t=\sqrt{5.99}\sigma ， σ \sigma 为估计不确定度，估计误差 ε \varepsilon 主要有两种度量方式：

① 代数误差 ε i = ∥ A i h t e s t ∥ \varepsilon _i=\left\| A_ih_{test} \right\| ，其中

A i = [ 0 0 0 − x − y − 1 y ′ x y ′ y y ′ − y ′ x − y ′ y − y ′ x ′ x x ′ y x ′ 0 0 0 ] A_i=\left[ \begin{matrix} 0& 0& 0& -x& -y& -1& y'x& y'y& y'\\ -y'x& -y'y& -y'& x'x& x'y& x'& 0& 0& 0\\\end{matrix} \right]

② 几何误差(二次投影误差) ε i = ∥ H X s r c , i , X d s t , i ∥ 2 2 + ∥ X s r c , i , H − 1 X d s t , i ∥ 2 2 \varepsilon _i=\left\| HX_{src,i}, X_{dst, i} \right\| _{2}^{2}+\left\| X_{src,i}, H^{-1}X_{dst, i} \right\| _{2}^{2} ，可以视作交叉检验。

(4) 若 S i n S_{in} 的大小小于阈值 T T ，则放弃该模型，重复(2)；若 S i n S_{in} 的大小大于阈值 t t ，则通过归一化DLT算法或Levenberg Marquardt等迭代优化算法，利用 S i n S_{in} 中的所有点重新估计模型 H t e s t ∗ H_{test}^{*} ；
(5) 计算当前模型 H t e s t ∗ H_{test}^{*} 下的内点率 ω = ∣ S i n ∣ ∣ S ∣ \omega =\frac{|S_{in}|}{|S|} ，根据 K = ln ⁡ ( 1 − p ) ln ⁡ ( 1 − ω n ) K=\frac{\ln \left( 1-p \right)}{\ln \left( 1-\omega ^n \right)} 更新迭代次数；
(6) 至此完成一次迭代，若 H t e s t ∗ H_{test}^{*} 下内点率为最大，则更新 H = H t e s t ∗ H=H_{test}^{*} ，重复(2) ~ (5)直至迭代次数满足要求。

🚀 计算机视觉基础教程说明

章号                                    内容
  0               色彩空间与数字成像
  1               计算机几何基础
  2               图像增强、滤波、金字塔
  3               图像特征提取
  4               图像特征描述
  5               图像特征匹配
  6               立体视觉
  7               项目实战

🔥 更多精彩专栏：

👇配套代码 · 优质体验 · 系统知识 请关注👇


【版权声明】本文为华为云社区用户原创内容，未经允许不得转载，如需转载请自行联系原作者进行授权。如果您发现本社区中有涉嫌抄袭的内容，欢迎发送邮件进行举报，并提供相关证据，一经查实，本社区将立刻删除涉嫌侵权内容，举报邮箱： cloudbbs@huaweicloud.com
人工智能 图像处理 机器视觉 模式识别
  • 点赞
  • 收藏
  • 关注作者

评论（0）

0/1000
抱歉，系统识别当前为高风险访问，暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称，即可参与社区互动！

*长度不超过10个汉字或20个英文字符，设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符，设置后3个月内不可修改。

确认 取消

加入云驻计划，成为创作者

  • 华为云周边好礼
  • 免费体验产品
  • 特殊身份标识
  • 线下官方门票
  • 内部专家零距离
  • 与10000+优质创作者共同成长
立即加入