这篇博客是建立在知道单应矩阵和相机参数相关概念的基础上进行讲解的,不了解的朋友可以参考我之前的博客:
单应矩阵:https://blog.csdn.net/qq_40369926/article/details/88918489
相机参数:https://blog.csdn.net/qq_40369926/article/details/89251296
在计算机视觉中,单张图像只能将真实的3D世界转化成2D图像,但是不能通过2D图像还原成真实的3D世界。但是,多张针对同一目标的图像就可以将2D图像还原成真实的3D模型,这也是我们所说的三维重建。基础矩阵描述了空间中的点在两个像平面中的坐标对应关系,这对于三维重建有非常重要,所以我们今天来学习一下基础矩阵。
一、外极几何
为了解释基础矩阵,我们要先了解一下外极几何的相关知识。
已知两个摄像头的光心和
,
为空间中的一点,
和
′是点
在两个摄像头成的像中的投影。
平面称为外极平面,显然
和
是
和
上的,所以该5点共面。外极平面
与两个相机的视平面相交于线
和
,这两条直线称为外极线。
在
上,
′在
上。
和
与相关视平面相交于点
和
,这两个点分别为
和
在对方视平面的投影,,称点
和
为外极点。如下图所示:
按照上面的定义,假设和
分别是空间中同一点在两个不同视平面上的像点,则
一定在
上,
一定在
上,这就是外极线约束。
如果我们已经知道相机的参数,那么我们需要解决的问题就是两张图像之间的关系,外极线约束就是对应特征点一定在外极线上,限制了对应点的搜索范围。
二、基础矩阵
有了外极几何的相关概念,我们就可以将讲基础矩阵了。
基础矩阵其实就是将上面所说的约束关系用代数表示出来。
假设左侧图像和右侧的图像之间的运动为和
。其中,
为旋转矩阵,
为相机中心位置的三维平移向量。
设在
,
坐标系中的相对坐标分别
,
,则有
设和
分别为两个相机的内参矩阵,根据约束关系可知:
即
即
又因为三线共面,所以
而叉积有可以写成矩阵相乘的形式,所以
定义矩阵为
将上式带入叉积公式中,得到:
即
定义本质矩阵,则:
根据前述,有:
带入得:
即
定义基础矩阵,则有:
基础矩阵描述了空间中的点在两个像平面中的坐标对应关系,是对积几何的代数表达方式,描述了图像中任意对应点 之间的约束关系。
基础矩阵有如下性质:
基础矩阵可以用于化简匹配和去除错配特征。
三、确定基础矩阵的方法
有相应的方程:
因为基础矩阵有7个自由度,所以确定基础矩阵最简单的方法就是8点算法,当有8对这样的点时,如下图所示:
则有如下方程:
令左边矩阵为,右边矩阵为
,即
优化方法一般使用最小二乘法,即优化
RANSAC算法可以用来消除错误匹配的的点,找到基础矩阵FF,算法思想如下:
(1)随机选择8个点;
(2)用这8个点估计基础矩阵;
(3)用这个进行验算,计算用
验算成功的点对数
;
重复多次,找到使最大的
作为基础矩阵。
8点算法的优缺点:
优点: 线性求解,容易实现,运行速度快
缺点:对噪声敏感
四、求解基础矩阵的实现
1.代码
# coding: utf-8
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
from PCV.geometry import homography, camera,sfm
from PCV.localdescriptors import sift
camera = reload(camera)
homography = reload(homography)
sfm = reload(sfm)
sift = reload(sift)
# 提取特征
im1 = array(Image.open('D:/test/test5/7.jpg'))
sift.process_image('D:/test/test5/7.jpg', 'im1.sift')
im2 = array(Image.open('D:/test/test5/8.jpg'))
sift.process_image('D:/test/test5/8.jpg', 'im2.sift')
l1, d1 = sift.read_features_from_file('im1.sift')
l2, d2 = sift.read_features_from_file('im2.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):
"""
使用RANSAC从点对应中稳健估计基本矩阵F.
(来自http://www.scipy.org/Cookbook/RANSAC的ransac.py)。
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
# 计算F并返回inlier索引
F, ransac_data = ransac.ransac(data.T, model,
8, maxiter, match_threshold, d, return_all=True)
return F, ransac_data['inliers']
# 通过RANSAC找到F.
model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-3)
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 is'
print F
# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)
# 绘制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]
# 读取特征
im3 = array(Image.open('D:/test/test5/9.jpg'))
sift.process_image('D:/test/test5/9.jpg', 'im3.sift')
l3, d3 = sift.read_features_from_file('im3.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 is'
print P1
print 'P2 is'
print P2
print 'P3 is'
print P3
2.实验结果
(1)外景
图像1和图像2的sift特征匹配结果图 图像1和图像2的投影效果图 图像1和图像3的投影效果图
基础矩阵为F,图像1、图像2和图像3的投影矩阵分别为P1、P2、P3(默认P1的投影矩阵为[[1 0 0 0] ,[0 1 0 0], [0 0 1 0]])
(2)落差大的外景
图像1和图像2的sift特征匹配结果图 图像1和图像2的投影效果图 图像1和图像3的投影效果图
基础矩阵为F,图像1、图像2和图像3的投影矩阵分别为P1、P2、P3(默认P1的投影矩阵为[[1 0 0 0] ,[0 1 0 0], [0 0 1 0]])
(3)内景
图像1和图像2的sift特征匹配结果图 图像1和图像2的投影效果图 图像1和图像3的投影效果图
基础矩阵为F,图像1、图像2和图像3的投影矩阵分别为P1、P2、P3(默认P1的投影矩阵为[[1 0 0 0] ,[0 1 0 0], [0 0 1 0]])
通过三组实验结果对比,我们可以知道,当图像之间落差较小时,sift算法能匹配出较多的特征点,基础矩阵计算出的投影效果较为理想。当图像之间落差较大时,sift算法能匹配出的特征点多集中在远处、变化不大的目标上,其范围较小,所以基础矩阵计算出来的投影效果不是很好,投影错误的点较多。当sift算法不能匹配出的特征点不够多时,基础矩阵计算出来的投影效果也不够好。