计算计视觉多视图几何、基本矩阵

一、单视几何：

对极几何基本概念：
**核点：**基线（baseline）与成像平面的交点。同时极点也可以理解为相邻影像成像中心在本影像上的像，因为基线是两个相邻影像成像中心的连线。
**核平面：**含有基线的平面，是一簇平面。可以看做是由基线与空间中任意一点构成的平面。
**核线：**核平面与成像平面的交线。可以看做是成像平面上的任意一（非核点）与核点所定义的直线。
计算计视觉多视图几何、基本矩阵
单视几何的目标是从单幅图像中恢复场景的全部或部分三维信息，运用摄影几何理论，探索利用单幅图像实现场景测量所需的图像信息以及场景信息，从而实现对场景中距离、面积、体积等的测量。
计算计视觉多视图几何、基本矩阵
用两幅或多幅图像对场景进行重建以后进行测量的方法以及摄影测量学的方法有很大的局限性，目前，国内外在此方面还没有系统的研究。

二、两视几何：

**外极几何：**外极几何是研究两幅图像之间存在的几何，它和场景结构无关，只依赖y于摄像机的内外参数，研究这种几何可以用在图像匹配、三维重建方面。包含的基本概念如下：
计算计视觉多视图几何、基本矩阵
**基线：**连接两个摄像机光心O的直线。
**外极点：**基线与像平面的交心。
**外极平面：**过基线的平面。
**外极线：**对极平面与图像平面的交线。
**基本矩阵F：**对应点之间的约束。
本质上两幅图之间的对极几何是图像平面与以基线为轴的平面束的交的几何，这种几何被广泛同于双目成像原理中
如下图所示，摄像机由相机中心C,C’以及各自的成像平面表示，对于任意一个空间中的点X,在两个像平面上的点分别为x,x’，第一幅图像上的点x反向投影成空间三维的一条射线，它由摄像机中心和x确定，这条射线向第二个图像平面上投影得到一条直线l’，显然x的投影x’必然在l’上
计算计视觉多视图几何、基本矩阵
如下图所示，我们把相机中心的连线叫做基线,基线与两幅图像平面交于对极点e和e’，任何一条包含基线的平面pi是一张对极平面，分别与图像平面相交于l和l’，实际上，当三维位置X变化时，对应的实际上是该对极平面绕基线”旋转”，这个旋转得到的平面簇叫做对极平面束，由所有对极平面和成像平面相交得到的对极限相交于对极点。
计算计视觉多视图几何、基本矩阵

三、基本矩阵F：

基本矩阵是从本质矩阵推导出来的，描述的是两个相机视点之间的三维投影变换关系，是对极约束的一种代数表示。常用F表示的一个3x3矩阵，该矩阵的秩为0。经典的8点算法就可以比较准确地计算出该关系。
单应矩阵一般描述的是两个平面之间的二维投影变换关系。常用H表示，一般也是3x3的矩阵。最少情况下用平面上的4个点对应就以求出该矩阵，上述基本矩阵和单应矩阵的求解对应点越多越好，以避免出现退化情况。
基本矩阵可以看做是将点投影（转换）为直线，将左影像上的一个点投影到右影像上形成一条核线。基本矩阵是表示两个三维空间对应点集之间的变换关系，单应矩阵一般是表示两个二维空间对应点集之间的变换关系，给出由一个平面引入的单应矩阵H可以计算出基本矩阵F=[e’]xH。另外，一般情况下，三维空间中也有单应矩阵，它也用H表示，却是个4x4的矩阵。
几何推导基本矩阵：
计算计视觉多视图几何、基本矩阵
*图像点x是空间点X在图像1中的投影点，那么，有
x=H1X⇒(1)
*同理，图像点x’是空间点X在图像1中的投影点，那么，有
x′=H2X⇒(2)
*根据式（1）可以得到
X=H1^-1x⇒(3)

将式（3）代入式（2）得到

*图像点x对应的基线l’是通过点x’和点e’的直线，那么，有

其中的××为叉乘；
*将式（4）代入式（5）可以得到

基本矩阵的归一化8点算法：
基本矩阵是由下述方程定义：

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

