计算机视觉学习九:对极几何与基础矩阵
对极几何
如果仅看一个相机,我们并不能知道深度信息,可如果有两个相机的话(就像人有两只眼睛)我们就能得到深度的信息,
上图O和O’是两个相机中心,P点是物体所在,如果我们只看左边图像 \pi 上的点p,我们不能知道物体到底是在哪,点P1、P2或其他地方,可有了右边图像 \pi’ 上的p’我们就能得到物体点P。
在上图,我们把两相机中心的连线OO’成为基线,把他们与观测物体的平面OO’P成为对极平面,对极平面与两相机图像的交线l和l’称为对极线,而OO’与两图像的交点e,e’就是对极点。
随着观测点P的上下移动,对极平面也会围绕基线旋转
我们可以看到在左图对极平面旋转时对极点是不变的,而在相机图像上所有对极线都会交于对极点,这个对极点就是另一个相机中心在其图像上的像,当然正如右图所示,对极点可以在图像外。
对极点= 基线与像平面相交点= 光心在另一幅图像中的投影
对极平面 = 包含基线的平面对极线 = 对极平面与像平面的交线
参考博客:https://www.cnblogs.com/yuanlibin/p/9462180.html
基础矩阵
基础矩阵:空间中同一点在两个不同平面图像坐标系下投影的关系矩阵。
我们已经知道对于一幅视图上的点x,在另一视图上有一条对应的极线l’,该点的对应点x’必然在l’上,因此存在映射: x-> l’,基本矩阵F实际上就是表示这样的一种点到直线的射影映射
如图所示,考虑空间中不通过任何两个摄像机中心的平面pi,过第一个相机中心C和x的射线与该平面pi交于X,这个点再投影到第二幅图像上的点x’,这个过程通过平面pi的平移。由于X在对应与x的射线上,所以投影点x’必然在对应于这条射线的图像即极限l’上,点x和x’都是在一张平面上的3D点X的像,在第一幅图像上的所有点xi和对应点x’i的射影实际上是等价的,因为他们都射影等价于共面点集Xi,所以存在一个2D映射H,把每一个xi映射到x’i
给定点x’,通过x’和对极点e’的对极线I’,可以记为I’=e’ * x’ = [e’]x’,又因为x’=Hx,所以必然有
I’ = [e’]Hx = Fx
•基础矩阵是对极几何的代数表达方式
• 基础矩阵描述了图像中任意对应点 x↔x’ 之间的
约束关系
F 为 3x3 矩阵,秩为2,对任意匹配点对 x↔x’ 均满足
8点算法估算基础矩阵F
实验
代码
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('../mydata/test1.jpg'))
sift.process_image('../mydata/test1.jpg', 'im1.sift')
l1, d1 = sift.read_features_from_file('im1.sift')
im2 = array(Image.open('../mydata/test2.jpg'))
sift.process_image('../mydata/test2.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()
结果
注:本实验图像依然拍摄于集美大学
室外图像对(一张近拍,一张远拍)
基础矩阵F:
相机矩阵(投影矩阵):
sift特征匹配:
蓝色为投影特征点,红色为原始特征点:
室内图像对
基础矩阵F:
相机矩阵(投影矩阵):