草庐IT

单应性(Homography)变换与单应性矩阵的求解

jiayuzhang128 2023-08-03 原文

单应性(Homography)变换

文章目录


单应性变换的严格数学定义请参考:

  • 《Multiple View Geometry in Computer Vision -2nd Edition》 by Richard Hartley, Andrew Zisserman
  • 第2.3节Projective transformations

1. 概念

  • 单应性变换又叫投影变换:应用在平面坐标变换中:

    • 平面投影变换是在三元素向量的齐次坐标下进行的线性变换,他由一个3×3的非奇异变换矩阵 H H H表示,具体如下:
    • x ′ x^{'} x x x x都是齐次坐标,一般的有: x 3 ′ = x 3 = 1 x_{3}^{'}=x_{3}=1 x3=x3=1
  • 单应性矩阵

    • 单应矩阵描述两个平面上的对应点之间的变换关系
    • 同一个平面在任意坐标系之间都可以建立单应性变换关系;
    • 如(a):plannar surface上的X点可以通过单应性矩阵 H 1 H_{1} H1 H 2 H_{2} H2变换到image1image2,(b)和(c)同理。

2. 在CV方面的应用

  • 图像校正、图像拼接、相机位姿估计、视觉SLAM等。

  • 图像校正

  • 视角变化

  • 图像拼接

  • AR

3. 求解单应性矩阵

3.1 假设

  • 首先,我们假设两张图像中的对应点对齐次坐标为 P x y = ( x , y , 1 ) T P_{xy}=(x,y,1)^{T} Pxy=(x,y,1)T P u v = ( u , v , 1 ) T P_{uv}=(u,v,1)^{T} Puv=(u,v,1)T,单应矩阵H定义为:
    • H = [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] H=\begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} H=h11h21h31h12h22h32h13h23h33
    • [ u v 1 ] = s [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] [ x y 1 ] ( 1 ) \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} =s \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix}\begin{bmatrix} x \\ y \\ 1 \end{bmatrix}\qquad\qquad\qquad\qquad\qquad\qquad(1) uv1=sh11h21h31h12h22h32h13h23h33xy1(1)

3.2 性质

  • 这里使用的是齐次坐标系,也就是说可以进行任意尺度的缩放(s为尺度因子),也就是说把H乘以任意一个非零常数k并不改变上式结果,无非就是尺度因子s有所改变。
    • 比如我们把H乘以 1 h 33 \frac{1}{h_{33}} h331可以得到:
      • H ′ = 1 h 33 H = 1 h 33 [ h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ] = [ h 11 ′ h 12 ′ h 13 ′ h 21 ′ h 22 ′ h 23 ′ h 31 ′ h 32 ′ 1 ] H'=\frac{1}{h_{33}}H=\frac{1}{h_{33}}\begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix}=\begin{bmatrix} h_{11}' & h_{12}' & h_{13}' \\ h_{21}' & h_{22}' & h_{23}' \\ h_{31}' & h_{32}' & 1 \end{bmatrix} H=h331H=h331h11h21h31h12h22h32h13h23h33=h11h21h31h12h22h32h13h231
      • 而上述等式依然成立
    • 同理H乘以 1 ∑ i = 1 3 ∑ j = 1 3 h i j 2 \frac{1}{\sqrt{\sum_{i=1}^3\sum_{j=1}^3h_{ij}^{2}}} i=13j=13hij2 1 可以得到约束:
      • h 11 ′ 2 + h 12 ′ 2 + h 13 ′ 2 + h 21 ′ 2 + h 22 ′ 2 + h 23 ′ 2 + h 31 ′ 2 + h 32 ′ 2 + h 33 ′ 2 = 1 h_{11}'^{2} + h_{12}'^{2} + h_{13}'^{2} + h_{21}'^{2} + h_{22}'^{2} + h_{23}'^{2} + h_{31}'^{2} + h_{32}'^{2} + h_{33}'^{2}=1 h112+h122+h132+h212+h222+h232+h312+h322+h332=1
      • ∑ i = 1 3 ∑ j = 1 3 h i j ′ 2 = 1 \sum_{i=1}^3\sum_{j=1}^3h_{ij}'^{2}=1 i=13j=13hij2=1
    • 由此我们可以得出单应性矩阵有8个自由度,并非9个自由度。

