详解立体匹配系列经典SGM: (5) 视差优化
点击上方“计算机视觉life”，选择“星标”
快速获得最新干货
作者李迎松授权发布，武汉大学 摄影测量与遥感专业 博士
https://ethanli.blog.csdn.net/article/details/105065660
详解立体匹配系列经典SGM: (1) 框架与类设计 详解立体匹配系列经典SGM: (2) 代价计算
详解立体匹配系列经典SGM: (3) 代价聚合
详解立体匹配系列经典SGM: (4) 视差计算、视差优化
代码已同步于Github开源项目： https://github.com/ethan-li-coding/SemiGlobalMatching
上回说到，路径聚合后，视差图英姿初现，我想初学者到这一步大多都会产生一种溢于言表的喜悦之情，这就是编程的魅力，你日以继夜的在键盘上敲打，与世隔绝甚至废寝忘食，当程序运行，屏幕出现一个还不错的结果时，你的体内会产生一种叫多巴胺的物质，它让你兴奋和开心，并激励着你不断前行。
让我们再回顾下上篇的结果：
本篇算是进入SGM编码教学的完结篇，即立体匹配的最后一步：视差优化（Disparity Refine）。优化的含义大家都懂，为了得到更好的视差图。
怎样算更好？如何能更好？本篇就是解答这两个问题。
我们先来看看，本篇我们将得到一个什么样的实验结果，大家自行判断是否有兴趣继续阅读。
客观的说，噪点剔除的很干净，遮挡区做了处理，我想你们应该会感点兴趣。
1 优化目的
1. 提高精度
2. 剔除错误
3. 弱纹理区优化
4. 填补空洞
提高精度：前面的视差计算步骤中，我们选择最小代价值对应的视差值，它是一个整数值（整数值我们才能有离散化的视差空间WHD
剔除错误：即剔除错误的视差值，比如a像素本应和b像素是同名点，而结果却是a像素和c像素是同名点，这就是错误的视差值。造成错误匹配的原因有遮挡、弱纹理等，所谓错误是永恒的，完美是不存在的。
弱纹理区优化：弱纹理区域是所有立体匹配都会面对的难题，极端情况下，一块白墙，怎么找同名点？SGM提出的是一种基于图像分割+平面拟合的处理弱纹理的方法。
填补空洞：剔除错误匹配后，被剔除的像素会造成无效值空洞，如何填补使视差图更加完整也是优化所研究的内容。
优化做的策略不少，科学家都是精益求精，限于篇幅，本篇只讲前两条，后两条我将在下一篇更新。
2 优化手段
1. 子像素拟合（Subpixel）
2. 一致性检查（Left/Right Consistency Check）
3. 唯一性约束（Uniqueness）（待续）
4. 剔除小连通区（Remove Peaks）（待续）
5. 中值滤波（Median Filter）
子像素拟合和一致性检查几乎是所有立体匹配算法必执行的策略。子像素拟合将整像素精度提高到子像素精度，而一致性检查可以说是剔除错误匹配的不二选择。
唯一性约束和剔除小连通区可以视情况而添加，比如在GPU实现的时候，找最小值是时间复杂度较高的操作，而唯一性约束要找两次最小值（一个最小一个次最小）；而区域跟踪用GPU也难以高效的实现，所以这两块一般也可省掉。
我将教大家编码实现这5个策略。
2.1 子像素拟合
看过我之前博客的同学们应该对这个图还有点印象，左边表示视差为18时代价最小，那么18就是我们得到的整像素视差值，而右边则是把最小代价左右两边的代价值也记录下来，3个代价值做一个一元二次曲线拟合，曲线的极值点横坐标就是视差值的子像素位置。
很简单对不对，一元二次拟合谁不会！
【一元二次求解图】
我们来看代码实现（只贴子像素部分）：
// 最优视差best_disparity前一个视差的代价值cost_1，后一个视差的代价值cost_2
const sint32 idx_1 = best_disparity - 1 - min_disparity;
const sint32 idx_2 = best_disparity + 1 - min_disparity;
const uint16 cost_1 = cost_local[idx_1];
const uint16 cost_2 = cost_local[idx_2];
// 解一元二次曲线极值
const uint16 denom = std::max(1, cost_1 + cost_2 - 2 * min_cost);
disparity[i * width + j] = best_disparity + (cost_1 - cost_2) / (denom * 2.0f);
嗯，确实很简单！子像素拟合代码位于视差计算函数体内，即在计算最优视差的同时完成子像素拟合。
2.2 一致性检查
我们先了解下一致性检查是什么东东。
我们在立体匹配里会区分左右影像，这是模拟人眼立体，左视图对应左眼，右视图对应右眼，对人来说左是左右是右，可不能搞反。但对计算机来说，无所谓，左可以是右，右可以是左，我才不管是否是合理，别让我死机就行。
确实，我们让左右对调，也可以形成双目立体，两个视图也有重叠区，只不过重叠区不在中间，在两边，但管他呢，只要告诉我怎么搜同名点就行。
一致性检查就是：把左右影像位置对调，再计算一个右影像视差图，对照两个视差图来看同名点对是否能够相互匹配成功。我这里有以下两个描述，你们看哪个能够理解。
对调前，左影像像素匹配右影像像素；则对调后，也匹配为一致，否则为不一致（比如对调后匹配c）。
对调前，左影像像素的视差为；则对调后右影像像素的视差为为一致，否则为不一致。
一致性检查的一般性操作步骤是：
获取左右视差图。
对左视差图的每个像素，计算出同名点在右视差图中的像素位置。
判断和的视差值之差的绝对值是否小于一个阈值（通常为1个像素）。
如果超过阈值，则一致性检查不通过，把对应位置的视差变为无效值。
借用下SGM作者老爷子的图，这里的就是左视图，是右视图，和为同名点对。
一致性检查有两种策略，一种是内部型，一种是外部型。
1. 内部型检查
2. 外部型检查
我会实现内部型（比较有意思），描述外部型（没技术含量）。
2.2.1 内部型检查
内部型就是直接通过左影像的代价数组，来推算右影像的代价数组，从而计算右影像的视差图。所以你只用代价聚合一次就可以做一致性检查。
听起来很爽，怎么办到的？
秘诀就是：右影像(i,j)视差为的代价 = 左影像(i,j+d)视差为的代价
能理解不，就是对于右影像的像素(i,j)，根据视差值可算出左影像的对应像素位置为(i,j+d)，然后把左影像同样视差值下的代价值取出来赋给。
似乎有点绕，但是大家读几遍应该能理解。不理解的话肯定是博主描述的不够清晰，欢迎和我交流。
根据此秘诀，我们可以将右影像每个像素的所有候选视差的代价值Cost(i,j,d)都得到，进而寻找最小代价值对应的视差，并做子像素优化，得到右影像视差图。
【右影像视差图计算代码】
// ---逐像素计算最优视差
// 通过左影像的代价，获取右影像的代价
// 右cost(xr,yr,d) = 左cost(xr+d,yl,d)
for (sint32 i = 0; i < height; i++) {
for (sint32 j = 0; j < width; j++) {
uint16 min_cost = UINT16_MAX;
sint32 best_disparity = 0;
// ---统计候选视差下的代价值
for (sint32 d = min_disparity; d < max_disparity; d++) {
const sint32 d_idx = d - min_disparity;
const sint32 col_left = j + d;
if (col_left >= 0 && col_left < width) {
const auto& cost = cost_local[d_idx] = cost_ptr[i * width * disp_range + col_left * disp_range + d_idx];
if (min_cost > cost) {
min_cost = cost;
best_disparity = d;
}
}
else {
cost_local[d_idx] = UINT16_MAX;
}
}
}
}
// ---子像素拟合
if (best_disparity == min_disparity || best_disparity == max_disparity - 1) {
disparity[i * width + j] = Invalid_Float;
continue;
}
// 最优视差前一个视差的代价值cost_1，后一个视差的代价值cost_2
const sint32 idx_1 = best_disparity - 1 - min_disparity;
const sint32 idx_2 = best_disparity + 1 - min_disparity;
const uint16 cost_1 = cost_local[idx_1];
const uint16 cost_2 = cost_local[idx_2];
// 解一元二次曲线极值
const uint16 denom = std::max(1, cost_1 + cost_2 - 2 * min_cost);
disparity[i * width + j] = static_cast<float32>(best_disparity) + static_cast<float32>(cost_1 - cost_2 ) / (denom * 2.0f);
这里只贴出了右视图视差计算的部分代码，大家可从Github/GemiGlobalMatching下载完整代码。
计算出右视差图后，执行一致性检查：根据左影像视差图可以算出像素(ileft,jleft)在右影像中的匹配像素是(ileft,jleft−d)，如果(ileft,jleft−d)的视差刚好也近似等于d，则满足一致性。
【一致性检查代码】
// ---左右一致性检查
void SemiGlobalMatching::LRCheck() const
{
const int width = width_;
const int height = height_;
const float32& threshold = option_.lrcheck_thres;
// ---左右一致性检查
for (int i = 0; i < height; i++) {
for (int j = 0; j < width; j++) {
// 左影像视差值
auto& disp = disp_left_[i * width + j];
// 根据视差值找到右影像上对应的同名像素
const auto col_right = static_cast<sint32>(j - disp + 0.5);
if(col_right >= 0 && col_right < width) {
// 右影像上同名像素的视差值
const auto& disp_r = disp_right_[i * width + col_right];
// 判断两个视差值是否一致（差值在阈值内）
if (abs(disp - disp_r) > threshold) {
// 左右不一致
disp = Invalid_Float;
}
}
else{
// 通过视差值在右影像上找不到同名像素（超出影像范围）
disp = Invalid_Float;
}
}
}
}
这部分代码还是比较简单的，理解起来比较容易。
2.2.2 外部型检查
相比内部型检查，外部型检查是比较笨的办法，就是在算法输入时把左右图像数据对调，再执行一次完整的立体匹配，得到右影像视差图，一致性检查则是采用同样的策略。
这里需要注意的是，左右对调后，视差在意义上和左影像是相反的，而立体匹配算法的设定是：视差=左减右，如果只做对调（就是简单的把右影像数据作为左影像数据传进算法），是得不到正确结果的，因为对调后重叠区在两边，不符合算法的设定，所以这里我一般会在对调后，把左右影像的像素来个水平镜像翻转，这样两张影像的重叠区到了中间，视差就等于左减右了。
外部型检查需要执行两次完整的匹配流程，所以时间效率不如内部型检查。
2.4 唯一性约束
唯一性约束的含义是：最优视差的代价值应该是所有候选视差中唯一的最小代价，换句话说它的代价值比其他视差的代价值足够小，这样它才能在众多候选视差中脱颖而出。如果有另外一个视差的代价值和它一样，那么它就不再是最好的那一个了，而是最好的两个之一，而我们不能接受有两个最好，只能忍痛舍弃它。
其实，这里面蕴含的另一层含义是：视差估计的可靠性！如果两个最小的代价值相差很小，比如一个是30，一个是31，因为视差估计是带有噪声影响的，所以其实不能肯定的说是30最好还是31最好，可能31才是正确的那个视差，因为两个值相差太小了，可能是由于一些轻微的噪声导致实际上最优的那个视差的代价值却是次最小的。所以就干脆把这种让人头疼的选择题给PASS掉，遇到这种情况直接给它一个NO！
用程序来处理时，我们会计算最小代价和次最小代价的相对差值，如果差值小于阈值，那就表明最小的两个代价值相差不明显，就给这个像素赋一个无效视差。
我们来看代码：
// 判断唯一性约束
// 若(min-sec)/min < min*(1-uniquness)，则为无效估计
if (sec_min_cost - min_cost <= static_cast<uint16>(min_cost * (1 - uniqueness_ratio))) {
disparity[i * width + j] = Invalid_Float;
continue;
}
代码很简单，次最小和最小代价的相对差值小于阈值，就给个无效值。这一步也在视差计算的函数体内完成，在遍历计算最小代价的时候，再遍历一次计算次最小代价，然后判断唯一性。
2.5 剔除小连通区
学过图像处理的同学对连通区应该并不陌生，它是图像处理中很常见的一个名词，它的含义是通过4-邻域或8-邻域连通在一起的像素集合，在SGM中，这一步用来剔除连通在一起的小块错误匹配像素，像这样的：
首先这些像素是一致性检查的漏网之鱼，看上去和周围的像素极不协调，而且连通成一定大小的块，所以中值滤波这类窗口类滤波算法也搞不定，只能通过区域跟踪，把它们跟踪成块，然后判断块的大小是否小于一定的阈值，如果是则剔除，即把整块都置为无效视差。思路不复杂，就是一个区域跟踪算法。
【剔除小连通区代码】
void sgm_util::RemoveSpeckles(float32* disparity_map, const sint32& width, const sint32& height,
const sint32& diff_insame, const uint32& min_speckle_aera, const float& invalid_val)
{
assert(width > 0 && height > 0);
if (width < 0 || height < 0) {
return;
}
// 定义标记像素是否访问的数组
std::vector<bool> visited(uint32(width*height),false);
for(sint32 i=0;i<height;i++) {
for(sint32 j=0;j<width;j++) {
if (visited[i * width + j] || disparity_map[i*width+j] == invalid_val) {
// 跳过已访问的像素及无效像素
continue;
}
// 广度优先遍历，区域跟踪
// 把连通域面积小于阈值的区域视差全设为无效值
std::vector<std::pair<sint32, sint32>> vec;
vec.emplace_back(i, j);
visited[i * width + j] = true;
uint32 cur = 0;
uint32 next = 0;
do {
// 广度优先遍历区域跟踪
next = vec.size();
for (uint32 k = cur; k < next; k++) {
const auto& pixel = vec[k];
const sint32 row = pixel.first;
const sint32 col = pixel.second;
const auto& disp_base = disparity_map[row * width + col];
// 8邻域遍历
for(int r=-1;r<=1;r++) {
for(int c=-1;c<=1;c++) {
if(r==0&&c==0) {
continue;
}
int rowr = row + r;
int colc = col + c;
if (rowr >= 0 && rowr < height && colc >= 0 && colc < width) {
if(!visited[rowr * width + colc] && abs(disparity_map[rowr * width + colc] - disp_base) <= diff_insame) {
vec.emplace_back(rowr, colc);
visited[rowr * width + colc] = true;
}
}
}
}
}
cur = next;
} while (next < vec.size());
// 把连通域面积小于阈值的区域视差全设为无效值
if(vec.size() < min_speckle_aera) {
for(auto& pix:vec) {
disparity_map[pix.first * width + pix.second] = invalid_val;
}
}
}
}
}
2.6 中值滤波
中值滤波在立体匹配中使用的还挺广泛的，作为一个平滑算法，它主要是用来剔除视差图中的一些孤立的离群外点，同时还能起到填补小洞的作用。这个部分我就不细说了，想必大家对中值滤波都不会太陌生，我会在实验中展示中值滤波的效果。
中值滤波的代码实现，我放到sgm_util工具集里，如下：
void sgm_util::MedianFilter(const float32* in, float32* out, const sint32& width, const sint32& height,
const sint32 wnd_size)
{
const sint32 radius = wnd_size / 2;
const sint32 size = wnd_size * wnd_size;
// 存储局部窗口内的数据
std::vector<float32> wnd_data;
wnd_data.reserve(size);
for (sint32 i = 0; i < height; i++) {
for (sint32 j = 0; j < width; j++) {
wnd_data.clear();
// 获取局部窗口数据
for (sint32 r = -radius; r <= radius; r++) {
for (sint32 c = -radius; c <= radius; c++) {
const sint32 row = i + r;
const sint32 col = j + c;
if (row >= 0 && row < height && col >= 0 && col < width) {
wnd_data.push_back(in[row * width + col]);
}
}
}
// 排序
std::sort(wnd_data.begin(), wnd_data.end());
// 取中值
out[i * width + j] = wnd_data[wnd_data.size() / 2];
}
}
}
3 实验
枯燥而或许有趣的理论部分介绍完毕，我们喜迎实验环节。
【实验1 子像素拟合对比图】
取视差图某一行的数据，上图为拟合前，下图为拟合后，可看到整像素的阶梯形被子像素拟合优化
【实验2 一致性检查效果图】
图中我们可以看到，非重叠区、遮挡区的匹配错误已经大部分都被剔除了，当然也让视差图有了很多黑色的无效区。
【实验3 唯一性约束+去小连通区&中值滤波效果】
唯一性约束和去小连通区可以进一步剔除不可靠的视差以及连通成块的错误视差区；最后通过中值滤波，进一步对噪点进行剔除，且填补一些小的无效区域，效果还不错。
夜已深，春困秋乏，博主也该歇息了，今天就介绍到这，不出意外，后面博主还会写视差填充及弱纹理区域优化的优化模块。
代码已同步于Github开源项目：Github/GemiGlobalMatching，大家可自行下载，点击项目右上角的star，有更新会实时通知到你的个人中心！
从0到1学习SLAM，戳↓
交流群
欢迎加入公众号读者群一起和同行交流，目前有SLAM、三维视觉、传感器、自动驾驶、计算摄影、检测、分割、识别、医学影像、GAN、算法竞赛等微信群（以后会逐渐细分），请扫描下面微信号加群，备注：”昵称+学校/公司+研究方向“，例如：”张三 + 上海交大 + 视觉SLAM“。请按照格式备注，否则不予通过。添加成功后会根据研究方向邀请进入相关微信群。请勿在群内发送广告，否则会请出群，谢谢理解~
投稿、合作也欢迎联系：simiter@126.com
长按关注计算机视觉life
相关内容
微信扫码咨询专知VIP会员