您的位置: 首页  >  文章  >  Python---基础矩阵

Python---基础矩阵

分类: 文章 • 2024-11-13 12:08:46

 先思考一个问题：用两个相机在不同的位置拍摄同一物体，如果两张照片中的景物有重叠的部分，就可以说这两张照片之间存在一定的对应关系，本文的任务就是描述它们之间的对应关系，描述工具是对极几何 ，和基础矩阵。

对极几何

1、基本概念

 它是图像平面与以基线为轴的平面束的交的几何（这里的基线是指连接摄像机中心的直线）。以下图为例：对极几何描述的是左右两幅图像(点x和x’对应的图像)与以CC’为轴的平面束的交的几何
Python---基础矩阵Python---基础矩阵
直线CC’为基线，以该基线为轴存在一个平面束，该平面束与两幅图像平面相交。
下图为该平面束的直观形象，可以看到，该平面束中不同平面与两幅图像相交于不同直线
Python---基础矩阵
过e, e′的平面 π ,其平面上所有点在两个像平面中的 投影分别为直线 l 与 l’
上图中的平面π，只是过基线的平面束中的一个平面。∀π 的集与像平面的交线集，分别相交于 e 与 e’

2、对极几何的一个重要约束—5点共面约束

在下图中，空间点X在两幅图像中的像分别为x和x’，这两个投影点之间是否存在关系？
Python---基础矩阵
 点x、x’与摄像机中心C和C’是共面的，并且与空间点X也是空面的，这5个点共面于平面π。5个点决定了一个平面π是最本质的一个约束
 从上述约束，可以推导出一个重要性质：由图像点x和x’反投影的射线共面，并且，在平面π上

3、 对极几何的几个相关概念

Python---基础矩阵

对极平面束：以基线为轴的平面束；上图给出了包含两个平面的对极平面束
对极点 = 摄像机基线与像平面相交点 = 光心在另一幅图像中的投影，如上图中的e和e’
对极平面 = 任何包含基线的平面，或者说是对极平面束中的平面，上图中的平面π就是一个对极平面
对极线 = 对极平面与像平面的交线，上图中的I和I’

4、 对应点的约束

假设只知道图像点x，那么，它的对应点x’如何约束呢？
Python---基础矩阵
 由于点x和x’一定位于平面π上，而平面π可以利用基线CC’和图像点x的反投影射线确定
 点x’又是右侧图像平面上的点，所以，点x’一定位于平面π与右侧图像平面的交线l’上
 直线l’为点x的对极线，也就是说，点x的对应点x’一定位于它的对极线上。

基础矩阵

基础矩阵的作用：如果已知基础矩阵F，和一个3D点在一个像平面上的像素坐标p，就可以求得在另一个像面上的像素坐标p’。
Python---基础矩阵
以C为原点，光轴方向为z轴，另外两个方向为x, y轴可以得到一个坐标系，在这个坐标系下，可以对x，p，x’得到三维坐标。同理，对C’也可以得到一个三维坐标，这两个坐标之间的转换矩阵为[R T]，即通过旋转Ｒ和平移Ｔ可以将C坐标系下的点X(x1, y1, z1)， 转换成C’坐标系下的X’(x2, y2, z2)。
则有 p=Rp’+T
根据三线共面可以得出：
Python---基础矩阵
其中，p, p’分别为X点的像点在两个坐标系下分别得到的坐标。Rp’为极面上一矢量，T为极面上一矢量，则两个矢量叉乘为极面的法向量， 这个法向量与极面上一矢量X一定是垂直的
（RTp’）T（T x p）=0     T x p=Sp
Python---基础矩阵
Python---基础矩阵
 E为本质矩阵（essential matrix），本质矩阵描述了空间中的点在两 个坐标系中的坐标对应关系
 本质矩阵p’TEp=0， K 和 K’ 分别为两个相机的内参矩阵，有：
Python---基础矩阵Python---基础矩阵
基础矩阵是对极几何的代数表达方式 ,描述了空间中的点在两 个像平面中的坐标对应关系
基础矩阵描述了图像中任意对应点 x↔x’ 之间的 约束关系
Python---基础矩阵
 基础矩阵的性质：
 一般情况下，3x3的矩阵，理论上9个自由度，但是需要符合以下两个约束
 (1)转置: 如果 F 是表述点对 (x, x’)之间的基础矩阵, 则 FT 是 表述点对 (x’,x)之间的基础矩阵；
 (2)秩为2
 所以减去两个自由度，F有7个自由度, i.e. 3x3-1(homogeneous)-1(rank2)
 对任意匹配点对 x↔x’ 均满足 xTFx’=0

 (3)对极线: F 可以将点 x 映射到对应像平面上一条线 l=Fx ，同理可得 l’=FTx
 (4) 对极点: 对于所有对极线, 有 eTFx’=0, 同理∀x’ =>eT=0 有 Fe’=0

