立体成像的基本几何就是对极几何。下图是最经典的对极几何示意图。

如果仅看一个相机，我们并不能知道深度信息，可如果有两个相机的话（就像人有两只眼睛）我们就能得到深度的信息。
在上图中O和O’是两个相机中心（也有可能是一个相机在不同时刻的位置），P点是物体所在位置，两个相对的白色平面是像面(严格按照光路应该是在O1 O2点的后方，与P点相反方向，CV中默认采用这种往前画的方式节省空间)。p1和p2是P点在像面上的对应点，e1 e2为像面和O1 O2的交点。O1O2为基线，也被称作相机的移动方向。

在对极几何中，e1和e2被称作极点，PO1O2平面为极面，p1e1为极线，同理p2e2也为极线。这是对极几何中重要的三个概念。

1.2 极点与像平面关系

在下面的实验中，会分析极点与像平面在三种不同情况下的结果，所以这里先对这三种情况做个介绍。
先来看极点的性质——
极点的性质：

基线与像平面相交点
光心在另一幅图像中的投影

三种情况——
情况一：极点位于像平面

如上图，两个相机相对放置（汇聚视角），汇聚摄像机的对极几何是比较一般的情形，如b,c图所示，视图之间的运动是一个平移加上旋转。这是因为拍摄过程中，相机1面向右边，相机2面向左边，可知极点位于b的右边，c的左边。图中花瓶上标示的横线即为极线。

情况二：基线平行像平面

平行于图像的平移运动是运动的特殊情况，基线与图像平面相交于无穷远，因此对极点是无限远点，对极线是平行线。

情况三：相机前后方位关系

上图b,c为相机前后方位拍摄的图像，极点在各自像平面上的位置相同，可以看到b,c中的极线呈现放射状。

二、基础矩阵

有了对极几何的相关概念，我们就可以将讲基础矩阵了。在前言中也提到过，基础矩阵就是对极几何的代数表示。下面对它进行介绍。

对基础矩阵的作用理解：
如果已知基础矩阵F，以及一个3D点在一个像面上的像素坐标p，则可以求得在另一个像面上的像素坐标p’。这个是基础矩阵的作用，可以表征两个相机的相对位置及相机内参数。

2.1 本质矩阵E的推导

如上图，设X在C，C’坐标系中的相对坐标分别p，p’.
我们可以建立坐标轴，方便理解，如上图，以C为原点，光轴方向为z轴，另外两个方向为x, y轴可以得到一个坐标系，在这个坐标系下，可以对X, p(即图中所标x), p’(即图中所标x’)得到三维坐标，同理，对C’也可以得到一个三维坐标，这两个坐标之间的转换矩阵为[R T]，（说明：R为旋转矩阵,T为相机中心位置的三维平移向量）即通过旋转Ｒ和平移Ｔ可以将C坐标系下的点p(x1, y1, z1)，转换成C’坐标系下的p’(x2, y2, z2)。

则有： $p'=Rp+T$ p′=Rp+T如果我们将其左叉乘一个T，即 $T\times p'=T\times Rp+T\times T=T\times Rp$ T×p′=T×Rp+T×T=T×Rp其中 $T\times p'$ T×p′ 表示对极平面的法线，若再左点乘一个p’得到 $p'(T\times p')=p'(T\times Rp)$ p′(T×p′)=p′(T×Rp)由于p’与法线 $T\times p'$ T×p′是垂直的，所以有 $0=p'(T\times Rp)$ 0=p′(T×Rp)我们知道两向量的叉乘可以转换为一向量的反对称矩阵与另一向量的点乘，于是定义矩阵S为 $S=\begin{bmatrix} 0 &-T_{z} &T_{y} \\ T_{z}& 0 & -T_{x}\\ -T_{y}& T_{x} & 0 \end{bmatrix}$ S=⎣⎡0Tz−Ty−Tz0TxTy−Tx0⎦⎤S表示T的反对称矩阵。叉积秩为2。
将上式带入叉积公式中，得到： $p'(T\times Rp）=p'SRp=0$ p′(T×Rp）=p′SRp=0我们让 E=SR ,那么 $p'^{T}Ep=0_{——（1）式}$ p′TEp=0——（1）式这个E就是本质矩阵.
本质矩阵的作用：本质矩阵描述了空间中的点在两个坐标系中的坐标对应关系。

对本质矩阵E的求解过程作个小结：以上内容中，p, p’分别为P点的像点在两个坐标系下分别得到的坐标（非二维像素坐标）。Rp为极面上一矢量，T为极面上一矢量，则两矢量一叉乘为极面的法向量，这个法向量与极面上一矢量p’一定是垂直的，所以上式一定成立。（这里采用转置是因为p’会表示为列向量的形式，此处需要为行向量）

2.2 基础矩阵F的推导

得到本质矩阵后，基础矩阵也就容易推导了
设K和K’ 分别为两个相机的内参矩阵，根据约束关系可知： $x=Kp即p=K^{-1}x$ x=Kp即p=K−1x $x'=K'p'即p'=K'^{-1}x'$ x′=K′p′即p′=K′−1x′代入（1）式得到： $(K'^{-1}x')^{T}E(K^{-1}x)=0$ (K′−1x′)TE(K−1x)=0
这样我们就得到了两个相机上的像素坐标和基础矩阵F之间的关系了。
基础矩阵的作用：基础矩阵F描述了空间中的点在两个像平面中的坐标对应关系。

2.3 基础矩阵F的性质

基础矩阵的性质：
简单说来，3x3的矩阵，理论上9个*度，但是需要符合以下两个约束：

a)如果F为基础矩阵，那么kF也为基础矩阵
b)秩为2

