计算机视觉学习（八）--基础矩阵（多视图几何）
一、原理
对极几何
对极几何，一定是对二幅图像而言，对极几何实际上是“两幅图像之间的对极几何”，它是图像平面与以基线为轴的平面束的交的几何（这里的基线是指连接摄像机中心的直线），以下图为例：对极几何描述的是左右两幅图像(点x和x’对应的图像)与以CC’为轴的平面束的交的几何！
上图中，点x、x’与摄像机中心C和C’是共面的，并且与空间点X也是空面的，这5个点共面于平面π！这是一个最本质的约束，即5个点决定了一个平面π。
由此可以推导出：由图像点x和x’反投影的射线共面，并且，在平面π上，在搜索点对应中，该性质非常重要。
对极平面束：以基线为轴的平面束；下图为包含两个平面的对极平面束：
直线CC’为基线，以该基线为轴存在一个平面束，该平面束与两幅图像平面相交，上图图给出了该平面束的直观形象，可以看到，该平面束中不同平面与两幅图像相交于不同直线；之前的灰色平面π，只是过基线的平面束中的一个平面（当然，这个平面才是平面束中最重要的、也是我们要研究的平面）。
对极点 = 基线与像平面相交点 = 光心在另一幅图像中的投影
对极平面 = 包含基线的平面
对极线 = 对极平面与像平面的交线
本质矩阵
本质矩阵描述了空间中的点在两个世界坐标系中的坐标对应关系。
如果已知基础矩阵F，以及一个3D点在一个像面上的像素坐标p，则可以求得在另一个像面上的像素坐标p‘。设X在C，C’坐标系中的相对坐标分别p，p’，则有：
根据三线共面：
其中的RS相乘得到的E即为本质矩阵，描述了空间中的点在两个坐标系中的坐标对应关系。
根据前述， K 和 K’ 分别为两个相机的内参矩阵：
基础矩阵是对极几何的代数表达方式，描述了图像中任意对应点 x↔x’ 之间的约束关系。F 为 3x3 矩阵，秩为2，对任意匹配点对 x↔x’ 均满足上式
参考博文：多视图几何-https://blog.csdn.net/weixin_43843780/article/details/89362546
对极几何-https://blog.csdn.net/tina_ttl/article/details/52749542
实现代码
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
# Read features
# 载入图像，并计算特征
im1 = array(Image.open('s1.jpg'))
sift.process_image('s1.jpg', 'im1.sift')
l1, d1 = sift.read_features_from_file('im1.sift')
im2 = array(Image.open('s2.jpg'))
sift.process_image('s2.jpg', 'im2.sift')
l2, d2 = sift.read_features_from_file('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.geometry 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()
实验结果
室内
两张图片的特征点匹配
得到的基础矩阵：
室外
室外匹配试了无数组图片后依旧是下面的这个问题
下图是我认为特征点匹配最成功的的一张了，但依旧报错（(╥╯^╰╥)）
三维点和照相机矩阵
算法使用的是八点算法，代码如下：
# coding: utf-8
# In[1]:
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import importlib
# In[2]:
from PCV.geometry import camera
from PCV.geometry import homography
from PCV.geometry import sfm
from PCV.localdescriptors import sift
camera = importlib.reload(camera)
homography = importlib.reload(homography)
sfm = importlib.reload(sfm)
sift =importlib.reload(sift)
# In[3]:
# Read features
im1 = array(Image.open('i5.jpg'))
sift.process_image('i5.jpg', 'im1.sift')
im2 = array(Image.open('i6.jpg'))
sift.process_image('image6.jpg', 'im2.sift')
# In[4]:
l1, d1 = sift.read_features_from_file('im1.sift')
l2, d2 = sift.read_features_from_file('im2.sift')
# In[5]:
matches = sift.match_twosided(d1, d2)
# In[6]:
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()
# In[7]:
figure(figsize=(16,16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()
# In[26]:
#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 = 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']
# In[27]:
# find F through RANSAC
model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-5)
print (F)
# In[28]:
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)
# In[29]:
print (P2)
print (F)
# In[30]:
# 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]]
# In[31]:
# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)
# In[32]:
# plot the projection of X
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)
x1p = cam1.project(X)
x2p = cam2.project(X)
# In[33]:
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()
# In[34]:
d1p = d1n[inliers]
d2p = d2n[inliers]
# In[35]:
# Read features
im3 = array(Image.open('i7.jpg'))
sift.process_image('i7.jpg', 'im3.sift')
l3, d3 = sift.read_features_from_file('im3.sift')
# In[36]:
matches13 = sift.match_twosided(d1p, d3)
# In[37]:
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)
# In[38]:
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()
# In[39]:
P3 = sfm.compute_P(x3_13, X[:, ndx_13])
# In[40]:
print (P3)
# In[41]:
print (P1)
print (P2)
print (P3)
# In[22]:
# Can't tell the camera position because there's no calibration matrix (K)