【码上实战】【立体匹配系列】经典AD-Census: （5）扫描线优化_综合

同学们国庆假期快乐呀！潇洒7天（带娃7天），难得坐下来更新博客。

下载完整源码，点击进入: https://github.com/ethan-li-coding/AD-Census
欢迎同学们在Github项目里讨论！

接上篇十字交叉域代价聚合，本篇的内容是AD-Census的扫描线优化步骤，实际上，该步的思路和SGM的代码聚合是基本一样的，只不过在P1/P2参数设置上做了一些修改。确实，SGM的P1、P2设置策略过于简单，优点是鲁棒性高，对大部分数据都能得到一个还不错的视差结果，但明显的弊端就是很难找到一组特别好的参数组合，使特定应用场景的数据达到比较完美的状态，P1/P2的设置对整体视差效果尤其是边缘处的视差很关键，所以AD-Census的改进方向是有实际意义的。

我们不妨直接先看下AD-Census扫描线优化的成果：

代价计算

代价聚合

扫描线优化

显然，扫描线优化后的视差图比代价聚合后的视差图更加完整，错误值更少。当然这并不能说明AD-Census的参数改进就是有效的，只能说明扫描线优化步骤是有效的。

我们来看编码介绍吧！

文章目录

- 算法
- 代码实现
- - 类设计
  - - 成员函数
    - 成员变量
  - 类实现
- 实验

算法

同样的，算法原理请看博文：

经典AD-Census: （3）扫描线优化（Scanline Optimization）

这里我就不再展开讲优化的原理了，和SGM（SemiGlobalMatching）的代价聚合策略确实是一模一样，看博主往期博客就行了，AD-Census采用4方向的扫描线优化，即上下左右4个方向。

AD-Census所做的修改在于 $P_1$ 和 $P_2$ 值的设定方式，在SGM中， $P_1$ 、 $P_2'$ 是预设的固定值，实际使用的 $P_2$ 是根据左视图相邻两个像素的亮度差值而实时调整的，调整公式为 $P_2=P_2'/(I_p-I_q)$ 。

而在Ad-Census中， $P_1$ 、 $P_2$ 不只是和左视图的相邻像素颜色差 $D_1=D_c(p,p-r)$ 有关，而且和右视图对应同名点的相邻像素颜色差 $D_2=D_c(pd,pd-r)$ 有关。

（注1：AD-Census算法默认输入彩色图，所以是算颜色差，如果是输入灰度图，则是亮度差，颜色差的定义是 $D_c(p_l,p)=max_{i=R,G,B}|I_i(p_l)-I_i(p)|$ ，即三个颜色分量差值的最大值）
（注2： $p d$ 实际就是像素 $p$ 通过视差 $d$ 找到的右视图上的同名点 $q = p ? d$ ）
（注3： $p ? r$ 代表聚合方向上的上一个像素，比如从左到右聚合，则 $p ? r$ 就是 $p ? 1$ ；从右到左聚合，则 $p ? r$ 就是 $p + 1$ ）

具体设定规则如下：

$P_1=Π_1,P_2=Π_2, if D_1<τ_{SO},D_2<τ_{SO}$
$P_1=Π_1/4,P_2=Π_2/4, if D_1<τ_{SO},D_2>τ_{SO}$
$P_1=Π_1/4,P_2=Π_2/4, if D_1>τ_{SO},D_2<τ_{SO}$
$P_1=Π_1/10,P_2=Π_2/10, if D_1>τ_{SO},D_2>τ_{SO}$

$Π_1,Π_2$ 是设定的固定阈值， $τ_{SO}$ 是设定的颜色差阈值。

代码实现

类设计

成员函数

同样，我们用一个扫描线优化器类ScanlineOptimizer来实现该功能。放在文件scanline_optimizer.h/scanline_optimizer.cpp中。

/*** \brief 扫描线优化器*/
class ScanlineOptimizer {
    
public:ScanlineOptimizer();~ScanlineOptimizer();
}

在公有成员函数的设计上，第一类接口是必不可少的 设置数据SetData 以及 设置参数SetParam ，完成算法的输入。第二类就是优化功能接口 Optimize 。