基本矩阵用途：
简化匹配
去除错配特征

极点性质

 前面说到，极点 = 摄像机基线与像平面相交点 = 光心在另一幅图像中的投影，如上图中的e和e’
（1）极点位于像平面
Python---基础矩阵
Python---基础矩阵
如图，两个相机相对放置， 相机1面向右边，相机2面向左边，可知极点位于1的右边，2的左边。图中花瓶上标示的横线即为平行于基线的线条。
为什么平行于基线的线条消失点为极点？上图中，平行于基线的线条所在极面与像平面必交于极点（基线必与像面交于极点），所以这些线条在像面上一定会交于极点。
Python---基础矩阵
（2）基线平行像平面
根据极点的定义，则极点位于无穷远处，极线与基线平行，如下图
Python---基础矩阵
Python---基础矩阵
这个时候，与基线平行的线条的在像面是一系列平行线，消失点在无穷远，和极点重合
（3）相机前后方位关系，当两个相机是前后放置且主点连线和像面垂直时
Python---基础矩阵Python---基础矩阵Python---基础矩阵
极点在各自像平面上的位置相同，且平行基线的线条在像面上的位置如上图右边所示。同样消失点为极点。

8点估算法

由于基础矩阵 F 定义为
Python---基础矩阵
任给两幅图像中的匹配点 x 与 x’ 。 令 x=(u,v,1)T ，x’=(u’,v’,1)T,
Python---基础矩阵
有相应方程：
Python---基础矩阵
Python---基础矩阵
在实际计算中，可以直接用ATA的分解来求解参数。 也可以用非线性优化，通过搜索f使得||Af||最小化, 同时满足||f||=1的约束。
上述求解后的F不一定能满足秩为2的约束，因此 还要在F基础上加以约束。
通过SVD分解可以解决上述问题，令 S=UƩVT则
Python---基础矩阵
则最终的解为 F’=UƩ’VT
优点: 线性求解，容易实现，运行速度快
缺点：对噪声敏感
由于矩阵各列的数据尺度差异太大， 最小二乘得到的结果精度一般很低，所以采用归一化8点算法

归一化8点算法

1.将图像坐标变换到合适的范围 ^Xi,=TXi ，^Xi’,=TXi’
2.根据变换后的坐标^Xi ，^Xi’ 计算归一化举出矩阵^F
3.还原原始基础矩阵 F=T’TF^T
Python---基础矩阵
归一化8点算法将图像坐标范围限定在 ~[-1,1]x[1,1]
实验表明可以得到更可靠的结果。
Python---基础矩阵

Python+opencv 实现

# coding: utf-8

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import sys
sys.path.append('C:/Users/WWW/Desktop/Chapter5/PCV')

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

# Read features
im1 = array(Image.open('C:/Users/WWW/Desktop/Chapter5/images/001.jpg'))
sift.process_image('C:/Users/WWW/Desktop/Chapter5/images/001.jpg', 'images_001.sift')

im2 = array(Image.open('C:/Users/WWW/Desktop/Chapter5/images/002.jpg'))
sift.process_image('C:/Users/WWW/Desktop/Chapter5/images/002.jpg', 'images_002.sift')

l1, d1 = sift.read_features_from_file('images_001.sift')
l2, d2 = sift.read_features_from_file('images_002.sift')


matches = sift.match_twosided(d1, d2)

# import cPickle as pkl
# fp=open("images_001_images_002_matches","w")
# pkl.dump(matches,fp)
# fp.close()


# import cPickle as pkl
# fp=open("images_001_images_002_matches","r")
# matches=pkl.load(fp)
# fp.close()


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):
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. """

    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']


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


#
# import cPickle as pkl
# fp=open("images_001_images_002_F_inliers","w")
# pkl.dump([F, inliers],fp)
# fp.close()

# import cPickle as pkl
# fp=open("images_001_images_002_F_inliers","r")
# F, inliers=pkl.load(fp)
# fp.close()

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)


# 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]]


# 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(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]

# Read features
im3 = array(Image.open('C:/Users/WWW/Desktop/Chapter5/images/003.jpg'))
sift.process_image('C:/Users/WWW/Desktop/Chapter5/images/003.jpg', 'images_003.sift')
l3, d3 = sift.read_features_from_file('images_003.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)
print (P2)
print (P3)

实验数据1：
Python---基础矩阵
结果：
Python---基础矩阵
Python---基础矩阵
Python---基础矩阵
实验数据2：
Python---基础矩阵
结果：
Python---基础矩阵
实验数据3：
Python---基础矩阵
结果：
Python---基础矩阵
从上面的实验数据可以看出，在室外拍摄并且光线较好的时候，从两幅图中找到对应点的效果较好，在室内进行实验时，由于光线的干扰，并且室内这组图相似度过高，匹配的结果一般。
 
 
 
 

相关推荐