密集匹配之置信度传播BP

       置信度传播是一种基于马尔科夫随机场理论的立体稠密匹配算法，马尔科夫随机场的具体理论这里不再详述，只对置信度传播立体匹配的实现原理做一定简述。

       成对的马尔科夫模型是BP的基础，成对的含义就是包含显式节点和隐含节点。假设我们观察到像素yi的一些信息，需要据此推断隐含场景xi的信息，可以假设xi与yi之间存在一些统计依赖关系，称之为似然函数。设，N为邻域系统，则可以用势函数来表示相邻节点之间的势能量。这样我们可以得到隐含场景xi和一个显式图像yi的联合概率：

其中Z是归一化常数。

四邻域的马尔科夫随机场如下图所示：

图中灰色点代表隐含节点，黑色点表示显示节点，实线表示似然函数，虚线表示势函数。

        置信度传播立体匹配使用四邻域马尔科夫随机场对图像建模，图像上每一个像素都可以作为随机场的节点，通过节点之间的概率进行迭代推理，能获得满足最大后验概率分布的隐含结果。

        置信度传播算法中，随机场网络中某一节点的概率求解过程是一种节点向节点传递消息的过程：假设节点xi传递给节点xj的消息为mij(xi,xj)，描述了两个节点之间的兼容关系，然后令为已知节点yi传递给节点xj的消息，最后定义bj(xj)为节点xj的置信度，即其联合概率分布函数。

       将变量mj(xj,yj)，mij(xi,xj)简化为mij(xj)和mj(xj)，那么标准置信度算法中迭代运算步骤如下所示：

（懒得打字了，从期刊上截图下来的，若理解不了请多看原版，另外这里是非营利性使用，若有侵权还望作者见谅啊！）

具体实现过程如下：

