基于python3计算机视觉编程(六)多视图几何

1. 简介

多视图几何是描述同一场景不同视角的多幅图像与物体之间投影关系的几何模型。为了描述这种关系首先讨论单个相机只能获得一幅图像。实际意义是三维空间在二维空间上的一个投影。准确的说是变形的投影，因为相机投影关系和三维空间在二维基上最小二乘映射并不完全相同。

2. 对极几何

• 对极点（epipole）：基线（baseline）与成像平面的交点。同时极点也可以理解为相邻影像成像中心在本影像上的像，因为基线是两个相邻影像成像中心的连线。
• 对极平面（epipolar plane）：含有基线的平面，是一簇平面。可以看做是由基线与空间中任意一点构成的平面。
• 对极线（epipolar line）：对极平面与成像平面的交线。可以看做是成像平面上的任意一点（非核点）与核点所定义的直线。

3. 基础矩阵F

基础矩阵可以看做是将点投影（转换）为直线，将左影像上的一个点投影到右影像上形成一条核线。

3.1 基础矩阵性质

• 转置对称性：如果$F$是一对影像 $(P,P^{\prime })$ 的基础矩阵（即 $x^{\prime }Fx=0$ ），反过来 ($P^{\prime }$,$P$) 的基础矩阵是 $F^T$。证明很简单，直接对 $x^{\prime }Fx=0$ 两侧分别转置，得到 $xTFT{x}^{\prime }=0$。
• 对极线：对于左影像上任意一点 $x$ ，其在右影像上的对极线为 $l^{\prime }=Fx$ 。
• 对极点：任何对极线都会经过核点，所以有对于左影像上任意一点$x$ ，$e^{\prime }{}^T{l}^{\prime }=e\prime {}^T(Fx)=0$，于是有$e^{\prime }{}^{T}F=0$。同理有 $Fe=0$。
• $F$ 具有$7$自由度：一个$3$x$3$ 的单应矩阵，具有$8$个自由度，而 $F$ 还满足 $detF=0$，所以$F$ 具有$7$个自由度。
• $F$ 是相关的：$F$ 将左影像上的一点 $x$ 投影到右影像上一条对极线$l^{\prime }$，投影本质上是将$x$ 与左极点的连线$l$ 投影到右影像上的对极线 $l^{\prime }$ ，所以右影像上的一条对极线 $l^{\prime }$ 对应的是左影像上的一条对极线$l$，这种点到线的投影不可逆。

3.2 基础矩阵原理

基础矩阵表示的是图像中的像点$p_1$到另一幅图像对极线$l_2$的映射，有如下公式：
$$\begin{array}{r}{l}_2=F{p}_1\end{array}$$
而和像点$p_1$匹配的另一个像点$p_2$必定在对极线$l_2$上，所以就有：
$$\begin{array}{r}{p}_2^T{l}_2=p_2^{T}F{p}_1=0\end{array}$$
这样仅通过匹配的点对，就可以计算出两视图的基础矩阵$F$。
基础矩阵F是一个$3$×$3$的矩阵，有$9$个未知元素。然而，上面的公式中使用的齐次坐标，齐次坐标在相差一个常数因子下是相等，$F$也就只有$8$个未知元素，也就是说，只需要$8$对匹配的点对就可以求解出两视图的基础矩阵$F$。下面介绍下$8$点法 Eight-Point-Algorithm计算基础矩阵的过程。

3.3 8点法估计基础矩阵

八点算法步骤：

1. 求线性解由系数矩阵$A$最小奇异值对应的奇异矢量$f$求的$F$。
2. 奇异性约束 是最小化Frobenius范数$||$ $F$ $-$ $F^{\prime }$ $||$ 的 $F^{\prime }$代替$F$

假设一对匹配的像点$p_1=[u_1,v_1,1{]}^T,p_2=[u_2,v_2,1{]}^T$，带入式子中，得到：

把基础矩阵F的各个元素当作一个向量处理

$$\begin{array}{r}f=[f_1,f_2,f_3,f_4,f_5,f_6,f_7,f_8,f_9]\end{array}$$
于是上面得式子可以写成：
$$\begin{array}{r}[u_1{u}_2,u_1{v}_2,u_1,v_1{u}_2,v_1{v}_2,v_1,u_2,v_2,1]\cdot f=0\end{array}$$
对于其他的点对也使用同样的表示方法。这样将得到的所有方程放到一起，得到一个线性方程组$(u^i,v^i)$表示第$i$个特征点:

求解上面的方程组就可以得到基础矩阵各个元素了。$8$点算法估算基础矩阵虽然具有线性求救，容易实现，运行速度快的优点。但求解出来的元素很容易被噪声，错误的匹配点影响，因此需要再该算法上进行改进，这里我们采用图像坐标归一化 Normalizing transformation。
归一算法步骤：
给定$n\ge 8$组对应点{$X_i<->X_i^{\prime }$}，确定基本矩阵$F$使得$x_i^{\prime T}F{x}_i=0$。

1. 归一化：根据$X_i=T{X}_i$,$X_i=T^{\prime }{X}_i^{\prime }$,变换图像坐标。其中$T$和$T^{\prime }$是有平移和缩放组成的归一变化。
2. 求解对应匹配的基本矩阵$F^{\prime }$。
   1. 求线性解:用由对应点集{$X_i<->x_i^{\prime }$}确定的系数矩阵$A$的最小奇异值的奇异矢量确定$F$
   2. 奇异性约束:用SVD对$F$进行分解，令其最小奇异值为$0$,得到$F^{\prime }$，使得$det{F}^{\prime }=0$。
