一、基礎矩陣 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