技術文章
視覺定位系統原理——基礎矩陣
閱讀:195 發布時間:2022-4-1一、視覺定位系統基礎矩陣的原理
介紹基礎矩陣我們先說一下對積幾和
對極幾何一定是對二幅圖像而言,對極幾何實際上是“兩幅圖像之間的對極幾何”,它是圖像平面與以基線為軸的平面束的交的幾何(這里的基線是指連接攝像機中心的直線),以下圖為例:對極幾何描述的是左右兩幅圖像(點x和x’對應的圖像)與以CC’為軸的平面束的交的幾何!
直線CC’為基線,以該基線為軸存在一個平面束,該平面束與兩幅圖像平面相交,下圖給出了該平面束的直觀形象,可以看到,該平面束中不同平面與兩幅圖像相交于不同直線;
上圖中的灰色平面π,只是過基線的平面束中的一個平面(當然,該平面才是平面束中最重要的、也是我們要研究的平面);
epipolarpoints極點
每一個相機的透鏡中心是不同的,會投射到另一個相機像面的不同點上。這兩個像點用eL和eR表示,被稱為epipolarpoints極點。兩個極點eL、eR分別與透鏡中心OL、OR在空間中位于一條直線上。
epipolarplane極面
將X、OL和OR三點形成的面稱為epipolarplane極面。
epipolarline極線
直線OL-X被左相機看做一個點,因為它和透鏡中心位于一條線上。然而,從右相機看直線OL-X,則是像面上的一條線直線eR-XR,被稱為epipolarline極線。從另一個角度看,極面X-OL-OR與相機像面相交形成極線。
極線是3D空間中點X的位置函數,隨X變化,兩幅圖像會生成一組極線。直線OL-X通過透鏡中心OL,右像面中對應的極線必然通過極點eR。一幅圖像中的所有極線包含了該圖像的所有極點。實際上,任意一條包含極點的線都是由空間中某一點X推導出的一條極線。
如果兩個相機位置已知,則:
1.如果投影點XL已知,則極線eR-XR已知,點X必定投影到右像面極線上的XR處。這就意味著,在一個圖像中觀察到的每個點,在已知的極線上觀察到該點的其他圖像。這就是Epipolarconstraint極線約束:X在右像面上的投影XR必然被約束在eR-XR極線上。對于OL-XL上的X,X1,X2,X3都受該約束。極線約束可以用于測試兩點是否對應同一3D點。極線約束也可以用兩相機間的基本矩陣來描述。
2.如果XL和XR已知,他們的投影線已知。如果兩幅圖像上的點對應同一點X,則投影線必然交于X。這就意味著X可以用兩個像點的坐標計算得到。
二、視覺定位系統基礎矩陣
如果已知基礎矩陣F,以及一個3D點在一個像面上的像素坐標p,則可以求得在另一個像面上的像素坐標p‘。這個是基礎矩陣的作用,可以表征兩個相機的相對位置及相機內參數。
面具體介紹基礎矩陣與像素坐標p和p’的關系。
以O1為原點,光軸方向為z軸,另外兩個方向為x,y軸可以得到一個坐標系,在這個坐標系下,可以對P,p1(即圖中所標p),p2(即圖中所標p‘)得到三維坐標,同理,對O2也可以得到一個三維坐標,這兩個坐標之間的轉換矩陣為[RT],即通過旋轉R和平移T可以將O1坐標系下的點P1(x1,y1,z1),轉換成O2坐標系下的P2(x2,y2,z2)。
則可知,P2=R(P1-T)(1)
采用簡單的立體幾何知識,可以知道

其中,p,p‘分別為P點的像點在兩個坐標系下分別得到的坐標(非二維像素坐標)。Rp’為極面上一矢量,T為極面上一矢量,則兩矢量一叉乘為極面的法向量,這個法向量與極面上一矢量p一定是垂直的,所以上式一定成立。(這里采用轉置是因為p會表示為列向量的形式,此處需要為行向量)
采用一個常用的叉乘轉矩陣的方法,

將我們的叉乘采用上面的轉換,會變成

紅框中所標即為本征矩陣E,他描述了三維像點p和p‘之間的關系

有了本征矩陣,我們的基礎矩陣也就容易推導了,
注意到將p和p‘換成P1和P2式(4)也是成立的,且有
q1=K1P1(6)
q2=K2P2(7)
上式中,K1K2為相機的校準矩陣,描述相機的內參數q1q2為相機的像素坐標代入式(4)中,得

上式中p->q1,p‘->q2
這樣我們就得到了兩個相機上的像素坐標和基礎矩陣F之間的關系了

二、不同圖像對的基礎矩陣
代碼
# coding: utf-8
# In[1]:
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
# In[2]:
from PCV.geometry import homography, camera, sfm
from PCV.localdescriptors import sift
# In[5]:
# Read features
im1 = array(Image.open('im1.jpg'))
sift.process_image('im1.jpg', 'im1.sift')
im2 = array(Image.open('im2.jpg'))
sift.process_image('im2.jpg.', 'im2.sift')
# In[6]:
l1, d1 = sift.read_features_from_file('im1.sift')
l2, d2 = sift.read_features_from_file('im2.sift')
# In[7]:
matches = sift.match_twosided(d1, d2)
# In[8]:
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()
# In[9]:
figure(figsize=(16,16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()
# In[10]:
#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
/Cookbook/RANSAC).
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
# 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']
# In[11]:
# find F through RANSAC
model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-3)
print(F)
# In[12]:
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = pute_P_from_fundamental(F)
# In[13]:
print(P2)
print(F)
# In[14]:
# P2, F (1e-4, d=20)
# [[ -1.e+00 1.e+01 5.e+02 4.e+03]
# [ 1.e+01 -9.e+01 -4.e+03 5.e+02]
# [ 2.e-02 3.e-03 -1.e+02 1.e+00]]
# [[ -1.e-07 4.e-06 -2.e-03]
# [ -1.e-06 6.e-07 2.e-02]
# [ 1.e-03 -2.e-02 1.e+00]]
# In[15]:
# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)
# In[16]:
# plot the projection of X
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)
x1p = cam1.project(X)
x2p = cam2.project(X)
# In[17]:
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]
plot([x1p[0][i], x2p[0][i]+cols1],[x1p[1][i], x2p[1][i]],'c')
axis('off')
show()
# In[18]:
d1p = d1n[inliers]
d2p = d2n[inliers]
# In[19]:
# Read features
im3 = array(Image.open('im3.jpg'))
sift.process_image('im3.jpg', 'im3.sift')
l3, d3 = sift.read_features_from_file('im3.sift')
# In[20]:
matches13 = sift.match_twosided(d1p, d3)
# In[21]:
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)
# In[22]:
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]
plot([x1_13[0][i], x3_13[0][i]+cols1],[x1_13[1][i], x3_13[1][i]],'c')
axis('off')
show()
# In[23]:
P3 = pute_P(x3_13, X[:, ndx_13])
# In[24]:
print(P3)
# In[25]:
print(P1)
print(P2)
print(P3)
# In[26]:
# Can't tell the camera position because there's no calibration matrix (K)
得到視覺對位系統基礎矩陣如下

室外

得到基礎矩陣如下
