python计算机视觉编程_计算机视觉编程系列-多视图几何

这篇文章主要讲解如何处理多个视图，以及如何利用多个视图的几何关系来恢复照相机位置信息和三维结构。通过在不同视点拍摄的图像，我们可以利用特征匹配来计算出三维场景点以及照相机位置。这篇文章会介绍一些基本的方法，展示一个三维重建的完整例子。

5.1 外极几何

多视图几何是利用在不同视点所拍摄图像间的关系，来研究照相机之间或者特征之间关系的一门科学。图像的特征通常是兴趣点，这篇文章使用的也是兴趣点特征。多视图几何中最重要的内容是双视图几何。

如果有一个场景的两个视图以及视图中的对应点图像点，那么根据照相机间的空间相对位置关系，照相机的性质以及三维场景点的位置，可以得到这些图像点的一些集合关关系约束。我们通过外极集合来描述这些几何关系。没有关于照相机的先验知识，会出现固有二义性，因为三维场景点X经过4x4的单应行矩阵H变换为HX后，则HX再想想及PH-1 里得到的图像点和X在照相机P里得到的图像点相同。利用照相机方程，可以将上述问题描述为：

因此，当我们分析双视图几何关系时，总是可以将照相机间的对应位置关系用单应性矩阵加以简化。这里的单应性矩阵通常是只做刚性变换，即只通过单应性矩阵变换了坐标系。一个比较好的做法是，将原点和坐标轴与第一个照相机对齐，则：

其中k1 和k2是标定矩阵，R是第二个照相机的旋转矩阵，t是第二个照相机的平移向量。利用这些照相机参数矩阵，我们可以找到X的投影点x1和x2的关系。这样，我们可以从寻找对应的图像出发，恢复照相机参数矩阵。同一个图像点经过不同的投影矩阵产生的不同投影点必须满足：

矩阵St为反对称矩阵：

基础矩阵可以由两个照相机的参数矩阵表示。由于det(F) = 0,所以基础矩阵的秩小于等于2.上面的公式表明，我们可以借助F来恢复出照相机参数，而F可以从对应的投影图像点计算出来。在不知道照相机内参数的情况下，仅能恢复照相机的投影变换矩阵。在知道照相机内参数的情况下，重建是基于度量的。所谓度量重建，即能够在三维重建中正确的表示距离和角度。利用上面的理论处理一些图像数据前，我们还需要了解一些几何学只是。给定图像中的一个点，例如第二个视图中的图像点x2，我们找到对应第一幅图像的一条直线：

其中

是第一幅图像中的一条直线。这条线称为对应于点x2的外极线，也就是说，点x2在第一副图像中的对应点一定在这条线上。所以基础矩阵可以将对应点的搜索限制在外极线上。

外极线都经过一点e，称为外极点。实际上，外极点是另一个照相机光心对应的图像点。外极点可以在我们看到的图像外，这取决于照相机间的相对方向。因为外极点在所有的外极线上，所以Fe1 = 0. 因此，我们可以通过计算F的零向量来计算出外极点。同理，另一个外极点可以通过计算eT2F = 0 得到。

5.1.1 一个简单的数据集

在接下来的内容中，我们需要一个带有图像点，三维点，和照相机参数矩阵的数据集。我们这里使用一个牛津多视图数据集；从http://www.robots.ox.ac.uk/~vgg/data/ 可以下载Merton1数据的压缩文件。下面的脚本可以加载Merton1的数据：

import 

上面的程序回家再前两个图像，三个视图中的所有特征点，对应不同视图图像重建后的三维点以及照相机参数矩阵。这里我们使用loadtxt()函数读取文本文件到Numpy数组中。因为并不是所有的点都可见，或都能够成功匹配到所有的视图，所以对应数据包含了确实的数据。加载对应数据时需要考虑到这一点。genfromtxt()函数通过将缺失的数值填充为-1 来解决这个问题。下面来可视化这些数据。将三维的点投影到一个视图，然后和观测道德图像点比较：

# make 3D points homogeneous and project
X = vstack( (points3D,ones(points3D.shape[1])) )
x = P[0].project(X)
# plotting the points in view 1
figure()
imshow(im1)
plot(points2D[0][0],points2D[0][1],'*')
axis('off')
figure()
imshow(im1)
plot(x[0],x[1],'r.')
axis('off')
show()