如果存在确定（非零）解，则系数矩阵AA的秩最多是8。由于FF是齐次矩阵，所以如果矩阵A的秩为8，则在差一个尺度因子的情况下解是唯一的。可以直接用线性算法解得。
　　如果由于点坐标存在噪声则矩阵A的秩可能大于8（也就是等于9，由于A是n×9的矩阵）。这时候就需要求最小二乘解，这里就可以用SVD来求解，ff的解就是系数矩阵A最小奇异值对应的奇异向量，也就是AA奇异值分解后A=UDVTA=UDVT中矩阵V的最后一列矢量，这是在解矢量ff在约束∥f∥下取‖Af‖最小的解。以上算法是解基本矩阵的基本方法，称为8点算法。
　　由于基本矩阵有一个重要的特点就是奇异性，F矩阵的秩是2。如果基本矩阵是非奇异的，那么所计算的对极线将不重合。所以在上述算法解得基本矩阵后，会增加一个奇异性约束。最简便的方法就是修正上述算法中求得的矩阵FF。设最终的解为F′，令detF′=0detF′=0下求得Frobenius范数（二范数）‖F−F′‖最小的F′。这种方法的实现还是使用了SVD分解，若F=UDVTF，此时的对角矩阵D=diag(r,s,t)，满足r≥s≥t，F′=Udiag(r,s,0)VT最小化范数‖F−F′‖，也就是最终的解。
　　求解基本矩阵时图像之间的特征点匹配运用sift特征点匹配方法，匹配效果为：
　　
　　上面两张图片分别为室内室外的两张图片进行特征点匹配，由于照片特征点明显，而且不同的地方极少，因此匹配出来的照片特征点很多，特此建议照片选取时选择容易找特征点的照片而且相差不要太大。
ransac算法：
它可以从一组包含“局外点”的观测数据集中，通过迭代方式估计数学模型的参数。它是一种不确定的算法——它有一定的概率得出一个合理的结果；为了提高概率必须提高迭代次数。
RANSAC的基本假设是：
（1）数据由“局内点”组成；例如：数据的分布可以用一些模型参数来解释；
（2）“局外点”是不能适应该模型的数据；
（3）除此之外的数据属于噪声。
如下所示：

输入外极点为：
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
输出基本矩阵为：
[[-1.44377513e-06 -9.20916708e-05 2.22578825e-02]
[ 9.03617043e-05 -3.01894182e-07 -3.48946109e-02]
[-2.10177283e-02 3.02720001e-02 1.00000000e+00]]
输出相机矩阵为：
[[ 5.29380016e+00 -8.29916666e+00 2.37804934e+02 3.27936270e+02]
[-7.29916840e+00 1.14432989e+01 -3.27957287e+02 2.37835206e+02]
[-2.98568184e-02 -2.15901966e-02 1.49260425e+01 1.00000000e+00]]
室内照片，去除错配后输出匹配结果为：

输出基本矩阵为：
[[ 1.62421639e-06 8.18366352e-05 -3.16876196e-02]
[-9.34339436e-05 -3.56829675e-06 1.64503881e-01]
[ 3.20181963e-02 -1.61522876e-01 1.00000000e+00]]
相机矩阵为：
[[-1.19551567e+01 6.20638710e+01 3.77440568e+02 1.99017371e+03]
[ 6.30638690e+01 -3.27391393e+02 -1.99014169e+03 3.77279046e+02]
[ 1.62256337e-01 2.81491387e-02 -3.33538376e+02 1.00000000e+00]]
8点算法的基本步骤：
1、求线性解由系数矩阵A最小奇异值对应的奇异矢量f求的F。
2、奇异性约束是最小化Frobenius范||F−F′||的F′代替F。
为了提高解的稳定性和精度，往往会对输入点集的坐标先进行归一化处理,使得各个点做平移缩放之后到坐标原点的均方根距离等于√2。对于归一化八点算法的总结如下：
给定n≥8组对应点{xi↔xi′}，确定基本矩阵F使得xi′TFxi=0
算法：

归一化：根据x′i^{=T′x′ixi′}=T′xi′，变换图像坐标。其中T和T′是有平移和缩放组成的归一化变换。、
求解对应匹配的基本矩阵F^′F′
求线性解：用由对应点集{xi^↔xi′}确定的系数矩阵A^A的最小奇异值的奇异矢量确定F^F。
奇异性约束：用SVD对F^F进行分解，令其最小奇异值为0，得到F^′F′，使得detF^′=0detF′=0。
解除归一化：令F=T′TF^′TF=T′TF′T。矩阵F就是数据xi↔xi′对应的基本矩阵。