而具体的优化子步骤，我们放在私有成员函数列表里，包括水平方向聚合 CostAggregateLeftRight 以及竖直方向聚合 CostAggregateUpDown。

同时，算法需要的一个小功能颜色距离计算函数 ColorDist，也放在私有函数中。

所有成员函数的声明代码如下：

public:ScanlineOptimizer();~ScanlineOptimizer();/*** \brief 设置数据* \param img_left // 左影像数据，三通道 * \param img_right // 右影像数据，三通道* \param cost_init // 初始代价数组* \param cost_aggr // 聚合代价数组*/void SetData(const uint8* img_left, const uint8* img_right, float32* cost_init, float32* cost_aggr);/*** \brief * \param width // 影像宽* \param height // 影像高* \param min_disparity // 最小视差* \param max_disparity // 最大视差* \param p1 // p1* \param p2 // p2* \param tso // tso*/void SetParam(const sint32& width,const sint32& height, const sint32& min_disparity, const sint32& max_disparity, const float32& p1, const float32& p2, const sint32& tso);/*** \brief 优化 */void Optimize();private:/*** \brief 左右路径聚合 → ←* \param cost_so_src 输入，SO前代价数据* \param cost_so_dst 输出，SO后代价数据* \param is_forward 输入，是否为正方向（正方向为从左到右，反方向为从右到左）*/void CostAggregateLeftRight(const float32* cost_so_src, float32* cost_so_dst, bool is_forward = true);/*** \brief 上下路径聚合 ↓ ↑* \param cost_so_src 输入，SO前代价数据* \param cost_so_dst 输出，SO后代价数据* \param is_forward 输入，是否为正方向（正方向为从上到下，反方向为从下到上）*/void CostAggregateUpDown(const float32* cost_so_src, float32* cost_so_dst, bool is_forward = true);/** \brief 计算颜色距离 */inline sint32 ColorDist(const ADColor& c1, const ADColor& c2) {
    return std::max(abs(c1.r - c2.r), std::max(abs(c1.g - c2.g), abs(c1.b - c2.b)));}

为每个函数都写了清晰的注释，便于快速理解。此外计算颜色距离的函数为内联函数，声明的同时也定义实现了它。

成员变量

成员变量全部都被设计为私有，仅在算法内部使用，他们是图像尺寸、影像数据、代价数据（初始/聚合）、算法参数等。

private:/** \brief 图像尺寸 */sint32	width_;sint32	height_;/** \brief 影像数据 */const uint8* img_left_;const uint8* img_right_;/** \brief 初始代价数组 */float32* cost_init_;/** \brief 聚合代价数组 */float32* cost_aggr_;/** \brief 最小视差值 */sint32 min_disparity_;/** \brief 最大视差值 */sint32 max_disparity_;/** \brief 初始的p1值 */float32 so_p1_;/** \brief 初始的p2值 */float32 so_p2_;/** \brief tso阈值 */sint32 so_tso_;

类实现

由于SetData和SetParam比较简单，代码量也很少，所以就不做介绍了，大家看代码就懂了。这里就介绍下扫描线优化的两个子步骤 CostAggregateLeftRight和 CostAggregateUpDown。

实际上，我是直接把SGM的代价聚合代码搬过来，修改 $P_1$ 和 $P_2$ 值的计算方式就行了。如下：