5.1.2 用matplotlib绘制三维数据

为了可视化三维重建结果，我们需要绘制出三维图像。matplotlib中的mplot3d 工具包可以方便的绘制出三维点，线，等轮廓线，表面以及其他基本图形组件，还可以通过图像窗口控件实现三维旋转和缩放。

可以通过在axes对象中的project = '3d' 关键字实现三维视图，如下：

from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection="3d")
# generate 3D sample data
X,Y,Z = axes3d.get_test_data(0.25)
# plot the points in 3D
ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'o')
show()

get_test_data()函数在x,y空间按照设定的空间间隔参数来产生均匀的采样点。压平这些网络，会产生散列数据点，然后我们可以将其输入plot()函数。这样我们就可以在立体表面上画出三维点。

现在通过画出merton样本数据来观察三维点的效果：

# plotting 3D points
from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection='3d')
ax.plot(points3D[0],points3D[1],points3D[2],'k.')

下图是三个不同视图中的三维图像点。图像窗口和控制界面外观效果像加上三维旋转工具的标准画图窗口。

5.1.3 计算F：八点法

八点法是通过对应点来计算基础矩阵的算法。

外极约束可以写成线性系统的形式：

其中f包含F的元素，基础矩阵中有9个元素，由于尺度是任意的，所以只需要8个方程。因为算法中需要8个对应点来计算基础矩阵F，所以该算法叫做8点法。新建一个文件sfm.py,写入下面八点法中最小化||Af|| 的函数：

def compute_fundamental(x1,x2):
""" Computes the fundamental matrix from corresponding points
(x1,x2 3*n arrays) using the normalized 8 point algorithm.
each row is constructed as
[x'*x, x'*y, x', y'*x, y'*y, y', x, y, 1] """
n = x1.shape[1]
if x2.shape[1] != n:
raise ValueError("Number of points don't match.")
# build matrix for equations
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] ]
# compute linear least square solution
U,S,V = linalg.svd(A)
F = V[-1].reshape(3,3)
104 | Chapter 5: Multiple View Geometry
# constrain F
# make rank 2 by zeroing out last singular value
U,S,V = linalg.svd(F)
S[2] = 0
F = dot(U,dot(diag(S),V))
return F

我们通常用SVD算法来计算最小二乘解。由于上面算法得出的解不为2（基础矩阵的秩小于2等于2），所以我们通过将最后一个奇异值置0来得到秩最接近2的基础矩阵。这是个很有用的技巧。上面的函数忽略了一个重要的步骤：对图像坐标进行归一化，这可能会带来数值问题。

5.1.4 外极点和外极线

由于外极点Fe1 = 0 ， 因此可以通过计算F的零空间来得到。把下面的函数添加到sfm.py中：

def compute_epipole(F):
""" Computes the (right) epipole from a
fundamental matrix F.
(Use with F.T for left epipole.) """
# return null space of F (Fx=0)
U,S,V = linalg.svd(F)
e = V[-1]
return e/e[2]

如果想获得另一幅图像的外极点，只需将F转置后输入上述函数即可。我们可以在之前的样本数据集的前两个视图上运行这两个函数：

import sfm
# index for points in first two views
ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
# get coordinates and make homogeneous
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])) )
# compute F
F = sfm.compute_fundamental(x1,x2)
# compute the epipole
e = sfm.compute_epipole(F)
# plotting
figure()
imshow(im1)
# plot each line individually, this gives nice colors
for i in range(5):
sfm.plot_epipolar_line(im1,F,x2[:,i],e,False)
axis('off')
figure()
imshow(im2)
# plot each point individually, this gives same colors as the lines
for i in range(5):
plot(x2[0,i],x2[1,i],'o')
axis('off')
show()

上面的函数将x轴的范围作为直线的参数，因此直线超出图像边界的部分会被截断。如果show_epipole为真，外极点也会被画出来(如果输入参数中没有外极点，外极点会在程序中计算获得)。程序结果如下图所示。在两幅图中，用不同的颜色将点和相应的外极线对应起来。

5.2 照相机和三维结构的计算

上一节讲述了视图和基础矩阵，外极线计算方法的关系。这一节我们将简单的介绍计算照相机参数和三维结构的工具。

5.2.1 三角剖分

给定照相机参数模型，图像点可以通过三角剖分来恢复这些点的三维位置。基本的算法思想如下。

