Variational-Bayes Optical Flow

Abstract

The Horn-Schunck (HS) optical flow method is widely employed to initialize many motion estimation algorithms. In this work, a variational Bayesian approach of the HS method is presented, where the motion vectors are considered to be spatially varying Student’s t-distributed unobserved random variables, i.e., the prior is a multivariate Student’s t-distribution, while the only observations available is the temporal and spatial image difference. The proposed model takes into account the residual resulting from the linearization of the brightness constancy constraint by Taylor series approximation, which is also assumed to be a spatially varying Student’s t-distributed observation noise. To infer the model variables and parameters we recur to variational inference methodology leading to an expectation-maximization (EM) framework with update equations analogous to the Horn-Schunck approach. This is accomplished in a principled probabilistic framework where all of the model parameters are estimated automatically from the data. Experimental results show the improvement obtained by the proposed model which may substitute the standard algorithm in the initialization of more sophisticated optical flow schemes.

Appendix

In what follows we present in detail the derivation of the update equations for the model variables and parameters.

In the fully Bayesian framework, the complete data likelihood, including the hidden variables and the parameters of the model, is given by:

$$\begin{aligned} p (\mathbf{d}, \mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta ) =& p ( \mathbf{d}|\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta ) p (\mathbf{u}| \tilde{\mathbf{A}}, \mathbf{b}; \theta ) \\ &{}\times p (\tilde{\mathbf{A}}; \theta )p ( \mathbf{b}; \theta ), \end{aligned}$$

(55)

where θ=[λ _noise ,λ _x ,λ _y ,μ,ν _x ,ν _y ] gathers the parameters of the model. Estimation of the model parameters could be obtained through maximization of the marginal distribution of the observations p(d;θ):

$$ \hat{\theta} = \mathop{\mathrm{arg\,max}}_{\theta}{\iiint p (\mathbf {d}, \mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta )d\mathbf{u} d \tilde{\mathbf{A}} d\mathbf{b}}. $$

(56)

However, in the present case, this marginalization is not possible, since the posterior of the latent variables given the observations $p(\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}|\mathbf{d})$ is not known explicitly and inference via the Expectation-Maximization (EM) algorithm may not be obtained. Thus, we resort to the variational methodology [5, 8] where we have to maximize a lower bound of $p (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b} )$:

$$\begin{aligned} &{ L (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta )} \\ &{\quad = \iiint q (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b} ) \log \frac{p (\mathbf{d}, \mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta )}{q (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b} )}d \mathbf{u} d\tilde{\mathbf{A}} d\mathbf{b}.} \end{aligned}$$

(57)

This involves finding approximations of the posterior distribution of the hidden variables, denoted by q(u), $q (\tilde{\mathbf{A}} )$, q(b) because there is no analytical form of the auxiliary function q for which the bound in (57) becomes equality. However, in the variational methodology, we employ the mean field approximation [8]:

$$ q (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b} ) = q (\mathbf{u} ) q (\tilde{\mathbf{A}} )q (\mathbf{b} ), $$

(58)

and (57) becomes:

$$ L (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta ) = \int _{\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}} q (\mathbf{u} ) q (\tilde{\mathbf{A}} )q (\mathbf{b} ) \log \frac{p (\mathbf{d}, \mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta )}{q (\mathbf{u} ) q (\tilde{\mathbf{A}} )q (\mathbf{b} )}. $$

(59)

In our case, in the E-step of the variational algorithm (VE-step), optimization of the functional $L (\mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta )$ is performed with respect to the auxiliary functions. Following the variational inference framework, the distributions q(u _k ), k∈{x, y}, are Normal:

$$ q (\mathbf{u} ) = \mathcal{N}\left (\left [ \begin{array}{c} \mathbf{m}_x \\ \mathbf{m}_y \\ \end{array} \right ], \left [ \begin{array}{c@{\quad}c} \mathbf{R}_x & \mathbf{0} \\ \mathbf{0} & \mathbf{R}_y \\ \end{array} \right ] \right ), $$

(60)

yielding

$$ q (\mathbf{u}_{x} ) = \mathcal {N} ( \mathbf{m}_x, \mathbf{R}_x ), $$

(61)

and

$$ q (\mathbf{u}_{y} ) = \mathcal {N} ( \mathbf{m}_y, \mathbf{R}_y ). $$

(62)

