Stereo Matching文献笔记之（六）：浅谈置信度传播算法（Belief-Propagation）在立体匹配中的应用~

这是我一个纠结过的问题，曾经反反复复的看相关的知识，Belief-Propagation是一个伴随着“马尔科夫随机场”提出的优化算法，我对优化算法情有独钟，一直觉得搞定了各种优化，机器学习剩下的也就是知识扩展而已，嘿嘿，我也不知道这么想是对是错，最近脑袋糊涂的厉害，请各位见谅。。。。

(转载请注明：http://blog.csdn.net/wsj998689aa/article/details/48417927, 作者：迷雾forest)

1. BP与stereo-matching

为什么要将Belief-Propagation和Stereo-Matching放在一起说？那是因为现有的全局立体匹配算法，很大一部分都是基于Belief-Propagation来进行视差求精，二者天然的纽带就是，一幅图加上对应像素标签就可以看成是一个马尔科夫场，只要有马尔科夫场，很自然地就可以利用BP算法进行优化求解。马尔科夫场是一种概率图模型，图中的节点就是像素，图中的隐节点就是标签，在Stereo-Matching里面，标签就是像素的视差值，视差值总是会有一个范围，BP就是在这个自定义空间中找到使得全局能量最小的那些视差值，一般全局能量函数的数学形式如下所示：

其中，D代表着像素本身的代价，又叫做一元势函数，W代表不同标签产生的代价，又叫做点对势函数。这个公式就是想告诉我们：一副图像的能量由每个节点的代价以及相邻节点各种标签下所产生的代价和所决定。如果我们想要这个能量函数最小化，一方面要尽可能压低各个节点的代价，另一方面，要考虑到各个节点的相互作用。一般来说，相邻节点标签越一致，代价越小。但是，如果全部像素都是一个标签，W值的和就为0，但是这个时候自身代价一定会非常高。反之，如果D的和想达到最小，那么只要每个像素取最小代价所对应的标签即可，但是这个时候往往发现，相邻像素往往不是一个标签了，这也不是最理想的结果。标签分布结果肯定是既兼顾了自身代价，又兼顾了相邻像素之间的标签差异，二者有点此消彼长的意义在里面。有趣。。。。。

如果我们把它放在Stereo-Matching里面考虑，意义就更加具体化了，D就代表着代价聚合值，标签就是视差值。一般的非全局算法，在D求出来之后便戛然而止，直接来一个WTA，就算是把视差求出来了，如果代价聚合准确还可以，如果不准，求出来的视差图“惨不忍睹”，全局算法就是在代价聚合的基础上加了一个双保险，不单单希望代价聚合值的和最小化，还要考虑图像的区域性，桌面肯定是一个区域，一个人肯定是一个区域等等等等，它们的视差应该差不多，所以直接利用这个数学模型来达到这种需求。

下面的问题就是如何来求解这个数学模型，答案就是：用Belief-Propagation算法，关于BP算法介绍很多，本文只对几个重要的数学公式进行解释。

2. 重要公式说明

1. 消息计算公式

这个公式代表的是像素p对其邻居像素q传递的消息向量，注意，一定是一个向量。这个公式曾经让我困惑，我在想，如果m是一个向量，那么D和W都应该是向量才对，那么问题来了，向量的最小值是什么意思？W中的两个标签向量相减，分量的对应关系又是啥啊？结果我一度怀疑这里的m不是一个向量，也曾经在其他博客中给博主留言，不知道大家看到这个公式之后有没有这样的问题。后来，我下载了几份BP的源代码，通过看代码弄清楚是怎么一回事，m的确是一个向量，D和W也都可以视作向量，这个公式的正确解释是：固定fq，fp取各个视差值（假设有n个），计算D、W以及m中的对应分量，共可以得到n个值，取最小的一个作为向量m中的对应分量值，循环往复，得到的m正是一个n维向量。

这个公式有什么意义呢？它想告诉我们像素p传递给q的消息是一个置信度向量，每个分量是p对q取各个视差值的支持力度。这个支持力度不单单由p本身的代价大小决定，还要有传递给p的消息决定，这样子就会形成一个传播的态势。一般都在第一步迭代中，将所有的message都设置为零向量。在HBP算法中，只在金字塔最顶层的图像上同样设定为零向量，其他层直接利用上一层的最新消息向量做为初始向量。

论文《Stereo Matching Using Belief Propagation》中认为，立体匹配要解决的主要问题有两个，一个是低纹理区域的视差估计，一个是深度不连续区域的视差估计，也就是遮挡区域，作者认为深度不连续区域往往就是颜色不连续区域（这点多篇论文都这样提到过，SegmenTree，DoubleBP等等）。而BP有两大优点，第一个是信息不对称，第二是BP的消息值是自适应的，第一点我也没有彻底弄明白原理，关于第二点，BP往往在低纹理区域中能够将消息传递到很远，但在深度不连续区域却马上就会停止传递，所以说它是自适应的。

注：在DoubleBP论文中的上述公式有个错误，那里面写成了argmin，这是不对的，argmin代表的是求出fp的值，而这个公式想要得到的是函数值。

2. 节点p的置信度向量

这是BP算法的最后一个公式，其不参与循环，当每个消息向量稳定之后，或者迭代步骤到达之后（因为有的BP算法根本不收敛），这时根据像素q的代价值以及周围像素对其的消息向量，确定q的各个视差值的可能性大小，这就是置信度向量。这个时候直接利用WTA，就能求出像素q的视差值了。

3. 代码

这部分的代码来源于CSBP算法，专门记述HBP部分，我添加了部分注释。

short*qx_csbp::disparity(unsigned char*left,unsigned char*right)
{
	short *disp=m_data_cost;
	memset(m_message,0,sizeof(short)*(m_max_nr_message>>1));
	memset(m_data_cost,0,sizeof(short)*m_h*m_w*m_max_nr_plane_pyramid[0]*2);
	m_data_cost_selected=&(m_data_cost[(m_h*m_w*m_max_nr_plane_pyramid[0])]);

    // 分层计算，每一层中迭代计算BP
    for(int i=m_nr_scale-1;i>=0;i--)
	{
        // 计算数据项
		if(i==(m_nr_scale-1)) 
		{
            // 每一层中的节点数据项相等，不论迭代多少次，更改的是消息
			compute_data_cost_init(left,right,m_h_pyramid[i],m_w_pyramid[i],
                                   i,m_max_nr_plane_pyramid[i],m_nr_plane,
                                   m_cost_max_data_term);
		}
		else
		{
            // 不同层中的节点数据项不等，需要更新
			compute_data_cost(left,right,m_h_pyramid[i],m_w_pyramid[i],
                              i,m_max_nr_plane_pyramid[i+1],
                              m_cost_max_data_term);

            // 初始化消息，每层基于上一层消息向量进行初始化，第一层为0向量
			init_message(i);
		}

        // 同一层中，不断的更新消息向量
		for(int j=0;j<m_iteration[i];j++) 
        {
            compute_message(m_h_pyramid[i],m_w_pyramid[i],m_max_nr_plane_pyramid[i],i);
        }
	}

    // 基于原始层计算视差图
	compute_disparity(disp,0);

    return(disp);
}

// 一次迭代中的消息向量更新
void qx_csbp::compute_message(int h,int w,int nr_plane,int scale)
{
	int i,y,x,yy=h-1,xx=w-1;
	short*c0,*p0,*p1,*p2,*p3,*p4;	
	short*d0,*d1,*d2,*d3,*d4;
	int count=0;
	int yshift=w*m_nr_neighbor*nr_plane;
	int xshift=m_nr_neighbor*nr_plane;
	int yshiftd=w*nr_plane;

    // 更新整幅图像中，每个像素对四个邻居像素的消息向量
	for(i=0;i<2;i++)        // 挑选像素方式没有研究？
	{
		for(y=1;y<yy;y++)   // 图像高
		{
			int yl=y*yshift;
			int yld=y*yshiftd;
			for(x=xx-1+(y+i)%2;x>=1;x-=2) //for(x=(y+i)%2+1;x<xx;x+=2)
			{
				int xl=x*xshift;
				int xld=x*nr_plane;
				c0=&(m_data_cost_selected[yld+xld]);
				p0=&(m_message[yl+xl]);
				p1=&(m_message[yl+xl-yshift+2*nr_plane]);
				p2=&(m_message[yl+xl-xshift+3*nr_plane]);
				p3=&(m_message[yl+xl+yshift]);
				p4=&(m_message[yl+xl+xshift+nr_plane]);
				d0=&(m_selected_disparity_pyramid[yld+xld]);
				d1=&(m_selected_disparity_pyramid[yld+xld-yshiftd]);
				d2=&(m_selected_disparity_pyramid[yld+xld-nr_plane]);
				d3=&(m_selected_disparity_pyramid[yld+xld+yshiftd]);
				d4=&(m_selected_disparity_pyramid[yld+xld+nr_plane]);

                // 计算当前像素p0对四个邻居像素的消息向量
				compute_message_per_pixel(c0,p0,p1,p2,p3,p4,d0,d1,d2,d3,d4,y,x,nr_plane,scale,count);
			}
		}
	}
}

// 计算当前像素对邻居像素的消息向量
void qx_csbp::compute_message_per_pixel(short*c0,short *p0,short *p1,short *p2,short *p3,short *p4,
                                        short*d0,short*d1,short*d2, short*d3,short*d4,
                                        int y,int x,int nr_plane,int scale,int &count)
{
	short minimum[4]={30000,30000,30000,30000};
	short *p0u=p0;
	short *p0l=&(p0[nr_plane]);
	short *p0d=&(p0[nr_plane+nr_plane]);
	short *p0r=&(p0[nr_plane+nr_plane+nr_plane]);
	count++;

    // 先计算jump cost
	for(int d=0;d<nr_plane;d++)
	{
        // 计算周围三个像素对当前像素的jump cost
		p0u[d]=c0[d]+p2[d]+p3[d]+p4[d];
		p0l[d]=c0[d]+p1[d]+p3[d]+p4[d];
		p0d[d]=c0[d]+p1[d]+p2[d]+p4[d];
		p0r[d]=c0[d]+p1[d]+p2[d]+p3[d];

        // 计算最小的jump cost值
		if(p0u[d]<minimum[0]) minimum[0]=p0u[d];
		if(p0l[d]<minimum[1]) minimum[1]=p0l[d];
		if(p0d[d]<minimum[2]) minimum[2]=p0d[d];
		if(p0r[d]<minimum[3]) minimum[3]=p0r[d];
	}

    // 当前像素传递给每个邻居像素消息向量
	compute_message_per_pixel_per_neighbor(p0u,minimum[0],d0,d1,nr_plane,scale);
	compute_message_per_pixel_per_neighbor(p0l,minimum[1],d0,d2,nr_plane,scale); 
	compute_message_per_pixel_per_neighbor(p0d,minimum[2],d0,d3,nr_plane,scale);
	compute_message_per_pixel_per_neighbor(p0r,minimum[3],d0,d4,nr_plane,scale);
}

// 计算当前像素对邻居像素的消息向量
void qx_csbp::compute_message_per_pixel_per_neighbor(short *comp_func_sub,short minimum,
                                                     short *disp_left,short *disp_right,
                                                     int nr_plane,int scale)
{
    // 计算当前像素对特定邻居像素的消息向量
	for(int d=0;d<nr_plane;d++)
	{
		short cost_min=minimum+m_cost_max_discontinuity;

        // 计算fq取每个分量下，以fp为自变量的消息最小值
		for(int i=0;i<nr_plane;i++)
		{
            // m_discontinuity_cost_single_jump*abs(disp_left[i]-disp_right[d])就是点对势函数
            // abs(disp_left[i]-disp_right[d]) - V(fp - fq)
            // 没有考虑到“discontimuity preserving”
			cost_min=min(cost_min,comp_func_sub[i]+m_discontinuity_cost_single_jump*abs(disp_left[i]-disp_right[d]));
		}

        // message(fq) = min(fp)
		m_temp[d]=cost_min;
	}
	memcpy(comp_func_sub,m_temp,sizeof(short)*nr_plane);

    // 对消息向量进行中心化处理
	bpstereo_normalize(comp_func_sub,nr_plane);
}

4. 总结

总结：本文主要说说BP算法与stereo-matching之间的关系，主要对两个公式进行解释。