所以减去两个*度，F有7个*度。

对第二条约束进行讨论（为什么秩为2？）：

我们知道，矩阵的秩有这么一个性质，矩阵相乘的秩不大于各矩阵的秩，那么，可以知道F其实是以下这些矩阵相乘的结果 $K'^{-T}SRK^{-1}$ K′−TSRK−1其中，S的具体形式在2.1中有表述 $S=\begin{bmatrix} 0 &-T_{z} &T_{y} \\ T_{z}& 0 & -T_{x}\\ -T_{y}& T_{x} & 0 \end{bmatrix}$ S=⎣⎡0Tz−Ty−Tz0TxTy−Tx0⎦⎤可以知道，第三列A2可以用第一列A0和第二列A1线性表示 $A2=-T_{y}*A1/T_{z}-T_{x}*A0/T_{z}$ A2=−Ty∗A1/Tz−Tx∗A0/Tz所以S秩为2，那么F秩为2也得证了。

2.4 基础矩阵F的求解

一般采用七点算法或八点算法对F进行求解，得知F后就可以对任意像面1上点找像面2上对应点了。下面介绍8点算法。

2.4.1 八点算法

由于基础矩阵F定义为： $x^{T}Fx'=0$ xTFx′=0其中 $x\leftrightarrow x'$ x↔x′是两幅图像的任意一对匹配点。由于每一组点的匹配提供了计算F系数的一个线性方程，当给定至少7个点（3×3的齐次矩阵减去一个尺度，以及一个秩为2的约束），方程就可以计算出未知的F。我们记点的坐标为 $x=(x,y,1),x'=(x',y',1)$ x=(x,y,1),x′=(x′,y′,1)，则对应的方程为

展开后有

把矩阵F写成列向量的形式，则有：

给定n组点的集合，我们有如下方程：

    如果存在确定（非零）解，则系数矩阵A的*度最多是8。由于F是齐次矩阵，所以如果矩阵A的*度为8，则在差一个尺度因子的情况下解是唯一的。可以直接用线性算法解得。
    如果由于点坐标存在噪声则矩阵A的*度可能大于8（也就是等于9，由于A是n×9的矩阵）。这时候就需要求最小二乘解，这里就可以用SVD来求解，f的解就是系数矩阵A最小奇异值对应的奇异向量，也就是A奇异值分解后 $A=UDV^{T}$ A=UDVT 中矩阵V的最后一列矢量，这是在解矢量f在约束||f||下取||Af||最小的解。以上算法是解基本矩阵的基本方法，称为8点算法。
    由于基本矩阵有一个重要的特点就是奇异性，F矩阵的秩是2。如果基本矩阵是非奇异的，那么所计算的对极线将不重合。所以在上述算法解得基本矩阵后，会增加一个奇异性约束。最简便的方法就是修正上述算法中求得的矩阵F。设最终的解为F′ ，令detF′=0下求得Frobenius范数（二范数）||F−F′|| 最小的F′。这种方法的实现还是使用了SVD分解，若 $F=UDV^{T}$ F=UDVT，此时的对角矩阵 $D=diag(r,s,t)$ D=diag(r,s,t)，满足 $r≥s≥t$ r≥s≥t，则 $F'=Udiag(r,s,0)V^{T}$ F′=Udiag(r,s,0)VT最小化范数||F−F′||，也就是最终的解。

