置信度传播是一种基于马尔科夫随机场理论的立体稠密匹配算法，马尔科夫随机场的具体理论这里不再详述，只对置信度传播立体匹配的实现原理做一定简述。
成对的马尔科夫模型是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)，那么标准置信度算法中迭代运算步骤如下所示：
（懒得打字了，从期刊上截图下来的，若理解不了请多看原版，另外这里是非营利性使用，若有侵权还望作者见谅啊！）
具体实现过程如下：
void BeliefPropagateMatch(const cv::Mat& img1, const cv::Mat& img2,
int min_disp, int max_disp, int smooth_d,
int P1, int P2, int m_size, cv::Mat& disp)
{
if (img1.empty() || img1.type() != CV_8UC1
|| img2.empty() || img2.type() != CV_8UC1
|| img1.size() != img2.size() || m_size < 1) return;
/******************************************建立似然概率空间*****************************************/
double*** simp = new double** [img1.rows];
cv::Vec2d*** up_msg = new cv::Vec2d** [img1.rows];
cv::Vec2d*** down_msg = new cv::Vec2d ** [img1.rows];
cv::Vec2d*** left_msg = new cv::Vec2d ** [img1.rows];
cv::Vec2d*** right_msg = new cv::Vec2d ** [img1.rows];
int disp_range = max_disp - min_disp + 1;
double init_msg = 1. / disp_range;
for (int r = 0; r < img1.rows; ++r)
{
printf("建立匹配代价空间...%d%%\r", (int)(r * 100.0 / img1.rows));
simp[r] = new double* [img1.cols];
up_msg[r] = new cv::Vec2d * [img1.cols];
down_msg[r] = new cv::Vec2d * [img1.cols];
left_msg[r] = new cv::Vec2d * [img1.cols];
right_msg[r] = new cv::Vec2d * [img1.cols];
for (int cl = 0; cl < img1.cols; ++cl)
{
simp[r][cl] = new double[disp_range];
up_msg[r][cl] = new cv::Vec2d[disp_range];
down_msg[r][cl] = new cv::Vec2d[disp_range];
left_msg[r][cl] = new cv::Vec2d[disp_range];
right_msg[r][cl] = new cv::Vec2d[disp_range];
for (int d = 0; d < disp_range; ++d)
{
up_msg[r][cl][d] = { init_msg ,init_msg };
down_msg[r][cl][d] = { init_msg ,init_msg };
left_msg[r][cl][d] = { init_msg ,init_msg };
right_msg[r][cl][d] = { init_msg ,init_msg };
int cr = cl - d - min_disp;///当前视差下左像素对应对的右像素坐标
if (cl >= 0)
{
double dif = SAD(img1, cl, r, img2, cr, r, 3);//计算匹配代价
simp[r][cl][d] = exp(-dif / smooth_d); //计算似然概率
}
else
{
simp[r][cl][d] = d > 0 ? simp[r][cl][d - 1] : 1;
}
}
}
}
printf("建立似然概率空间...ok \n");
/******************************************四邻域消息迭代*****************************************/
int iteration_num = 0;///迭代次数
disp = cv::Mat(img1.size(), CV_8UC1, cv::Scalar::all(0));
double img_belief = 0;///总的置信度
double img_belief_pre = 0;
do
{
iteration_num++;
int i_cur = iteration_num % 2;
int i_pre = i_cur == 0 ? 1 : 0;
printf("置信度消息迭代第%d次\r", iteration_num);
img_belief_pre = img_belief;
img_belief = 0;
for (int r = 1; r < img1.rows - 1; ++r)//更新左、右邻域传递的消息
{
for (int cl = 1; cl < img1.cols - 1; ++cl)
{
int cr = img1.cols - 1 - cl;
double sum_pl = 0;///左领域归一化分母
double sum_pr = 0;///右领域归一化分母
for (int d = 0; d < disp_range; ++d)
{
double max_pl = 0;//左邻域当前状态下的最大消息值
double max_pr = 0;//右邻域当前状态下的最大消息值
for (int pre_d = 0; pre_d < disp_range; ++pre_d)
{
int deta_d = abs(d - pre_d);
double smooth = deta_d > 1 ? P2 : (deta_d == 0 ? 0 : P1);///平滑值
double pri_p = exp(-smooth / smooth_d); //先验概率
double tmp_pl = pri_p * simp[r][cl][d] * up_msg[r - 1][cl][pre_d].val[i_pre]
* down_msg[r + 1][cl][pre_d].val[i_pre] * left_msg[r][cl - 1][pre_d].val[i_pre]; //计算左邻域消息
double tmp_pr = pri_p * simp[r][cr][d] * up_msg[r - 1][cr][pre_d].val[i_pre]
* down_msg[r + 1][cr][pre_d].val[i_pre] * right_msg[r][cr + 1][pre_d].val[i_pre]; //计算右邻域消息
if (tmp_pl > max_pl)
{
max_pl = tmp_pl;
}
if (tmp_pr > max_pr)
{
max_pr = tmp_pr;
}
}
left_msg[r][cl][d].val[i_cur] = max_pl;//左邻域更新消息
sum_pl += max_pl;
right_msg[r][cr][d].val[i_cur] = max_pr;//右邻域更新消息
sum_pr += max_pr;
}
for (int d = 0; d < disp_range; ++d)
{
left_msg[r][cl][d].val[i_cur] /= sum_pl; //左邻域消息归一化
right_msg[r][cr][d].val[i_cur] /= sum_pr; //右邻域消息归一化
}
}
}
for (int c = 1; c < img1.cols - 1; ++c)//更新上、下邻域传递的消息
{
for (int ru = 1; ru < img1.rows - 1; ++ru)
{
int rd = img1.rows - 1 - ru;
double sum_pu = 0;///归一化分母
double sum_pd = 0;///归一化分母
for (int d = 0; d < disp_range; d++)
{
double max_pu = 0;//当前状态下的最大消息值
double max_pd = 0;//当前状态下的最大消息值
for (int pre_d = 0; pre_d < disp_range; ++pre_d)
{
int deta_d = abs(d - pre_d);
double smooth = deta_d > 1 ? P2 : (deta_d == 0 ? 0 : P1);///平滑值
double pri_p = exp(-smooth / smooth_d); //先验概率
double tmp_pu = pri_p * simp[ru][c][d] * up_msg[ru - 1][c][pre_d].val[i_pre]
* left_msg[ru][c - 1][pre_d].val[i_pre] * right_msg[ru][c + 1][pre_d].val[i_pre]; //计算消息
double tmp_pd = pri_p * simp[rd][c][d] * down_msg[rd + 1][c][pre_d].val[i_pre]
* left_msg[rd][c - 1][pre_d].val[i_pre] * right_msg[rd][c + 1][pre_d].val[i_pre]; //计算消息
if (tmp_pu > max_pu)
{
max_pu = tmp_pu;
}
if (tmp_pd > max_pd)
{
max_pd = tmp_pd;
}
}
up_msg[ru][c][d].val[i_cur] = max_pu;//上邻域更新消息
sum_pu += max_pu;
down_msg[rd][c][d].val[i_cur] = max_pd;//下邻域更新消息
sum_pd += max_pd;
}
for (int d = 0; d <disp_range; d++)
{
up_msg[ru][c][d].val[i_cur] /= sum_pu; //上邻域归一化
down_msg[rd][c][d].val[i_cur] /= sum_pd; //下邻域归一化
}
}
}
for (int i = 1; i < img1.rows - 1; i++)
{
for (int j = 1; j < img1.cols - 1; j++)
{
double max_belief = 0;
int max_bp = NULL;
for (int d = 0; d < disp_range; d++)
{
double belief = down_msg[i + 1][j][d].val[i_cur] * up_msg[i - 1][j][d].val[i_cur]
* left_msg[i][j - 1][d].val[i_cur] * right_msg[i][j + 1][d].val[i_cur] * simp[i][j][d];
if (belief > max_belief)
{
max_belief = belief;
max_bp = d;
}
}
disp.ptr<uchar>(i)[j] = (max_bp + min_disp);
img_belief = img_belief + max_belief;
}
disp.ptr<uchar>(i)[0] = disp.ptr<uchar>(i)[1];
disp.ptr<uchar>(i)[img1.cols - 1] = disp.ptr<uchar>(i)[img1.cols - 2];
}
for (int i = 0; i < img1.cols; ++i)
{
disp.ptr<uchar>(0)[i] = disp.ptr<uchar>(1)[i];
disp.ptr<uchar>(img1.rows - 1)[i] = disp.ptr<uchar>(img1.rows - 2)[i];
}
} while (abs(img_belief - img_belief_pre) > 0.01);
printf("置信度消息迭代 %d 次：ok \n", iteration_num);
for (int r = 0; r < img1.rows; ++r)
{
for (int c = 0; c < img1.cols; ++c)
{
delete[] simp[r][c];
delete[] up_msg[r][c];
delete[] down_msg[r][c];
delete[] left_msg[r][c];
delete[] right_msg[r][c];
}
delete[] simp[r];
delete[] up_msg[r];
delete[] down_msg[r];
delete[] left_msg[r];
delete[] right_msg[r];
}
delete[] simp;
delete[] up_msg;
delete[] down_msg;
delete[] left_msg;
delete[] right_msg;
}
void BeliefPropagationTest()
{
cv::Mat img1 = cv::imread("cones\\im2.png");
cv::Mat img2 = cv::imread("cones\\im6.png");
cv::imshow("img1", img1);
cv::imshow("img2", img2);
cv::Mat gray1, gray2;
cv::cvtColor(img1, gray1, cv::COLOR_BGR2GRAY);
cv::cvtColor(img2, gray2, cv::COLOR_BGR2GRAY);
cv::Mat disp;
dense_matching::BeliefPropagateMatch(gray1, gray2, 0, 100,100, 1, 4, 1, disp);
cv::Mat disp_show;
cv::normalize(disp, disp_show, 0, 255, cv::NORM_MINMAX);
cv::imshow("disp", disp_show);
cv::waitKey(0);
}代码中将左、右方向和上下方向的消息传递进行了合并，以提高迭代速度。
本人水平有限，如有错误，还望不吝指正，代码有一定删减，没有重新编译，如有错误，请自行调试，有问题请邮箱联系1299771369@qq.com。
扫一扫
875
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
