A Study in Vision: Loopy BP, MRF

今天读到一篇介绍BP和MRF的blog, 在stereo vision中，如果只考虑单独的像素做匹配，而忽略像素周围的pattern，就容易导致不连贯的disparity map。实际上我们希望disparity map中相邻像素有相同的值。这个时候就可以用到MRF。MRF是无向图，但是跟Bayesian不同的是它可以有cycle和loop。

蓝色的点表示观测到的变量，就是像素值。粉色的点表示隐藏变量，就是disparity。通常隐藏变量叫做label。node之间的link表示dependency，比如中间粉色的点只跟周围四个点和上面的蓝色的点有关。这个某点只跟周围点有关的假设就是Markov假设。这个假设使我们能够高效的求解隐藏变量。

如果用MRF来表达stereo vision，它的energy function就是

$energy(Y,X) = \sum\limits_{i} DataCost\left(y_i, x_i\right) \; + \sum\limits_{j = \mbox{neighbours of i}} SmoothnessCost\left(x_i,x_j\right)$

Y表示观测变量，X表示隐藏变量。i是pixel的index，j是xi相邻的node。给定一个图像Y和一些label X，这个能量方程求得了每个link的cost的和。我们的目标是找到一个label X，比如disparity map，使得这个能量方程最小化。接下来我们分开来看data cost和smoothness cost。

Datacost主要指把label xi赋值给data yi造成的cost。对于正确的匹配，datacost很低。对错误的匹配datacost就很高。常用的衡量datacost的有差值绝对值的和，SSD等。

Smoothness cost确保相邻像素有相同的label。我们需要一个函数来惩罚相邻像素有不同label的情况。常用的函数有如下几种

$f\left(n\right) = \begin{cases} 0 & \mbox{if } n = 0 \\ \lambda & \mbox{otherwise} \end{cases}$ Also known as the Potts model.

$f\left(n\right) = \lambda\times\mbox{min}\left(\left|n\right|, K\right)$ Truncated linear model.

$f\left(n\right) = \lambda\times\mbox{min}\left(n^2, K\right)$ Truncated quadratic model.

如何选取合适的函数和参数看起来就像黑魔法，大多数paper并没有提到应该怎样做。

Loopy Belief Propagation

因为图像中有很多像素，disparity value也有很多可能，所以很难找到MRF的精确解。LBP提供了一种方法来寻找近似解，类似的方法还有graph cut, ICM.不过LBP不保证convergence。

LBP是中用来传输信息的方法，当一个node收到了所有信息的时候，它就发给相邻node一个信息。下图展示了从x1传送到x2的过程。

x1首先需要从A,B,C,D接收到信息，然后才会给x2传输信息。x2不会返回给x1信息。准确来说信息的定义是 $msg_{i \rightarrow j}\left(l\right)$ ，表示从node i发送label l的信息给node j。换言之就是node i对node j属于label l的belief。这些信息只在隐藏变量之间传递。一个完整的信息包含所有可能的label。比如node i会给node j发送如下信息

hey node j，我认为你是label 0，概率是s0

hey node j，我认为你是label 1，概率是s1

。。。

Node i记载了所有关于node j的可能性。概率的计算取决于MRF。

LBP的第一步是初始化信息。因为node要等到所有相邻node都发送信息，这就变成了一个鸡生蛋蛋生鸡的问题，因为所有node都会等待其他node发送信息，实际上谁也没有发送任何东西。为了解决这个问题，我们把所有信息都初始化成一个常数，通常是0或1.

LBP主体算法是iterative的。如同其他iterative的算法，我们可以在一定循环次数后结束，或者到energy的变化小于一个阈值。在每个iteration，信息在MRF中传递。信息传递的次序是随机的。一旦这个过程结束，我们就可以根据每个node的belief计算这个node的label。

接下来我们一个个来看信息更新，初始化，和belief的步骤，和三个不同算法sum product，max product, min sum。