所以8点算法由下面两个步骤组成：

1.求线性解由系数矩阵A最小奇异值对应的奇异矢量f 求的F。
2.奇异性约束是最小化Frobenius范数||F−F’||的F’代替F。

2.4.2 归一化八点算法

上面8点算法估算基础矩阵F的优缺点如下：

优点: 线性求解，容易实现，运行速度快
缺点：对噪声敏感

8点算法是计算基本矩阵的最简单的方法，但它存在问题，如：如果矩阵各列的数据尺度差异太大，最小二乘得到的结果精度一般很低。为了提高解的稳定性和精度，往往会对输入点集的坐标先进行归一化处理。对于归一化八点算法的总结如下：

给定 $n≥8$ n≥8 组对应点{ $x_{i}\leftrightarrow x_{i}'$ xi↔xi′}，确定基本矩阵F使得 $x_{i}'^{T}Fx_{i}=0$ xi′TFxi=0
算法：
1.归一化：根据 $\widehat{x_{i}}=Tx_{i},\widehat{x_{i}'}=T'x_{i}'$ xi=Txi,xi′=T′xi′变换图像坐标。其中T和T’是有平移和缩放组成的归一化变换。
2.求解对应匹配的基本矩阵 $\widehat{F'}$ F′
①求线性解：用由对应点集{ $\widehat{x_{i}}\leftrightarrow \widehat{x_{i}'}$ xi↔xi′}确定的系数矩阵 $\widehat{A}$ A的最小奇异值的奇异矢量确定 $\widehat{F}$ F
②奇异性约束：用SVD对 $\widehat{F}$ F进行分解，令其最小奇异值为0，得到 $\widehat{F'}$ F′，使得 $det\widehat{F'}=0$ detF′=0
③解除归一化：令 $F=T'^{T}\widehat{F}'T$ F=T′TF′T。矩阵F就是数据 $x_{i}\leftrightarrow x_{i}'$ xi↔xi′对应的基本矩阵。

归一化8点算法将图像坐标范围限定在 ~[-1,1]x[- 1,1]，实验表明可以得到更可靠的结果。

三、实验内容

3.1 实验目的与要求

分别用七点、八点、十点（匹配点），计算基础矩阵
实验图片包含三种情况，即（1）左右拍摄，极点位于图像平面上（2）像平面接*行，极点位于无穷远（3）图像拍摄位置位于前后
针对上述情况，画出极点和极线，其中点坐标要均匀分布于各行

3.2 代码实现

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt


def drawlines(img1, img2, lines, pts1, pts2):
    ''' img1 - image on which we draw the epilines for the points in img2
        lines - corresponding epilines '''
    r, c = img1.shape
    img1 = cv.cvtColor(img1, cv.COLOR_GRAY2BGR)
    img2 = cv.cvtColor(img2, cv.COLOR_GRAY2BGR)
    for r, pt1, pt2 in zip(lines, pts1, pts2):
        color = tuple(np.random.randint(0, 255, 3).tolist())
        x0, y0 = map(int, [0, -r[2] / r[1]])
        x1, y1 = map(int, [c, -(r[2] + r[0] * c) / r[1]])
        img1 = cv.line(img1, (x0, y0), (x1, y1), color, 1)
        img1 = cv.circle(img1, tuple(pt1), 5, color, -1)
        img2 = cv.circle(img2, tuple(pt2), 5, color, -1)
    return img1, img2


img1 = cv.imread('D:/Alike/python_dw/Code_a/admin_code/data/data_five/p58.jpg', 0)  # queryimage # left image
img2 = cv.imread('D:/Alike/python_dw/Code_a/admin_code/data/data_five/p59.jpg', 0)  # trainimage # right image
sift = cv.xfeatures2d.SIFT_create()
# find the keypoints and descriptors with SIFT
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)
# FLANN parameters
FLANN_INDEX_KDTREE = 1
index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks=50)
flann = cv.FlannBasedMatcher(index_params, search_params)
matches = flann.knnMatch(des1, des2, k=2)
good = []
pts1 = []
pts2 = []
# ratio test as per Lowe's paper
for i, (m, n) in enumerate(matches):
    if m.distance < 0.8 * n.distance:
        good.append(m)
        pts2.append(kp2[m.trainIdx].pt)
        pts1.append(kp1[m.queryIdx].pt)

pts1 = np.int32(pts1)
pts2 = np.int32(pts2)
F, mask = cv.findFundamentalMat(pts1, pts2, cv.FM_LMEDS,cv.FM_8POINT)
print F
# We select only inlier points
pts1 = pts1[mask.ravel() == 1]
pts2 = pts2[mask.ravel() == 1]

# Find epilines corresponding to points in right image (second image) and
# drawing its lines on left image
lines1 = cv.computeCorrespondEpilines(pts2.reshape(-1, 1, 2), 2, F)
lines1 = lines1.reshape(-1, 3)
img5, img6 = drawlines(img1, img2, lines1, pts1, pts2)
# Find epilines corresponding to points in left image (first image) and
# drawing its lines on right image
lines2 = cv.computeCorrespondEpilines(pts1.reshape(-1, 1, 2), 1, F)
lines2 = lines2.reshape(-1, 3)
img3, img4 = drawlines(img2, img1, lines2, pts2, pts1)
plt.subplot(121), plt.imshow(img5)
plt.subplot(122), plt.imshow(img3)
plt.show()

若要改变匹配点数量（7点、8点、10点），修改代码：

F, mask = cv.findFundamentalMat(pts1, pts2, cv.FM_LMEDS,cv.FM_8POINT)

3.3 实验结果

分三种情况讨论，每种情况输出极点极线图，并分别用7点、8点、10点（匹配点）计算基础矩阵。
（1）左右拍摄，极点位于图像平面上
极点极线图：

基础矩阵F:

（2）像平面接*行
极点极线图：

基础矩阵F:

（3）图像拍摄位置位于前后
极点极线图：

基础矩阵F:

3.4 实验分析

首先观察三种不同情况的结果图：
极点位于图像平面上的情况下，两幅图像中的极线各自有序排列，但延长线能够汇聚，图中的亮点表示的是对应点。
像平面接*行的情况下，极线呈现的是几乎平行的状态，这里对其展开解释:

如上图所示，当图像平面互相平行时，由 O1 、O2 构成的基线平行于图像平面，故极点 e 和 e’ 位于无穷远处，对极线是平行线。
回到上述实验，由于拍摄过程中相机是平移的关系，所以包包的像面互相平行，因此得到平行的对极线。
图像拍摄位置位于前后情况下，极线呈现放射状。在这组图片的拍摄过程中，相机前后移动，从结果图中可以看到所有极线的汇聚点为对面房间的桌子方位附近，极点在各自像平面上的位置相同。
观察图中的对应点，发现绝大多数对应点落在极线上，但也有部分对应点散落在图像上。
观察基础矩阵F，发现当选取的算法匹配点数不同时，得到的结果互不相同，一般情况下我们选取的是8点法求解基础矩阵，当然还可以选取7点、10点等，经查阅资料，要评估误差，可以采用3种方式：代数误差、几何误差和重投影误差。当然因为拍摄过程中位移与旋转角度不同，即R、T矩阵不同（或者外界光线等原因），看到不同情况下得到的是互不相同的F矩阵。

四、小结

实验内容小结：
通过本次实验了解了对极几何的概念以及基础矩阵的求导过程，对极几何解决的是已知两幅图像中两点是对应关系，如何求解两相机的相对位置和姿态的问题，而基础矩阵则是它的代数表示。通过三组不同情况的实验进一步加深对本章内容的理解，发现三组实验中，相机左右拍摄与前后拍摄得到的图像中，能看到极线明显地相交或者有相交趋势，而在像平面接近于平行的情况下，极线却是呈现平行的状态，这时极点在无穷远处。
实验过程小结：

本次实验在取材（拍摄图片）上我花了不少时间，要注意的是图片如果不压缩，可能需要很长的运行时间，但压缩了会导致找到的匹配点数量减少，影响sift特征匹配和Ransac算法的结果，所以要把握好图片像素大小。
遇到错误:pycharm中的array突然报错，检查开头已经导入numpy包，可是无法使用。
解决办法：根据系统的提示，在开头处添加from array import array即可。但是仍然存有疑惑，后续再做检查更正。

码农公寓

前言

一、对级几何

1.1 基本概念