一、基础矩阵 F
如果已知基础矩阵F，以及一个3D点在一个像面上的像素坐标p，则可以求得在另一个像面上的像素坐标p‘。这个是基础矩阵的作用，可以表征两个相机的相对位置及相机内参数。
下面具体介绍基础矩阵与像素坐标p和p’的关系。
以O1为原点，光轴方向为z轴，另外两个方向为x，y轴可以得到一个坐标系，在这个坐标系下，可以对P，p1（即图中所标p），p2（即图中所标p‘）得到三维坐标，同理，对O2也可以得到一个三维坐标，这两个坐标之间的转换矩阵为［RT］，即通过旋转Ｒ和平移Ｔ可以将O1坐标系下的点P1（x1，y1，z1），转换成O2坐标系下的P2（x2，y2，z2）。
则可知，P2=R（P1-T）（1）
采用简单的立体几何知识，可以知: PT .[Tx(R p')]=0 (2)
其中，p，p‘分别为P点的像点在两个坐标系下分别得到的坐标（非二维像素坐标）。Rp’为极面上一矢量，T为极面上一矢量，则两矢量一叉乘为极面的法向量，这个法向量与极面上一矢量p一定是垂直的，所以上式一定成立。（这里采用转置是因为p会表示为列向量的形式，此处需要为行向量）
采用一个常用的叉乘转矩阵的方法
(3)
将我们的叉乘采用上面的转换，会变成:PT*[Tx]* Rp'=0 (4)
红框中所标即为本征矩阵E，他描述了三维像点p和p‘之间的关系: E=[Tx]R (5)
有了本征矩阵，我们的基础矩阵也就容易推导了
注意到将p和p‘换成P1和P2式（4）也是成立的，且有
q1=K1P1（6）
q2=K2P2（7）
上式中，K1K2为相机的校准矩阵，描述相机的内参数q1q2为相机的像素坐标代入式（4）中，得
(8)
上式中p->q1，p‘->q2
这样我们就得到了两个相机上的像素坐标和基础矩阵F之间的关系了 : q1TFq2=0 (9)
基础矩阵的用途：
1、简化匹配
2、去除错误特征
二、基础矩阵8点算法：
（1）基本矩阵是由下述方程定义：x'TFx=0 ，其中x↔x′是两幅图像的任意一对匹配点。由于每一组点的匹配提供了计算F系数的一个线性方程，
当给定至少7个点（3×3的齐次矩阵减去一个尺度，以及一个秩为2的约束），方程就可以计算出未知的F。我们记点的坐标为x=(x,y,1)T x′=(x′,y′,1)T ，
则对应的方程为：
展开后有：
，
把矩阵F写成列向量的形式，则有：
给定n组点的集合，我们有如下方程：
(2)分析：
如果存在确定（非零）解，则系数矩阵A的自由度最多是8。由于F是齐次矩阵，所以如果矩阵A的自由度为8，则在差一个尺度因子的情况下解是唯一的。可以直接用线性算法解得。
如果由于点坐标存在噪声则矩阵A的自由度可能大于8（也就是等于9，由于A是n×9的矩阵）。这时候就需要求最小二乘解，这里就可以用SVD来求解，f的解就是系数矩阵A最小奇异值对应的奇异向量，也就是A奇异值分解后A=UDVT中矩阵V的最后一列矢量，这是在解矢量f在约束下取∥Af∥最小的解。以上算法是解基本矩阵的基本方法，称为8点算法。
由于基本矩阵有一个重要的特点就是奇异性，F矩阵的秩是2。如果基本矩阵是非奇异的，那么所计算的对极线将不重合。所以在上述算法解得基本矩阵后，会增加一个奇异性约束。最简便的方法就是修正上述算法中求得的矩阵F。设最终的解为F′，令detF′=0下求得Frobenius范数（二范数）∥F−F′∥最小的F。这种方法的实现还是使用了SVD分解，若F=UDVT，此时的对角矩阵D=diag(r,s,t)D=diag(r,s,t) D=diag(r,s,t)，满足r≥s≥t，则F′=Udiag(r,s,0)VT最小化范数∥F−F′∥，也就是最终的解。
所以8点算法由下面两个步骤组成：
如果由于点坐标存在噪声则矩阵A的自由度可能大于8（也就是等于9，由于A是n×9的矩阵）。这时候就需要求最小二乘解，这里就可以用SVD来求解，f的解就是系数矩阵A最小奇异值对应的奇异向量，也就是A奇异值分解后A=UDVT中矩阵V的最后一列矢量，这是在解矢量f在约束下取∥Af∥最小的解。以上算法是解基本矩阵的基本方法，称为8点算法。
由于基本矩阵有一个重要的特点就是奇异性，F矩阵的秩是2。如果基本矩阵是非奇异的，那么所计算的对极线将不重合。所以在上述算法解得基本矩阵后，会增加一个奇异性约束。最简便的方法就是修正上述算法中求得的矩阵F。设最终的解为F′，令detF′=0下求得Frobenius范数（二范数）∥F−F′∥最小的F。这种方法的实现还是使用了SVD分解，若F=UDVT，此时的对角矩阵D=diag(r,s,t)D=diag(r,s,t) D=diag(r,s,t)，满足r≥s≥t，则F′=Udiag(r,s,0)VT最小化范数∥F−F′∥，也就是最终的解。
所以8点算法由下面两个步骤组成：
1.求线性解由系数矩阵A最小奇异值对应的奇异矢量f求的F。
2.奇异性约束是最小化Frobenius范数∥F−F′∥的F′代替F。
2.奇异性约束是最小化Frobenius范数∥F−F′∥的F′代替F。
(3) 算法：
1.归一化：根据xiˆ=T，变换图像坐标。其中T和T′是有平移和缩放组成的归一化变换。、
2.求解对应匹配的基本矩阵Fˆ
2.求解对应匹配的基本矩阵Fˆ
2.1求线性解：用由对应点集{xiˆ↔x′iˆ}确定的系数矩阵Aˆ的最小奇异值的奇异矢量确定Fˆ。
2.2奇异性约束：用SVD对Fˆ进行分解，令其最小奇异值为0，得到Fˆ′，使得detFˆ′=0。
2.2奇异性约束：用SVD对Fˆ进行分解，令其最小奇异值为0，得到Fˆ′，使得detFˆ′=0。
3.解除归一化：令F=T'。矩阵F就是数据xi↔x′i对应的基本矩阵。
代码：
# coding: utf-8
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import sys
sys.path.append('C:/Users/WWW/Desktop/Chapter5/PCV')
import camera
import homography
import sfm
from PCV.localdescriptors import sift
# Read features
im1 = array(Image.open('C:/Users/WWW/Desktop/Chapter5/images/001.jpg'))
sift.process_image('C:/Users/WWW/Desktop/Chapter5/images/001.jpg', 'images_001.sift')
im2 = array(Image.open('C:/Users/WWW/Desktop/Chapter5/images/002.jpg'))
sift.process_image('C:/Users/WWW/Desktop/Chapter5/images/002.jpg', 'images_002.sift')
l1, d1 = sift.read_features_from_file('images_001.sift')
l2, d2 = sift.read_features_from_file('images_002.sift')
matches = sift.match_twosided(d1, d2)
ndx = matches.nonzero()[0]
x1 = homography.make_homog(l1[ndx, :2].T)
ndx2 = [int(matches[i]) for i in ndx]
x2 = homography.make_homog(l2[ndx2, :2].T)
d1n = d1[ndx]
d2n = d2[ndx2]
x1n = x1.copy()
x2n = x2.copy()
figure(figsize=(16,16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()
#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. """
import ransac
data = np.vstack((x1, x2))
d = 10 # 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 F through RANSAC
model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-5)
print (F)
#
# import cPickle as pkl
# fp=open("images_001_images_002_F_inliers","w")
# pkl.dump([F, inliers],fp)
# fp.close()
# import cPickle as pkl
# fp=open("images_001_images_002_F_inliers","r")
# F, inliers=pkl.load(fp)
# fp.close()
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)
print (P2)
print (F)
# P2, F (1e-4, d=20)
# [[ -1.48067422e+00 1.14802177e+01 5.62878044e+02 4.74418238e+03]
# [ 1.24802182e+01 -9.67640761e+01 -4.74418113e+03 5.62856097e+02]
# [ 2.16588305e-02 3.69220292e-03 -1.04831621e+02 1.00000000e+00]]
# [[ -1.14890281e-07 4.55171451e-06 -2.63063628e-03]
# [ -1.26569570e-06 6.28095242e-07 2.03963649e-02]
# [ 1.25746499e-03 -2.19476910e-02 1.00000000e+00]]
# 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(figsize=(16, 16))
imj = sift.appendimages(im1, im2)
imj = vstack((imj, imj))
imshow(imj)
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()
d1p = d1n[inliers]
d2p = d2n[inliers]
# Read features
im3 = array(Image.open('C:/Users/WWW/Desktop/Chapter5/images/003.jpg'))
sift.process_image('C:/Users/WWW/Desktop/Chapter5/images/003.jpg', 'images_003.sift')
l3, d3 = sift.read_features_from_file('images_003.sift')
matches13 = sift.match_twosided(d1p, d3)
ndx_13 = matches13.nonzero()[0]
x1_13 = homography.make_homog(x1p[:, ndx_13])
ndx2_13 = [int(matches13[i]) for i in ndx_13]
x3_13 = homography.make_homog(l3[ndx2_13, :2].T)
figure(figsize=(16, 16))
imj = sift.appendimages(im1, im3)
imj = vstack((imj, imj))
imshow(imj)
cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1_13[0])):
if (0<= x1_13[0][i]<cols1) and (0<= x3_13[0][i]<cols1) and (0<=x1_13[1][i]<rows1) and (0<=x3_13[1][i]<rows1):
plot([x1_13[0][i], x3_13[0][i]+cols1],[x1_13[1][i], x3_13[1][i]],'c')
axis('off')
show()
P3 = sfm.compute_P(x3_13, X[:, ndx_13])
print (P3)
print (P1)
print (P2)
print (P3)
三、运行结果：
数据集：
（1）左右拍摄，极点位于图像平面上
（2）像平面接近平行，极点位于无穷远
（3）图像拍摄位置位于前后
四、实验总结
（1）室外场景中，远近不同的图片组经过算法匹配效果相对于远近相同的图片组效果要差一些，甚至出现了明显的错误匹配，室内场景中，第一组的图片匹配点较多，但是优化后所剩下的较少，第二张匹配点相对第一组较少，但是优化过后明显特征的匹配点仍被保留。
（2）运行代码时遇到最多的问题，原因都出在图片本身上，因为图片的匹配点不够多，不符合匹配标准，但是即使是两张相似度很高的图片也有可能会出现这个错误，这时候可以去尝试调整阈值，或对换图片比对顺序，或者不要压缩图片。
（3）sift算法本身有时候会因为匹配点附近太过相似而出现错误匹配，这个算法便能在这个基础上相对优化，避开错误的匹配，但是以上的结果都出现了一个情况，优化后的结果仍存在错误的匹配点，整体效果并不好，由于图像远近以及光线的影响，呈现的效果也不够理想，而相对的，远近相同的测试图片效果会比远近相差较大的图片要好，更容易匹配。
（2）运行代码时遇到最多的问题，原因都出在图片本身上，因为图片的匹配点不够多，不符合匹配标准，但是即使是两张相似度很高的图片也有可能会出现这个错误，这时候可以去尝试调整阈值，或对换图片比对顺序，或者不要压缩图片。
（3）sift算法本身有时候会因为匹配点附近太过相似而出现错误匹配，这个算法便能在这个基础上相对优化，避开错误的匹配，但是以上的结果都出现了一个情况，优化后的结果仍存在错误的匹配点，整体效果并不好，由于图像远近以及光线的影响，呈现的效果也不够理想，而相对的，远近相同的测试图片效果会比远近相差较大的图片要好，更容易匹配。
（4）由于矩阵各列的数据尺度差异太大， 最小二乘得到的结果精度一般很低，所以采用归一化8点算法，从上面的实验数据可以看出，在室外拍摄并且光线较好的时候，从两幅图中找到对应点的效果较好，在室内进行实验时，由于光线的干扰，并且室内这组图相似度过高，匹配的结果一般。
五、外极点和外极线
（1）如果只是用一台摄像机我们不可能知道 3D 空间中的 X 点到图像平面的距离，因为 OX 连线上的每个点投影到图像平面上的点都是相同的。但是如果我 们也考虑上右侧图像的话，直线 OX 上的点将投影到右侧图像上的不同位置。 所以根据这两幅图像，我们就可以使用三角测量计算出 3D 空间中的点到摄像 机的距离（深度）。这就是整个思路。
直线 OX 上的不同点投射到右侧图像上形成的线 l′ 被称为与 x 点对应的极线。也就是说，我们可以在右侧图像中沿着这条极线找到 x 点。它可能在这条 直线上某个位置（这意味着对两幅图像间匹配特征的二维搜索就转变成了沿着极线的一维搜索。这不仅节省了大量的计算，还允许我们排除许多导致虚假匹配的点）。这被称为对极约束。与此相同，所有的点在其他图像中都有与之对应 的极线。平面 XOO' 被称为对极平面。
O 和 O' 是摄像机的中心。从上面的示意图可以看出，右侧摄像机的中心O' 投影到左侧图像平面的 e 点，这个点就被称为极点。极点就是摄像机中心连 线与图像平面的交点。因此点 e' 是左侧摄像机的极点。有些情况下，我们可能 不会在图像中找到极点，它们可能落在了图像之外（这说明这两个摄像机不能拍摄到彼此）。
所有的极线都要经过极点。所以为了找到极点的位置，我们可以先找到多条极线，这些极线的交点就是极点。 本节我们的重点就是找到极线和极点。为了找到它们，我们还需要两个元素，本征矩阵（E）和基础矩阵（F）。本征矩阵包含了物理空间中两个摄像机 相关的旋转和平移信息。如下图所示（本图来源自：Learning OpenCV）
基础矩阵 F 除了包含 E 的信息外还包含了两个摄像机的内参数。由于 F 包含了这些内参数，因此它可以它在像素坐标系将两台摄像机关联起来。（如果 使用是校正之后的图像并通过除以焦距进行了归一化，F=E）。简单来说，基 础矩阵 F 将一副图像中的点映射到另一幅图像中的线（极线）上。这是通过匹 配两幅图像上的点来实现的。要计算基础矩阵至少需要 8 个点（使用 8 点算 法）。点越多越好，可以使用 RANSAC 算法得到更加稳定的结果。
（2）代码
为了得到基础矩阵我们应该在两幅图像中找到尽量多的匹配点。我们可以使用 SIFT 描述符，FLANN 匹配器和比值检测。
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
img1 = cv.imread('myleft.jpg', 0) # queryimage # left image
img2 = cv.imread('myright.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)
现在得到了一个匹配点列表，我们就可以使用它来计算基础矩阵了。
下一步我们要找到极线。我们会得到一个包含很多线的数组。所以我们要 定义一个新的函数将这些线绘制到图像中。
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
现在计算并绘制极线。
# 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()
（3）运行结果：
六、实验分析与总结：
(1)选择两幅图像的对应点，然后将它们转换为齐次坐标，这里的对应点是从一个文本文件中读取得到的；而实际上，可以按照sift提取图像特征的方式，然后通过匹配来找到它们。由于缺失的数据在对应列表corr中为-1，所以程序中有可能选取这些点。因此，上面的程序通过数组操作符&只选取了索引大于等于0的点。借助plot_epipolar_line()函数，这个函数将x轴的范围作为直线的参数，所以直线超出了图像边界的部分会被截断。如果show_epipole为真，外极点会被画出来（如果输入参数中没有外极点，外极点会在程序中计算获得）。用不同的颜色将点和对应的外极线对应起来。
(1)选择两幅图像的对应点，然后将它们转换为齐次坐标，这里的对应点是从一个文本文件中读取得到的；而实际上，可以按照sift提取图像特征的方式，然后通过匹配来找到它们。由于缺失的数据在对应列表corr中为-1，所以程序中有可能选取这些点。因此，上面的程序通过数组操作符&只选取了索引大于等于0的点。借助plot_epipolar_line()函数，这个函数将x轴的范围作为直线的参数，所以直线超出了图像边界的部分会被截断。如果show_epipole为真，外极点会被画出来（如果输入参数中没有外极点，外极点会在程序中计算获得）。用不同的颜色将点和对应的外极线对应起来。
(2)RANSAC 的稳健估计对于远近不同的图像集可能会最终遇到无法找到适配点对的情况，需要在调用时多次尝试对阈值（match_threshold=1e-3）进行调整以适用当前图像集；
(3)对于较远的目标适配点对数量大，越容易增加错配无法去除的点对数量，对于较近的目标则相反；
(4)错配情况时常发生且没有被剔除，甚至保留了错配很明显的点对而将较好的适配点剔除，这种现象与几乎完全相同的物体结构外观相对应，其中还有因追求时间效率将高分辨率图像进行了压缩处理的影响。
(4)错配情况时常发生且没有被剔除，甚至保留了错配很明显的点对而将较好的适配点剔除，这种现象与几乎完全相同的物体结构外观相对应，其中还有因追求时间效率将高分辨率图像进行了压缩处理的影响。
————————————————
参考链接:
原文链接：https://blog.csdn.net/weixin_42424674/java/article/details/89361267
原文链接：https://blog.csdn.net/kokerf/java/article/details/72630863
原文链接：https://blog.csdn.net/mary_0830/article/details/91492422
原文链接：https://blog.csdn.net/yukinoai/article/details/89813463
来源:CSDN