3.3 求解

  • 由3.1中假设公式(1)可得:
    • { u = h 11 x + h 12 y + h 13 h 31 x + h 32 y + h 33 v = h 21 x + h 22 y + h 23 h 31 x + h 32 y + h 33 \left\{\begin{array}{ll} u=\frac{h_{11}x + h_{12}y + h_{13}}{h_{31}x + h_{32}y + h_{33}} \\ v=\frac{h_{21}x + h_{22}y + h_{23}}{h_{31}x + h_{32}y + h_{33}}\end{array} \right. {u=h31x+h32y+h33h11x+h12y+h13v=h31x+h32y+h33h21x+h22y+h23,进一步变换得:
    • { u ( h 31 x + h 32 y + h 33 ) = h 11 x + h 12 y + h 13 v ( h 31 x + h 32 y + h 33 ) = h 21 x + h 22 y + h 23 \left\{\begin{array}{ll} u(h_{31}x + h_{32}y + h_{33})=h_{11}x + h_{12}y + h_{13} \\ v(h_{31}x + h_{32}y + h_{33})=h_{21}x + h_{22}y + h_{23}\end{array} \right. {u(h31x+h32y+h33)=h11x+h12y+h13v(h31x+h32y+h33)=h21x+h22y+h23,进一步得到:
    • { x h 11 + y h 12 + h 13 − u x h 31 − u y h 32 − u h 33 = 0 x h 21 + y h 22 + h 23 − v x h 31 − v y h 32 − v h 33 = 0 \left\{\begin{array}{ll} xh_{11} + yh_{12} + h_{13} - uxh_{31} - uyh_{32} - uh_{33}=0\\ xh_{21} + yh_{22} + h_{23} - vxh_{31} - vyh_{32} - vh_{33}=0\end{array} \right. {xh11+yh12+h13uxh31uyh32uh33=0xh21+yh22+h23vxh31vyh32vh33=0,化成矩阵形式有:
    • [ x y 1 0 0 0 − u x − u y − u 0 0 0 x y 1 − v x − v y − v . . . . . . . . . . . . . . . . . . . . . ] [ h 11 h 12 h 113 h 21 h 22 h 23 h 31 h 32 h 33 ] = 0 \begin{bmatrix} x & y& 1& 0&0&0& - ux& -uy& -u \\ 0&0&0&x & y& 1& - vx& -vy& -v \\ ...\\...\\...\\...\\...\\...\\... \end{bmatrix}\begin{bmatrix} h_{11} \\ h_{12} \\ h_{113}\\h_{21}\\h_{22}\\h_{23}\\h_{31}\\h_{32}\\h_{33} \end{bmatrix}=0 x0.....................y0100x0y01uxvxuyvyuvh11h12h113h21h22h23h31h32h33=0,更进一步抽象:
    • [ P x y T 0 ⃗ T − u P x y T 0 ⃗ T P x y T − v P x y T ] h = 0 \begin{bmatrix} P_{xy}^{T} & \vec{0}^{T} & - uP_{xy}^{T} \\ \vec{0}^{T} &P_{xy}^{T} &- vP_{xy}^{T} \end{bmatrix}h=0 [PxyT0 T0 TPxyTuPxyTvPxyT]h=0
      • 其中 h = [ h 11 h 12 h 113 h 21 h 22 h 23 h 31 h 32 h 33 ] T h=\begin{bmatrix} h_{11} & h_{12}& h_{113}&h_{21}&h_{22}&h_{23}&h_{31}&h_{32}&h_{33} \end{bmatrix}^{T} h=[h11h12h113h21h22h23h31h32h33]T
      • 0 ⃗ T = [ 0 0 0 ] \vec{0}^{T}=\begin{bmatrix} 0&0&0 \end{bmatrix} 0 T=[000]
    • 我们发现一对点可以提供两个方程
    • 由于单应矩阵H包含了 ∣ ∣ H ∣ ∣ F = ∑ i = 1 3 ∑ j = 1 3 h i j ′ 2 = 1 ||H||_{F}=\sum_{i=1}^3\sum_{j=1}^3h_{ij}'^{2}=1 HF=i=13j=13hij2=1或者 h 33 = 1 h_{33}=1 h33=1约束,因此根据上式的线性方程组8自由度的H我们至少需要4对点才能计算出单应矩阵。

4. 优化

但是,以上只是理论推导,在真实的应用场景中,我们计算的点对中都会包含噪声。比如点的位置偏差几个像素,甚至出现特征点对误匹配的现象,如果只使用4个点对来计算单应矩阵,那会出现很大的误差。因此,为了使得计算更精确,一般都会使用远大于4个点对来计算单应矩阵。另外上述方程组采用直接线性解法通常很难得到最优解,所以实际使用中一般会用其他优化方法如奇异值分解、Levenberg-Marquarat(LM)算法等进行求解,具体的优化过程可以参考张正友相机标定中得优化方法

5. 推荐阅读

(超详细)张氏标定法原理及公式推导

有关单应性(Homography)变换与单应性矩阵的求解的更多相关文章

  1. 旋转矩阵的几何意义 - 2

    点向量坐标矩阵的几何意义介绍旋转矩阵的几何含义之前,先介绍一下点向量坐标矩阵的几何含义点:在一维空间下就是一个标量,如同一条直线上,以任意某一个位置为0点,以一定的尺度间隔为1,2,3...,相反方向为-1,-2,-3...;如此就形成了一维坐标系,这时候任何一个点都可以用一个数值表示,如点p1=5,即即从原点出发沿着x轴正方向移动5个尺度;点p2=-3,负方向移动3个尺度;     在一维坐标系上过原点做垂直于一维坐标系的直线,则形成了二维坐标系,此时描述一个点需要两个数值来表示点p3=(3,2),即从原点出发沿着x轴正方向移动3个尺度,在此基础上沿着y轴正方向移动两个尺度的位置就是点p3。

  2. 华为OD机试真题 C++ 实现【带传送阵的矩阵游离】【2023 Q2 | 200分】 - 2

            所有题目均有五种语言实现。C实现目录、C++实现目录、Python实现目录、Java实现目录、JavaScript实现目录题目n行m列的矩阵,每个位置上有一个元素你可以上下左右行走,代价是前后两个位置元素值差的绝对值.另外,你最多可以使用一次传送阵(只能从一个数跳到另外一个相同的数)求从走上角走到右下角最少需要多少时间。输入描述:第一行两个整数n,m,分别代表矩阵的行和列。后面n行,每行m个整数,分别代表矩阵中的元素。输出描述:一个整数,表示最少需要多少时间。

  3. 欧拉角表示的姿态矩阵(313和312转序) - 2

    一、习惯约定图片来自PSINS(高精度捷联惯导算法)PSINS工具箱入门与详解.pptx二、基本旋转矩阵绕x轴逆时钟旋转α\alphaα角度Rx(α)=[ 1000cos⁡αsin⁡α0−sin⁡αcos⁡α]R_x(\alpha)=\begin{bmatrix}\1&0&0\\0&\cos\alpha&\sin\alpha\\0&-\sin\alpha&\cos\alpha\end{bmatrix}Rx​(α)=​ 100​0cosα−sinα​0sinαcosα​​绕y轴逆时钟旋转α\alphaα角度Ry(α)=[ cos⁡α0−sin⁡α010sin⁡α0cos⁡α]R_y(\alpha

  4. 欧拉角、旋转矩阵及四元数 - 2

    欧拉角、旋转矩阵及四元数1.简介2.欧拉角2.1欧拉角定义2.2右手系和左手系2.3转换流程3.旋转矩阵4.四元数4.1四元数与欧拉角和旋转矩阵之间等效变换4.2测试Matlab代码5.总结1.简介常用姿态参数表达方式包括方向余弦矩阵、欧拉轴/角参数、欧拉角、四元数以及罗德里格参数等。高分辨率光学遥感卫星主要采用欧拉角与四元数对姿态参数进行描述。这里着重讲解欧拉角、旋转矩阵和四元数。2.欧拉角2.1欧拉角定义欧拉角是表征刚体旋转的一种方法之一,由莱昂哈德·欧拉引入的三个角度,用于描述刚体相对于固定坐标系的方向。在摄影测量、空间科学或其它技术领域,一般用一组(三个)欧拉角描述两个空间坐标之间的旋

  5. ruby - 如何修改矩阵(Ruby std-lib Matrix 类)? - 2

    我理解RubystdlibMatrix是不可修改的,也就是说,例如。m=Matrix.zero(3,4)不会写m[0,1]=7但我非常想做...我可以用笨拙的编程来做,比如defmodify_value_in_a_matrix(matrix,row,col,newval)ary=(0...m.row_size).map{|i|m.rowi}.map(&:to_a)ary[row][col]=newvalMatrix[*ary]end...或者作弊,比如Matrix.send:[]=,0,1,7但我想知道,这一定是人们一直遇到的问题。有没有一些标准的、习惯的方法可以做到这一点,而不必使用

  6. 线性代数让我想想:快速求三阶矩阵的逆矩阵 - 2

    快速求三阶矩阵的逆矩阵前言一般情况下,我们求解伴随矩阵是要注意符号问题和位置问题的(如下所示)A−1=1[  ][−[  ]−[  ]−[  ]  −[  ]]=A−1=1[  ][   M11−[M12]   M13−[M21]   M22−[M23]     M31−[M32]   M33]⊤\begin{aligned}&A^{-1}=\frac{1}{[\\]}\left[\begin{array}{cccccc}&-[\\]&\\-[\\]&&-[\\]\\\\&-[\\]&\\\end{array}\right]=\\\\&A^{-1}=\frac{1}{[\\]}\left[\b

  7. 相机校准—外参矩阵 - 2

    在本文中,我们将探讨摄影机的外参,并通过Python中的一个实践示例来加强我们的理解。相机外参摄像头可以位于世界任何地方,并且可以指向任何方向。我们想从摄像机的角度来观察世界上的物体,这种从世界坐标系到摄像机坐标系的转换被称为摄像机外参。那么,我们怎样才能找到相机外参呢?一旦我们弄清楚相机是如何变换的,我们就可以找到从世界坐标系到相机坐标系的基变换的变化。我们将详细探讨这个想法。具体来说,我们需要知道相机是如何定位的,以及它在世界空间中的位置,有两种转换可以帮助我们:有助于确定摄影机方向的旋转变换。有助于移动相机的平移变换。让我们详细看看每一个。旋转通过旋转改变坐标让我们看一下将点旋转一个角度

  8. ruby - Ruby 中的有限矩阵 - 2

    为什么Matrix类没有方法来编辑它的向量和组件?似乎矩阵中的所有内容都可以读取但不能写入。我错了吗?是否有一些类似于Matrix的第三方优雅类允许我删除行并有意地编辑它们?如果没有这样的类(class),请通知我——我将停止搜索。 最佳答案 Matrix类的设计者一定是不可变数据结构和函数式编程的爱好者。是的,你是对的。无论如何,总有一个简单的解决方案可以满足您的需求。使用Matrix它可以做的事情,然后,只需使用.to_a来获得一个真正的数组。>>Matrix.identity(2).to_a=>[[1,0],[0,1]]另见N

  9. ruby - 从文件中提取快速傅里叶变换数据 - 2

    我正在构建一个应该在服务器上运行并分析声音文件的工具。我想在Ruby中执行此操作,因为我的所有其他工具也是用Ruby编写的。但我很难找到完成此任务的好方法。我发现的很多例子都是在做可视化和图形化的东西。我只需要FFT数据,仅此而已。我既需要获取音频数据,又需要对其进行FFT。我的最终目标是计算一些东西,例如所有频率(加权幅度)的均值/中值/众数、第25个百分位数和第75个百分位数、BPM,也许还有其他一些好的特性,以便以后能够将相似的声音聚集在一起.首先,我尝试使用ruby-audio和fftw3,但我从未将两者真正结合使用。文档也不好,所以我真的不知道有什么数据被洗牌了。接下来,我尝

  10. ruby - 在 Ruby 中打印可读矩阵 - 2

    在Ruby中是否有内置的打印可读矩阵的方法?例如require'matrix'm1=Matrix[[1,2],[3,4]]printm1让它显示=>1234在REPL中代替:=>Matrix[[1,2][3,4]]matrix的Ruby文档让它看起来像应该显示的那样,但这不是我所看到的。我知道编写一个函数来执行此操作是微不足道的,但如果有“正确”的方法,我宁愿学习! 最佳答案 您可以将其转换为数组:m1.to_a.each{|r|putsr.inspect}=>[1,2][3,4]编辑:这是一个“无积分”版本:putsm1.to_a

随机推荐