对于两个照相机P1 和P2 的视图，三维实物点X的投影点维x1和x2，照相机方程定义了下列关系：

由于图像噪声，照相机参数误差和其他系统误差，上面的方程可能没有精确解。我们可以通过SVD算法来得到三维点的最小二乘估值。

下面的函数用于计算一个点对的最小二乘三角剖分，把它添加到sfm.py中：

def triangulate_point(x1,x2,P1,P2):
""" Point pair triangulation from
least squares solution. """
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]

最后一个特征向量的前4个值就是齐次坐标系下的对应三维坐标。我们可以增加下面的函数来实现多个点的三角剖分：

def triangulate(x1,x2,P1,P2):
""" Two-view triangulation of points in
x1,x2 (3*n homog. coordinates). """
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
# index for points in first two views
ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
# get coordinates and make homogeneous
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])) )
# check first 3 points
Xest = sfm.triangulate(x1,x2,P[0].P,P[1].P)
print Xest[:,:3]
print Xtrue[:,:3]
# plotting
from mpl_toolkits.mplot3d import axes3d
fig = 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()

上面的代码首先利用前两个视图的信息来对图像点进行三角剖分，然后把前三个图像点的齐次坐标输出到控制台，最后绘制出恢复的最接近三维图像点。输出到控制台的信息如下：

算法估计出的三维图像点和实际图像点很接近。如下图所示，估计点和实际点可以很好的匹配。

5.2.2 由三维点计算照相机矩阵

如果已经知道了一些三维点及其图像投影， 我们可以使用直接线性变换的方法来计算照相机矩阵P。本质上，这是三角剖分方法的逆问题，有时我们将其称为照相机反切法。利用该方法恢复照相机矩阵同样也是一个最小二乘问题。

我们从照相机方程可以得出，每个三维点Xi按照λixi = PXi，投影到图像点xi = [xi , yi, 1]，相应的点满足下面的关系：

其中P1，P2 和P3是矩阵P的三行。上面的式子可以改写得更紧凑，如下所示：Mv = 0.

然后，我们可以使用SVD分解估计得出照相机矩阵。利用上面讲述的矩阵操作，我们可以直接写出相应的代码。将下面的函数添加到sfm.py文件中：

def compute_P(x,X):
""" Compute camera matrix from pairs of
2D-3D correspondences (in homog. coordinates). """
n = x.shape[1]
if X.shape[1] != n:
raise ValueError("Number of points don't match.")
# create matrix for DLT solution
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个元素是照相机矩阵的元素，经过重新排列成矩阵形状后返回。

下面，在我们的样本数据集上测试算法的性能。下面的代码会选出第一个视图中的一些可见点，将它们转化为齐次坐标表示，然后估计照相机矩阵：

import sfm, camera
corr = corr[:,0] # view 1
ndx3D = where(corr>=0)[0] # missing values are -1
ndx2D = corr[ndx3D]

# select visible points and make homogeneous
x = points2D[0][:,ndx2D] # view 1
x = vstack( (x,ones(x.shape[1])) )
X = points3D[:,ndx3D]
X = vstack( (X,ones(X.shape[1])) )
# estimate P
Pest = camera.Camera(sfm.compute_P(x,X))
# compare!
print Pest.P / Pest.P[2,3]
print P[0].P / P[0].P[2,3]
xest = Pest.project(X)
# plotting
figure()
imshow(im1)
plot(x[0],x[1],'bo')
plot(xest[0],xest[1],'r.')
axis('off')
show()

为了 检查照相机矩阵的正确性，将他们以归一化的格式打印到控制台。输出如下所示：

上面是估计出的照相机矩阵，下面是数据集到创建者计算出的照相机矩阵。可以看到他们的元素几乎完全相同。组后使用估计出的照相机矩阵投影这些三维点，然后绘制出投影后的结果。真实点用圆圈表示，估计出的照相机投影点用点表示。

5.2.3 由基础矩阵计算照相机矩阵

在两个视图的场景中，照相机矩阵可以由基础矩阵恢复出来。假设第一个照相机矩阵归一化为P1 = [I/0] ,现在我们需要计算出第二个照相机矩阵P2。 研究分为两类，未标定的情况和已经标定的情况。

1. 未标定的情况 -- 投影重建

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变化恢复出来。也就是说如果利用照相机的信息来重建三维点，那么该重建只能由射影变换来计算出来。在这里，我们不考虑角度和距离。