3. 解除归一化： 令$F=T^{\prime T}{F}^{\prime T}$。矩阵$F$就是数据$x_i<->x_i^{\prime }$对应的基本矩阵。

4. 实验

4.1 代码

# -*- coding: utf-8 -*-
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
from PCV.geometry import camera
from PCV.geometry import homography
from PCV.geometry import sfm
from PCV.localdescriptors import sift
# -*- coding: utf-8 -*-

# Read features
# 载入图像，并计算特征
im1 = array(Image.open('E:\\test\\2.jpg'))
sift.process_image('E:\\test\\2.jpg', 'E:\\test\\im1.sift')
l1, d1 = sift.read_features_from_file('E:\\test\\im1.sift')

im2 = array(Image.open('E:\\test\\3.jpg'))
sift.process_image('E:\\test\\3.jpg', 'E:\\test\\im2.sift')
l2, d2 = sift.read_features_from_file('E:\\test\\im2.sift')

# 匹配特征
matches = sift.match_twosided(d1, d2)
ndx = matches.nonzero()[0]

# 使用齐次坐标表示，并使用 inv(K) 归一化
x1 = homography.make_homog(l1[ndx, :2].T)
ndx2 = [int(matches[i]) for i in ndx]
x2 = homography.make_homog(l2[ndx2, :2].T)

x1n = x1.copy()
x2n = x2.copy()
print(len(ndx))

figure(figsize=(16,16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()

# Don't use K1, and K2

#def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    """ Robust estimation of a fundamental matrix F from point
    correspondences using RANSAC (ransac.py from
    http://www.scipy.org/Cookbook/RANSAC).

    input: x1, x2 (3*n arrays) points in hom. coordinates. """

    from PCV.tools import ransac
    data = np.vstack((x1, x2))
    d = 20 # 20 is the original
    # compute F and return with inlier index
    F, ransac_data = ransac.ransac(data.T, model,
                                   8, maxiter, match_threshold, d, return_all=True)
    return F, ransac_data['inliers']

# find E through RANSAC
# 使用 RANSAC 方法估计 E
model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-4)

print(len(x1n[0]))
print(len(inliers))

# 计算照相机矩阵（P2 是 4 个解的列表）
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)

# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)

# plot the projection of X
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)
x1p = cam1.project(X)
x2p = cam2.project(X)

figure()
imshow(im1)
gray()
plot(x1p[0], x1p[1], 'o')
#plot(x1[0], x1[1], 'r.')
axis('off')

figure()
imshow(im2)
gray()
plot(x2p[0], x2p[1], 'o')
#plot(x2[0], x2[1], 'r.')
axis('off')
show()

figure(figsize=(16, 16))
im3 = sift.appendimages(im1, im2)
im3 = vstack((im3, im3))

imshow(im3)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1p[0])):
    if (0<= x1p[0][i]<cols1) and (0<= x2p[0][i]<cols1) and (0<=x1p[1][i]<rows1) and (0<=x2p[1][i]<rows1):
        plot([x1p[0][i], x2p[0][i]+cols1],[x1p[1][i], x2p[1][i]],'c')
axis('off')
show()

print(F)

x1e = []
x2e = []
ers = []
for i,m in enumerate(matches):
    if m>0: #plot([locs1[i][0],locs2[m][0]+cols1],[locs1[i][1],locs2[m][1]],'c')
        x1=int(l1[i][0])
        y1=int(l1[i][1])
        x2=int(l2[int(m)][0])
        y2=int(l2[int(m)][1])
        # p1 = array([l1[i][0], l1[i][1], 1])
        # p2 = array([l2[m][0], l2[m][1], 1])
        p1 = array([x1, y1, 1])
        p2 = array([x2, y2, 1])
        # Use Sampson distance as error
        Fx1 = dot(F, p1)
        Fx2 = dot(F, p2)
        denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2
        e = (dot(p1.T, dot(F, p2)))**2 / denom
        x1e.append([p1[0], p1[1]])
        x2e.append([p2[0], p2[1]])
        ers.append(e)
x1e = array(x1e)
x2e = array(x2e)
ers = array(ers)

indices = np.argsort(ers)
x1s = x1e[indices]
x2s = x2e[indices]
ers = ers[indices]
x1s = x1s[:20]
x2s = x2s[:20]

figure(figsize=(16, 16))
im3 = sift.appendimages(im1, im2)
im3 = vstack((im3, im3))

imshow(im3)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1s)):
    if (0<= x1s[i][0]<cols1) and (0<= x2s[i][0]<cols1) and (0<=x1s[i][1]<rows1) and (0<=x2s[i][1]<rows1):
        plot([x1s[i][0], x2s[i][0]+cols1],[x1s[i][1], x2s[i][1]],'c')
axis('off')
show()

实验分为3组:角度变换图像、平移变换图像、远近变换图像
数据集：

第一组，平移变换(7点，8点，10点)
7点

8点

10点


第二组，角度变换(7点，8点，10点)
7点

8点

10点

第三组，远近变换(7点，8点，10点)
7点

8点

10点

4.2 实验总结

• 8点算法是计算基本矩阵的最简单的方法，运行速度相对于7点10点较快。
• 基础矩阵FFF描述的是3维场景的对极约束关系，和三维场景的结构无关，只依赖于相机的内参数以及外参数其描述的是点和线的约束关系，当已知图像上的点x1和基础矩阵F时是无法获得对应的像素点下x2的，只能获得x1对应的极线。
• 图像拍摄的环境会影响结果，匹配的点数太少会得到以下错误