计算机视觉python--基础矩阵和极点极线

1 外极几何

1.1 基本概念

在这里插入图片描述
在这里插入图片描述

1.2 基础矩阵原理和介绍

在计算机视觉中,基础矩阵(Fundamental matrix)F是一个3×3的矩阵,表达了立体像对的像点之间的对应关系。在对极几何中,对于立体像对中的一对同名点,它们的齐次化图像坐标分别为p与 p’, 表示一条必定经过p’的直线(极线)。这意味着立体像对的所有同名点对都满足: 在这里插入图片描述
F矩阵中蕴含了立体像对的两幅图像在拍摄时相互之间的空间几何关系(外参数)以及相机检校参数(内参数),包括旋转、位移、像主点坐标和焦距。因为F矩阵的秩为2,并且可以自由缩放(尺度化),所以只需7对同名点即可估算出F的值。
基础矩阵这一概念由Q. T. Luong在他那篇很有影响力的博士毕业论文中提出。 [1] Faugeras则是在1992年发表的著作 [2] 中以上面的关系式给出了F矩阵的定义。尽管Longuet-Higgins提出的本质矩阵也满足类似的关系式,但本质矩阵中并不蕴含相机检校参数。本质矩阵与基础矩阵之间的关系可由下式表达:在这里插入图片描述
其中K和K’分别为两个相机的内参数矩阵。
在这里插入图片描述

2 基本矩阵之归一化8点算法

基本矩阵是由下述方程定义:
在这里插入图片描述
其中x↔x’是两幅图像的任意一对匹配点。由于每一组点的匹配提供了计算F系数的一个线性方程,当给定至少7个点(3×3的齐次矩阵减去一个尺度,以及一个秩为2 的约束),方程就可以计算出未知的F FF。我们记点的坐标为x=(x,y,1)T,则对应的方程为在这里插入图片描述
展开后有
在这里插入图片描述
把矩阵F FF写成列向量的形式,则有:在这里插入图片描述
给定n nn组点的集合,我们有如下方程:在这里插入图片描述
如果存在确定(非零)解,则系数矩阵A的秩最多是8。由于F是齐次矩阵,所以如果矩阵A的秩为8,则在差一个尺度因子的情况下解是唯一的。可以直接用线性算法解得。

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

由于基本矩阵有一个重要的特点就是奇异性,F矩阵的秩是2。如果基本矩阵是非奇异的,那么所计算的对极线将不重合。所以在上述算法解得基本矩阵后,会增加一个奇异性约束。最简便的方法就是修正上述算法中求得的矩阵F。设最终的解为F′ ,令detF′=0 下求得Frobenius范数(二范数)∥F−F′∥ 最小的F′ 。这种方法的实现还是使用了SVD分解,
若F=(UDV)T,此时的对角矩阵D=diag(r,s,t) ,满足r≥s≥t r≥s≥tr≥s≥t,则F′=Udiag(r,s,0)VT
最小化范数∥F−F′∥ 也就是最终的解。
所以8点算法由下面两个步骤组成:

1.求线性解 由系数矩阵A最小奇异值对应的奇异矢量f 求的F。
2.奇异性约束 是最小化Frobenius范数∥F−F′∥ 代替F 。

8点算法是计算基本矩阵的最简单的方法。为了提高解的稳定性和精度,往往会对输入点集的坐标先进行归一化处理。对于归一化八点算法的总结如下:

在这里插入图片描述

3 实验及代码

3.1 基础矩阵之八点算法

# -*- coding: utf-8 -*-
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
# -*- coding: utf-8 -*-

# Read features
# 载入图像,并计算特征
im1 = array(Image.open('D:\\picture\\007\\29.jpg'))
sift.process_image('D:\\picture\\007\\29.jpg', 'D:\\picture\\007\\im1.sift')
l1, d1 = sift.read_features_from_file('D:\\picture\\007\\im1.sift')

