1. 单应性矩阵的基本概念
什么是单应性矩阵?单应性变换包含什么样的射影组合(projective transformation)?
单应性关系的前提条件?
单应性与极几何的联系?
2. 单应性矩阵的计算
本质矩阵和基础矩阵的性质,上一篇博文有详细介绍,所以此处只讲计算方法了。
3. 基础矩阵的计算及计算前提条件
4.本质矩阵的计算
5.给定本质矩阵,如何提取相机朝向R和相对位移关系(基线方向T)
1. 单应性矩阵的基本概念
什么是单应性矩阵?
单应性是几何中的一个概念。单应性是一个从实射影平面到射影平面的可逆变换,直线在该变换下仍映射为直线。它是一对透视投影的组合射影变换群。它描述了当观察者视角改变时,被观察物体的感知位置会发生何种变化。射影变换并不保持大小和角度,但会保持重合关系和交比。
在计算机视觉领域中,空间中同一平面的任意两幅图像通过单应性关联在一起(假定是针孔相机)。比如,一个物体可以通过旋转相机镜头获取两张不同的照片(这两张照片的内容不一定要完全对应,部分对应即可),我们可以把单应性设为一个二维矩阵M,那么照片1乘以M就是照片2. 这有着很多实际应用,比如图像校正、图像对齐或两幅图像之间的相机运动计算(旋转和平移)等。一旦旋转和平移从所估计的单应性矩阵中提取出来,那么该信息将可被用来导航或是把3D物体模型插入到图像或视频中,使其可根据正确的透视来渲染,并且成为原始场景的一部分。
假设图A中有一个点x的坐标为(u,v,1),图B中对应匹配点x’坐标为(u’,v’,1),那么单应性矩阵可以将两个坐标联系起来:
注意,这种单应性联系不依赖于观察场景结构,即无论相机拍摄什么画面,两个画面的单应性关系保持不变。
单应性变换包含什么样的射影组合(projective transformation)?
首先看等距变换Isometries transformation,它相当于是平移变换和旋转变换的复合。
上式也可以写作:
上式中ε等于正负1,当ε=1时,等距变换与原结构保持相同方向,即是欧氏变换,表达刚性物体运动;当ε=-1时,等距变换是逆向的。另外,R是正交矩阵。
等距变换有三个自由度,一个旋转,两个平移。等距变换后目标结构的长度、角度、面积与原结构相同。
再看相似变换(Similarity Transformation),它相当于是等距变换和均匀缩放的一个复合。
相似变换可以从两组匹配点计算出来,它有四个自由度,一个尺度,一个旋转,一个平移。相似变换后目标结构内部的角度、长度比、面积比保持不变。R是正交矩阵
再看仿射变换(Affine Transformation),它相当于一个平移变换和一个非均匀变换的复合。
仿射矩阵A是2阶非奇异矩阵,A通过奇异值分解为旋转rotation和非均匀缩放deformation:
上式中U和V都是正交矩阵,对A简单描述下:λ1和λ2代表旋转后的目标在x和y方向的缩放比例,所以是非均匀的。Φ代表缩放的方向,θ代表旋转的方向。
仿射变换中A有六个自由度,2个尺度,2个旋转,2个平移。经过仿射变换后,原来平行的线还是平行的,并且平行线段的长度之比保持不变,面积之比也是不变的。另外仿射变换是否保向,取决于A的行列式值的正负。
最后看射影变换projectivetransformation,也就是单应性矩阵,他表达齐次坐标的非奇异线性变换。
v是射影变换后平行直线可能不再平行的原因。
由于齐次坐标对于尺度具有不变性,可以从H中提取该尺度(比如v),因此单应性矩阵有8个自由度。两个平面之间的射影关系,可以通过4组匹配点计算得到(注意任意3点不可共线)。经过射影变换后,目标结构的交比cross ratio保持不变(两个平面上分别有4个点,平面1上线段AB与线段CD之比,等于平面2上线段A’B’与线段C’D’之比)。另外变换前后共点,共线,相切,拐点,切线的不连续性和岐点保持不变。
射影变换可以看作是一些列变换的组合:
上三角矩阵K的行列式值为1,v不等于0,如果s的正负确定,则H是唯一的。
HP (2 dof) moves theline at infinity; HA (2 dof) affects the affine properties, but does not movethe line at infinity; and finally, HS is a general similarity transformation (4dof) which does not affect the affine or projective properties.
总结:
- 射影变换后,面积的缩放与位置有关,直线变换方向也与面积有关。
- 射影变换是不可逆的。
- 仿射变换位于射影变换和相似变换之间,仿射变换推广相似变换使得夹角不再保持不变,造成物体形状在变换后产生歪斜,他在变换中,直线只与方向有关,与位置无关,不会像射影变换那样,越远越小
下图中参数数量与上面不符,因为他针对的是3D点,上面是2D点
单应性关系的前提条件?
如果两幅图像之间的相机运动只有旋转而没有平移,那么这两幅图像通过单应性关联在一起(假定是针孔相机)。也就是说(if and only if):
a. Both images are viewing thesame plane from a different angle
b. Both images are taken from thesame camera but from a different angle.( Camera is rotated about its center ofprojection without any translation)
也就是说:
a.当两个相机所代表的视平面(image plane)相互平行,stereo normal case,单应性矩阵才是有效的。(baseline不为0,但是极点在无穷远)
b.当且仅当相机是绕相机中心旋转,没有偏移,baseline什么的都是零,单应性矩阵才是有效的。(baseline为0)
单应性与极几何的联系?
下图为极几何的经典视图,相机A的中心为O1,相机B的中心为O2,它们共同看到空间一点P。根据极几何的性质,空间中所有落在直线O1P上的点,最终都会落在在Image Plane B的极线上。
根据单应性的前提条件,单应性与极几何的区别在于两个相机主点之间的距离为0时,
总结起来:
当相机之间存在偏移,那么该用极几何关系描述点的对应关系,此时两个图片之间点的匹配依赖于物点的深度,即物点到相机的距离。
It maps a point in one image to a line in the other image. Actualmatching point depends on the depth of that point,that is its distance from thecamera
当相机之间不存在偏移时,极几何关系不再适用,该用单应性描述点的对应关系,此时两个图片之间点的匹配不依赖于物点的深度。(此时三角方程不成立)
It maps a point in one image to a point in the other image. Themapping does not depend on the depth (distance) of point to the camera. Thismakes sense because baseline z zero means ordinary stereo triangulationequations fail.
从下图看出相机拍摄时并未有偏移,故通过单应性做的图像拼接。
2. 单应性矩阵的算法
从单应性矩阵的定义出发:
除此之外,我们可以通过归一化坐标,使得h的结果是numerical stable。但是即使这样,这种做法对于噪声非常敏感,有几个outlier就完蛋了。
而第二种方法就很好,
为了使上式为0,咱只要对左侧ATA进行奇异值分解,挑选奇异值前九个最小的v的列向量作为h即可,SVD:
Let h be the column ofV (unit eigenvector) associated with the smallest eigenvalue in D.(if only 4points, that eigenvalue will be 0)
Matlab代码:
- function H = find_H(P1, P2)
- % H = find_H(P1, P2)
- %
- % P1 are coordinates from source frame ([n*2], n>= 4)
- % P2 are coordinates from destination frame ([n*2], n>= 4)
- % H is the 3*3 homography matrix from P1 to P2
- % [P2; 1] ~ H [P1; 1]
- %
- x = P1(:, 1);
- y = P1(:, 2);
- X = P2(:, 1);
- Y = P2(:, 2);
- A = zeros(length(x(:))*2,9);
- for i = 1:length(x(:)),
- a = [x(i),y(i),1];
- b = [0 0 0];
- c = [X(i);Y(i)];
- d = -c*a;
- A((i-1)*2+1:(i-1)*2+2,1:9) = [[a b;b a] d]; %拼接矩阵A,自己对照截图,一次两行,总共2N行
- end
- [U S V] = svd(A);
- h = V(:,9);
- H = reshape(h,3,3)';
3. 基础矩阵F的计算及前提条件
分别取[xn’,yn’,1],[xn’’,yn’’,1]与F的第一个元素,得到展开式
按照下图定义an和f,得到线性方程
我们同样使用奇异值分解去求Af=0这个方程。由于基础矩阵F是齐次的,矩阵A的秩不大于8,因此我们只要8对匹配点就可以解该方程。实际场景中有可能有不止8对匹配点,而且这些匹配点中存在噪声,矩阵A可能是满秩的(should not be)。此时,我们把A中最小的奇异值对应的奇异向量作为基础矩阵F的解。
通过下图可以清晰得看到SVD全貌,即V的最后一列(对应D的最小的九个奇异值)是F的解。但是这还不够,现在并不能保证F的秩等于2(否则Af全为0了,不符合噪声实际情况)。
现在的目标是使r(F)=2,并且F应该与我们的估计F’无限接近。怎么做?
我们对F本身进行第二次奇异值分解,从F’的奇异值中挑最小的,并将它设置为0.也就是说,根据上一步得到的F’,我们强制修改F’的奇异值,而不改变U’,V’,使得rank(F’)=2。
Matlab代码:
- function F = F_from_point_pairs(xs, xss)
- %xs, xss: Nx3 homologous point coordinates, N>7
- %F: 3x3 fundamental matrix
- %coefficient matrix
- for n = 1 : size(xs, 1)
- A(n, :) = kron(xss(n, :), xs(n, :));
- end
- %singular value decomposition
- [U, D, V] = svd(A);
- %approximate F, possibly regular
- Fa = reshape(V(:,9), 3, 3)';
- % svd decomposition of F
- [Ua, Da, Va] = svd(Fa);
- %algebraically best F, singular
- F = Ua * diag([Da(1,1), Da(2,2), 0]) * Va';
注:
在求单应性矩阵时扯到过坐标归一化,防止numerically unstable,这里详细讲下。
数值不稳定是指,比如对于一副3000*4000的图片,由于其图片像素坐标过大,使得计算结果不稳定。我们按照下图所示,将原点移动到图像中心,将像素点坐标约束到(-1,1)。
坐标归一化,会改变真实F的值,那么坐标归一化前后基础矩阵的对应关系为:
设存在变换T使得坐标x归一化
哪些情况下,基础矩阵的F会很不稳定?
a.如果所有匹配点对应的物点均在同一平面上,即rank(A)<8,此时Af=0无解。就好比很难对墙面上的画做三维重建
b.如果两个相机的光心重合,即相机没有移动,纯粹旋转,则skew-symmetric矩阵全为0,此时需要考虑使用单应性关系。
4. 本质矩阵E的计算
我们可以模仿之前求基础矩阵的方法求本质矩阵,不过此时E的约束关系除了秩必须为2,还要求E的奇异值均为1。
Matlab代码:
- function E = E_from_point_pairs(xs, xss)
- %xs, xss: Nx3 homologous point coordinates, N>7
- %E: 3x3 fundamental matrix
- %coefficient matrix
- for n = 1 : size(xs, 1)
- A(n, :) = kron(xss(n, :), xs(n, :));
- end
- %singular value decomposition
- [U, D, V] = svd(A);
- %approximate F, possibly regular
- Ea = reshape(V(:,9), 3, 3)';
- % svd decomposition of E
- [Ua, Da, Va] = svd(Ea);
- %algebraically best F, singular
- E = Ua * diag([1, 1, 0]) * Va';
注意:
上面代码函数的输入参数,不是图像坐标,而是乘以相机内参后的图像坐标,xs=K*xs
坐标归一化同样使用此处
5. 给定本质矩阵E,如何计算相对位移关系及相机朝向
给定本质矩阵,我们可以提取R和基线b方向T,但是解并不唯一,实际上有四种!!
只有左上角这种情况是我们需要的(右上角是绕b旋转了180度,不对的哈),衡量正确与否的标准是:根据R和T复原的点,必须同时位于两个相机的前方。接下来具体展开如何求R和T。
根据E的奇异值分解:
上述做法是从E的两种定义上套用的,现在很轻松的写出基线方向T和旋转矩阵R的表达式。然后这还没完,因为Z和W有尼玛4种组合:
因此E有四种结果,而基线方向T和旋转矩阵R分别有两种情况:
我们测试所有基线方向T和旋转矩阵R组合,还原场景3D坐标(并不是真实尺度),判断所有点是否都位于两个相机的前方,确定唯一的T和R