void ScanlineOptimizer::CostAggregateLeftRight(const float32* cost_so_src, float32* cost_so_dst, bool is_forward)
{
    const auto width = width_;const auto height = height_;const auto min_disparity = min_disparity_;const auto max_disparity = max_disparity_;const auto p1 = so_p1_;const auto p2 = so_p2_;const auto tso = so_tso_;assert(width > 0 && height > 0 && max_disparity > min_disparity);// 视差范围const sint32 disp_range = max_disparity - min_disparity;// 正向(左->右) ：is_forward = true ; direction = 1// 反向(右->左) ：is_forward = false; direction = -1;const sint32 direction = is_forward ? 1 : -1;// 聚合for (sint32 y = 0u; y < height; y++) {
    // 路径头为每一行的首(尾,dir=-1)列像素auto cost_init_row = (is_forward) ? (cost_so_src + y * width * disp_range) : (cost_so_src + y * width * disp_range + (width - 1) * disp_range);auto cost_aggr_row = (is_forward) ? (cost_so_dst + y * width * disp_range) : (cost_so_dst + y * width * disp_range + (width - 1) * disp_range);auto img_row = (is_forward) ? (img_left_ + y * width * 3) : (img_left_ + y * width * 3 + 3 * (width - 1));const auto img_row_r = img_right_ + y * width * 3;sint32 x = (is_forward) ? 0 : width - 1;// 路径上当前颜色值和上一个颜色值ADColor color(img_row[0], img_row[1], img_row[2]);ADColor color_last = color;// 路径上上个像素的代价数组，多两个元素是为了避免边界溢出（首尾各多一个）std::vector<float32> cost_last_path(disp_range + 2, Large_Float);// 初始化：第一个像素的聚合代价值等于初始代价值memcpy(cost_aggr_row, cost_init_row, disp_range * sizeof(float32));memcpy(&cost_last_path[1], cost_aggr_row, disp_range * sizeof(float32));cost_init_row += direction * disp_range;cost_aggr_row += direction * disp_range;img_row += direction * 3;x += direction;// 路径上上个像素的最小代价值float32 mincost_last_path = Large_Float;for (auto cost : cost_last_path) {
    mincost_last_path = std::min(mincost_last_path, cost);}// 自方向上第2个像素开始按顺序聚合for (sint32 j = 0; j < width - 1; j++) {
    color = ADColor(img_row[0], img_row[1], img_row[2]);const uint8 d1 = ColorDist(color, color_last);uint8 d2 = d1;float32 min_cost = Large_Float;for (sint32 d = 0; d < disp_range; d++) {
    const sint32 xr = x - d;if (xr > 0 && xr < width - 1) {
    const ADColor color_r = ADColor(img_row_r[3 * xr], img_row_r[3 * xr + 1], img_row_r[3 * xr + 2]);const ADColor color_last_r = ADColor(img_row_r[3 * (xr - direction)],img_row_r[3 * (xr - direction) + 1],img_row_r[3 * (xr - direction) + 2]);d2 = ColorDist(color_r, color_last_r);}// 计算P1和P2float32 P1(0.0f), P2(0.0f);if (d1 < tso && d2 < tso) {
    P1 = p1; P2 = p2;}else if (d1 < tso && d2 >= tso) {
    P1 = p1 / 4; P2 = p2 / 4;}else if (d1 >= tso && d2 < tso) {
    P1 = p1 / 4; P2 = p2 / 4;}else if (d1 >= tso && d2 >= tso) {
    P1 = p1 / 10; P2 = p2 / 10;}// Lr(p,d) = C(p,d) + min( Lr(p-r,d), Lr(p-r,d-1) + P1, Lr(p-r,d+1) + P1, min(Lr(p-r))+P2 ) - min(Lr(p-r))const float32  cost = cost_init_row[d];const float32 l1 = cost_last_path[d + 1];const float32 l2 = cost_last_path[d] + P1;const float32 l3 = cost_last_path[d + 2] + P1;const float32 l4 = mincost_last_path + P2;float32 cost_s = cost + static_cast<float32>(std::min(std::min(l1, l2), std::min(l3, l4)));cost_s /= 2;cost_aggr_row[d] = cost_s;min_cost = std::min(min_cost, cost_s);}// 重置上个像素的最小代价值和代价数组mincost_last_path = min_cost;memcpy(&cost_last_path[1], cost_aggr_row, disp_range * sizeof(float32));// 下一个像素cost_init_row += direction * disp_range;cost_aggr_row += direction * disp_range;img_row += direction * 3;x += direction;// 像素值重新赋值color_last = color;}}
}

如果不了解聚合代码，可以看我此前博客：

编码实现经典SGM：（3）代价聚合

本篇我们重点看下P1和P2的计算方式：

