计算机视觉编程——多视图几何
文章目錄
- 多視圖幾何
- 1 外極幾何
- 1.1 一個簡單的數據集
- 1.2 用Matplotlib繪制三維數據
- 1.3 計算F：八點法
- 1.4 外極點和外極線
- 2 照相機和三維結構的計算
- 2.1 三角剖分
- 2.2 由三維點計算照相機矩陣
- 2.3 由基礎矩陣計算照相機矩陣
- 2.3.1 投影重建
- 2.3.1 度量重建
- 3 多視圖重建
- 3.1 穩健估計基礎矩陣
- 3.2 三維重建示例
- 4 立體圖像
多視圖幾何
多視圖幾何是利用在不同視點所拍攝的圖像間的關系，來研究照相機之間或者特征之間關系的一門學科。圖像的特征通常就是興趣點。
1 外極幾何
如果一個場景的兩個視圖以及視圖中的對應圖像點，那么根據照相機間的空間相對位置關系、照相機的性質以及三維場景點的位置，可以得到對這些圖像點的幾何約束關系，即外極幾何。
1.1 一個簡單的數據集
使用下面的程序加載兩個圖像、三個視圖中的所有圖像特征點、對應不同視圖圖像點重建后的三維點以及照相機參數矩陣。
import my_cameraim1 = array(Image.open('001.jpg')) im1 = array(Image.open('002.jpg'))points2D = [loadtxt('2D/00' + str(i + 1) + '.corners').T for i in range(3)]points3D = loadtxt('3D/p3d').Tcorr = genfromtxt('2D/nview-corners', dtype = 'int', missing = '*')P = [my_camera.Camera(loadtxt('2D/00' + str(i + 1) + '.P')) for i in range(3)]使用loadtxt()函數讀取文本文件到NumPy數組中，因為并不是所有的點都可見，或都能夠成功匹配到所有的視圖，所以對應數據里包含了缺失的數據。genfromtxt()函數通過將缺失的數值填充為-1來解決這個問題。
將上述代碼保存到load_vggdata.py，然后使用命令excefile()可以運行該腳本，從而獲取數據：
execfile('load_vggdata.py')結果報錯：NameError: name ‘execfile’ is not defined
這是由于execfile()在python3中已被廢除，可以使用代替函數： exec(open(filename).read())來完成。
下面來可視化這些數據，將三維的點投影到一個視圖上，然后和觀測到的圖像點比較：
X = vstack((points3D, ones(points3D.shape[1]))) x = P[0].project(X)figure() imshow(im1) plot(points2D[0][0], points2D[0][1], '*') axis('off')figure() imshow(im1) plot(x[0], x[1], 'r.') axis('off')show() 圖1 圖像點與三維投影點上邊的代碼繪制出第一個視圖以及該視圖中的圖像點。仔細觀察發現第二幅圖比第一幅圖多一些點，這些多出來的點是從視圖2和視圖3重建出來的，而不是在視圖1中。
1.2 用Matplotlib繪制三維數據
為了可視化三維重建結果，需要繪制出三維圖像。Matplotlib中的mplot3d工具包可以方便地繪制出三維點、線、等輪廓線、表面以及其他基本圖形組件，還可以通過圖像窗口控件實現三維旋轉和縮放。
可以通過在axes對象中加上projection=“3d”關鍵字實現三維繪圖，如下：
from mpl_toolkits.mplot3d import axes3dfig = figure() ax = fig.gca(projection = "3d")X, Y, Z = axes3d.get_test_data(0.25)ax.plot(X.flatten(), Y.flatten(), Z.flatten(), 'o')show()get_test_data()函數在x，y空間按照設定的空間間隔參數來產生均勻的采樣點。壓平這些網格會產生三列數據點，然后可以將其輸入plot()函數，就可以在立體表面畫出這些三維點。
通過畫出Merton樣本數據來觀察三維點的效果：
fig = figure() ax = fig.gca(projection = "3d") ax.plot(points3D[0], points3D[1], points3D[2], 'k.')show() 圖3 建筑墻體和屋頂上的三維點1.3 計算F：八點法
八點法是通過對應點來計算基礎矩陣的算法。
下面的代碼寫入八點法中最小化||Af||的函數：
def compute_fundamental(x1, x2):n = x1.shape[1]if x2.shape[1] != n:raise ValueError("Number of points don't match.")A = zeros((n, 9))for i in range(n):A[i] = [x1[0, i] * x2[0, i], x1[0, i] * x2[1, i], x1[0, i] * x2[2, i],x1[1, i] * x2[0, i], x1[1, i] * x2[1, i], x1[1, i] * x2[2, i],x1[2, i] * x2[0, i], x1[2, i] * x2[q, i], x1[2, i] * x2[0, i]]U, S, V = linalg.svd(A)F = V[-1].reshape(3, 3)U, S, V = linalg.svd(F)S[2] = 0F = dot(U, dot(diag(S), V))return F通常使用SVD算法來計算最小二乘解。由于上面算法得出的解可能秩不為2，所以通過將最后一個奇異值置為0來得到秩最接近2的基礎矩陣。
1.4 外極點和外極線
如果要獲得一幅圖像的外極點，只需要將F轉置后輸入下述函數：
def compute_epipole(F):U, S, V = linalg.svd(F)e = V[-1]return e / e[2] def plot_epipolar_line(im, F, x, epipole = None, show_epipole = True):m, n = im.shape[:2]line = dot(F, x)t = linspace(0, n, 100)lt = array([(line[2] + line[0] * tt) / (-line[1]) for tt in t])ndx = (lt >= 0) & (lt < m)plot(t[ndx], lt[ndx], linewidth = 2)if show_epipole:if epipole is None:epipole = compute_epipole(F)plot(epipole[0] / epipole[2], epipole[1] / epipole[2], 'r*')上面的函數將x軸的范圍作為直線的參數，因此直線超出圖像邊界的部分會被截斷。如果show_epipole為真，外極點也會被畫出來（如果輸入參數沒有外極點，則外極點會在程序中計算出來）。
可以在之前樣本數據集的前兩個視圖上運行這兩個函數：
首先選擇兩幅圖像的對應點，然后將他們轉換為齊次坐標。然后畫出了對應匹配點。在下圖中用不同顏色將點和相應的外極線對應起來。
2 照相機和三維結構的計算
2.1 三角剖分
給定照相機參數模型，圖像點可以通過三角剖分來恢復出這些點的三維位置。基本算法思想如下：
對于兩個照相機P1，P2的視圖，三維實物點X的投影點x1和x2，照相機方程定義了下列關系：
由于圖像噪聲、照相機參數誤差和其他系統誤差，上面的方程可能沒有精確解。需要用到SVD算法來得到三維點的最小二乘估值。
下面的函數用于計算一個點對的最小二乘三角剖分：
def triangulate_points(x1, x2, P1, P2)：M = zeros((6, 6))M[:3, :4] = P1M[3:, :4] = P2M[:3, 4] = -x1M[3:, 5] = -x2U, S, V = linalg.svd(M)X = V[-1, :4]return X / X[3]最后一個特征向量的前4個值就是齊次坐標系下的對應三維坐標，可以增加下面的函數來實現多個點的三角剖分：
def triangulate(x1, x2, P1, P2):n = x1.shape[1]if x2.shape[1] != n:raise ValueError("Number of points don't match.")X = [triangulate_points(x1[:, i], x2[:, i], P1, P2) for i in range(n)]return array(X).T這個函數的輸入是兩個圖像點數組，輸出為一個三維坐標數組。
可以利用下面的代碼來實現三角剖分：
ndx = (corr[:, 0] >= 0) & (corr[:, 1] >= 0)x1 = points2D[0][:, corr[ndx, 0]] x1 = vstack((x1, ones(x1.shape[1]))) x2 = points2D[1][:, corr[ndx, 1]] x2 = vstack((x2, ones(x2.shape[1])))Xtrue = points3D[:, ndx] Xtrue = vstack((Xtrue, ones(Xtrue.shape[1])))Xest = sfm.triangulate(x1, x2, P[0].P, P[1],P) print(Xest[:, :3]) print(Xtrue[:, :3])from mpl_toolkits.mplot3d import axes3dfig = figure() ax = fig.gca(projection = '3d') ax.plot(Xest[0], Xest[1], Xest[2], 'ko') ax.plot(Xtrue[0], Xtrue[1], Xtrue[2], 'r.')axis('equal')show()上面的代碼首先利用前兩個視圖的信息來對三維點進行三角剖分，然后把前三個圖像點的齊次坐標輸出到控制臺，最后繪制出灰度的最接近三維圖像點。輸出到控制臺的信息如下：
[[ 1.03743725 1.56125273 1.40720017]
[-0.57574987 -0.55504127 -0.46523952]
[ 3.44173797 3.44249282 7.53176488]
[ 1. 1. 1. ]]
[[ 1.0378863 1.5606923 1.4071907 ]
[-0.54627892 -0.5211711 -0.46371818]
[ 3.4601538 3.4636809 7.5323397 ]
[ 1. 1. 1. ]]
算法估計出的三維圖像點和實際圖像點很接近。
2.2 由三維點計算照相機矩陣
如果已經知道了一些三維點及其圖像投影，可以使用直接線性變換的方法來計算照相機矩陣P。本質上這也是三角剖分方法的逆問題，也稱照相機反切法。
相應代碼如下：
def compute_P(x, X):n = x.shape[1]if X.shape[1] != n:raise ValueError("Number of points don't match.")M = zeros((3 * n, 12 + n))for i in range(n):M[3 * i, 0: 4] = X[:, i]M[3 * i + 1, 4: 8] = X[:, i]M[3 * i + 2, 8: 12] = X[:, i]M[3 * i: 3 * i + 3, i + 12] = -x[:, i]U, S, V = linalg.svd(M)return V[-1, :12].reshape((3, 4))該函數的輸入參數為圖像點和三維點，構造出上述所示的M矩陣。最后一個特征向量的前12個元素是照相機矩陣的元素，經過重新排列成矩陣形狀后返回。
下面，在樣本數據集上測試算法的性能。
corr = corr[:, 0] ndx3D = where(corr >= 0)[0] ndx2D = corr[ndx3D]x = points2D[0][:, ndx2D] x = vstack((x, ones(x.shape[1]))) X = points3D[:, ndx3D] X = vstack((X, ones(X.shape[1])))Pest = my_camera.Camera(sfm.compute_P(x, X))print(Pest.P / Pest.P[2, 3]) print(P[0].P / P[0].P[2, 3])xest = Pest.project(X)figure() imshow(im1) plot(x[0], x[1], 'bo') plot(xest[0], xest[1], 'r.') axis('off')show()上面的代碼選出第一個視圖中的一些可見點，將它們轉換為齊次坐標表示，然后估計照相機矩陣。為了檢查照相機矩陣的正確性，將他們歸一化的格式（除以最有一個元素）打印到控制臺。輸出如下所示：
[[ 1.06520794e+00 -5.23431275e+01 2.06902749e+01 5.08729305e+02]
[-5.05773115e+01 -1.33243276e+01 -1.47388537e+01 4.79178838e+02]
[ 3.05121915e-03 -3.19264684e-02 -3.43703738e-02 1.00000000e+00]]
[[ 1.06774679e+00 -5.23448212e+01 2.06926980e+01 5.08764487e+02]
[-5.05834364e+01 -1.33201976e+01 -1.47406641e+01 4.79228998e+02]
[ 3.06792659e-03 -3.19008054e-02 -3.43665129e-02 1.00000000e+00]]
可以看出估計出的照相機矩陣和數據集計算出的照相機矩陣，它們的元素幾乎完全相同。最后使用估計出的照相機矩陣投影這些三維點，然后繪制出投影后的結果，如上圖所示，真實點用圓圈表示，估計出的照相機投影點用點表示。
2.3 由基礎矩陣計算照相機矩陣
在兩個視圖場景中，照相機矩陣可以由基礎矩陣恢復出來。研究分為兩類：未標定的情況和已標定的情況。
2.3.1 投影重建
在沒有任何照相機內參數知識的情況下，照相機矩陣只能通過射影變換恢復出來。也就是說，如果利用照相機信息來重建三維點，那么該重建只能由射影變換計算出來。（在這里不考慮角度和距離）
具體代碼如下：
2.3.1 度量重建
在已標定的情況下，重建會保持歐式空間中的一些度量特性（除了全局的尺度參數）。從本質矩陣恢復出的照相機矩陣中存在度量關系，但有四個可能解。因為只有一個解產生位于兩個照相機前的場景，所以可以輕松選出。
具體代碼如下：
3 多視圖重建
由于照相機的運動給我們提供了三維結構，所以這樣計算三維重建的方法通常稱為SfM。
假設照相機已經標定，計算重建可以分為下面4個步驟：
3.1 穩健估計基礎矩陣
這部分類似于穩健計算單應性矩陣，當存在噪聲和不正確的匹配時，需要估計基礎矩陣。使用RANSAC方法并結合八點法來完成。
class RansacModel(object):def __init__(self, debug = False):self.debug = debugdef fit(self, data):data = data.Tx1 = data[:3, :8]x2 = data[3:, :8]F = compute_fundamental_normalized(x1, x2)return Fdef get_error(self, data, F):data = data.Tx1 = data[:3]x2 = data[3:]Fx1 = dot(F, x1)Fx2 = dot(F, x2)denom = Fx1[0] ** 2 + Fx[1] ** 2 + Fx2[0] ** 2 + Fx2[1] ** 2err = (diag(dot(x1.T, dot(F, x2)))) ** 2 / denomreturn err def compute_fundamental_normalized(x1, x2):n = x1.shape[1]if x2.shape[1] != n:raise ValueError("Number of points don't match.")x1 = x1 / x1[2]mean_1 = mean(x1[:2], axis = 1)S1 = sqrt(2) / std(x1[:2])T1 = array([[S1, 0, -S1 * mean_1[0]], [0, S1, -S1 * mean_1[1]], [0, 0, 1]])x1 = dot(T1, x1)x2 = x2 / x2[2]mean_2 = mean(x2[:2], axis = 1)S2 = sqrt(2) / std(x2[:2])T2 = array([[S2, 0, -S2 * mean_2[0]], [0, S2, -S2 * mean_2[1]], [0, 0, 1]])x2 = dot(T2, x2)F = compute_fundamental(x1, x2)F = dot(T1.T, dot(F, T2))return F / F[2, 2] def F_from_ransac(x1, x2, model, maxiter = 5000, match_theshold = 1e-6):import ransacdata = vstack((x1, x2))F, ransac_data = ransac.ransac(data.T, model, 8, maxiter, match_theshold, 20, return_all = True)return F, ransac_data['inliers']這里返回最佳基礎矩陣F，以及正確點的索引，可以知道那些匹配和F矩陣是一致的。與單應性矩陣估計相比，增加了默認的最大迭代次數，改變了匹配的閾值。
3.2 三維重建示例
在接下來的處理中，將該示例的處理代碼分成若干塊，使得代碼更容易理解。首先，提取、匹配特征，然后估計基礎矩陣和照相機矩陣：
K = array([[2394, 0, 932], [0, 2398, 628], [0, 0, 1]])im1 = array(Image.open('alcatraz1.jpg')) my_sift.process_image('alcatraz1.jpg', 'im1.sift') l1, d1 = my_sift.read_features_from_file('im1.sift')im2 = array(Image.open('alcatraz2.jpg')) my_sift.process_image('alcatraz2.jpg', 'im2.sift') l2, d2 = my_sift.read_features_from_file('im2.sift')matches = my_sift.match_twosided(d1, d2) ndx = matches.nonzero()[0]x1 = my_homography.make_homog(l1[ndx, :2].T) ndx2 = [int(matches[i]) for i in ndx] x2 = my_homography.make_homog(l2[ndx2, :2].T)x1n = dot(inv(K), x1) x2n = dot(inv(K), x2)model = sfm.RansacModel() E, inliers = sfm.F_from_ransac(x1n, x2n, model)P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]]) P2 = sfm.compute_P_from_essential(E)在獲得了標定矩陣后，對矩陣K進行硬編碼。使用K的逆矩陣來對其進行歸一化，并使用歸一化的八點算法來運行Ransac估計，返回一個本質矩陣并保存正確匹配點的索引。從本質矩陣出發，可以計算出第二個照相機矩陣的四個可能解。
從照相機矩陣的列表中，挑選出經過三角剖分后，在兩個照相機前郡含有場景點的照相機矩陣：
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 = iinfront = (d1 > 0) & (d2 > 0)X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2[ind])X = X[:, infront]循環遍歷這四個解，每次對對應于正確點的三維點進行三角剖分。將三角剖分后的圖像投影回圖像后，深度的符號由每個圖像點的第三個數值給出。我們保存了正向最大深度的索引；對于和最優解一致的點，用相應的布爾變量保存了信息，這樣可以取出真正在照相機前面的點。因為所有估計中都存在噪聲和誤差，所以即便使用正確的照相機矩陣，也存在一些點仍位于某個照相機后面的風險。首先獲得正確的解，然后對這些正確的點進行三角剖分，最后保留位于照相機前方的點。
現在可以繪制出該三維重建：
cam1 = my_camera.Camera(P1) cam2 = my_camera.Camera(P2)x1p = cam1.project(X) x2p = cam2.project(X)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()將這些三維點投影后，可以通過乘以標定矩陣的方式來彌補初始歸一化對點的影響。
程序在運行中出現一個問題：
最后得到的結果如下圖所示：
圖8 使用圖像匹配從一對圖像中計算三維重建
從圖中可以看到，二次投影后的點和原始特征位置不完全匹配，但是相當接近。可以進一步調整照相機矩陣來提高重建和二次投影的性能。
4 立體圖像
一個多視圖成像的特殊例子是立體視覺，即使用兩臺只有水平偏移的照相機觀測同一場景。當照相機的位置如上設置，兩幅圖像具有相同的圖像平面，圖像的行是垂直對齊的，那么稱圖像是經過矯正的。該設置在機器人學中很常見，被稱為立體平臺。
立體重建就是恢復深度圖，圖像中每個像素的深度都需要計算出來。在計算視差圖這一立體重建算法中，要對每個像素嘗試不同的偏移，并按照局部圖像周圍歸一化的互相關值選擇具有最好分數的偏移，然后記錄下該最佳偏移。因為每個偏移在某種程度上對應于一個平面，所以該過程稱為掃平面法。我們使每個像素周圍的圖像塊來計算歸一化的互相關，然后對該像素周圍局部像素塊中的像素求和。這一部分使用圖像濾波器來快速計算。
總結
以上是生活随笔為你收集整理的计算机视觉编程——多视图几何的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 计算机视觉编程——增强现实基础
- 下一篇： 计算机视觉编程——图像聚类