因此，在无标定情况下，第二个照相机矩阵可以使用一个3x3 的射影变换得出。一个简单的方法是：P2 = [SeF | e],其中e是左极点，满足eT F =0。下面是具体的代码：

def 

这里，我们使用skew()函数定义如下：

def skew(a):
""" Skew matrix A such that a x v = Av for any v. """
return array([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])

2. 已标定的情况 --度量重建

在已经标定的情况下，重建会保持欧氏空间中的一些度量特性。在重建三维场景中，已经标定的例子更为有趣。

给定标定矩阵K，我们可以将它的逆作用于图像点，因此，在新的图像坐标下，照相机方程变为：

在新的图像坐标系下，点同样满足之前的基础矩阵方程：

在标定归一化的坐标系下，基础矩阵称为本质矩阵。为了区别为标定后的情况，以及归一化了的图像坐标，我们通常将其标记为E，而非F。

从本质矩阵恢复出的照相机矩阵中存在度量关系，但有四个可能解。因为只有一个解产生位于两个照相机前的场景，所以我么可以轻松地从中选出来。

下面是计算这4个解地算法。

def compute_P_from_essential(E):
""" Computes the second camera matrix (assuming P1 = [I 0])
from an essential matrix. Output is a list of four
possible camera matrices. """
# make sure E is rank 2
U,S,V = svd(E)
if det(dot(U,V))<0:
V = -V
E = dot(U,dot(diag([1,1,0]),V))
# create matrices (Hartley p 258)
Z = skew([0,0,-1])
W = array([[0,-1,0],[1,0,0],[0,0,1]])
# return all four solutions
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.然后，计算出这4个解。

5.3 多视图重建

由于照相机的运动给我们提供了三维结构，所以这样计算三维重建的方法通常称为SFM，

Structure from Motion, 从运动中恢复 结构。

假设照相机已经标定，计算重建可以分为下面4个步骤：

1.检测特征点，然后在两幅图像间匹配

2.由匹配计算基础矩阵

3.由基础矩阵计算照相机矩阵

4.三角剖分这些三维点。

5.3.1 文件估计基础矩阵

类似于稳健计算单应行矩阵，当存在噪声和不正确的匹配时，我们需要估计基础矩阵。和前面的方法一样，我们使用RANSAC 方法，只不过这次结合了八点算法。

class 

和之前一样，我们需要fit() 和 get_error() 方法。这里采用的错误衡量方法是SAMPSON距离。fit() 方法会选择8个点，然后使用归一化的八点算法：

def compute_fundamental_normalized(x1,x2):
""" Computes the fundamental matrix from corresponding points
(x1,x2 3*n arrays) using the normalized 8 point algorithm. """
n = x1.shape[1]
if x2.shape[1] != n:
raise ValueError("Number of points don't match.")
# normalize image coordinates
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)
# compute F with the normalized coordinates
F = compute_fundamental(x1,x2)
# reverse normalization
F = dot(T1.T,dot(F,T2))
return F/F[2,2]

这里，我们返回最佳矩阵F，以及正确点的索引，这样可以知道哪些匹配和F矩阵一致。与单应性矩阵估计相比，我们增加了默认的最大迭代次数，改变了匹配的阙值。

5.3.2 三维重建示例

接下来，我们将观察一个重建三维场景的完整例子。我们使用由已知标定矩阵的照相机拍摄两幅图像。该图像是著名的恶魔岛，如下图所示：

我们将该处理代码分成若干块，使得代码更容易理解。首先提取，匹配特征，然后估计基础矩阵和照相机矩阵：

import homography
import sfm
import sift
# calibration
K = array([[2394,0,932],[0,2398,628],[0,0,1]])
# load images and compute features
im1 = array(Image.open('alcatraz1.jpg'))
sift.process_image('alcatraz1.jpg','im1.sift')
l1,d1 = sift.read_features_from_file('im1.sift')
im2 = array(Image.open('alcatraz2.jpg'))
sift.process_image('alcatraz2.jpg','im2.sift')
l2,d2 = sift.read_features_from_file('im2.sift')
# match features
matches = sift.match_twosided(d1,d2)
ndx = matches.nonzero()[0]
# make homogeneous and normalize with 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)
# estimate E with RANSAC
model = sfm.RansacModel()
E,inliers = sfm.F_from_ransac(x1n,x2n,model)
# compute camera matrices (P2 will be list of four solutions)
P1 = array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])
P2 = sfm.compute_P_from_essential(E)

