局部立体匹配算法介绍及代码实现
时间：2022-08-18 16:30:03
作者Idulingwen@CSDN
编辑I3D视觉开发者社区
01
局部匹配算法是什么？优点是什么？
局部(Local)三维匹配相对于半全局和全局(Non-Local)就立体匹配算法而言，它不构建能量函数，而是使用某种成本函数（或相似度量），只通过比较左右视图中相同大小的图像块来确定视差。其基本过程一般为成本计算、成本聚合、视差计算和视差细化。虽然非局部立体匹配算法在性能上可能优于局部算法，但也有很多困难，并不是所有情况下的最佳选择，比如:
全局算法或半全局算法可能需要相当多的计算量，特别是对于1080等高分辨率图像x720、1920x即使使用硬件加速，也很难实时实现1080等。然而，局部算法占用的计算资源很少，计算效率也很高。滑动窗口的数据再利用策略和并行编程可以实时实现。如果能实现硬件，就能充分发挥其无与伦比的速度优势，特别适合计算能力低的计算平台。
对于包含结构光发射器的双目或单目传感器，与局部匹配算法的匹配效果相比，全局匹配算法的差距可能不大。当对视差图的精度和完整性要求不高时，选择运行效率较高的局部匹配算法可能是更好的选择。
02
代价函数
匹配成本的常用方法如下：SAD/SSD/ZSAD/BT、NCC/ZNCC、Census/StarCensus、rank这里就不详细描述这些成本方法的原理了。下面通过OpenCV为了实现局部三维匹配算法，为了简化代码的复杂性，成本聚合部分将被忽略，视差细化部分仅包括亚像素插入值、唯一检测、左右一致性检测和中值滤波。假设输入图像的算法是通过极线校正的。
03
代码实现
让我们把代码放在这里供您参考，核心匹配部分的代码约为100行。完整的代码可以在这里下载：
https://download.csdn.net/download/dulingwen/20449594
1、头文件LocalStereoMatch.h
#include typedef char int8; typedef unsigned char uint8; typedef short int16; typedef int int32; typedef float float32; enum { LSM_SAD = 1, LSM_NCC = 2, LSM_CENSUS = 3, LSM_RANK = 4, }; enum { LSM_LINEAR = 1, LSM_PARABOLA = 2, LSM_EQUALIZE = 3, }; /** * @brief 局部立体匹配算法 * @param left 左图像 * @param right 右图像 * @param[out] disp 视差图 * @param window_size 匹配窗口尺寸 * @param max_disparities 最大搜索视差 * @param mode 匹配成本类型 * @param ip_mode 亚像素插值方式：线性插值、抛物线插值、直方图均衡线性插值 * @param uinqueness_ratio 唯一检测比率 * @param lrc_threshold 左右一致性检测阈值 */ void LocalStereoMatching(const cv::Mat& left, const cv::Mat& right, cv::Mat& disp, intwindow_size=11,intmax_disparities=64,intmode=LSM_SAD,intip_mode=LSM_PARABOLA,intuinqueness_ratio=5,intlrc_threshold=1); 2.源文件LocalStereoMatch.cpp
#include "LocalStereoMatch.h" // ---------------------------------------------------------------------------------------------------------------------------------- template inline float32 parabola_ip(T c1, T c2, T c, float32 d) { float32 denom = std::max((float32)1, (float32)(c1 - c c2 - c)); float32 diff = (float32)(c1 - c2) / (denom * 2.f); float32 ds = d diff; return ds; } template inline float32 linear_ip(T c1, T c2, T c, float32 d) { float32 ldiff = c1 - c; float32 rdiff = c2 - c; float32 diff = 0; if (ldiff <= rdiff) diff = -0.5f 0.5f * (ldiff / rdiff); else diff = 0.5f - 0.5f * (rdiff / ldiff); float32 ds = d diff; return ds; } template inline float32 equiangular_equalize_ip(T c1, T c2, T c, float32 d) { float32 ldiff = c1 - c; float32 rdiff = c2 - c; float32 diff = 0; if (ldiff <= rdiff) { float32 x = ldiff / rdiff; diff = -0.5f 0.25f * (x * x x); } else { float32 x = rdiff / ldiff; diff = 0.5f - 0.25f * (x * x x); } float32 ds = d diff; return ds; } // ---------------------------------------------------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------------------------------------------------- static int32 compute_cost_sad(const uint8* lptr, const uint8* rptr, int x, int y, int w, int rad, int d) { int32 sad = 0; for (int i = x - rad; i < x rad; i ) { for (int j = y - rad; j < y rad; j ) { uint8 val1 = lptr[i * w j]; uint8 val2 = rptr[i * w std::max(j - d, 0)]; sad = std::abs(val1 - val2); } } return sad; } static float32 compute_cost_ncc(const uint8* lptr, const uint8* rptr, int x, int y, int w, int rad, int d) { int32 a = 0, b = 0, c = 0; for (int i = x - rad; i < x rad; i ) { for (int j = y - rad; j < y rad; j ) { uint8 val1 = lptr[i * w j]; uint8 val2 = rptr[i * w std::max(j - d, 0)]; a = val1 * val1; b = val2 * val2; c = val1 * val2; } } float32 ncc = FLT_MAX; if (c != 0) { ncc = (float32)(std::sqrt(a) * std::sqrt(b)) / (float32)c; } return ncc; } static int16 compute_cot_census(const uint8* lptr, const uint8* rptr, int x, int y, int w, int rad, int d)
{
int16 hamming = 0;
const uint8& cval1 = lptr[x * w + y];
const uint8& cval2 = rptr[x * w + std::max(y - d, 0)];
for (int i = x - rad; i < x + rad; i++)
{
for (int j = y - rad; j < y + rad; j++)
{
if (i == x && j == y) continue;
uint8 val1 = lptr[i * w + j];
uint8 val2 = rptr[i * w + std::max(j - d, 0)];
if (val1 > cval1 && val2 <= cval2)
hamming += 1;
else if (val1 <= cval1 && val2 > cval2)
hamming += 1;
}
}
return hamming;
}
static int16 compute_cost_rank(const uint8* lptr, const uint8* rptr, int x, int y, int w, int rad, int d)
{
int16 F = 0;
const uint8& cval1 = lptr[x * w + y];
const uint8& cval2 = rptr[x * w + std::max(y - d, 0)];
int16 u = 2;
int16 v = 9;
int16 iu = -u;
int16 iv = -v;
int8 rank1 = 0;
int8 rank2 = 0;
for (int i = x - rad; i < x + rad; i++)
{
for (int j = y - rad; j < y + rad; j++)
{
int16 val1 = lptr[i * w + j];
int16 val2 = rptr[i * w + std::max(j - d, 0)];
int16 diff1 = val1 - cval1;
int16 diff2 = val2 - cval2;
if (diff1 < iv) rank1 = -2;
else if (diff1 >= iv && diff1 < iu) rank1 = -1;
else if (diff1 >= iu && diff1 <= u) rank1 = 0;
else if (diff1 > u && diff1 <= v) rank1 = 1;
else if (diff1 > v) rank1 = 2;
if (diff2 < iv) rank2 = -2;
else if (diff2 >= iv && diff2 < iu) rank2 = -1;
else if (diff2 >= iu && diff2 <= u) rank2 = 0;
else if (diff2 > u && diff2 <= v) rank2 = 1;
else if (diff2 > v) rank2 = 2;
F += std::abs(rank1 - rank2);
}
}
return F;
}
// ----------------------------------------------------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------------------------------------------------
template
void pseudo_lrc_check(CostType* cost_mat, float32* disp, const int& w, const int& h, const double& lrc_thr, float32 invalid_disp)
{
if (!cost_mat || lrc_thr <= 0)
return;
CostType* disp2cost = (CostType*)malloc(w * sizeof(CostType));
float32* disp2buf = (float32*)malloc(w * sizeof(float32));
if (!disp2cost || !disp2buf)
return;
CostType max_cost = (CostType)0;
if (std::is_same::value)
max_cost = FLT_MAX;
else if (std::is_same::value)
max_cost = INT_MAX;
else if (std::is_same::value)
max_cost = SHRT_MAX;
int y = 0;
for (int x = 0; x < h; x++)
{
float32* dptr = disp + x * w;
for (y = 0; y < w; y++)
{
disp2buf[y] = invalid_disp;
disp2cost[y] = max_cost;
}
const CostType* cptr = cost_mat + x * w;
for (y = 0; y < w; y++)
{
float32 d = dptr[y];
CostType c = cptr[y];
if (d == invalid_disp)
continue;
int y2 = (int)(y - d);
if (y2 < 0)
continue;
if (disp2cost[y2] > c)
{
disp2cost[y2] = c;
disp2buf[y2] = d;
}
}
for (y = 0; y < w; y++)
{
float32 d0 = dptr[y];
if (d0 == invalid_disp)
continue;
float32 d1 = d0 + 1.f;
int y0 = (int)(y - d0);
int y1 = (int)(y - d1);
if ((0 <= y0 && y0 < w && disp2buf[y0] > invalid_disp && std::fabs(disp2buf[y0] - d0) > lrc_thr) &&
(0 <= y1 && y1 < w && disp2buf[y1] > invalid_disp && std::fabs(disp2buf[y1] - d0) > lrc_thr))
dptr[y] = invalid_disp;
}
}
free(disp2cost);
free(disp2buf);
}
// ----------------------------------------------------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------------------------------------------------
template
void matchImpl(const uint8* lptr, const uint8* rptr, float32* dptr, const int& w, const int& h, int& ws, int& ndisp, int mode, int ip_mode, int uniq, int lrc_thr)
{
int num = ws * ws;
int rad = ws / 2;
CostType* cost0 = (CostType*)malloc((ndisp + 2) * sizeof(CostType));
memset(cost0, 0, (ndisp + 2) * sizeof(CostType));
CostType* cost = cost0 + 1;
if (std::is_same::value)
uniq = 0; // ncc代价不是在任何情况下都适合用下面的唯一性检测，故设为0
float32 invalid_disp = -1.f;
CostType* cost_mat = nullptr;
if (lrc_thr > 0)
{
cost_mat = (CostType*)malloc(w * h * sizeof(CostType));
}
for (int i = rad; i < h - rad; i++)
{
for (int j = rad; j < w - rad; j++)
{
for (int k = 0; k < ndisp + 2; k++)
{
cost0[k] = (CostType)0;
}
// 代价计算
for (int d = 0; d < ndisp; d++)
{
if (mode == 1)
cost[d] = compute_cost_sad(lptr, rptr, i, j, w, rad, d);
if (mode == 2)
cost[d] = compute_cost_ncc(lptr, rptr, i, j, w, rad, d);
if (mode == 3)
cost[d] = compute_cost_census(lptr, rptr, i, j, w, rad, d);
if (mode == 4)
cost[d] = compute_cost_rank(lptr, rptr, i, j, w, rad, d);
}
// 寻找最优视差
CostType mincost = 0;
int32 bestd = 0;
if (std::is_same::value) mincost = FLT_MAX;
else if (std::is_same::value) mincost = INT_MAX;
else if (std::is_same::value) mincost = SHRT_MAX;
for (int d = 0; d < ndisp; d++)
{
const CostType& currcost = cost[d];
if (currcost < mincost)
{
mincost = currcost;
bestd = d;
}
}
if (cost_mat)
cost_mat[i * w + j] = mincost;
// 唯一性检测
if (uniq > 0)
{
CostType thresh = mincost + mincost * ((CostType)(uniq)) / ((CostType)(100));
int d = 0;
for (d = 0; d < ndisp; d++)
{
const auto& costp = cost[d];
if ((d < bestd - 1 || d > bestd + 1) && costp <= thresh)
break;
}
if (d < ndisp)
{
dptr[i * w + j] = invalid_disp;
continue;
}
}
// 亚像素插值
const CostType cost1 = cost[bestd - 1];
const CostType cost2 = cost[bestd + 1];
if (((bestd - 1) == -1) || ((bestd + 1) == ndisp))
{
dptr[i * w + j] = (float32)bestd;
continue;
}
if (ip_mode == 1)
dptr[i * w + j] = parabola_ip(cost1, cost2, mincost, bestd);
else if (ip_mode == 2)
dptr[i * w + j] = linear_ip(cost1, cost2, mincost, bestd);
else if (ip_mode == 3)
dptr[i * w + j] = equiangular_equalize_ip(cost1, cost2, mincost, bestd);
else
dptr[i * w + j] = (float32)bestd;
}
}
// 伪左右一致性检测
pseudo_lrc_check(cost_mat, dptr, w, h, lrc_thr, invalid_disp);
if (cost_mat)
free(cost_mat);
cost_mat = nullptr;
}
void LocalStereoMatching(const cv::Mat& left, const cv::Mat& right, cv::Mat& disp, int window_size, int max_disparities,
int mode, int ip_mode, int uinqueness_ratio, int lrc_threshold)
{
assert(!left.empty());
assert(!right.empty());
assert(left.size() == right.size());
const int width = left.cols;
const int height = left.rows;
int ws = window_size;
int rad = ws / 2;
int num = ws * ws;
int ndisp = max_disparities;
const uint8* lptr = left.data;
const uint8* rptr = right.data;
disp.release();
disp.create(cv::Size(width, height), CV_32FC1);
float32* dptr = disp.ptr();
if (mode == 1 /*sad*/)
matchImpl(lptr, rptr, dptr, width, height, ws, ndisp, mode, ip_mode, uinqueness_ratio, lrc_threshold);
else if (mode == 2 /*ncc*/)
matchImpl(lptr, rptr, dptr, width, height, ws, ndisp, mode, ip_mode, uinqueness_ratio, lrc_threshold);
else if (mode == 3 /*census*/)
matchImpl(lptr, rptr, dptr, width, height, ws, ndisp, mode, ip_mode, uinqueness_ratio, lrc_threshold);
else if (mode == 4 /*rank*/)
matchImpl(lptr, rptr, dptr, width, height, ws, ndisp, mode, ip_mode, uinqueness_ratio, lrc_threshold);
else
{
printf("Error: unsupported alogrithm.\n");
return;
}
cv::medianBlur(disp, disp, 3);
}
// ---------------------------------------------------------------------------------------------------------------------------------- 04
结果及运算效率分析
测试图像采用了MiddleBurry网站的双目图像，各个匹配算法采用了统一的参数设置：window_size = 17, max_disparities = 128~145, uniqueness_ratio = 5, lrc_threshold = 1，结果如下图，从左到右使用的匹配算法依次为SAD、NCC、Census、Rank。从中可以看出NCC的匹配效果是最好的，具有良好的抗噪性能，生成的视差图完整度最高，但同时噪声也最多，SAD次之，然后是rank和census。其中SAD的运算量最低，因此运算速度最快，只需十几秒便可以计算一张图像。
05
如何提升局部立体匹配算法性能？
匹配代价至关重要，最常用的改善方式是使用自适应权重对匹配代价进行加权，这通常适用于pixel-wise的匹配代价。还有的是使用多个不同的窗口尺寸进行匹配，从中选择置信度最高的结果作为最终结果，但这种方法的计算量会比较高，实际使用起来并不划算。随着深度学习的发展，深度学习主导的特征提取方式大大超越了传统的手工特征设计，因此一些基于卷积神经网络的度量学习开始用来替代传统的匹配代价，比较典型的有MC-CNN、HashMatch等。
代价聚合能够很好的提升匹配精度与鲁棒性，至今为止已发展出多种代价聚合方式：普通的邻域求和、基于引导滤波的代价聚合、基于双边滤波的代价聚合、基于十字交叉臂的代价聚合以及跨尺度代价聚合等等，另外SGM中的动态规划部分通常也被认为是一种代价聚合方式。
版权声明：本文为作者授权转载，由3D视觉开发者社区编辑整理发布，仅做学术分享，未经授权请勿二次传播，版权归原作者所有，图片来源于网络，若涉及侵权内容请联系删文。
本文仅做学术分享，如有侵权，请联系删文。
3D视觉精品课程推荐：
1.面向自动驾驶领域的多传感器数据融合技术
2.面向自动驾驶领域的3D点云目标检测全栈学习路线！(单模态+多模态/数据+代码)
3.彻底搞透视觉三维重建：原理剖析、代码讲解、及优化改进
4.国内首个面向工业级实战的点云处理课程
5.激光-视觉-IMU-GPS融合SLAM算法梳理和代码讲解
6.彻底搞懂视觉-惯性SLAM：基于VINS-Fusion正式开课啦
7.彻底搞懂基于LOAM框架的3D激光SLAM: 源码剖析到算法优化
8.彻底剖析室内、室外激光SLAM关键算法原理、代码和实战(cartographer+LOAM +LIO-SAM)
9.从零搭建一套结构光3D重建系统[理论+源码+实践]
10.单目深度估计方法：算法梳理与代码实现
11.自动驾驶中的深度学习模型部署实战
12.相机模型与标定(单目+双目+鱼眼）
重磅！3DCVer-学术论文写作投稿 交流群已成立
扫码添加小助手微信，可申请加入3D视觉工坊-学术论文写作与投稿 微信交流群，旨在交流顶会、顶刊、SCI、EI等写作与投稿事宜。
同时也可申请加入我们的细分方向交流群，目前主要有3D视觉、CV&深度学习、SLAM、三维重建、点云后处理、自动驾驶、多传感器融合、CV入门、三维测量、VR/AR、3D人脸识别、医疗影像、缺陷检测、行人重识别、目标跟踪、视觉产品落地、视觉竞赛、车牌识别、硬件选型、学术交流、求职交流、ORB-SLAM系列源码交流、深度估计等微信群。
一定要备注：研究方向+学校/公司+昵称，例如：”3D视觉 + 上海交大 + 静静“。请按照格式备注，可快速被通过且邀请进群。原创投稿也请联系。
▲长按加微信群或投稿
▲长按关注公众号
3D视觉从入门到精通知识星球：针对3D视觉领域的视频课程（三维重建系列、三维点云系列、结构光系列、手眼标定、相机标定、激光/视觉SLAM、自动驾驶等）、知识点汇总、入门进阶学习路线、最新paper分享、疑问解答五个方面进行深耕，更有各类大厂的算法工程人员进行技术指导。与此同时，星球将联合知名企业发布3D视觉相关算法开发岗位以及项目对接信息，打造成集技术与就业为一体的铁杆粉丝聚集区，近4000星球成员为创造更好的AI世界共同进步，知识星球入口：
学习3D视觉核心技术，扫描查看介绍，3天内无条件退款
圈里有高质量教程资料、可答疑解惑、助你高效解决问题
觉得有用，麻烦给个赞和在看~
