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之间的关系,主要对两个公式进行解释。
Stereo Matching文献笔记之(六):浅谈置信度传播算法(Belief-Propagation)在立体匹配中的应用~相关推荐
- 浅谈标签传播算法LPA
研究生期间第一次对相关内容做了一个汇报,查找了大量文献,发现很多的介绍对于新手来说都看不懂,这里采用最简单的方法来浅谈一下,如有错误,欢迎指正. 标签传播算法是一种基于图的半监督学习方法,其 ...
- 浅谈标签传播算法:LPA
标签传播算法:LPA 1.半监督学习 让学习器不依赖外界交互,自动的利用未标记样本来提升学习性能,这种就是半监督学习,主要用来处理现实中有标记数据少.未标记数据多的问题,要利用未标记的数据,必须要做一 ...
- 浅谈Base64编码算法
一.什么是编码解码 编码:利用特定的算法,对原始内容进行处理,生成运算后的内容,形成另一种数据的表现形式,可以根据算法,再还原回来,这种操作称之为编码. 解码:利用编码使用的算法的逆运算,对经过编码的 ...
- 浅谈列控系统的阶梯式分级速度控制中专有名词的表述问题
浅谈列控系统的阶梯式分级速度控制中专有名词的表述问题 一.问题提出 个人在学习列车运行控制系统的过程中发现,不同的文献中,在介绍和论述列车运行控制系统的阶梯式分级速度控制方式时,总会提到下面专有名词中 ...
- 用matlab2012b计算自控原理的稳态误差,浅谈用终值定理计算自控原理中的稳态误差...
·115· 哈尔滨职业技术学院学报 2013年第2 期 Journal of Harbin Vocational & Technical College 在自动控制原理中,控制系统的稳态误差是 ...
- 样条线怎么挤出平面_浅谈“沿样条线挤出”在多边形建模中的应用
龙源期刊网 http://www.qikan.com.cn 浅谈 " 沿样条线挤出 " 在多边形建模中的应用 作者:谢巍 来源:<硅谷> 2011 年第 04 期 摘要 ...
- 计算机在财务核算中的应用,浅谈计算机在财务核算和财务管理工作中的辅助应用...
浅谈计算机在财务核算和财务管理工作中的辅助应用 计算机在财务管理中的应用日益广泛,已成为企业财务管理的必要手段.计算机的应用改善了企业财务管理环境,提高了财 (本文共1页) 阅读全文>> ...
- 浅谈流处理算法 (1) – 蓄水池采样
转载自 浅谈流处理算法 (1) – 蓄水池采样 前言 现如今,"大数据 "已经不是什么新概念,"一千个人眼中有一千个大数据".社交网络,智能穿戴设备,智能家居 ...
- 微型计算机在机械设计中的应用,浅谈计算机技术在机械设计制造及自动化中的应用.docx...
浅谈计算机技术在机械设计制造及自动化中的应用 当前科学技术与机械制造与自动化技术相互融合,将多种学科中的复合型技术加以整合,形成综合性的机械设计制造自动化学科.作为机械制造的核心内容,自动化在人们的生 ...
最新文章
- oozie调度中的重试和手工rerun一个workflow
- 解决MyBatis中 Could not set property ~ o f ~异常
- WinCE启动失败的原因与解决办法分析
- 塞尔达传说gba_《塞尔达传说缩小帽》是系列一年级生?,回忆众多玩友的启蒙之作...
- 自带的jvm监控不准_如何实时监控 Flink 集群和作业?
- Spread for Windows Forms高级主题(6)---数据绑定管理
- 天猫苛费猛如虎,天猫抽检潜“坑爹”
- msfconsole启动失败并报错`not_after=‘: bignum too big to convert into `long‘的解决方法
- BZOJ 1717: [Usaco2006 Dec]Milk Patterns 产奶的模式( 二分答案 + 后缀数组 )
- c#使用正则表达式获取TR中的多个TD_Linux之正则表达式
- keepalived做nginx的高可用,企业版简单介绍。
- 自定义 feign 调用实现 hystrix 超时、异常熔断
- Tomcat 设置系统默认文件编码
- 储能系统下垂控制,蓄电池通过双向dc/dc变换器并联负载,变换器输出电流按虚拟电阻比例分配
- 推荐一款免费,不限流量的内网穿透软件
- ESP8266的AT指令集(基础 Wi-Fi)
- 华为的PBC个人绩效评价模板
- 小心 transmittable-thread-local 的这个坑
- x4提示你的产品已经被禁用_win10系统注册表已被管理员禁用的解决方法
- vue2和vue3关闭语法检查