前言
立体匹配主要是通过找出每对图像间的对应关系,根据三角测量原理,得到视差图;在获得了视差信息后,根据投影模型很容易地可以得到原始图像的深度信息和三维信息。立体视觉匹配在立体视觉研究中是比较核心的问题。立体视觉应用:车导航,3D场景重建等。本文主要介绍窗口匹配方法及NCC视差匹配。
一、立体视觉原理
1.1 几何原理
一个多视图成像的特殊例子是主体视觉(或者立体成像),即使用两台只有水平(向一侧) 偏移的照相机观测同一场景。当照相机的位置如上设置,两幅图像具有相同的图像平面,图像的行是垂直对齐的,那么称图像对是经过矫正的。
假设两幅图像经过了矫正,那么对应点的寻找限制在图像的同一行上。一旦找到对应点,由于深度是和偏移成正比的,那么深度(Z坐标)可以直接由水平偏移来计算。那么Z该如何计算呢?先看下图:
通过几何知识运算,得到
T
+
x
l
−
x
r
T
=
Z
−
f
Z
\frac{T+x_{l}-x_{r}}{T}=\frac{Z-f}{Z}
TT+xl−xr=ZZ−f进一步化简,可得
Z
=
f
T
x
r
−
x
l
Z=f\frac{T}{x_{r}-x_{l}}
Z=fxr−xlT其中,f是经过矫正图像的焦距,T是两个照相机中心之间的距离,
x
l
和
x
r
x_{l}和x_{r}
xl和xr是左右两幅图像中对应点的x坐标。分开照相机中心的距离称为基线。
1.2 立体匹配的步骤
1)匹配代价计算:
计算匹配代价,即计算参考图像上每个像素点 I R ( P ) IR(P) IR(P) ,以所有视差可能性去匹配目标图像上对应点 I T ( p d ) IT(pd) IT(pd)的代价值,因此计算得到的代价值可以存储在一个 h ∗ w ∗ d ( M A X ) h*w*d(MAX) h∗w∗d(MAX)的三维数组中,通常称这个三维数组为视差空间图。匹配代价时立体匹配的基础,设计抗噪声干扰、对光照变化不敏感的匹配代价,能提高立体匹配的精度。因此,匹配代价的设计在全局算法和局部算法中都是研究的重点。
2)代价聚合:
通常全局算法不需要代价聚合,而局部算法需要通过求和、求均值或其他方法对一个支持窗口内的匹配代价进行聚合而得到参考图像上一点p在视差d处的累积代价 C A ( p , d ) CA(p,d) CA(p,d),这一过程称为代价聚合。通过匹配代价聚合,可以降低异常点的影响,提高信噪比进而提高匹配精度。代价聚合策略通常是局部匹配算法的核心,策略的好坏直接关系到最终视差图(Disparity maps)的质量。
3)视差计算:
局部立体匹配算法的思想,在支持窗口内聚合完匹配代价后,获取视差的过程通常采用‘胜者为王’策略(WTA,Winner Take All),即在视差搜索范围内选择累积代价最优的点作为对应匹配点,与之对应的视差即为所求的视差。即P点的视差为 d = a r g m i n C A ( p , d ) d=arg minCA(p,d) d=argminCA(p,d)但是这种方法可能导致噪声大,因为如果直接用点的像素值坐标匹配,风险太大,找到的匹配点不一定是正确的,这时可以采用窗口匹配方法,关于这种方法在下面会做详细介绍。
4)后处理:
一般的,分别以左右两图为参考图像,完成上述三个步骤后可以得到左右两幅视差图像。但所得的视差图还存在一些问题,如遮挡点视差不准确、噪声点、误匹配点等存在,因此还需要对视差图进行优化,采用进一步执行后处理步骤对视差图进行修正。常用的方法有插值(Interpolation)、亚像素增强(Subpixel Enhancement)、精细化(Refinement)、图像滤波(Image Filtering)等操作。
1.3 窗口匹配方法
在上面提到过,如果使用点匹配方法的话,可能导致噪声大,因为直接找最匹配的点误差太大,如下图,比如我们要在(T)图中找(R)图标示的眼睛位置匹配点,由于两边眼睛的像素值相近,很可能找错,造成误差。
这时,我们就可以考虑改进的方法:窗口匹配方法。它的主要思路是:以当前点为中心,切出一个图像块,每个图像块依次做对比,找到最接近的图像块。这时就不只是拿单个点做对比了,可以有效减小误差。举例见下图:
匹配代价
相似性测度函数用于度量参考图像中的匹配基元和目标图像中的匹配基元的相似性,即判断参考图像和目标图像中两点为对应匹配点的可能性,也称为匹配代价,这里记为
C
(
x
,
y
,
d
)
C(x,y,d)
C(x,y,d),表示像素点
(
x
,
y
)
(x,y)
(x,y)在视差为d情况下的匹配误差。
下面介绍最常见的三种匹配代价:
(1)绝对差值和:(Sum of Absolute Differences,SAD)
C
(
x
,
y
,
d
)
=
∑
x
ε
S
∣
I
R
(
x
,
y
)
−
I
T
(
x
+
d
,
y
)
∣
C(x,y,d)=\sum_{x\varepsilon S}^{}\left | I_{R}(x,y)-I_{T}(x+d,y)\right |
C(x,y,d)=xεS∑∣IR(x,y)−IT(x+d,y)∣(2)差值平方和(Sum of squared Differences,SSD)
C
(
x
,
y
,
d
)
=
∑
x
ε
S
(
I
R
(
x
,
y
)
−
I
T
(
x
+
d
,
y
)
)
2
C(x,y,d)=\sum_{x\varepsilon S}^{}(I_{R}(x,y)-I_{T}(x+d,y))^{2}
C(x,y,d)=xεS∑(IR(x,y)−IT(x+d,y))2(3)截断绝对差值和(Sum of Truncated Absolute Differences,STAD)
C
(
x
,
y
,
d
)
=
∑
x
ε
S
m
i
n
{
∣
I
R
(
x
,
y
)
−
I
T
(
x
+
d
,
y
)
∣
,
T
}
C(x,y,d)=\sum_{x\varepsilon S}^{}min\left \{\left | I_{R}(x,y)-I_{T}(x+d,y)\right |,T\right \}
C(x,y,d)=xεS∑min{∣IR(x,y)−IT(x+d,y)∣,T}另外常用的匹配代价还有归一化互相关NCC(Normalized Cross Correlation),在下面做详细介绍。
1.4 归一化互相关(NCC)
原理:
对于原始的图像内任意一个像素点 ( p x , p y ) (p_{x},p_{y}) (px,py)构建一个n x n的邻域作为匹配窗口。然后对于目标相素位置 ( p x + d , p y ) (p_{x}+d,p_{y}) (px+d,py)同样构建一个n x n大小的匹配窗口,对两个窗口进行相似度度量,注意这里的d有一个取值范围。对于两幅图像来说,在进行NCC计算之前要对图像处理,也就是将两帧图像校正到水平位置,即光心处于同一水平线上,此时极线是水平的,否则匹配过程只能在倾斜的极线方向上完成,这将消耗更多的计算资源。
NCC计算公式:
其中
N
C
C
(
p
,
d
)
NCC(p,d)
NCC(p,d)得到的值得范围将在[−1,1]之间。
- W p W_{p} Wp为之前提到的匹配窗口。
- I 1 ( x , y ) I_{1}(x,y) I1(x,y)为原始图像的像素值。
- I 1 ˉ ( p x , p y ) \bar{I_{1}}(p_{x},p_{y}) I1ˉ(px,py)为原始窗口内像素的均值。
- I 2 ( x + d , y ) I_{2}(x+d,y) I2(x+d,y)为原始图像在目标图像上对应点位置在x方向上偏移d后的像素值。
- I 2 ˉ ( p x + d , p y ) \bar{I_{2}}(p_{x}+d,p_{y}) I2ˉ(px+d,py)为目标图像匹配窗口像素均值。
若NCC=−1,则表示两个匹配窗口完全不相关,相反,若NCC=1时,表示两个匹配窗口相关程度非常高。
匹配流程:
-
采集图像:通过标定好的双目相机采集图像,当然也可以用两个单目相机来组合成双目相机。
-
极线校正:校正的目的是使两帧图像极线处于水平方向,或者说是使两帧图像的光心处于同一水平线上。通过校正极线可以方便后续的NCC操作。
由标定得到的内参中畸变信息中可以对图像去除畸变。
通过校正函数校正以后得到相机的矫正变换R和新的投影矩阵P,接下来是要对左右视图进行去畸变,并得到重映射矩阵。 -
特征匹配:这里便是我们利用NCC做匹配的步骤啦,匹配方法如上所述,右视图中与左视图待测像素同一水平线上相关性最高的即为最优匹配。完成匹配后,我们需要记录其视差d,即待测像素水平方向xl与匹配像素水平方向xr之间的差值 d = x r − x l d=x_{r}-x_{l} d=xr−xl,最终我们可以得到一个与原始图像尺寸相同的视差图D。
-
深度恢复:通过上述匹配结果得到的视差图D,我们可以很简单的利用相似三角形反推出以左视图为参考系的深度图。计算原理同上述1.1中介绍的几何原理(这里不做重述),即通过公式: Z = f T x r − x l Z=f\frac{T}{x_{r}-x_{l}} Z=fxr−xlT可以简单得到以左视图为参考系的深度图了。
二、实验内容
2.1 实验目的与要求
- 实现NCC 视差匹配方法,即给定左右两张视图,根据NCC计算视差图
- 分析不同窗口值对匹配结果的影响,重点考查那些点(或者哪些类型的点)在不同窗口大小下的匹配精度影响
2.2 代码实现
# -*- coding: utf-8
from pylab import *
from numpy import *
np.seterr(divide='ignore', invalid='ignore')
from PIL import Image
import numpy as np
import math
from scipy import ndimage
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))
# 计算图像块的平均值
ndimage.filters.uniform_filter(im_l, wid, mean_l)
ndimage.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):
# 将左边图像移动到右边,计算加和
ndimage.filters.uniform_filter(np.roll(norm_l, -displ - start) * norm_r, wid, s) # 和归一化
ndimage.filters.uniform_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, s_l)
ndimage.filters.uniform_filter(norm_r * norm_r, wid, s_r) # 和反归一化
# 保存 ncc 的分数
dmaps[:, :, displ] = s / np.sqrt(s_l * s_r)
# 为每个像素选取最佳深度
return np.argmax(dmaps, axis=2)
def plane_sweep_gauss(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))
# 计算平均值
ndimage.filters.gaussian_filter(im_l, wid, 0, mean_l)
ndimage.filters.gaussian_filter(im_r, wid, 0, mean_r)
# 归一化图像
norm_l = im_l - mean_l
norm_r = im_r - mean_r
# 尝试不同的视差
for displ in range(steps):
# 将左边图像移动到右边,计算加和
ndimage.filters.gaussian_filter(np.roll(norm_l, -displ - start) * norm_r, wid, 0, s) # 和归一化
ndimage.filters.gaussian_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, 0, s_l)
ndimage.filters.gaussian_filter(norm_r * norm_r, wid, 0, s_r) # 和反归一化
# 保存 ncc 的分数
dmaps[:, :, displ] = s / np.sqrt(s_l * s_r)
# 为每个像素选取最佳深度
return np.argmax(dmaps, axis=2)
"""载入图像,并使用该函数计算偏移图"""
im_l = np.array(Image.open('../data/data_six/p1.png').convert('L'), 'f')
im_r = np.array(Image.open('../data/data_six/p2.png').convert('L'), 'f')
figure()
gray()
subplot(221)
imshow(im_l)
subplot(222)
imshow(im_r)
# 开始偏移,并设置步长
steps = 12
start = 4
# ncc 的宽度
wid = 13
res = plane_sweep_ncc(im_l, im_r, start, steps, wid)
wid2 = 5
gauss = plane_sweep_gauss(im_l, im_r, start, steps, wid2)
subplot(223)
imshow(res)
subplot(224)
imshow(gauss)
show()
import scipy.misc
scipy.misc.imsave('depth.png', res)
2.3 实验结果与分析
图片来源:
下载链接为: 视差匹配实验图片集.
NCC版本与高斯版本的视差图:
分析:如上图所示,在标准版本中设置wid=9,高斯版本中设置wid=3。上面一行图像所示为图像对,左下方图像是标准的NCC扫平面法重建的视差图,右下方是高斯版本的视差图。可以看到,与标准版本相比,高斯版本具有较少的噪声,但缺少很多细节信息。
NCC视差匹配及不同窗口值对匹配结果的影响:
为了实现一次性输出不同wid值下的视差图,修改代码:
wid = np.array([1,3,5,7,9,11])
res=[]
for i in range(0, 6):
res.append(plane_sweep_ncc(im_l, im_r, start, steps, wid[i]))
a = np.array([231,232,233,234,235,236])
color_d = ['ncc wid=1', 'ncc wid=3', 'ncc wid=5','ncc wid=7','ncc wid=9','ncc wid=11']
for i in range(0, 6):
subplot(a[i])
imshow(res[i])
title(color_d[i])
axis('off')
原图:
NCC视差匹配(不同窗口值):
分别为wid=1、3、5、7、9、11的视差图:
分析:
- 如上图所示,首先对单张NCC视差匹配的结果做分析,比如观察当wid=9时的结果(如下图):
可以看到离相机越近的物体亮度越高,比如靠近相机的笔,它们整体呈现的更亮,而后面的背景则更暗。这是因为前景的位移比背景的位移更多,越靠近摄像机的目标,它在左右视图位移的距离会更多,这是由近大远小导致。结合原理来分析,我们可以观察公式: Z = f T x r − x l Z=f\frac{T}{x_{r}-x_{l}} Z=fxr−xlT发现 Z 和 x r − x l Z和x_{r}-x_{l} Z和xr−xl呈反比关系,因此视差越大,Z越小,即目标越近,视差越大,深度越小。 - 观察不同窗口值(wid)的运行结果,初步观察可得知:随着wid值的增大,图像深度信息的图像的轮廓逐渐明显,可以看到离照相机最近的目标是笔,随着wid增大,笔的轮廓逐渐清晰,但较远处的物体,比如墙壁以及画布,则越来越难以辨认。
- 进一步观察结果,这次拿出wid=3与wid=11的结果作对比,如下图:
可以看到wid越小,能得到更多的细节,比如见wid=3时的结果图,虽然图像整体比较模糊,但是可以看到我们能够得到更多笔的细节信息,从原理分析,当wid值较小,进行滤波时考虑的因素少,只受某像素点自身影响以及周围少量像素点的影响,所以保留下来更多的自身特征。但同时噪声也大,见wid=3时的背景墙,上面出现了许多噪声点。
- 而当窗口值增大时,鲁棒性增强,噪声明显减少,但我们可以发现图像的细节也在逐渐减少,可物体的轮廓逐渐清晰,尤其是近景目标,如wid=11的视差图,我们能够更清晰地看出笔的轮廓,但是由于细节减少,远处的物体就比较不易判断。在wid值较大情况下,NCC匹配下充分考虑了周围像素点的影响,能更好的观察出视差。
三、总结
实验内容总结:
经过本次实验,主要实现了NCC视差匹配,发现越靠近摄像机的目标,视差越大,因此视差图中近景目标的图像更亮。对于不同wid值,wid越小,能得到更多的细节,但是噪声大;wid越大,具有更好的鲁棒性,但是得到的细节少。
实验过程总结:
①问题:在使用高斯滤波器做视差匹配时报错:
解决办法:改变sqrt()函数调用的库,即将原本的从math库调用改为从numpy库调用,修改代码:
dmaps[:, :, displ] = s / np.sqrt(s_l * s_r)
②问题:输出结果模糊不清,无法看到轮廓。
解决办法:调整图片的像素值,而且注意左右两张图拍摄时位移不能太大。