现在，我们已经获得了标定矩阵，所以只对矩阵K进行硬编码。与之前的例子一样，我们挑选属于匹配关系的特定点。然后使用K-1 来对其归一化，并使用归一化的八点算法来运行RANSAC估计。因为该点已经归一化，所以这里会返回一个本质矩阵。从本质矩阵出发，我们可以计算出第二个照相机矩阵的四个可能解。

从照相机矩阵的列表中，挑选出经过三角剖分后，在两个照相机前均含有最多场景点的照相机矩阵：

# pick the solution with points in front of cameras
ind = 0
maxres = 0
for i in range(4):
# triangulate inliers and compute depth for each camera
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)
# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[ind])
X = X[:,infront]

循环遍历这4个解，每次对于正确点的三维点进行三角剖分。将三角剖分后的图像投影回图像后，深度的符号由每个图像点的第三个数值给出。我们保存了正向最大深度的索引；对于和最优解一致的点，用相应的布尔值保存了信息，这样可以取出真正在照相机前面的点。因为所有估计中都存在噪声和误差，所以即便使用正确的照相机矩阵，也存在一些仍位于某个照相机后面的风险。首先获得正确的解，然后对这些正确的点进行三角剖分，然后保留位于照相机前方的点。现在可以绘制出该三维重建：

# 3D plot
from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection='3d')
ax.plot(-X[0],X[1],X[2],'k.')
axis('off')

和我们的坐标系相比，使用mplot3绘制出三维图像需要将第一个坐标值取相反数，所以这里改变其符号。

然后，可以在每个视图中绘制出二次投影：

# plot the projection of X
import camera
# project 3D points
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2[ind])
x1p = cam1.project(X)
x2p = cam2.project(X)
# reverse K normalization
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()

将这些三维点投影后，可以通过乘以标定矩阵的方式来弥补初始归一化对点的影响。结果输出如下图所示，可以看出，二次投影后的和原始特征位置不完全匹配，但是相当接近。当然，可以进一步调整照相机矩阵来提高重建和二次投影的性能。

5.4 立体图像

一个多视图成像的特殊例子是立体视觉，或者立体成像，即使用两台只有水平偏移的照相机观测同一场景。当照相机的位置如上设置，两幅图像具有相同的图像平面，图像的行为是垂直对齐的，那么称图像对是经过矫正的。该设置在机器人学中很常见，常被称为，立体平台。

通过将图像扭曲到公共平面上，使外极线位于图像行上，任何立体照相机设置都能得到矫正。结社两幅图像经过了矫正，那么对应点的寻找限制在图像的同一行上。一旦找到对应点，由于深度是和偏移成正比的，那么深度可以直接由水平偏移来计算，

其中，f是经过矫正图像的焦距，b是两个照相机中心之间的距离，x1和x2是左右两幅图像中对应点的x坐标。分开照相机中心的距离称为基线。矫正后的立体照相机设置如下图所示。

立体重建，有时候称为致密深度重建，就是恢复深度图，图像中每个像素的深度，都需要计算出来。这是计算机视觉中的经典问题。

计算视差图

在该立体重建算法中，我们将对于每个图像尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面，所以该过程有时称为扫平面法。虽然该方法并不是立体重建中最好的方法，但是非常简单，通常会得出令人满意的结果。

当密集地应用在图像中时，归一化的互相关值可以很快的计算出来。这和我们在第2章中应用于稀疏点对应的不同。我们使用每个像素周围的图像块来计算归一化的互相关。对于这里的情形，我们可以在像素周围重新写出公式中的NCC，如下所示：

我们在前面跳过了归一化常数，然后对该像素周围局部像素块中的像素求和。

现在，我们要将图像中的每个像素进行该操作。这三个求和操作是在局部图像块区域上进行的，我们可以使用图像滤波器来快速计算，这与我们在模糊和求导操作中使用的技巧一样。ndimage.filters 模块中uniform_filter()函数可以在一个矩形图像块中计算相加。

下面是扫平面法的具体实现代码，该函数返回每个像素的最佳视差。

