文章目錄

• 多視圖幾何
• • 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()函數，就可以在立體表面畫出這些三維點。

圖2 立體表面的三維點

通過畫出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為真，外極點也會被畫出來（如果輸入參數沒有外極點，則外極點會在程序中計算出來）。
可以在之前樣本數據集的前兩個視圖上運行這兩個函數：

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])))F = sfm.compute_fundamental(x1, x2)e = sfm.compute_epipole(F)figure() imshow(im1)for i in range(5):sfm.plot_epipolar_line(im1, F, x2[:, i], e, False) axis('off')figure() imshow(im2)for i in range(5):plot(x2[0, i], x2[1, i], 'o') axis('off')show()

首先選擇兩幅圖像的對應點，然后將他們轉換為齊次坐標。然后畫出了對應匹配點。在下圖中用不同顏色將點和相應的外極線對應起來。

圖4 數據視圖中的5個外極點

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. ]]
算法估計出的三維圖像點和實際圖像點很接近。

圖5 使用照相機矩陣和點的對應關系三角剖分數據點

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]]

圖6 使用估計出的照相機矩陣計算視圖1中的投影點
可以看出估計出的照相機矩陣和數據集計算出的照相機矩陣，它們的元素幾乎完全相同。最后使用估計出的照相機矩陣投影這些三維點，然后繪制出投影后的結果，如上圖所示，真實點用圓圈表示，估計出的照相機投影點用點表示。

2.3 由基礎矩陣計算照相機矩陣

在兩個視圖場景中，照相機矩陣可以由基礎矩陣恢復出來。研究分為兩類：未標定的情況和已標定的情況。

2.3.1 投影重建

在沒有任何照相機內參數知識的情況下，照相機矩陣只能通過射影變換恢復出來。也就是說，如果利用照相機信息來重建三維點，那么該重建只能由射影變換計算出來。（在這里不考慮角度和距離）
具體代碼如下：

def skew(a):return array([[0, -a[2], a[1]], [a[2], 0, -a[0]], [-a[1], a[0], 0]])def compute_P_from_fundamental(F):e = compute_epipole(F.T)Te = skew(e)return vstack((dot(Te, F.T).T, e)).T

2.3.1 度量重建

在已標定的情況下，重建會保持歐式空間中的一些度量特性（除了全局的尺度參數）。從本質矩陣恢復出的照相機矩陣中存在度量關系，但有四個可能解。因為只有一個解產生位于兩個照相機前的場景，所以可以輕松選出。
具體代碼如下：

def compute_P_from_essential(E):U, S, V = svd(E)if det(dot(U, V)) < 0:V = -VE = dot(U, dot(diag([1, 1, 0]), V))Z = skew([0, 0, -1])W = array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])P2 = [vstack((dot(U, dot(W, V)).T, U[:, 2])).T,vstack((dot(U, dot(W, V)).T, -U[:, 2])).T,vstack((dot(U, dot(W.T, V)).T, U[:, 2])).T,vstack((dot(U, dot(W.T, V)).T, -U[:, 2])).T]return P2

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()

將這些三維點投影后，可以通過乘以標定矩陣的方式來彌補初始歸一化對點的影響。
程序在運行中出現一個問題：

圖7 三維重建中遇到的問題
最后得到的結果如下圖所示：
圖8 使用圖像匹配從一對圖像中計算三維重建
從圖中可以看到，二次投影后的點和原始特征位置不完全匹配，但是相當接近。可以進一步調整照相機矩陣來提高重建和二次投影的性能。

4 立體圖像

一個多視圖成像的特殊例子是立體視覺，即使用兩臺只有水平偏移的照相機觀測同一場景。當照相機的位置如上設置，兩幅圖像具有相同的圖像平面，圖像的行是垂直對齊的，那么稱圖像是經過矯正的。該設置在機器人學中很常見，被稱為立體平臺。

立體重建就是恢復深度圖，圖像中每個像素的深度都需要計算出來。在計算視差圖這一立體重建算法中，要對每個像素嘗試不同的偏移，并按照局部圖像周圍歸一化的互相關值選擇具有最好分數的偏移，然后記錄下該最佳偏移。因為每個偏移在某種程度上對應于一個平面，所以該過程稱為掃平面法。我們使每個像素周圍的圖像塊來計算歸一化的互相關，然后對該像素周圍局部像素塊中的像素求和。這一部分使用圖像濾波器來快速計算。

總結

以上是生活随笔為你收集整理的计算机视觉编程——多视图几何的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。