4 Estimation – 2D Projective Transformations
本章主要估计这么几种2D投影矩阵:
- 2D齐次矩阵,就是从一个图像中的点到另外一个图像中的点的转换,由于点的表示都是齐次的,所以叫齐次矩阵
- 3D到2D的摄像机矩阵
- 基本矩阵
- 三视图之间的转换矩阵
形式化的表示,我们已知 x i ↔ x i ′ x_i \leftrightarrow x'_i xi↔xi′,我们需要计算一个 3 × 3 3 \times 3 3×3 的 H H H,满足 x i ′ = H x i x'_i=H x_i xi′=Hxi。
需要多少对点 第一个问题就是,计算这个projective transformation H矩阵需要多少对对应点。 H H H是 3 × 3 3 \times 3 3×3 有9个数,但是因为尺度等价性,只需要解出8组方程,都用最后1个数字表示就可以了,因此是8个自由度。每个点对点对应关系都占两个约束,因为对于第一个图像中的每个点,第二个图像中的点必须对应于映射的点。 2D点具有与X和Y相对应的两个自由度。 另外,该点是齐次向量,有3个值,由于尺度是任意的,因此它也具有两个自由度。所以最少需要4对点。
估计算法的种类 我们考虑到图像中的对应点都是由噪声的,那么只有4个点其实估计的结果就不准确了,如果我们找出多于4对点,那么我们就需要一个损失函数,并且把它最小。下文介绍的主要分成2类,一类是代数误差,一类是几何误差。
黄金标准算法 黄金标准算法不是一个具体的算法。它是一个评判标准。通俗来讲:假设我们有n个损失函数,哪一个是最好的?让 H H H能达到最小的损失函数就是最好的。那么哪个损失函数让 H H H最小?在“估计2D投影矩阵”这个问题里,最好的损失函数就是书中4.8式。优化最优损失函数的算法就叫黄金标准算法,其他损失函数给出的结果都是和4.8式的结果来比较的。
文章目录
- 4 Estimation – 2D Projective Transformations
- 4.1 The Direct Linear Transformation (DLT) algorithm
- 4.1.1 Over-determined solution
- 4.1.2 Inhomogeneous solution
- 4.1.3 Degenerate configurations
- 4.1.4 Solutions from lines and other entities
- 4.2 Different cost functions
- 4.2.1 Algebraic distance
- 4.2.2 Geometric distance
- 4.2.3 Reprojection error – both images
- 4.2.4 Comparison of geometric and algebraic distance
- 4.2.5 Geometric interpretation of reprojection error
- 4.2.6 Sampson error
- 4.2.7 Another geometric interpretation
- 4.3 Statistical cost functions and Maximum Likelihood estimation
- 4.4 Transformation invariance and normalization
- 4.4.1 Invariance to image coordinate transformations
- 4.4.3 Invariance of geometric error
- 4.4.4 Normalizing transformations
- 4.5 Iterative minimization methods
- 4.6 Experimental comparison of the algorithms
- 4.7 Robust estimation
- 4.7.1 RANSAC
- 4.7.2 Robust Maximum Likelihood estimation
- 4.7.3 Other robust algorithms
- 4.8 Automatic computation of a homography
4.1 The Direct Linear Transformation (DLT) algorithm
原理很简单,既然 x i ′ = H x i x'_i=Hx_i xi′=Hxi,那么就可以建立方程 x i ′ × H x i = 0 x'_i \times Hx_i=0 xi′×Hxi=0,把这个式子改写成 A i h = 0 A_{i}h=0 Aih=0,找4对点解方程就行。我们想要找到一个 h h h的非零解(显然零解对我们来说是没有意义的)。
4.1.1 Over-determined solution
假设A依旧是rank为8且在1维null-space,但是我们找出了多于4对点,那么怎么利用上多出来的信息?事实上,我们实际得到的点是有噪声的,我们可以考虑获得一个 h h h的近似解。那么我们需要最优化一个目标函数。问题来了:我们需要一个啥样的目标函数呢?首先,为了得到一个非零解,我们需要额外加一个约束,那就是满足 ∣ ∣ h ∣ ∣ = 1 ||h||=1 ∣∣h∣∣=1。然后,对于目标函数,最直观的想法就是优化 A h Ah Ah。这个解是 A T A A^{T}A ATA的(unit) eigenvector有最小特征值的那个。等效地,也是对应于 A A A的最小单数值的单位奇异矢量。
4.1.2 Inhomogeneous solution
除了使用齐次向量直接求解H之外,一种替代方法是将一组方程(4.3)转换为一组不齐次的线性方程组,通过对某些 h h h施加条件 h j = 1 h_{j} = 1 hj=1。这个解法不稳定,因为是非齐次坐标,丢失了尺度这个因子。
4.1.3 Degenerate configurations
关于退化的情况。如果 x 1 , x 2 , x 3 x_1,x_2,x_3 x1,x2,x3共线且 x 1 ′ , x 2 ′ , x 3 ′ x'_1,x'_2,x'_3 x1′,x2′,x3′也共线,那么 H H H的解不是唯一的。如果 x 1 , x 2 , x 3 x_1,x_2,x_3 x1,x2,x3共线但是 x 1 ′ , x 2 ′ , x 3 ′ x'_1,x'_2,x'_3 x1′,x2′,x3′不共线,这种情况下 H H H不存在。但是我们用上述方程确实可以解出来一个 H H H,那么这个 H H H是什么?代表什么呢?书中给出的答案是: H H H是一个奇异矩阵,根本不代表任何映射。
4.1.4 Solutions from lines and other entities
除了点对应,我们也可以找线对应。比如我们可以用线对应来构造一个方程 l i = H T l i ′ l_i = H^T l'_i li=HTli′。那么如果用线对应,需要找多少对对应线呢?一般规则是,约束的数量必须等于或超过转换自由度的数量。所以我们需要的约束应该多于 H T H^T HT的自由度。所以4对对应线足够了(每一对是2个约束, 4 × 2 = 8 4 \times 2 =8 4×2=8)。如果我们考虑3D空间,有15个自由度,因为一对对应点或平面提供3个约束,那么求解3D的 H H H用5对就够。
如果我们把点和线混合起来使用,就要小心了。比如说,2对对应点+2对对应线是不能确定
H
H
H的,但是3对对应线+1对对应点或者1对对应线+3对对应点。这是因为,三线和一个点的情况在几何上等同于四个点,因为这三线定义了三角形,三角形的顶点独特地定义了三个点。同样,三个点和一条线的情况相当于四个线,而四个线的对应关系唯一地决定了H。 但是,正如图4.1,两个点和两条线的情况相当于五条线,有四个相交,或五个有四个共线。 如上一节所示,此配置是退化的。
4.2 Different cost functions
4.2.1 Algebraic distance
既然 x ′ = H x x'=Hx x′=Hx,那么我们可以优化 x ′ x' x′和 H x Hx Hx的距离。DLT算法可最大程度地减少 A h Ah Ah。 向量 ε = A h ε= Ah ε=Ah称为residual vector,它的norm也就是待最小化的误差。
4.2.2 Geometric distance
接下来,我们将根据图像中几何距离的测量,以及测得的图像坐标和估计图像坐标之间的差异,讨论替代的误差函数。几何距离考虑的是理论点与实际点的误差。
符号说明 我们使用向量 x x x表示测得的图像坐标, x ^ \hat{x} x^表示点的估计值, x ˉ \bar{x} xˉ表示点的真实值。
我们先考虑只有第二个图像中有error的情况,假设第一个图像上的点都是测得准确的。例如,估计校准模式或世界平面之间的投射转换,该点可以测量到非常高的精度及其图像。 要最小化的是transfer error。 这是第二个图像中的在测得的点
x
′
x'
x′和点
H
x
ˉ
H\bar{x}
Hxˉ之间的欧几里得图像距离,其中相应点
x
ˉ
\bar{x}
xˉ从第一个图像映射而来。那我们就可以来优化:
Σ
d
(
x
i
′
,
H
x
ˉ
)
2
\Sigma d(x'_i,H\bar{x})^2
Σd(xi′,Hxˉ)2
当两张图片当中都有误差的时候,我们把这个式子推广到第一幅图像,考虑前向 H H H和反向 H − 1 H^{-1} H−1,就可以得到:
Σ d ( x i , H − 1 x i ′ ) 2 + d ( x i ′ , H x i ) 2 \Sigma d(x_i,H^{-1}x'_i)^2 + d(x'_i,Hx_i)^2 Σd(xi,H−1xi′)2+d(xi′,Hxi)2
此总和中的第一个项是第一个图像中的误差,第二项是第二个图像中的误差。也就是第1幅图像的点变换到第二幅图像,它们需要尽可能靠近。反过来第2幅图像的点变换到第1幅图像,它们也需要尽可能靠近。
4.2.3 Reprojection error – both images
理论上完美的点,书中用 x ^ \hat{x} x^表示。首先我们要知道理论完美点不是 x ˉ \bar{x} xˉ, x ˉ \bar{x} xˉ就是机器学习里的Ground Truth,是不存在误差的。 x x x是我们已经知道的点,是有噪声的。那么我们用 x , x ′ x,x' x,x′估计出一个 H H H, H x Hx Hx肯定不等于 x ′ x' x′。但是我们可以找出一对 x ^ , x ^ ′ \hat{x},\hat{x}' x^,x^′,它们可以完美满足 x ^ = H x ^ ′ \hat{x}=H\hat{x}' x^=Hx^′。注意, H H H是从有噪声的数据里估计出来的,而不是从 x ˉ \bar{x} xˉ里估计出来的。
根据这个思想,我们可以去寻找一个 H ^ \hat{H} H^和一对完美的匹配点 x ^ , x ^ ′ \hat{x},\hat{x}' x^,x^′,它们之间满足 x ^ = H x ^ ′ \hat{x}=H\hat{x}' x^=Hx^′,然后让 x , x ′ x,x' x,x′靠近 x ^ , x ^ ′ \hat{x},\hat{x}' x^,x^′。这其实是先从 x , x ′ x,x' x,x′估计出一个三维空间点,然后再重新往图像上投影,所以叫重投影。
4.2.4 Comparison of geometric and algebraic distance
它们之间差一个常数因子,该常数因子由 x ′ x' x′的第三维坐标和 x ^ ′ \hat{x}' x^′相乘得到。
4.2.5 Geometric interpretation of reprojection error
重投影误差的几何含义到底是什么呢?一句话概括:重投影误差相当于在4维空间中拟合一个平面。
通俗解释,我们可以把
x
,
x
′
x,x'
x,x′的坐标拼在一起,得到
X
=
(
x
i
,
y
i
,
x
i
′
,
y
i
′
)
X=(x_i,y_i,x'_i,y'_i)
X=(xi,yi,xi′,yi′),然后把
x
^
,
x
^
′
\hat{x},\hat{x}'
x^,x^′的坐标拼在一起,得到
X
^
=
(
x
i
^
,
y
i
^
,
x
i
^
′
,
y
i
^
′
)
\hat{X}=(\hat{x_i},\hat{y_i},\hat{x_i}',\hat{y_i}')
X^=(xi^,yi^,xi^′,yi^′)。
那么 ∣ ∣ X − X ^ ∣ ∣ 2 ||X-\hat{X}||^2 ∣∣X−X^∣∣2就是重投影误差。
以上思想还可以被用来拟合圆锥,假设说我们要拟合圆锥 C C C,找5个点直接解圆锥方程很很麻烦,因为不是线性的。那么就求点到 C C C的距离,让它们加起来最短就可以了。
4.2.6 Sampson error
主要思想是用泰勒公式把要拟合的几何体展开到1阶,让它余项等于0。
4.2.7 Another geometric interpretation
上文说了求解 H H H相当于在4维空间中拟合一个平面。那么本节给出在n维度量空间中的几何解释。求解H相当于在n维度量空间中找出一个和 X X X最接近的向量。可以想象4维空间的平面在n维就成了点。
4.3 Statistical cost functions and Maximum Likelihood estimation
我们知道图像中都会包含噪声,如果我们认为噪声服从高斯分布,那么我们就可以用高斯分布然后求解最大后验概率。
比如,我们可以构造以下式子:
P
(
x
)
=
(
1
2
π
σ
2
)
e
x
p
(
d
(
x
,
x
ˉ
)
2
)
/
2
σ
2
P(x)=(\frac{1}{2 \pi \sigma^2}) exp(d(x,\bar{x})^2)/2 \sigma^2
P(x)=(2πσ21)exp(d(x,xˉ)2)/2σ2
其中 d d d可以换成任何你喜欢的距离,然后求解最大后验概率就可以。
4.4 Transformation invariance and normalization
现在,我们开始讨论第4.1节的DLT算法的属性和性能,以及与算法最小化几何误差的对比。 第一个主题是算法对图像中坐标的不同选择的不变性。 显然,对于算法而言,通常不希望取决于图像中坐标系的原点,比例甚至方向等任意选择。
4.4.1 Invariance to image coordinate transformations
图像原点有时候在左下角,有时候在右上角,那么不同的原点对DLT算法有啥影响吗?答案是有的,对几何损失函数是没有影响。
那么我们有这样的结论,假设有一个相似变换 T T T,变换前的点 x x x计算出的 H H H,和变换后计算出的 H ′ H' H′之间,相差一个倍数,也就是矩阵的向量之间相差一个倍数。
4.4.3 Invariance of geometric error
几何损失函数在相似变换下保持不变。
4.4.4 Normalizing transformations
首先我们明确为什么要对点坐标归一化,因为有的点坐标很大,有的很小,那么用SVD求解出来的答案会有很大的误差。
那么怎么归一化?先把所有点的质心旋转到原点,然后让缩放所有点,使所有点到原点的平均距离等于 2 \sqrt{2} 2,也就是让所有点的平均值等于 ( 1 , 1 , 1 ) (1,1,1) (1,1,1)。
需要明确的是,归一化是必须的,不是可选的。
4.5 Iterative minimization methods
如何用迭代的方法优化几何损失函数?
书中114页给出了一个案例,我们来一起看一下。
已知多于4对点
x
↔
x
′
x \leftrightarrow x'
x↔x′,求解重投影误差。
Σ
d
(
x
i
,
x
i
^
)
2
+
d
(
x
i
′
,
x
i
^
′
)
2
x
i
^
′
=
H
^
x
i
^
\Sigma d(x_i,\hat{x_i})^2 + d(x'_i,\hat{x_i}')^2 \\ \hat{x_i}'=\hat{H} \hat{x_i}
Σd(xi,xi^)2+d(xi′,xi^′)2xi^′=H^xi^
- 先初始化 H ^ \hat{H} H^ 用DLT解出一个初始值或者用RANSAC找出一个初始值
- 使用这个初始值来计算重投影误差的sampson余项
或者用初始值以及 x x x找出 x ^ , x ^ ′ \hat{x},\hat{x}' x^,x^′,根据 H ^ 和 x ^ \hat{H}和\hat{x} H^和x^最优化重投影误差。
4.6 Experimental comparison of the algorithms
DLT的效果是最接近黄金标准算法的。
4.7 Robust estimation
4.7.1 RANSAC
我们计算 H H H需要匹配点,但其实很多点都是误匹配的,所以我们先要找出合适的匹配点,就用RANSAC找。
4.7.2 Robust Maximum Likelihood estimation
4.7.3 Other robust algorithms
4.8 Automatic computation of a homography
- 先选4对点,找出 H H H
- 用 H H H估计出那些点是内点,也就是与 H x Hx Hx差距小于阈值的点
- 内点数量如果够多就转下步,否则再重复1到3
- 用所有的内点重新估计一次 H H H,怎么估计呢?用最大后验概率,因为现在的点比4个多,直接解DLT肯定不合适
- 知道了 H H H那么对应点就确定下来了
剩下的就是一些小问题。比如对应点怎么找?ORB、SIFT、FAST等等。点与点之间距离怎么衡量?ORB等等找出来的都是描述子,是一个向量,直接算向量之间距离就可以。怎么采样?应该对整个图像均匀采样。