目录
1 外极几何
多视图几何是利用在不同视点所拍摄图像间的关系,来研究照相机之间或者特征之间关系的一门科学。最重要的内容是双视图几何。如果有一个场景的两个视图以及视图中的对应图像点,那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置,可以得到对这些图像点的一些几何关系约束。我们通过外极几何来描述这些几何关系。
同一图像点经过不同的投影矩阵产生不同投影点必须满足:
x
2
T
F
x
1
=
0
\mathbf{x}_2^TF\mathbf{x}_1=0
x2TFx1=0
其中:
F
=
K
2
−
T
S
t
R
K
1
−
1
F=K_2^{-T}S_tR\:K_1^{-1}
F=K2−TStRK1−1
矩阵
S
t
S_\mathrm{t}
St为反对称矩阵:
S
1
=
[
0
−
t
3
t
2
t
3
0
−
t
1
−
t
2
t
1
0
]
\boldsymbol{S}_1=\begin{bmatrix}0&-t_3&t_2\\[2ex]t_3&0&-t_1\\[2ex]-t_2&t_1&0\end{bmatrix}
S1=
0t3−t2−t30t1t2−t10
矩阵F为基础矩阵。基础矩阵可以由两照相机的参数矩阵。
可以借助F来恢复出照相机参数,F可以从对应的投影图像点计算出来。在不知道照相机内参数的情况下,仅能恢复照相机的投影变换矩阵。在知道照相机内参数的情况下,重建是基于度量的。
基础矩阵可以将对应点的搜索限制在外极线上。外极线经过点e,称为外极点。
1.1 简单数据集
从相应网址上下载一个牛津多视图数据集,加载Merton1数据并且可视化数据,将三维的点投影到一个视图。
import camera
import numpy as np
from PIL import Image
from pylab import *
im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('002.jpg'))
points2D = [loadtxt('2D/00'+str(i+1)+'.corners').T for i in range(3)]
points3D = loadtxt('3D/p3d').T
corr = genfromtxt('2D/nview-corners',dtype='int',missing_values='*')
P = [camera.Camera(loadtxt('2D/00'+str(i+1)+'.P')) for i in range(3)]
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 用Matplpotlib绘制三维数据
Matplotlib中的mplot3d工具包可以方便地绘制出三维点,线,等轮廓线,表面以及其他基本图形组件,还可以通过图像窗口控件实现三维旋转和缩放。
通过再axes对象中加上projection="3d"关键字实现三维绘图,代码如下:
import camera
import numpy as np
from PIL import Image
from pylab import *
from mpl_toolkits.mplot3d import axes3d
im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('002.jpg'))
points2D = [loadtxt('2D/00'+str(i+1)+'.corners').T for i in range(3)]
points3D = loadtxt('3D/p3d').T
corr = genfromtxt('2D/nview-corners',dtype='int',missing_values='*')
P = [camera.Camera(loadtxt('2D/00'+str(i+1)+'.P')) for i in range(3)]
fig = figure()
ax = fig.add_subplot(projection = '3d')
ax.plot(points3D[0],points3D[1],points3D[2],'k.')
X,Y,Z = axes3d.get_test_data(0.25)
ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'o')
show()
get_test_data()函数再x,y空间按照设定的空间间隔参数来产生均匀的采样点。
得到结果
注意:出现错误,FigureBase.gca() got an unexpected keyword argument ‘projection’。解决方法:使用 figure() 函数来创建一个新的figure,然后使用 add_subplot() 或者 add_axes() 函数来添加一个Axes。
1.3 计算F:八点法
八点法是通过对应点来计算基础矩阵的算法。
外极约束可以写成线性系统的形式:
[
x
2
1
x
1
1
x
2
1
y
1
1
x
2
1
w
1
1
⋯
w
2
1
w
1
1
x
2
2
x
1
2
x
2
2
y
1
2
x
2
2
w
1
2
⋯
w
2
2
w
1
2
⋮
⋮
⋮
⋱
⋮
x
2
n
x
1
n
x
2
n
y
1
n
x
2
n
w
1
n
⋯
w
2
n
w
1
n
]
[
F
11
F
12
F
13
⋮
F
33
]
=
A
f
=
0
\begin{bmatrix}x_2^1x_1^1&x_2^1y_1^1&x_2^1w_1^1&\cdots&w_2^1w_1^1\\x_2^2x_1^2&x_2^2y_1^2&x_2^2w_1^2&\cdots&w_2^2w_1^2\\\vdots&\vdots&\vdots&\ddots&\vdots\\x_2^nx_1^n&x_2^ny_1^n&x_2^nw_1^n&\cdots&w_2^nw_1^n\end{bmatrix}\begin{bmatrix}F_{11}\\F_{12}\\F_{13}\\\vdots\\F_{33}\end{bmatrix}=A\boldsymbol{f}=0
x21x11x22x12⋮x2nx1nx21y11x22y12⋮x2ny1nx21w11x22w12⋮x2nw1n⋯⋯⋱⋯w21w11w22w12⋮w2nw1n
F11F12F13⋮F33
=Af=0
其中,
f
f
f包含
F
F
F的元素,
x
1
i
=
[
x
1
i
,
y
1
i
,
w
1
i
]
x_1^i=\left[x_1^i,y_1^i,w_1^i\right]
x1i=[x1i,y1i,w1i]和
x
2
i
=
[
x
2
i
,
y
2
i
,
w
2
i
]
\mathbf{x}_2^i=\left[x_2^i,y_2^i,w_2^i\right]
x2i=[x2i,y2i,w2i]是一对图像点,共有
n
n
n 对对应点。基础矩阵中有 9 个元素,由于尺度是任意的,所以只需要 8 个方程。因为算法中需要 8 个对应点来计算基础矩阵
F
F
F,所以该算法叫做八点法。
八点法中最小化 ∣ ∣ A f ∣ ∣ ||Af|| ∣∣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[1, i], x1[2, i] * x2[2, i]]
U, S, V = linalg.svd(A)
F = V[-1].reshape(3, 3)
U, S, V = linalg.svd(F)
S[2] = 0
F = dot(U, dot(diag(S), V))
return F
通常使用SVD算法来计算最小二乘解。由于上面算法得出的解可能秩不为2,所以通过将最后一个奇异值置0来得到秩最接近2的基础矩阵。
1.4 外极点和外极线
外极点满足 F e 1 = 0 Fe_1=0 Fe1=0,可以通过计算F的零空间来得到。将函数添加到sfm.py中。
def compute_epipole(F):
U, S, V = linalg.svd(F)
e = V[-1]
return e / e[2]
在之前样本数据集运行代码如下:
import sfm
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()
在第一个视图中画出了前五个外极线,在第二个视图中画出来对应匹配点。
2 照相机和三维结构的计算
2.1 三角剖分
给定照相机参数模型,图像点可以通过三角剖分来恢复出这些点的三维位置。思想如下:
对于两个照相机
P
1
P_1
P1和
P
2
P_2
P2的视图,三维实物点
X
\mathbf{X}
X的投影点为
x
1
\mathbf{x}_1
x1和
x
2
\mathbf{x}_2
x2 (这里用齐次坐
标表示),照相机方程定义了下列关系:
[
P
1
−
x
1
0
P
2
0
−
x
2
]
[
X
λ
1
λ
2
]
=
0
\begin{bmatrix}P_1&-\mathbf{x}_1&0\\P_2&0&-\mathbf{x}_2\end{bmatrix}\begin{bmatrix}\mathbf{X}\\\lambda_1\\\lambda_2\end{bmatrix}=0
[P1P2−x100−x2]
Xλ1λ2
=0
由于图像噪声、照相机参数误差和其他系统误差,上面的方程可能没有精确解。可以通过 SVD 算法来得到三维点的最小二乘估值。
通过如下函数来计算一个点对的最小二乘三角剖分:
def triangulate_point(x1, x2, P1, P2):
M = zeros((6, 6))
M[:3, :4] = P1
M[3:, :4] = P2
M[:3, 4] = -x1
M[3:, 5] = -x2
U, S, V = linalg.svd(M)
X = V[-1, :4]
return X / X[3]
可使用下面函数实现多个点的三角剖分(输入两个图像点数组,输出为一个三维坐标数组):
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_point(x1[:, i], x2[:, i], P1, P2) for i in range(n)]
return array(X).T
可以使用下列代码实现Merton1 数据集上的三角剖分:
import sfm
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])
fig = figure()
ax = fig.add_subplot(projection = '3d')
ax.plot(Xest[0],Xest[1],Xest[2],'ko')
ax.plot(Xtrue[0],Xtrue[1],Xtrue[2],'r.')
axis('equal')
show()
首先利用前两个视图的信息来对图像点进行三角剖分,将前三个图像点的齐次坐标输出到控制台,绘制出恢复的最接近三维图像点。
打印和绘制信息如下:
估计点和实际点可以很好匹配。
2.2 由三维点计算照相机矩阵
当知道一些三维点及其图像投影,可以使用直接线性变换来计算照相机矩阵P。称为照相机反切法。利用该方法恢复照相机矩阵同样也是最小二乘问题。每个三维点
X
i
\mathbf{X}_i
Xi(齐次坐标系下)按照
λ
i
x
i
=
P
X
i
\lambda_i\mathbf{x}_i=\boldsymbol{P}\mathbf{X}_i
λixi=PXi
投影到图像点
x
i
=
[
x
i
,
y
i
,
1
]
\mathbf{x}_i=[x_i,y_i,1]
xi=[xi,yi,1],相应的点满足下面的关系:
[
X
1
T
0
0
−
x
1
0
0
⋯
0
X
1
T
0
−
y
1
0
0
⋯
0
0
X
1
T
−
1
0
0
⋯
X
2
T
0
0
0
−
x
2
0
⋯
0
X
2
T
0
0
−
y
2
0
⋯
0
0
X
2
T
0
−
1
0
⋯
⋮
⋮
⋮
⋮
⋮
⋮
⋮
]
[
p
1
T
p
2
T
p
3
T
λ
1
λ
2
⋮
]
=
0
\begin{bmatrix}\mathbf{X}_1^T&0&0&-x_1&0&0&\cdots\\0&\mathbf{X}_1^T&0&-y_1&0&0&\cdots\\0&0&\mathbf{X}_1^T&-1&0&0&\cdots\\\mathbf{X}_2^T&0&0&0&-x_2&0&\cdots\\0&\mathbf{X}_2^T&0&0&-y_2&0&\cdots\\0&0&\mathbf{X}_2^T&0&-1&0&\cdots\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}\mathbf{p}_1^T\\\mathbf{p}_2^T\\\mathbf{p}_3^T\\\lambda_1\\\lambda_2\\\vdots\end{bmatrix}=0
X1T00X2T00⋮0X1T00X2T0⋮00X1T00X2T⋮−x1−y1−1000⋮000−x2−y2−1⋮000000⋮⋯⋯⋯⋯⋯⋯⋮
p1Tp2Tp3Tλ1λ2⋮
=0
其中
p
1
,
p
2
\mathbf{p}_{1},\mathbf{p}_{2}
p1,p2和
p
3
\mathbf{p}_{3}
p3是矩阵
P
\boldsymbol{P}
P的三行。
可以使用SVD分解估计出照相机矩阵。
将下面函数添加到sfm.py文件中(输入参数为图像点和三维点):
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))
在样本数据集上测试算法性能,选出第一个视图中的一些可见点,将它们转换为齐次坐标表示,然后估计照相机矩阵:
import sfm
import camera
import numpy as np
from PIL import Image
from pylab import *
from mpl_toolkits.mplot3d import axes3d
im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('002.jpg'))
points2D = [loadtxt('2D/00'+str(i+1)+'.corners').T for i in range(3)]
points3D = loadtxt('3D/p3d').T
corr = genfromtxt('2D/nview-corners',dtype='int',missing_values='*')
P = [camera.Camera(loadtxt('2D/00'+str(i+1)+'.P')) for i in range(3)]
corr = corr[:,0]
ndx3D = where(corr>=0)[0]
ndx2D = corr[ndx3D]
x = points2D[0][:,ndx2D] # 视图 1
x = vstack( (x,ones(x.shape[1])) )
X = points3D[:,ndx3D]
X = vstack( (X,ones(X.shape[1])) )
Pest = 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()
下图为照相机矩阵投影后的结果
2.3 由基础矩阵计算照相机矩阵
1 未标定的情况——投影重建
在没有任何照相机内参数知识的情况下,照相机矩阵只能通过射影变换恢复出来,在无标定的情况下,第二个照相机矩阵可以使用一个3*3的射影变换得出。
P
2
=
[
S
e
F
∣
e
]
其中,e 是左极点,满足 e
T
F
=
0
,
S
e
是如公式所示的反对称矩阵。
P_{2}=[S_{e}F|\mathbf{e}]\\\text{其中,e 是左极点,满足 e}^{T}F=0 , S_{e} 是如公式所示的反对称矩阵。
P2=[SeF∣e]其中,e 是左极点,满足 eTF=0,Se是如公式所示的反对称矩阵。
代码如下:
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 已标定的情况——度量重建
已经标定情况下,重建会保持欧式空间中的一些度量特性。
给定标定矩阵
K
K
K,我们可以将它的逆
K
−
1
K^{-1}
K−1作用于图像点
x
κ
=
K
−
1
x
\mathbf{x}_\kappa=K^{-1}\mathbf{x}
xκ=K−1x,因此,在新的图像
坐标系下,照相机方程变为:
x
K
=
K
−
1
K
[
R
∣
t
]
X
=
[
R
∣
t
]
X
\mathbf{x}_{K}=K^{-1}K[R|\mathbf{t}]\mathbf{X}=[R|\mathbf{t}]\mathbf{X}
xK=K−1K[R∣t]X=[R∣t]X
在新的图像坐标系下,点同一满足基础矩阵方程:
x
K
2
T
F
x
K
1
=
0
\mathbf{x}_{K_2}^TF\mathbf{x}_{K_1}=0
xK2TFxK1=0
从本质矩阵恢复出的照相机矩阵中存在度量关系,但有四个可能解。
下面是计算这4个解的算法:
def compute_P_from_essential(E):
U, S, V = svd(E)
if det(dot(U, V)) < 0:
V = -V
E = 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
该函数确保本质矩阵秩为2。
3 多视图重建
假设照相机已经标定,步骤如下:
(1)检测特征点,然后在两幅图像间匹配$,
(2) 由匹配计算基础矩阵,
(3) 由基础矩阵计算照相机矩阵,
(4)三角剖分这些三维点。
当图像间的点对应包含不正确匹配时,需要一个稳健的方法来计算基础矩阵。
3.1 稳健估计基础矩阵
使用RANSAC方法,这次结合了八点算法(八点算法在平面场景中会失效)。
首先添加类到sfm.py文件中
class RansacModel(object):
def __init__(self, debug=False):
self.debug = debug
def fit(self, data):
data = data.T
x1 = data[:3, :8]
x2 = data[3:, :8]
F = compute_fundamental_normalized(x1, x2)
return F
def get_error(self, data, F):
""" Compute x^T F x for all correspondences,
return error for each transformed point. """
data = data.T
x1 = data[:3]
x2 = data[3:]
Fx1 = dot(F, x1)
Fx2 = dot(F, x2)
denom = Fx1[0] ** 2 + Fx1[1] ** 2 + Fx2[0] ** 2 + Fx2[1] ** 2
err = (diag(dot(x1.T, dot(F, x2)))) ** 2 / denom
return err
fit()方法选择8个点,使用归一化的八点算法
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]
该函数将图像点归一化为零均值固定方差。
将下列函数添加到sfm.py文件中:
def F_from_ransac(x1, x2, model, maxiter=5000, match_theshold=1e-6):
import ransac
data = vstack((x1, x2))
# compute F and return with inlier index
F, ransac_data = ransac.ransac(data.T, model, 8, maxiter, match_theshold, 20, return_all=True)
return F, ransac_data['inliers']
3.2 三维重建
代码如下:
import sfm
import camera
import numpy as np
from PIL import Image
from pylab import *
import sift
import homography
from mpl_toolkits.mplot3d import axes3d
im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('002.jpg'))
sift.process_image('001.jpg', 'im1.sift')
l0, d0 = sift.read_features_from_file('im1.sift')
sift.process_image('002.jpg', 'im2.sift')
l1, d1 = sift.read_features_from_file('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 = dot(inv(K), x1)
x2n = dot(inv(K), x2)
# 使用RANSAC方法估计E
model = sfm.RansacModel()
E, inliers = sfm.F_from_ransac(x1n, x2n, model)
# 计算照相机矩阵(P2是4个解的列表)
P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_essential(E)
# 选取点在照相机前的解
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 = i
infront = (d1 > 0) & (d2 > 0)
# 三角剖分正确点,并移除不在所有照相机前面的点
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2[ind])
X = X[:, infront]
# 绘制三维图像
from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection='3d')
ax.plot(-X[0], X[1], X[2], 'k.')
axis('off')
# 绘制X的投影 import camera
# 绘制三维点
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2[ind])
x1p = cam1.project(X)
x2p = cam2.project(X)
# 反K归一化
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()
3.3 多视图的扩展示例
1 多视图
当由同一场景的多个视图时,三维重建会变得更准确,包含更多细节信息。
对于视频序列,可以使用时序关系,在连续的帧对中匹配特征。相对方位需要从每个帧对增量地加入下一个帧对。
对于静止的图像,一个办法是找到一个中央参考视图,然后计算与之有关的所有其他照相机矩阵,另一个办法是计算一个图像对的照相机矩阵和三维重建,然后增量地加入新的图像和三维点。
2 光束法平差
多视图重建的最后一步,通常是通过优化三维点的位置和照相机参数来减少二次投影误差。
3 自标定
在未标定照相机的情形中,有时可以从图像特征来计算标定矩阵。该过程称为自标定。
4 立体图像
一个多视图成像的特殊例子是立体视觉(或者立体成像),即使用两台只有水平(向一侧)偏移的照相机观测同一场景。当照相机的位置如上设置,两幅图像具有相同的图像平面,图像的行是垂直对齐的,那么称图像对是经过矫正的。该设置在机器人学中很常见,常被称为立体平台。
通过将图像扭曲到公共的平面上,使外极线位于图像行上,任何立体照相机设置都能得到矫正(我们通常构建立体平台来产生经过矫正的图像对)。
立体重建 (有时称为致密深度重建)就是恢复深度图 (或者相反,视差图),图像中每个像素的深度 (或者视差)都需要计算出来。
- 计算视差图
在该立体重建算法中,我们将对于每个像素尝试不同的偏移,并按照局部图像周围归一化的互相关值,选择具有最好分数的偏移,然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面,所有该过程又是称为扫平面法。下面是扫平面法代码,该函数返回每个像素的最佳视差。
def plane_sweep_ncc(im_l,im_r,start,steps,wid):
m,n = im_l.shape
mean_l = np.zeros((m,n))
mean_r = np.zeros((m,n))
s = np.zeros((m,n))
s_l = np.zeros((m,n))
s_r = np.zeros((m,n))
dmaps = np.zeros((m,n,steps))
filters.uniform_filter(im_l,wid,mean_l)
filters.uniform_filter(im_r,wid,mean_r)
norm_l = im_l - mean_l
norm_r = im_r - mean_r
for displ in range(steps):
filters.uniform_filter(np.roll(norm_l,-displ-start)*norm_r,wid,s)
filters.uniform_filter(np.roll(norm_l,-displ-start)*np.roll(norm_l,-displ-start),wid,s_l)
filters.uniform_filter(norm_r*norm_r,wid,s_r)
dmaps[:,:,displ] = s/np.sqrt(s_l*s_r)
return np.argmax(dmaps,axis=2)