多视几何
多視圖幾何是主要研究用幾何的方法,通過(guò)若干幅多視圖幾何二維圖像,來(lái)恢復(fù)三維物體。間言之就是研究三維重構(gòu)。
1.對(duì)極幾何基礎(chǔ)概念
核點(diǎn)(epipole):基線(baseline)與成像平面的交點(diǎn)。同時(shí)極點(diǎn)也可以理解為相鄰影像成像中心在本影像上的像,因?yàn)榛€是兩個(gè)相鄰影像成像中心的連線。
核平面(epipolar plane):含有基線的平面,是一簇平面。可以看做是由基線與空間中任意一點(diǎn)構(gòu)成的平面。
核線(epipolar line):核平面與成像平面的交線。可以看做是成像平面上的任意一點(diǎn)(非核點(diǎn))與核點(diǎn)所定義的直線。
2. 基礎(chǔ)矩陣 F
基礎(chǔ)矩陣可以看做是將點(diǎn)投影(轉(zhuǎn)換)為直線,將左影像上的一個(gè)點(diǎn)投影到右影像上形成一條核線。
基礎(chǔ)矩陣表示的是圖像中的像點(diǎn)p1到另一幅圖像對(duì)極線l2的映射,有如下公式:
而和像點(diǎn)P1P1匹配的另一個(gè)像點(diǎn)p2必定在對(duì)極線l2上,所以就有:
這樣僅通過(guò)匹配的點(diǎn)對(duì),就可以計(jì)算出兩視圖的基礎(chǔ)矩陣F。
基礎(chǔ)矩陣FF是一個(gè)3×3的矩陣,有9個(gè)未知元素,由于尺度是任意的,所以只需要8個(gè)方程。因?yàn)樗惴ㄖ行枰?個(gè)對(duì)應(yīng)點(diǎn)來(lái)計(jì)算基礎(chǔ)矩陣F,所以該算法叫做八點(diǎn)法。
3.8點(diǎn)法估算基礎(chǔ)矩陣F
由于基礎(chǔ)矩陣 F 定義為
任給兩幅圖像中的匹配點(diǎn) x 與 x’,令
有相應(yīng)方程:
最后得到AF=0.
4.實(shí)驗(yàn)過(guò)程
#!/usr/bin/env python # coding: utf-8# In[1]: from PIL import Image from numpy import * from pylab import * import numpy as np from imp import reload# In[2]: import importlib from PCV.geometry import camera from PCV.geometry import homography from PCV.geometry import sfm from PCV.localdescriptors import siftcamera = importlib.reload(camera) homography = importlib.reload(homography) sfm = importlib.reload(sfm) sift =importlib.reload(sift)# In[3]:# Read features im1 = array(Image.open('image5/1.jpg')) sift.process_image('image5/1.jpg', 'im1.sift') l1, d1 = sift.read_features_from_file('im1.sift')im2 = array(Image.open('image5/2.jpg')) sift.process_image('image5/2.jpg', 'im2.sift') l2, d2 = sift.read_features_from_file('im2.sift')# In[9]:matches = sift.match_twosided(d1, d2)# In[10]: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)x1n = x1.copy() x2n = x2.copy()# In[11]:print (len(ndx))# In[12]:figure(figsize=(16,16)) sift.plot_matches(im1, im2, l1, l2, matches, True) show()# In[13]:# Chapter 5 Exercise 1 # 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 pointcorrespondences using RANSAC (ransac.py fromhttp://www.scipy.org/Cookbook/RANSAC).input: x1, x2 (3*n arrays) points in hom. coordinates. """from PCV.tools import ransacdata = np.vstack((x1, x2))d = 20 # 20 is the original# compute F and return with inlier indexF, ransac_data = ransac.ransac(data.T, model,8, maxiter, match_threshold, d, return_all=True)return F, ransac_data['inliers']# In[15]:# find E through RANSAC model = sfm.RansacModel() F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-3)# In[16]:print (len(x1n[0])) print (len(inliers))# In[17]:P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]]) P2 = sfm.compute_P_from_fundamental(F)# In[18]:# triangulate inliers and remove points not in front of both cameras X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)# In[19]:# plot the projection of X cam1 = camera.Camera(P1) cam2 = camera.Camera(P2) x1p = cam1.project(X) x2p = cam2.project(X)# In[20]: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()# In[21]: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()# In[22]:print (F)# In[23]:# Chapter 5 Exercise 2x1e = [] 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')p1 = array([l1[i][0], l1[i][1], 1])p2 = array([l2[m][0], l2[m][1], 1])# Use Sampson distance as errorFx1 = dot(F, p1)Fx2 = dot(F, p2)denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2e = (dot(p1.T, dot(F, p2)))**2 / denomx1e.append([p1[0], p1[1]])x2e.append([p2[0], p2[1]])ers.append(e) x1e = array(x1e) x2e = array(x2e) ers = array(ers)# In[24]:indices = np.argsort(ers) x1s = x1e[indices] x2s = x2e[indices] ers = ers[indices] x1s = x1s[:20] x2s = x2s[:20]# In[25]: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()# In[ ]:室外結(jié)果如下:
注:拍攝照片時(shí)要注意圖片的景深,否則運(yùn)行中很容易出現(xiàn)問(wèn)題。
室外結(jié)果如下圖
得到的基礎(chǔ)矩陣為
參考博客:https://www.cnblogs.com/wangguchangqing/p/8214032.html
https://www.cnblogs.com/JingeTU/p/6390915.html
總結(jié)
- 上一篇: 关于Java单例模式的思考
- 下一篇: IoT原型开发利用现成的单板设计---凯