def plane_sweep_ncc(im_l,im_r,start,steps,wid):
""" Find disparity image using normalized cross-correlation. """
m,n = im_l.shape
# arrays to hold the different sums
mean_l = zeros((m,n))
mean_r = zeros((m,n))
s = zeros((m,n))
s_l = zeros((m,n))
s_r = zeros((m,n))
# array to hold depth planes
dmaps = zeros((m,n,steps))
# compute mean of patch
filters.uniform_filter(im_l,wid,mean_l)
filters.uniform_filter(im_r,wid,mean_r)
# normalized images
norm_l = im_l - mean_l
norm_r = im_r - mean_r
# try different disparities
for displ in range(steps):
# move left image to the right, compute sums
filters.uniform_filter(roll(norm_l,-displ-start)*norm_r,wid,s) # sum nominator
filters.uniform_filter(roll(norm_l,-displ-start)*roll(norm_l,-displ-start),wid,
s_l)
filters.uniform_filter(norm_r*norm_r,wid,s_r) # sum denominator
# store ncc scores
dmaps[:,:,displ] = s/sqrt(s_l*s_r)
# pick best depth for each pixel
return argmax(dmaps,axis=2)

首先，因为uniform_filter()函数的输入参数为数组，我们需要创建用于保存滤波结果的一些数组。然后，我们创建一个数组来保存每个平面，从而能够在最后一个维度上使用argmax()函数找到对于每个像素的最佳深度。该函数从start偏移出发，在所有的steps偏移上迭代。使用roll()函数平移一副图像，然后使用滤波计算NCC的三个求和操作。下面是载入图像，并使用该函数计算偏移图的完整例子：

import stereo
im_l = array(Image.open('scene1.row3.col3.ppm').convert('L'),'f')
im_r = array(Image.open('scene1.row3.col4.ppm').convert('L'),'f')
# starting displacement and steps
steps = 12
start = 4
# width for ncc
wid = 9
res = stereo.plane_sweep_ncc(im_l,im_r,start,steps,wid)
import scipy.misc
scipy.misc.imsave('depth.png',res)

这里首先从经典的“tsukuba”数据集中载入一对图像，然后将其灰度化。接下来，设置扫平面函数所需的参数，包括尝试偏移的数目，初始值和NCC路径的宽度。你会发现，该方法非常快，至少快于使用NCC进行特征匹配。这也是使用滤波器来计算的原因。该方法同样使用于其他滤波器。均匀滤波器给正方形中所有像素相同的权值。但是，在一些例子中，其他滤波器对于NCC的计算可能更为适用。

def plane_sweep_gauss(im_l,im_r,start,steps,wid):
""" Find disparity image using normalized cross-correlation
with Gaussian weighted neigborhoods. """
m,n = im_l.shape
# arrays to hold the different sums
mean_l = zeros((m,n))
mean_r = zeros((m,n))
s = zeros((m,n))
s_l = zeros((m,n))
s_r = zeros((m,n))
# array to hold depth planes
dmaps = zeros((m,n,steps))
# compute mean
filters.gaussian_filter(im_l,wid,0,mean_l)
filters.gaussian_filter(im_r,wid,0,mean_r)
# normalized images
norm_l = im_l - mean_l
norm_r = im_r - mean_r
# try different disparities
for displ in range(steps):
# move left image to the right, compute sums
filters.gaussian_filter(roll(norm_l,-displ-start)*norm_r,wid,0,s) # sum nominator
filters.gaussian_filter(roll(norm_l,-displ-start)*roll(norm_l,-displ-start),wid,
0,s_l)
filters.gaussian_filter(norm_r*norm_r,wid,0,s_r) # sum denominator
# store ncc scores
dmaps[:,:,displ] = s/sqrt(s_l*s_r)
# pick best depth for each pixel
return argmax(dmaps,axis=2)

除了在滤波器中使用了额外的参数，该代码和均匀滤波的代码相同。我们需要在gaussian_filter()函数中传入零参数来表示我们使用的是标准高斯函数，而不是其他任何导数。

你可以像前面的扫平面函数一样来使用该函数。以下两个图像所示，为两个扫平面实现操作在一些标准立体基准图像上的结果。这里我们使用 tsukuba 和cones图像，在标准版本中设置wid为9，高斯版本中设置wid为3.上面一行图像所示为图像对，坐下图像是标准的NCC扫平面法重建的视差图。右下方是高斯版本的视差图。可以看到，于标准版本相比，高斯版本具有较少的噪声，但缺少很多细节信息。