用于信息更新的sum product
$msg_{i \rightarrow j}\left( l \right) = \sum\limits_{l' \in \mbox{all labels}} \left[ \begin{array}{c} exp\left(-DataCost\left(y_i,l'\right)\right) exp\left(-SmoothnessCost\left(l,l'\right)\right) \times \\ \prod\limits_{k=\left( \begin{array}{c} \mbox{neighbours of i} \\ \mbox{except j} \end{array} \right)} msg_{k\rightarrow i}\left(l'\right) \end{array} \right]$
等式左边表示从node i发到node j，关于label l的信息。右边的y表示图像像素。这里我们遍历所有的label，在disparity map中共有15种。因为有加和，内积的计算，所以叫sum product。这个算法用于概率的计算，所以要用exp函数把data cost, smoothness cost，转换到在0到1之间的概率，这个概率越大越好。在中括号里面的是data cost, smoothness cost对于label l的所有信息的joint probability. 中括号外面的加和是对概率在变量l上的marginalization.

一个完整的信息是一个矢量
$msg_{i \rightarrow j}=\left[ \begin{array}{c} msg_{i \rightarrow j}\left(0\right)\\ msg_{i \rightarrow j}\left(1\right)\\ msg_{i \rightarrow j}\left(2\right)\\{..} \end{array} \right]$
所以对于每个label都要遍历所有可能，复杂度就是O（L^2）.

连续对概率做乘积的时候，很快就会接近0.为了避免这个情况，我们要把信息向量normalize
$msg_{i \rightarrow j} = \frac{msg_{i \rightarrow j}}{\sum\limits_{l}msg_{i \rightarrow j}\left(l\right)}$
进行初始化的时候，所有信息的概率都设为1.每个node的belief是所有信息的乘积。
$Belief \left(x_i = l\right) = exp\left(-DataCost\left(y_i, l\right)\right) \times \prod\limits_{k=\mbox{neighbours of i}} msg_{k \rightarrow i}\left(l\right)$
这是node i对于label l的belief。为了找到最合适的label，需要遍历所有label然后找到最高的belief。

用于信息更新的max product
sum product可以找到每个node的最佳label。但是总体来说并不一定是最优解。举例来说，假设有两个变量x，y

P(x,y)	x=0	x=1
y=0	0.5	0.4	P(y=0) = 0.9
y=1	0.1	0.3	P(y=1) = 0.4
	P(x=0) = 0.6	P(x=1) = 0.7

表格外边的是变量的marginal。如果用单独的marginal计算，我们会选择x=1, y = 0,得到P(x=1,y=0) = 0.4。但是最佳的解是p(x=0,y=0) = 0.5。我们最关心的是Joint probability。此类问题经常会在maximum a posteriori (MAP)求解中出现，因为这时我们想找到全局的最优解。max product在sum product的基础上做了一点点改变
$msg_{i \rightarrow j}\left( l \right) = \max\limits_{l' \in \mbox{all labels}} \left[ \begin{array}{c} exp\left(-DataCost\left(y_i,l'\right)\right) exp\left(-SmoothnessCost\left(l,l'\right)\right) \times \\ \prod\limits_{k=\left(\begin{array}{c} \mbox{neighours of i} \\ \mbox{except j} \end{array} \right)} msg_{k\rightarrow i}\left(l'\right) \end{array} \right]$
现在不再求和，而是计算marginal probability的最大值。

用来更新信息的min sum
和max sum相似，min sum也是计算每个node的max marginal，不过是在log space中。
$msg_{i \rightarrow j}\left( l \right) = \min\limits_{l' \in \mbox{all labels}} \left[ \begin{array}{c} DataCost\left(y_i,l'\right) + SmoothnessCost\left(l,l'\right) + \\ \sum\limits_{k=\mbox{neighours of i except j}} msg_{k\rightarrow i}\left(l'\right) \end{array} \right]$
这是个求解最小值的问题。在初始化的时候所有的数值都是0. 这时的belief是

$Belief \left(x_i = l\right) = DataCost\left(y_i, l\right) + \sum\limits_{k=neighbours of x_i} msg_{k \rightarrow i}\left(l\right)$

不过因为我们其实在找最小值，称它为cost更合适。

在这些方法中，min sum是最方便实现的，它没有exp函数，只有加和。如果用sum product的话，就要在exp里面加上scaling来避免underflow。eg. exp(-DataCost(…)*scaling) * exp(-SmoothnessCost(..)*scaling), scaling是 0 到1之间的数.

A Study in Vision

Monday, 19 October 2015

Loopy BP, MRF

No comments:

Post a Comment