我们首先在轮到每个像素时，计算了左视图上它与上一个像素的颜色距离（颜色差） $d_1$ ：

const uint8 d1 = ColorDist(color, color_last);

然后在遍历像素每个视差时，计算右视图对应像素与其上一个像素的颜色距离 $d_2$ 。

const sint32 xr = x - d;
if (xr > 0 && xr < width - 1) {
    const ADColor color_r = ADColor(img_row_r[3 * xr], img_row_r[3 * xr + 1], img_row_r[3 * xr + 2]);const ADColor color_last_r = ADColor(img_row_r[3 * (xr - direction)],img_row_r[3 * (xr - direction) + 1],img_row_r[3 * (xr - direction) + 2]);d2 = ColorDist(color_r, color_last_r);
}

接下来根据 $d_1$ 和 $d_2$ 与阈值的比较情况，判定为四种情况中的某一种，计算P1和P2的值。

// 计算P1和P2
float32 P1(0.0f), P2(0.0f);
if (d1 < tso && d2 < tso) {
    P1 = p1; P2 = p2;
}
else if (d1 < tso && d2 >= tso) {
    P1 = p1 / 4; P2 = p2 / 4;
}
else if (d1 >= tso && d2 < tso) {
    P1 = p1 / 4; P2 = p2 / 4;
}
else if (d1 >= tso && d2 >= tso) {
    P1 = p1 / 10; P2 = p2 / 10;
}

其中，小写的p1、p2，以及tso都是输入的算法参数。

const auto p1 = so_p1_;
const auto p2 = so_p2_;
const auto tso = so_tso_;

竖直方向的代码我就不贴了，除了方向不同，和水平方向并无其他区别，照葫芦画瓢。

在公有的优化接口 Optimize 内，只需要依次调用四个方向的优化函数就行了。

void ScanlineOptimizer::Optimize()
{
    if (width_ <= 0 || height_ <= 0 ||img_left_ == nullptr || img_right_ == nullptr ||cost_init_ == nullptr || cost_aggr_ == nullptr) {
    return;}// 4方向扫描线优化// 模块的首次输入是上一步代价聚合后的数据，也就是cost_aggr_// 我们把四个方向的优化按次序进行，并利用cost_init_及cost_aggr_间次保存临时数据，这样不用开辟额外的内存来存储中间结果// 模块的最终输出也是cost_aggr_// left to rightCostAggregateLeftRight(cost_aggr_, cost_init_, true);// right to leftCostAggregateLeftRight(cost_init_, cost_aggr_, false);// up to downCostAggregateUpDown(cost_aggr_, cost_init_, true);// down to upCostAggregateUpDown(cost_init_, cost_aggr_, false);
}

这里用了一个小技巧，即交替使用cost_aggr和cost_init，不用额外开辟四个方向的代价数组，只用两个代价数据即完成整个优化操作。

实验

我们做了三组实验，一组是只做左右水平方向的扫描线优化，一组是只做上下竖直方向的扫描线优化，剩下一组是做四个方向的优化。我们来看看效果。

代价聚合

水平方向优化

竖直方向优化

4方向优化

看上去，只做水平或者竖直优化，视差图已有明显的改进，但单方向的优化会存在方向条纹效应，而4方向的优化结果则能够消除这一现象，达到更佳的状态。

最后，我们再贴一下文章开头的实验图：

代价计算

代价聚合

扫描线优化

好了，本篇到此结束，下一篇将为大家带来的是后处理部分。感谢观看！

下载AD-Census完整源码，点击进入: https://github.com/ethan-li-coding/AD-Census
欢迎同学们在Github项目里讨论，如果觉得博主代码质量不错，右上角给颗星！感谢！

博主简介：
Ethan Li 李迎松（知乎：李迎松）
武汉大学摄影测量与遥感专业博士
主方向立体匹配、三维重建
2019年获测绘科技进步一等奖（省部级）

爱三维，爱分享，爱开源
GitHub： https://github.com/ethan-li-coding （欢迎follow和star）

个人微信：

欢迎交流！

关注博主不迷路，感谢！
博客主页：https://ethanli.blog.csdn.net