求解基础矩阵

关注

原创最新推荐文章于 2026-07-12 14:24:26 发布 6k 阅读 AI 写同款· GEO 优化 ›

基础矩阵文章目录一、基础矩阵1、简介2、F的几何推导3、F的代数推导![在这里插入图片描述](https://img-blog.csdnimg.cn/2020042115321268.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2xwMjAxN... 阅读详情

文章目录

一、对极几何
- 1.1概念
- 1.2基本模型
二、基础矩阵
- 2.1基本矩阵原理
- 2.2八点估算法
三、实验内容
四、总结

一、对极几何

1.1概念

对极几何（Epipolar Geometry）是Structure from Motion问题中，在两个相机位置产生的两幅图像的之间存在的一种特殊几何关系，是sfm问题中2D-2D求解两帧间相机姿态的基本模型。

1.2基本模型

在这里插入图片描述
其中c0、c1为两个相机中心，p为空间中一点，p在c0、c1对应像平面上的投影分别为x0、x1。c0、c1连线与像平面的交点e0、e1称为极点（Epipoles），l0、l1称为极线（Epipolar Lines），c0、c1、p三点组成的平面称为极平面（Epipolar Plane）。

二、基础矩阵

通过对极几何一副图像上的点可以确定另外一幅图像上的一条直线，这种情况用基础矩阵来表示。通过一种映射，一幅图像上的点可以确定另外一副图像上的一个点，这种情况用单应矩阵。
本质矩阵是基础矩阵的一种特殊情况，是在归一化图像坐标下的基础矩阵。
本质矩阵
存在一个不经过两个相机光心的的平面π，光心C与x的射线与平面π相交与一点X。该点X又投影到第二幅图像平面上的点x′。这个称为点x通过平面π的转移。点x，x′是平面ππ上的3D点X在两个相机平面上的像。对应每一个3D点X都存在一个2D的单应R把每一个x映射到x′。
在这里插入图片描述
如上图所示，给定一个目标点P，以左摄像头光心Ol为原点。点P相对于光心Ol的观察位置为Pl，相对于光心Or的观察位置为Pr。点P在左摄像头成像平面上的位置为pl，在右摄像头成像平面上的位置为pr。
现在我们要寻找由点P、Ol和Or确定的对极平面的表达式。注意到平面上任意一点x与点a的连线垂直于平面法向量n，即向量 (x-a) 与向量 n 的点积为0：(x-a)·n = 0。在Ol坐标系中，光心Or的位置为T，则P、Ol和Or确定的对极平面可由下式表示：
在这里插入图片描述
由Pr = R(Pl-T) 得
另一方面，向量的叉积又可表示为矩阵与向量的乘积，记向量T的矩阵表示为S，得：

也就可以得到
就可以得到本质矩阵E=RS。
通过矩阵E我们知道Pl和Pr的关系满足：

根据相似三角形定理，pl = flPl/Zl 和 pr = frPr/Zr 我们可以得到点P在左右两个摄像机坐标系中的观察点 pl 和 pr 应满足的极线约束关系为：
在这里插入图片描述
注意到 E 是不满秩的，它的秩为2，那么

表示的是一条直线，也就是对极线。

2.1基本矩阵原理

为了描述对极几何，引入基础矩阵F。对于一幅图像上的点x（图上p1），在另一幅图像上存在对极线l’，并且在第二幅图像上，与x匹配的点x’（点p2）必然在l’上。为了表示x与l’关系，我们令基础矩阵F定义为： l′=Fx。

假设第一幅图到第二幅图之间存在单应性变换H，则 x′=Hx。

又因为l’是表示过对极点e2和图像点x’的直线，可以表示为： l′=e2×x′=[e2]xx′=[e2]xHx。

所以可得 F=[e2]xH。
其中e′=[a1,a2,a3]Te′=[a1,a2,a3]T，则它的反对称阵定义为：在这里插入图片描述
在已知两个摄像机的射影矩阵P和P’时，对于图像上的一点x，可以得到其反投影射线方程为： X(λ)=P+ x+λC

其中P+ 是P的伪逆，即PP+ = I，C为摄像机中心。射线上的两点：P+ x（当λ=0）、C（当λ=∞）在第二个摄像机P’拍摄下，在第二幅视图上的点分别为P′P+ 和P′C，这两点过对极线l’，即 l′=(P′C)×(P′P+ x)=[e′]x(P′P+ )x

可以推得 F=[e′]xP′P+

2.2八点估算法

基本矩阵是由该方程定义的： x′TFx=0

其中x↔x′是两幅图像的任意一对匹配点。由于每一组点的匹配提供了计算F系数的一个线性方程，当给定至少7个点（3×3的齐次矩阵减去一个尺度，以及一个秩为2的约束），方程就可以计算出未知的F。我们记点的坐标为x=(x,y,1)T，x′=(x′,y′,1)T

又因为F为：
在这里插入图片描述
所以可得到方程：

即相应方程式为

给定n组点的集合，我们有如下方程：

如果存在确定（非零）解，则系数矩阵A的秩最多是8。由于F是齐次矩阵，所以如果矩阵A的秩为8，则在差一个尺度因子的情况下解是唯一的。可以直接用线性算法解得。

如果由于点坐标存在噪声则矩阵AA的秩可能大于8（也就是等于9，由于A是n×9的矩阵）。这时候就需要求最小二乘解，这里就可以用SVD来求解，f的解就是系数矩阵A最小奇异值对应的奇异向量，也就是A奇异值分解后A=UDVT中矩阵V的最后一列矢量，这是在解矢量f在约束∥f∥下取∥Af∥最小的解。以上算法是解基本矩阵的基本方法，称为8点算法。

上述求解后的F不一定能满足秩为2的约束，因此还要在F的基础上加以约束。通过SVD分解可以解决，令F=UΣVT，则
在这里插入图片描述
因为要秩为2，所以取最后一个元素设置为0，则

最终的解

三、实验内容

3.1实验图片

在这里插入图片描述

3.2实验结果

前后拍摄场景

特征匹配
在这里插入图片描述
极点极线

基础矩阵

平行

特征匹配
在这里插入图片描述
极点极线

基础矩阵

图片左右拍摄

匹配
在这里插入图片描述
极点极线

在这里插入图片描述
基础矩阵

3.3实验代码

特征匹配

# 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(‘C:/Users/jxtx/计算机视觉/5/1.jpg’))
sift.process_image('C:/Users/jxtx/计算机视觉/5/1.jpg', 'im1.sift')

im2 = array(Image.open('C:/Users/jxtx/计算机视觉/5/2.jpg'))
sift.process_image('C:/Users/jxtx/计算机视觉/5/2.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

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]

绘制极点极线

# 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('C:/Users/jxtx/计算机视觉/5/1.jpg'))
sift.process_image('C:/Users/jxtx/计算机视觉/5/1.jpg', 'im1.sift')

im2 = array(Image.open('C:/Users/jxtx/计算机视觉/5/2.jpg'))
sift.process_image('C:/Users/jxtx/计算机视觉/5/2.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()

#import sfm111
### 计算 F
#F = sfm111.compute_fundamental(x1n,x2n)
#print(F)


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

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]

def compute_epipole(F):
 """ 从基础矩阵 F 中计算右极点（可以使用 F.T 获得左极点）"""
 # 返回 F 的零空间（Fx=0）
 U,S,V = linalg.svd(F)
 e = V[-1]
 return e/e[2]



def plot_epipolar_line(im,F,x,epipole=None,show_epipole=True):
 """ 在图像中，绘制外极点和外极线 F×x=0。F 是基础矩阵，x 是另一幅图像中的点 """
 m,n = im.shape[:2]
 line = dot(F,x)
 # 外极线参数和值
 t = linspace(0,n,100)
 lt = array([(line[2]+line[0]*tt)/(-line[1]) for tt in t])
 # 仅仅处理位于图像内部的点和线
 ndx = (lt>=0) & (lt<m)
 plot(t[ndx],lt[ndx],linewidth=2)
 if show_epipole:
     if epipole is None:
         epipole = compute_epipole(F)
     plot(epipole[0]/epipole[2],epipole[1]/epipole[2],'r*')

e = compute_epipole(F)
for i in range(5):
    plot_epipolar_line(im1,F,x2[:,i],e,False)
axis('off')
figure()
imshow(im2)
# 分别绘制每个点，这样会绘制出和线同样的颜色
for i in range(5):
    plot(x2[0,i],x2[1,i],'o')
axis('off')
show()

d1p = d1n[inliers]
d2p = d2n[inliers]