多视图几何 - 爱码网

1、对极几何（Epipolar Geometry）

对极几何，用来描述的是如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。
图片来源及相关参考：https://blog.csdn.net/u012936940/article/details/80723609
多视图几何
极点e和e’：e是左边相机中心在右图像平面上的投影，e’右相机中心在左像平面上的投影
极平面：两个相机中心和点p形成的平面
极线l和l’：极平面分别和两个像平面的交线

多视图几何
若其中一相机中心e分别与像平面交于ee’,o、o’、p已知，即不知道p点的深度信息的情况下，射线op是p点可能的空间位置，o’p的连线同另一个像平面的交点即为p对应在该平面上的位置，可知P点对应在另一个像平面上的位置总是会在一条直线上，即在极线l’上

2、本质矩阵（Essential Matrix）、基础矩阵（fundamental matrix）

本质矩阵描述了空间中的点在两个坐标系中的坐标对应关系。
多视图几何
设X在C，C’坐标系中的相对坐标分别p，p’，则有 p = Rp’+T
x =Kp
p= K ^-1 x
x’= K’p
p’ =K’^-1 x’

基础矩阵描述了空间中的点在两个像平面中的坐标对应关系。
基础矩阵是一个7个自由度，秩为2 的 3*3矩阵。

多视图几何

8点算法估算基础矩阵F

基础矩阵 F 定义为
多视图几何
任给两幅图像中的匹配点 x 与 x’，

则有有相应方程：

在实际计算中，可以直接用A^TA的分解来求解参数。也可以用非线性优化，通过搜索f使得||Af||最小化,
同时满足||f||=1的约束。
上述求解后的F不一定能满足秩为2的约束，因此还要在F基础上加以约束。
多视图几何
矩阵各列的数据尺度差异太大，最小二乘得到的结果精度一般很低，所以应对矩阵各个列向量进行归一化操作！

•归一化8点算法
多视图几何

3、实验结果

实验代码：
代码来源：https://github.com/moizumi99/CVBookExercise/tree/master/Chapter-5
打开ipynb后缀文件方法参考：https://blog.csdn.net/zt384347827/article/details/56485253

from PIL import Image
from numpy import *
from pylab import *
import numpy as np


# In[2]:

from PCV.geometry import camera
from PCV.geometry import homography
from PCV.geometry import sfm
from PCV.localdescriptors import sift



# In[3]:

# Read features
im1 = array(Image.open('G:/PCV/0001.jpg'))
sift.process_image('G:/PCV/0001.jpg', 'im1.sift')
l1, d1 = sift.read_features_from_file('im1.sift')

im2 = array(Image.open('G:/PCV/0002.jpg'))
sift.process_image('G:/PCV/0002.jpg', 'im2.sift')
l2, d2 = sift.read_features_from_file('im2.sift')


# In[9]:

matches = sift.match_twosided(d1, d2)


# In[10]:

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)

x1n = x1.copy()
x2n = x2.copy()


# In[11]:

print (len(ndx))


# In[12]:

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


# In[13]:

# Chapter 5 Exercise 1
# 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-3):
    """ 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']


# In[15]:

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


# In[16]:

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


# In[17]:

P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)


# In[18]:

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

# In[19]:

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


# In[20]:

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()


# In[21]:

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()


# In[22]:

print (F)


# In[23]:

# Chapter 5 Exercise 2

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')
        p1 = array([int(l1[i][0]), int(l1[i][1]), 1])
        p2 = array([int(l2[int(m)][0]), int(l2[int(m)][1]), 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)


# In[24]:

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


# In[25]:

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()


# In[ ]:

①室外图像组（近景+远景）
sift特征匹配：多视图几何

多次实验发现，用室外图像组进行实验时，若景深差距相对较大，sift特征匹配效果不好，具有局限性
②室内图像组
sift特征匹配：