///lImg、rImg分别为左右影像，MinDeta、MaxDeta分别为最小最大视差值，SmoothD为平滑控制参数，P1、P2为plots模型参数void BP::Matcher(Mat lImg, Mat rImg,int MinDeta,int MaxDeta,int SmoothD,int P1,int P2){if (lImg.empty() || rImg.empty() || (lImg.type() != CV_8UC3&&lImg.type() != CV_8UC1) || (rImg.type() != CV_8UC3&&rImg.type() != CV_8UC1)) return;Mat lGrayImg = lImg.type == CV_8UC3 ? cvtColor(lImg, lGrayImg, COLOR_BGR2GRAY) : lImg;Mat rGrayImg = rImg.type == CV_8UC3 ? cvtColor(rImg, rGrayImg, COLOR_BGR2GRAY) : rImg;/******************************************建立似然概率空间*****************************************/float ***simp = new float**[lGrayImg.rows];float ***upmsg = new float**[lGrayImg.rows];float ***downmsg = new float**[lGrayImg.rows];float ***leftmsg = new float**[lGrayImg.rows];float ***rightmsg = new float**[lGrayImg.rows];for (int i = 0; i<lGrayImg.rows; i++){printf("建立匹配代价空间...%d%%\r", (int)(i*100.0 / lGrayImg.rows));simp[i] = new float*[lGrayImg.cols];upmsg[i] = new float*[lGrayImg.cols];downmsg[i] = new float*[lGrayImg.cols];leftmsg[i] = new float*[lGrayImg.cols];rightmsg[i] = new float*[lGrayImg.cols];for (int j = 0; j<lGrayImg.cols; j++){int deta = MaxDeta - MinDeta + 1;simp[i][j] = new float[deta];upmsg[i][j] = new float[deta];downmsg[i][j] = new float[deta];leftmsg[i][j] = new float[deta];rightmsg[i][j] = new float[deta];for (int d = MinDeta; d <= MaxDeta; d++){int dd = d - MinDeta;///标号upmsg[i][j][dd] = 1.0f / deta;downmsg[i][j][dd] = 1.0f / deta;leftmsg[i][j][dd] = 1.0f / deta;rightmsg[i][j][dd] = 1.0f / deta;int rc = j - d - MinDeta;///当前视差下左像素对应对的右像素坐标if (rc >= 0){CMatchCost mc;float dif =SAD(lGrayImg, rGrayImg, j, i, rc, i, 0);//////计算匹配代价simp[i][j][dd] = exp(-dif / SmoothD);////计算似然概率}else{simp[i][j][dd] = d>0?simp[i][j][dd - 1]:1;}}}}printf("建立似然概率空间...ok    \n");/******************************************四邻域消息迭代*****************************************/int iteNum = 0;///迭代次数Mat dispMap(Size(lGrayImg.cols, lGrayImg.rows), CV_8UC1, Scalar::all(0));float Imgbelief = 0;///总的置信度float Imgbeliefpre = 0;do{iteNum++;Imgbeliefpre = Imgbelief;Imgbelief = 0;printf("置信度消息迭代第%d次\r", iteNum);for (int r = 1; r<lGrayImg.rows - 1; r++)////更新左、右邻域传递的消息{for (int cl = 1; cl<lGrayImg.cols - 1; cl++){int cr = lGrayImg.cols - 1 - cl;float sumpl = 0;///左领域归一化分母float sumpr = 0;///右领域归一化分母for (int d = MinDeta; d <= MaxDeta; d++){int dd1 = d - MinDeta;///当前阶段标号标号float maxpl = 0;//////左邻域当前状态下的最大消息值float maxpr = 0;//////右邻域当前状态下的最大消息值for (int pred = MinDeta; pred <= MaxDeta; pred++){int dd2 = pred - MinDeta;///上一阶段标号标号float smooth = abs(d - pred)>1 ? P2 : (abs(d - pred) == 0 ? 0 : P1);///平滑值float prip = exp(-smooth / SmoothD);////先验概率float tmpPl = prip*simp[r][cl][dd1] * upmsg[r - 1][cl][dd2] * downmsg[r + 1][cl][dd2] * leftmsg[r][cl - 1][dd2];////计算左邻域消息float tmpPr = prip*simp[r][cr][dd1] * upmsg[r - 1][cr][dd2] * downmsg[r + 1][cr][dd2] * rightmsg[r][cr + 1][dd2];////计算右邻域消息if (tmpPl>maxpl){maxpl = tmpPl;}if (tmpPr>maxpr){maxpr = tmpPr;}}leftmsg[r][cl][dd1] = maxpl;//左邻域更新消息sumpl = sumpl + maxpl;rightmsg[r][cr][dd1] = maxpr;//右邻域更新消息sumpr = sumpr + maxpr;}for (int d = MinDeta; d <= MaxDeta; d++){int dd = d - MinDeta;///标号leftmsg[r][cl][dd] = leftmsg[r][cl][dd] / sumpl;////左邻域消息归一化rightmsg[r][cr][dd] = rightmsg[r][cr][dd] / sumpr;////右邻域消息归一化}}}for (int c = 1; c<lGrayImg.cols - 1; c++)////更新上、下邻域传递的消息{for (int ru = 1; ru<lGrayImg.rows - 1; ru++){int rd = lGrayImg.rows - 1 - ru;float sumpu = 0;///归一化分母float sumpd = 0;///归一化分母for (int d = MinDeta; d <= MaxDeta; d++){int dd1 = d - MinDeta;float maxpu = 0;//////当前状态下的最大消息值float maxpd = 0;//////当前状态下的最大消息值for (int pred = MinDeta; pred <= MaxDeta; pred++){int dd2 = pred - MinDeta;float smooth = abs(d - pred)>1 ? P2 : (abs(d - pred) == 0 ? 0 : P1);///平滑值float prip = exp(-smooth / SmoothD);////先验概率float tmpPu = prip*simp[ru][c][dd1] * upmsg[ru - 1][c][dd2] * leftmsg[ru][c - 1][dd2] * rightmsg[ru][c + 1][dd2];////计算消息float tmpPd = prip*simp[rd][c][dd1] * downmsg[rd + 1][c][dd2] * leftmsg[rd][c - 1][dd2] * rightmsg[rd][c + 1][dd2];////计算消息if (tmpPu>maxpu){maxpu = tmpPu;}if (tmpPd>maxpd){maxpd = tmpPd;}}upmsg[ru][c][dd1] = maxpu;//上邻域更新消息sumpu = sumpu + maxpu;downmsg[rd][c][dd1] = maxpd;//下邻域更新消息sumpd = sumpd + maxpd;}for (int d = MinDeta; d <= MaxDeta; d++){int dd = d - MinDeta;upmsg[ru][c][dd] = upmsg[ru][c][dd] / sumpu;////上邻域归一化downmsg[rd][c][dd] = downmsg[rd][c][dd] / sumpd;////下邻域归一化}}}for (int i = 1; i<lGrayImg.rows - 1; i++){for (int j = 1; j<lGrayImg.cols - 1; j++){float maxbelief = 0;int maxbp = NULL;for (int d = 0; d<MaxDeta - MinDeta + 1; d++){float belief = downmsg[i + 1][j][d] * upmsg[i - 1][j][d] * leftmsg[i][j - 1][d] * rightmsg[i][j + 1][d] * simp[i][j][d];if (belief>maxbelief){maxbelief = belief;maxbp = d;}}dispMap.ptr<uchar>(i)[j] = (maxbp + MinDeta);Imgbelief = Imgbelief + maxbelief;}}} while (abs(Imgbelief - Imgbeliefpre)>0.01);printf("置信度消息迭代：ok    \n");cout << "迭代了" << iteNum << "次！" << endl;/***********************************计算每点的置信度获取视差图*****************************/}

代码中将左、右方向和上下方向的消息传递进行了合并，以提高迭代速度。

本人水平有限，如有错误，还望不吝指正，代码有一定删减，没有重新编译，如有错误，请自行调试，有问题请邮箱联系。