im2 = array(Image.open('D:\\picture\\007\\30.jpg'))
sift.process_image('D:\\picture\\007\\30.jpg', 'D:\\picture\\007\\im2.sift')
l2, d2 = sift.read_features_from_file('D:\\picture\\007\\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.tools 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()

a.左右拍摄,极点位于图像平面上

sift特征匹配:
在这里插入图片描述

ransac算法:
在这里插入图片描述

8点算法:
在这里插入图片描述
7点算法:
在这里插入图片描述
10点算法:
在这里插入图片描述

基础矩阵:

8点算法
在这里插入图片描述
7点算法
在这里插入图片描述
10点算法
在这里插入图片描述
b.像平面接近平行,极点位于无穷远

sift特征匹配:
在这里插入图片描述

ransac算法:
在这里插入图片描述

8点算法:
在这里插入图片描述
基础矩阵:
在这里插入图片描述

图像拍摄位置位于前后
sift特征匹配:
在这里插入图片描述
ransac算法:
在这里插入图片描述
8点算法:
在这里插入图片描述
基础矩阵:
在这里插入图片描述

3.2 极线和极点

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt


def drawlines(img1, img2, lines, pts1, pts2):
    ''' img1 - image on which we draw the epilines for the points in img2
        lines - corresponding epilines '''
    r, c = img1.shape
    img1 = cv.cvtColor(img1, cv.COLOR_GRAY2BGR)
    img2 = cv.cvtColor(img2, cv.COLOR_GRAY2BGR)
    for r, pt1, pt2 in zip(lines, pts1, pts2):
        color = tuple(np.random.randint(0, 255, 3).tolist())
        x0, y0 = map(int, [0, -r[2] / r[1]])
        x1, y1 = map(int, [c, -(r[2] + r[0] * c) / r[1]])
        img1 = cv.line(img1, (x0, y0), (x1, y1), color, 1)
        img1 = cv.circle(img1, tuple(pt1), 5, color, -1)
        img2 = cv.circle(img2, tuple(pt2), 5, color, -1)
    return img1, img2


img1 = cv.imread('D:\\picture\\007\\29.jpg', 0)  # queryimage # left image
img2 = cv.imread('D:\\picture\\007\\30.jpg', 0)  # trainimage # right image
sift = cv.xfeatures2d.SIFT_create()
# find the keypoints and descriptors with SIFT
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)
# FLANN parameters
FLANN_INDEX_KDTREE = 1
index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks=50)
flann = cv.FlannBasedMatcher(index_params, search_params)
matches = flann.knnMatch(des1, des2, k=2)
good = []
pts1 = []
pts2 = []
# ratio test as per Lowe's paper
for i, (m, n) in enumerate(matches):
    if m.distance < 0.8 * n.distance:
        good.append(m)
        pts2.append(kp2[m.trainIdx].pt)
        pts1.append(kp1[m.queryIdx].pt)

pts1 = np.int32(pts1)
pts2 = np.int32(pts2)
F, mask = cv.findFundamentalMat(pts1, pts2, cv.FM_LMEDS)
# We select only inlier points
pts1 = pts1[mask.ravel() == 1]
pts2 = pts2[mask.ravel() == 1]

# Find epilines corresponding to points in right image (second image) and
# drawing its lines on left image
lines1 = cv.computeCorrespondEpilines(pts2.reshape(-1, 1, 2), 2, F)
lines1 = lines1.reshape(-1, 3)
img5, img6 = drawlines(img1, img2, lines1, pts1, pts2)
# Find epilines corresponding to points in left image (first image) and
# drawing its lines on right image
lines2 = cv.computeCorrespondEpilines(pts1.reshape(-1, 1, 2), 1, F)
lines2 = lines2.reshape(-1, 3)
img3, img4 = drawlines(img2, img1, lines2, pts2, pts1)
plt.subplot(121), plt.imshow(img5)
plt.subplot(122), plt.imshow(img3)
plt.show()

素材准备
在这里插入图片描述
角度不同
在这里插入图片描述
平移相机
在这里插入图片描述
远近视觉
在这里插入图片描述

在这里插入图片描述

4 实验总结

1 遇到的问题
在计算基础矩阵的代码中所拍摄的大部分图片运行时都遇到一个问题
在这里插入图片描述
意思是 不符合合适的接受标准
通过不断地尝试和运行发现这个问题是因为sift特征匹配的匹配点太少了, 有的只有几十个 ,导致无法进行后面的代码运算 ,而且匹配点太少的话, 后面你实现ransac算法的时候效果也会特别差.可以通过修改匹配特征的阙值来改善,但还是选择好的参照物来拍摄会使实验更加清晰.

2

F, ransac_data = ransac.ransac(data.T, model, 9, maxiter, match_threshold, d, return_all=True)

通过修改第三个参数来实现基础矩阵8\9\10点算法的实现

通过对比三种算法所的出来的图片的匹配点数量可知,点数越大匹配点越多,也就是说删除的匹配点越少
不同的算法算出来的基础矩阵也不一样

3 极线代码的实现可以看出效果还是不错的
但是对于我做第一组实验的图片进行极线代码运行后三组图片全部都是极线向一个中心发射的形状,也就是都和远近不同的两张图片跑出来的效果一样.我猜测这可能是匹配点过多导致的

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值