8点算法结果为：
计算计视觉多视图几何、基本矩阵
投影的基础矩阵为：
[[-7.56314977e-04 -3.37888902e-03 8.96512424e-01 1.32028992e-02]
[-6.42999160e-04 -1.32247386e-03 4.42705734e-01 8.86308301e-03]
[-3.38317487e-06 -1.03347672e-05 3.14943593e-03 3.49340981e-05]]
由于两幅图像在匹配的时候错误的匹配较少，所以计算的基本矩阵误差较小。
室内照片匹配结果为：
计算计视觉多视图几何、基本矩阵
投影的基础矩阵为：
[[ 5.70397637e-02 -4.65621354e-02 -1.22135945e-01 4.98519366e-02]
[ 1.76046606e-01 5.70680848e-03 2.27340319e-01 9.45193098e-01]
[ 4.60312787e-04 -1.07180588e-04 3.27134481e-02 2.41179585e-03]]

四、三维建模：

三维模型是物体的多边形表示，通常用计算机或者其它视频设备进行显示。显示的物体可以是现实世界的实体，也可以是虚构的物体，任何物理自然界存在的东西都可以用三维模型表示。
实现代码为：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import camera
import homography
import sfm
import sift
# 标定矩阵
K = array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像，并计算特征


im1 = array(Image.open('D:/qita/cai_view/picture/7.jpg'))
sift.process_image('D:/qita/cai_view/picture/7.jpg','im1.sift')
l1,d1 = sift.read_features_from_file('im1.sift')

im2 = array(Image.open('D:/qita/cai_view/picture/8.jpg'))
sift.process_image('D:/qita/cai_view/picture/8.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 = dot(inv(K),x1)
x2n = dot(inv(K),x2)

# 使用 RANSAC 方法估计 E
model = sfm.RansacModel()
E,inliers = sfm.F_from_ransac(x1n,x2n,model)
# 计算照相机矩阵（P2 是 4 个解的列表）
P1 = array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])
P2 = sfm.compute_P_from_essential(E)

# 选取点在照相机前的解
ind = 0
maxres = 0
for i in range(4):
    # 三角剖分正确点，并计算每个照相机的深度
    X = sfm.triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[i])
    d1 = dot(P1,X)[2]
    d2 = dot(P2[i],X)[2]

    if sum(d1>0)+sum(d2>0) > maxres:
        maxres = sum(d1>0)+sum(d2>0)
        ind = i
        infront = (d1>0) & (d2>0)
    # 三角剖分正确点，并移除不在所有照相机前面的点
    X = sfm.triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[ind])
    X = X[:,infront]

    # 绘制三维图像
    from mpl_toolkits.mplot3d import axes3d
    fig = figure()
    ax = fig.gca(projection='3d')
    ax.plot(-X[0], X[1], X[2], 'k.')
    axis('off')
    # 绘制 X 的投影 import camera
    # 绘制三维点
    cam1 = camera.Camera(P1)
    cam2 = camera.Camera(P2[ind])
    x1p = cam1.project(X)
    x2p = cam2.project(X)

    # 反 K 归一化
    x1p = dot(K, x1p)
    x2p = dot(K, x2p)
    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()

实现结果为：
计算计视觉多视图几何、基本矩阵
上图显示为ransac算法算出的模型参数，下图二里面、蓝色部分为局内点，而红色部分就是局外点，而这个算法要算出的就是蓝色部分那个模型的参数。实验呈现的结果为呈现一次函数的模型参数。
计算计视觉多视图几何、基本矩阵
RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点。

上面两幅图显示的结果蓝色为投影特征点，红色为原始特征点，由于照片特征点比较明显而且照片差异不大，所以投影特征点较多。
多次迭代标定特征点，此处省略其他特征点标定的结果。
最后附上运行代码：

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import camera
import homography
import sfm
import sift
#camera = reload(camera)
#homography = reload(homography)
#sfm = reload(sfm)
#sift = reload(sift)
# Read features
im1 = array(Image.open('D:/qita/cai_view/picture/13.jpg'))
sift.process_image('D:/qita/cai_view/picture/13.jpg', 'im1.sift')

im2 = array(Image.open('D:/qita/cai_view/picture/14.jpg'))
sift.process_image('D:/qita/cai_view/picture/14.jpg', 'im2.sift')
#processed tmp.pgm to im1.sift
#processed tmp.pgm to 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):
    """ 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')
print(F)
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)

print ('P2:')
print (P2)
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)

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:/qita/cai_view/picture/15.jpg'))
sift.process_image('D:/qita/cai_view/picture/15.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 (P3)
print('P1:')
print( P1)
print('P2:')
print (P2)
print('P3:')
print(P3)

计算计视觉多视图几何、基本矩阵

一、单视几何：

二、两视几何：

三、基本矩阵F：

四、三维建模：

相关推荐