Therefore, this bound is actually a function of the parameters R _k and m _k , k∈{x,y} and a functional with respect to the auxiliary functions q(a _k ),q(b). Using (58), the variational bound in our problem becomes:

$$\begin{aligned} &{ L \bigl(q(\mathbf{u}_{x}), q(\mathbf{u}_{y}), q( \mathbf{a}_x), q(\mathbf{a}_y), q(\mathbf{b}), \theta_1, \theta_2 \bigr) } \\ &{\quad =\iiint \biggl(\prod_{k\in\{x, y\}} q( \mathbf{u}_{k}; \theta_1) q(\mathbf{a}_k) \biggr) q(\mathbf{b}) } \\ &{\qquad{}\times\log p(\mathbf{d}, \mathbf{u}, \tilde{\mathbf{A}}, \mathbf{b}; \theta _2) d\mathbf{u}d\tilde{\mathbf{A}}d\mathbf{b} } \\ &{\qquad{}-\iiint \biggl(\prod_{k\in\{x, y\}}q( \mathbf{u}_{k}; \theta_1) q(\mathbf{a}_k) \biggr) q(\mathbf{b}) } \\ &{\qquad{}\times\log \biggl( \biggl( \prod_{k\in\{x, y\}} p(\mathbf{u}_{k}; \theta_1) q(\mathbf{a}_k) \biggr) q(\mathbf{b}) \biggr) d\mathbf{u}d\tilde{\mathbf{A}}d\mathbf{b}} \end{aligned}$$

(63)

where we have separated the parameters into two sets:

$$ \theta_1 = \{\mathbf{R}_x, \mathbf{R}_x, \mathbf{m}_x, \mathbf{m}_y \}, $$

(64)

and

$$ \theta_2 = \{\mathbf{a}_x, \mathbf{a}_y, \mathbf{b}, \lambda_x, \lambda_y, \nu_x, \nu_y \}. $$

(65)

Thus, in the VE-step of the variational EM algorithm the bound must be optimized with respect to R _k , m _k , q(a _k ) and q(b).

Taking the derivative of (63) with respect to m _k , R _k , q(α _k ) and q(b) and setting the result equal to zero, we obtain the following update equations:

$$ \mathbf{m}_x^{(t+1)} = \lambda_\mathit{noise}^{(t)} \mathbf{R}_x^{(t)}\hat{\mathbf {B}}^{(t)} \mathbf{G}_x \bigl(\mathbf{d}-\mathbf{G}_y \mathbf{u}_{y}^{(t)} \bigr), $$

(66)

and

$$ \mathbf{m}_y^{(t+1)} = \lambda_\mathit{noise}^{(t)} \mathbf{R}_y^{(t)}\hat{\mathbf {B}}^{(t)} \mathbf{G}_y \bigl(\mathbf{d}-\mathbf{G}_x \mathbf{u}_{x}^{(t)} \bigr), $$

(67)

where

$$ \mathbf{R}_x^{(t+1)} = \bigl( \lambda_\mathit{noise}^{(t)}\mathbf{G}_x^T\hat{ \mathbf{B}}^{(t)}\mathbf{G}_x + \lambda_x^{(t)} \mathbf{Q}^T \hat{\mathbf{A}}_x^{(t)} \mathbf{Q} \bigr)^{-1}, $$

(68)

and

$$ \mathbf{R}_y^{(t+1)} = \bigl( \lambda_\mathit{noise}^{(t)}\mathbf{G}_y^T\hat{ \mathbf{B}}^{(t)}\mathbf{G}_y + \lambda_y^{(t)} \mathbf{Q}^T \hat{\mathbf{A}}_y^{(t)} \mathbf{Q} \bigr)^{-1}. $$

(69)

Notice that the final estimates for u _x , u _y are m _x and m _y , in (37) and (38), respectively.

After some manipulation, we obtain the update equations for the model parameters which maximize (63) with respect to q(a _k ), q(b). The form of all q approximating-to-the-posterior functions will remain the same as the corresponding prior (due to the conjugate priors we employ) namely q(a _k ), q(b) which approximate p(a _k |u _k ,λ _k ,C _k ;ν _k ), p(b|u,λ _noise ,F;μ) will follow Gamma distributions, ∀i=1,…,N,∀k∈{x,y}:

$$\begin{aligned} &{ q^{(t+1)}\bigl(\boldsymbol{\alpha}_k(i)\bigr)} \\ &{\quad= \mathrm{Gamma} \biggl(\frac{\nu_k^{(t)}}{2} + \frac{1}{2},} \\ &{\phantom{\quad= \mathrm{Gamma} \biggl(}\frac{\nu_k^{(t)}}{2} + \frac{1}{2}\lambda_k^{(t)} \bigl( \bigl[\mathbf{Q}\mathbf{u}_{k}^{(t)} \bigr]_i^{2} + \mathbf{C}_k^{(t)} (i, i ) \bigr) \biggr),} \end{aligned}$$

(70)

and

$$\begin{aligned} &{q^{(t+1)}\bigl(\mathbf{b}(i)\bigr) } \\ &{\quad= \mathrm{Gamma} \biggl( \frac{\mu^{(t)}}{2} + \frac{1}{2},} \\ &{\phantom{\quad= \mathrm{Gamma} \biggl(} \frac{\mu^{(t)}}{2} + \frac{1}{2} \lambda_\mathit{noise}^{(t)} \bigl( \bigl[\mathbf{G}\mathbf{u}^{(t)} - \mathbf{d} \bigr]_i^2 + \mathbf{F}^{(t)} (i, i ) \bigr) \biggr),} \end{aligned}$$

(71)

where the N×N matrix

$$ \mathbf{C}_k^{(t)} = \mathbf{Q} \mathbf{R}_k^{(t)} \mathbf{Q}^T, $$

(72)

the N×N matrix

$$ \mathbf{F}^{(t)} = \mathbf{G}_x \mathbf{R}_x^{(t)}\mathbf{G}_x^T+ \mathbf{G}_y \mathbf{R}_y^{(t)} \mathbf{G}_y^T, $$

(73)

and [⋅] _i denotes the i-th element of the vector inside the brackets.

The size of matrices R _x , R _y and consequently C _x , C _y and F makes their direct calculation prohibitive. In order to overcome this difficulty, we employ the iterative Lanczos method [29] for their calculation. For matrices C _x , C _y and F only the diagonal elements are needed in (70) and (71) and they are obtained as a byproduct of the Lanczos method.

Note that as we can see from (66) and (67), there is a dependency between u _x and u _y , as it is the case in the standard Horn-Schunck method.

Notice also that since each q ^(t+1)(α _k (i)) is a Gamma pdf, it is easy to derive its expected value:

$$ \bigl\langle \boldsymbol{\alpha}_k(i) \bigr\rangle _{q^{(t+1)}(\boldsymbol {\alpha}_k(i))} = \frac{\nu_k^{(t)} + 1}{\nu_k^{(t)} + \lambda_k^{(t)} ( [\mathbf{Q}\mathbf{u}_{k}^{(t)} ]_i^{2} + \mathbf{C}_k^{(t)}(i, i) )}, $$

(74)

and the same stands for the expected value of b(i):

$$ \bigl\langle \mathbf{b}(i) \bigr\rangle _{q^{(t+1)}(\mathbf{b}(i))} = \frac{\mu^{(t)} + 1}{\mu^{(t)} + \lambda_\mathit{noise}^{(t)} ( [\mathbf{G}\mathbf{u}^{(t)} - \mathbf{d} ]_i^{2} + \mathbf{F}^{(t)}(i, i) )}, $$

(75)

where 〈.〉 _q(.) denotes the expectation with respect to an arbitrary distribution q(⋅). These estimates are used in (66), (67), (68) and (69), where $\hat{\mathbf{A}}_{k}^{(t)}$ and $\hat{\mathbf{B}}^{(t)}$ are diagonal matrices with elements:

$$\hat{\mathbf{A}}_k^{(t)}(i, i) = \bigl\langle \boldsymbol{ \alpha}_k(i) \bigr\rangle _{q^{(t)}(\boldsymbol{\alpha}_k(i))}, $$

and

$$\hat{\mathbf{B}}^{(t)}(i, i) = \bigl\langle \mathbf{b}(i) \bigr\rangle _{q^{(t)}(\mathbf{b}(i))}, $$

for i=1,…,N.

At the variational M-step, the bound is maximized with respect to the model parameters:

$$\begin{aligned} &{\theta_2^{(t+1)} = \mathop{\mathrm{arg\,max}}_{\theta_2} L \bigl(q^{(t+1)} (\mathbf{u}_{k} ), q^{(t+1)} (\hat{ \mathbf{A}}_k ),} \\ &{\phantom{\theta_2^{(t+1)} = \mathop{\mathrm{arg\,max}}_{\theta_2} L \bigl(} q^{(t+1)} (\hat{\mathbf{B}} ), \theta_1^{(t+1)}, \theta_2 \bigr),} \end{aligned}$$

(76)

where

$$\begin{aligned} &{L \bigl(q^{(t+1)} (\mathbf{u}_{k} ), q^{(t+1)} (\hat { \mathbf{A}}_k ), q^{(t+1)} (\hat{\mathbf{B}} ), \theta _1^{(t+1)}, \theta_2 \bigr) }\\ &{\quad{}\propto\bigl\langle \log p (\mathbf{d}, \mathbf{u}, \hat{\mathbf{A}}_k, \hat{\mathbf{B}}; \theta_2 ) \bigr\rangle _{q (\mathbf{u}_{k}; \theta_1^{(t+1)} ), q^{(t+1)} (\hat{\mathbf{A}}_k ), q^{(t+1)} (\hat{\mathbf{B}} )}} \end{aligned}$$

is calculated using the results from (66)–(69).

The update for λ _noise is obtained after taking the derivative of $L (q^{(t+1)} (\mathbf{u}_{k} ), q^{(t+1)} (\hat{\mathbf{A}}_{k} ), q^{(t+1)} (\hat{\mathbf{B}} ), \theta_{1}^{(t+1)}, \theta_{2} )$ in (63) with respect to it and setting it to zero:

$$\begin{aligned} \lambda_\mathit{noise}^{(t+1)} = \frac{N}{\sum_{i=1}^N \mathbf{b}^{(t+1)}(i) ( [\mathbf{G}\mathbf{u}^{(t+1)} - \mathbf{d} ]^2_i + \mathbf{F}^{(t+1)}(i, i) )}. \end{aligned}$$

(77)

By the same means we obtain the estimates for λ _x and λ _y :

$$ \lambda_k^{(t+1)} = \frac{N}{\sum_{i=1}^N \boldsymbol{\alpha}_k^{(t+1)}(i) ( [\mathbf{Q}\mathbf{u}_{k}^{(t+1)} ]^2_i + \mathbf{C}_k^{(t+1)}(i, i) )}, $$

(78)

with k∈{x,y}.

The degrees of freedom parameters ν _k of the Student’s t-distributions are also computed accordingly through the roots of the following equation:

$$\begin{aligned} &{ \frac{1}{N} \Biggl(\sum_{i=1}^N \log\bigl\langle \boldsymbol{\alpha}_k(i) \bigr\rangle _{q^{(t+1)} (\mathbf{A}_k )} - \sum_{i=1}^N \bigl\langle \boldsymbol{ \alpha}_k(i) \bigr\rangle _{q^{(t+1)} (\mathbf{A}_k )} \Biggr)} \\ &{\quad{} + \digamma \biggl( \frac{\nu_k^{(t)}}{2} + \frac{1}{2} \biggr) - \log \biggl(\frac{\nu_k^{(t)}}{2} + \frac{1}{2} \biggr)} \\ &{\quad{} - \digamma \biggl(\frac{\nu_k}{2} \biggr) + \log \biggl(\frac{\nu_k}{2} \biggr) + 1= 0,} \end{aligned}$$

(79)

for ν _k , k∈{x,y}, where Ϝ(x) is the digamma function (derivative of the logarithm of the Gamma function) and $\nu_{k}^{(t)}$ is the value of ν _k at the previous iteration.

Finally, by the same procedure we obtain estimates for the parameter μ of the noise distribution:

$$\begin{aligned} &{ \frac{1}{N} \Biggl(\sum_{i=1}^N \log\bigl\langle \mathbf{b}(i) \bigr\rangle _{q^{(t+1)} (\mathbf{b}(i) )} - \sum _{i=1}^N \bigl\langle \mathbf{b}(i) \bigr\rangle _{q^{(t+1)} (\mathbf{b}(i) )} \Biggr) } \\ &{\quad{}+ \digamma \biggl( \frac{\mu^{(t)}}{2} + \frac{1}{2} \biggr) - \log \biggl(\frac{\mu^{(t)}}{2} + \frac{1}{2} \biggr) - \digamma \biggl(\frac{\mu}{2} \biggr) } \\ &{\quad{}+ \log \biggl(\frac{\mu}{2} \biggr) + 1 = 0. } \end{aligned}$$

(80)

In our implementation Eqs. (79) and (80) are solved by the bisection method, as also proposed in [23].