文章目录
一、以窗口代价计算视差的原理
1.概念
视差: 同一个场景在两个相机下成像的像素的位置偏差
深度: 指场景中每个点离相机的距离
当已知两个相机的相关参数,可由视差计算深度。
为什么需要双目(两个相机)才能得到视差图?
如图所示,从物理原理上,我们看到红色线条上三个不同远近的黑色的点在下方相机上投影在同一个位置,因此单目相机无法分辨成的像到底是远的那个点还是近的那个点,但是它们在上方相机的投影却位于三个不同位置,因此通过两个相机的观察可以确定到底是哪一个点。
2.如何得到深度图
2.1公式推导
图像俯视图如下:
立体视图如下:
如图,
X
X
X是空间中的一个3D点,投影在矫正后的左右视图上为
x
l
x_{l}
xl、
x
r
x_{r}
xr,两相机中心的连线为
T
T
T,相机焦距为
f
f
f,图像宽度
W
W
W,图像深度
Z
Z
Z,实际的坐标计算根据相似三角形的原理有得到如下公式
从公式可以看出,视差
x
l
−
x
r
x_{l}-x_{r}
xl−xr和深度
Z
Z
Z成反比关系。视差越大,可以探测的深度越小。
由上可知,要得到图像深度
Z
Z
Z,需要有以下参数
(1)相机焦距
f
f
f,左右相机基线
T
T
T。这些参数可以通过先验信息或者相机标定得到。
(2)视差 x l − x r x_{l}-x_{r} xl−xr。需要知道左相机的每个像素点( x l x_{l} xl, y l y_{l} yl )和右相机中对应点( x r x_{r} xr, y r y_{r} yr)的对应关系
2.2获得视差图流程
- 首先需要对双目相机进行标定,得到两个相机的内外参数、单应矩阵。
- 根据标定结果对原始图像校正,校正后的两张图像位于同一平面且互相平行。
- 对校正后的两张图像进行像素点匹配。
- 根据匹配结果计算每个像素的深度,从而获得视差图。
2.3对流程的解析
2.3.1极线约束——确定两图像中像素的位置
对于左图中的一个像素点,如何确定该点在右图中的位置?
极线约束:在一幅图像中上的 p 1 p_{1} p1点在另一幅图像中的对应点 p 2 p_{2} p2一定在基线 l ′ l' l′上
因为有极线约束,对于左图的一点
p
1
p_{1}
p1点,寻找它在右图中的对应点
p
2
p_{2}
p2,这样就能确定P点的空间位置,即空间物体和相机的距离(深度)。
2.3.2图像矫正
图像矫正是通过分别对两张图片用单应(homography)矩阵变换(可以通过标定获得)得到的,目的就是把两个不同方向的图像平面(下图中灰色平面)重新投影到同一个平面且光轴互相平行(下图中黄色平面),使得两个图像的极线变为水平(如图2)。
2.3.3两张图像进行像素点匹配——归一化相关(NCC)
NCC,就是用于归一化待匹配目标之间的相关程度,注意这里比较的是原始像素。通过在待匹配像素位置
p
(
p
x
,
p
y
)
p(p_{x},p_{y})
p(px,py)构建3*3邻域匹配窗口,与目标像素位置
p
′
(
p
x
+
d
,
p
y
)
p'(p_{x}+d,p_{y})
p′(px+d,py)同样构建邻域匹配窗口的方式建立目标函数来对匹配窗口进行度量相关性(前提图像极线已矫正为水平),公式如下:
公式变量解释:
-
p p p点表示图像I1待匹配像素坐标 p x − p y p_{x}-p_{y} px−py,d表示在图像I2被查询像素位置在水平方向上与 p x p_{x} px的距离。如下图所示:
-
左边为图像 I 1 I_{1} I1,右边为图像 I 2 I_{2} I2。图像 I 1 I_{1} I1,蓝色方框表示待匹配像素坐标 ( p x − p y ) (p_{x}-p_{y}) (px−py),图像 I 2 I_{2} I2蓝色方框表示坐标位置为 ( p x , p y ) (p_{x},p_{y}) (px,py),红色方框表示坐标位置 ( p x + d , p y ) (p_{x}+d,p_{y}) (px+d,py)。
-
W x W_{x} Wx表示以待匹配像素坐标为中心的匹配窗口,通常为 3 ∗ 3 3*3 3∗3匹配窗口。
-
没有上划线的 I 1 I_{1} I1表示匹配窗口中某个像素位置的像素值,带上划线的 I 1 I_{1} I1表示匹配窗口所有像素的均值。 I 2 I_{2} I2同理。
上述公式表示度量两个匹配窗口之间的相关性,通过归一化将匹配结果限制在 [-1,1]的范围内,可以非常方便得到判断匹配窗口相关程度:若NCC = -1,则表示两个匹配窗口完全不相关,相反,若NCC = 1时,表示两个匹配窗口相关程度非常高。
二、代码
stereo.py
功能: 对图像中的每个像素进行归一化相关(NCC)操作
实现: 因为 uniform_filter() 函数的输入参数为数组,我们需要创建用于保存滤波结果的一些数组。然后,我们创建一个数组来保存每个平面,从而能够在最后一个纬度上使用 argmax() 函数找到对于每个像素的最佳深度。该函数从 start 偏移出发,在所有的 steps 偏移上迭代。使用 roll() 函数平移一幅图像,然后使用滤波计算NCC 的三个求和操作。
import numpy
from scipy.ndimage import filters
#ncc算法
def plane_sweep_ncc(im_l, im_r, start, steps, wid):
'''Find disparity image using normalized cross-correlation.'''
m, n = im_l.shape # Must match im_r.shape.
mean_l = numpy.zeros(im_l.shape)
mean_r = numpy.zeros(im_l.shape)
s = numpy.zeros(im_l.shape)
s_l = numpy.zeros(im_l.shape)
s_r = numpy.zeros(im_l.shape)
dmaps = numpy.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 disp in range(steps):
filters.uniform_filter(numpy.roll(norm_l, -disp - start) * norm_r, wid, s)
filters.uniform_filter(numpy.roll(norm_l, -disp - start) *
numpy.roll(norm_l, -disp - start), wid, s_l)
filters.uniform_filter(norm_r * norm_r, wid, s_r)
dmaps[:, :, disp] = s / numpy.sqrt(s_l * s_r)
return numpy.argmax(dmaps, axis=2)
def plane_sweep_ssd(im_l, im_r, start, steps, wid):
'''Find disparity image using sum of squared differences.'''
m, n = im_l.shape # Must match im_r.shape.
s = numpy.zeros(im_l.shape)
dmaps = numpy.zeros((m, n, steps))
for disp in range(steps):
filters.uniform_filter((numpy.roll(im_l, -disp - start) - im_r) ** 2,
wid, s)
dmaps[:, :, disp] = s
return numpy.argmin(dmaps, axis=2)
#高斯滤波器
def plane_sweep_gauss(im_l, im_r, start, steps, wid):
'''Find disparity image using normalized cross-correlation with Gaussian
weighted neighborhoods.'''
m, n = im_l.shape # Must match im_r.shape.
mean_l = numpy.zeros(im_l.shape)
mean_r = numpy.zeros(im_l.shape)
s = numpy.zeros(im_l.shape)
s_l = numpy.zeros(im_l.shape)
s_r = numpy.zeros(im_l.shape)
dmaps = numpy.zeros((m, n, steps))
filters.gaussian_filter(im_l, wid, 0, mean_l)
filters.gaussian_filter(im_r, wid, 0, mean_r)
norm_l = im_l - mean_l
norm_r = im_r - mean_r
for disp in range(steps):
filters.gaussian_filter(numpy.roll(norm_l, -disp - start) *
norm_r, wid, 0, s)
filters.gaussian_filter(numpy.roll(norm_l, -disp - start) *
numpy.roll(norm_l, -disp - start), wid, 0, s_l)
filters.gaussian_filter(norm_r * norm_r, wid, 0, s_r)
dmaps[:, :, disp] = s / numpy.sqrt(s_l * s_r)
return numpy.argmax(dmaps, axis=2)
def plane_sweep_gauss_ssd(im_l, im_r, start, steps, wid):
'''Find disparity image using sum of squared differences with Gaussian
weighted neighborhoods.'''
m, n = im_l.shape # Must match im_r.shape.
s = numpy.zeros(im_l.shape)
dmaps = numpy.zeros((m, n, steps))
for disp in range(steps):
filters.gaussian_filter((numpy.roll(im_l, -disp - start) - im_r) ** 2,
wid, 0, s)
dmaps[:, :, disp] = s
return numpy.argmin(dmaps, axis=2)
主代码
实现: 首先从经典“tsukuba”数据集中载入一对图像,然后将其灰度化。接下来,设置扫平面函数所需的参数,包括尝试偏移的数目、初始值和 NCC 路径的宽度。
from PIL import Image
import numpy
import stereo
im_l = numpy.array(Image.open('data/scene1.row3.col3.ppm').convert('L'),'f')
im_r = numpy.array(Image.open('data/scene1.row3.col4.ppm').convert('L'),'f')
steps = 12
start = 4
wid = 5
res = stereo.plane_sweep_ncc(im_l, im_r, start, steps, wid)
import scipy.misc
scipy.misc.imsave('out_depth_5.png', res)
三、结果与分析
3.1使用均匀滤波器(plane_sweep_ncc)
不同窗口值下视差图的像素深度值:
分析:
- wid=1时,由于窗口值过小,窗口均值即为当前像素点的灰度值,归一化计算相减后得到的像素值为0,故左右两图片没有匹配出结果。
- wid=2时,由图片结果与原图对比可看出,图片中噪声点很多,最清晰的区域时近处的橙色台灯,在视差图中呈白色,其次是在台灯后方的雕塑,虽然和后方景物有一定色差,肉眼细看可以辨别出雕塑轮廓,最远处的书架仅有部分轮廓,桌子及其上方的物体与后方书架,难以辨别出清晰的轮廓,稍远处的摄像机不仔细看,几乎看不到轮廓。
- wid=3时,图片中噪声点依旧多,最清晰的还是近处的台灯,雕塑和后方背景色差大了一些,稍远处的摄像机可以依稀看见轮廓信息,桌子及其上方的物体与后方书架轮廓也比wid=2时明显了一些。
- wid=4时,图片中噪声点相比wid=2和wid=3时明显减少,各个景物的轮廓也比前两张图更加清晰可见。
- 当wid较小时,存在许多与正确匹配点像素值相近的点,因而匹配时找到的一些点并非实际匹配点,导致视差图噪点很多。当wid值变大一些,当前窗口值小,窗口像素平均值(与更小的窗口值比)大,归一化计算的分子值增大,则结果NCC值增大,窗口匹配度高,因而各物体区域变清晰。近处视差大,表现为视差图中物体越亮(如视差图图中台灯是白色的)。
分析:
- 这4张视差图,各物体轮廓逐渐更加清晰,远处的书架和书桌细节逐渐消失,摄像机和雕塑几乎只有形状,摄像机形状区域已经形变,雕塑形状区域较清晰,近处台灯仍可以看见部分细节。
- 由近处视差大,图像深度小,因而近处台灯形状较清晰。远处视差小,图像深度大,因而远处书架区域灰度值近似,匹配后呈现在视差图中是灰度值较小(趋于黑色)的一片。
- 窗口值小时可以看见的细节点也清晰,窗口值增大,细节变得模糊,窗口值增大,匹配的区域增大,因而细节逐渐消失,变成区域与区域相匹配,匹配出的结果更多是物体的形状。
不同窗口值下视差图的像素深度值:
分析:
- 这4张图,所有物体细节已经看不清了,并且随着窗口值增大,视差图中各物体轮廓模糊,甚至看不出物体原本的形状,只能大致看出是几块灰度值不同的区域。
- 由视差图的像素深度值,也可以看出,在wid=9和wid=13时,图片按左上、右上、左下,右下区域深度值是一致的,窗口值增大,物体区域灰度值与其周围灰度值区域趋于相同,像素点匹配后, 视差图中各物体轮廓模糊。
不同窗口值下视差图的像素深度值:
分析:
- 当wid继续增大,取值15/17/20/25,得到的视差图相似度很高,可以看到视差图明显的大值分为5块区域,远处的书架完全为黑色,近处的台灯则是白色区域,雕塑、摄像机和桌子呈由浅到深的灰度,每个物体的轮廓信息已经不是物体本身的轮廓区域,但各区域边缘较清晰。
- 分析原因在于,窗口值增大后,在像素点匹配时,存在很多灰度值相似的点,归一化计算匹配后,误匹配了很多实际不匹配的像素点,导致灰度值相同的区域增大,从而在视差图中显示物体没有了原本的轮廓。
3.2使用高斯滤波器(plane_sweep_gauss)(去除噪声)
结果截取了不同窗口值下,使用高斯滤波器与均匀滤波器的NCC算法中使用普通滤波器的视差图,对比如下
- wid=2时,高斯滤波器的视差图几乎看不到噪点但也已经没有细节信息,近处台灯和雕塑轮廓清晰,较远的摄像机和桌面轮廓模糊,只能看到大致区域而使用均匀滤波器的NCC算法噪点明显,但可以看到物体细节信息。
- wid=9时,高斯滤波器的视差图由于窗口值增大,匹配的窗口像素区域已经不再只是物体本身形状的区域,所以完全失去物体轮廓信息,而均匀滤波器的NCC算法台灯的雕塑的轮廓信息还是很清晰的。
- 高斯滤波器是一种线性平滑滤波器,能够有效的抑制噪声,其滤波器的模板是对二维高斯函数离散得到。查阅资料,得知高斯滤波器的作用原理和均值滤波器类似,都是取滤波器窗口内的像素的均值作为输出。其窗口模板的系数和均值滤波器不同,均值滤波器的模板系数都是相同的为1。由于高斯模板的中心值最大,四周逐渐减小,其滤波后的结果相对于均值滤波器来说更平滑。
- 高斯版本具有较少的噪声,视差图更加平滑,但缺少很多细节信息。NCC算法在窗口较小时,噪声较多,随着窗口值增大,看到的噪声变少,物体细节也减少。
四、小结
4.1关于计算视差图的结论
-
近处视差大,图像深度小,远处视差小,图片深度大,表现为在视差图中近处物体,轮廓更加清晰
-
窗口值小视差图对噪声点敏感,但精度更高、细节更丰富;窗口大,对噪声点比较不敏感,精度不高、细节不够,只能看到物体的轮廓
-
窗口值不是越小越好,窗口值太小,视差图会存在很多噪点;窗口值也不是越大越好,窗口值增大,窗口匹配区域增大,窗口值太大,会失去物体轮廓信息。
4.2 高斯滤波器与均匀滤波器对比结论
- 高斯版本的NCC算法具有较少的噪声,视差图更加平滑,但缺少很多细节信息。均匀滤波器版本的NCC算法在窗口较小时,噪声较多,随着窗口值增大,噪声变少,物体细节也减少。
4.3关于NCC算法的结论
优点:
- 双目立体视觉对环境光照非常敏感,拍摄的两张图片亮度有差别,曝光程度不同,对计算视差图会有较大影响,而用NCC算法,采用当前像素值与窗口平均值做差的方法,可以一定程度减小两张图片的曝光差异。
缺点:
- NCC算法是基于窗口的匹配方法计算得到视差图,但这种方法效果并不理想,窗口值太小会有明显噪声,太大又会失去物体轮廓信息,适中得到的视差图每个物体的轮廓也无法十分清晰;由于要逐点进行滑动窗口匹配,计算效率也不是很高。
展望:
本文是用NCC算法,也可以是用SAD(绝对差值和)算法、SSD(差值平方和)算法、STAD(截断绝对差值和)算法等计算匹配代价。
匹配代价:判断参考图像和目标图像中两点为对应匹配点的可能性。