FALDOI: A New Minimization Strategy for Large Displacement Variational Optical Flow

Abstract

We propose a large displacement optical flow method that introduces a new strategy to compute a good local minimum of any optical flow energy functional. The method requires a given set of discrete matches, which can be extremely sparse, and an energy functional which locally guides the interpolation from those matches. In particular, the matches are used to guide a structured coordinate descent of the energy functional around these keypoints. It results in a two-step minimization method at the finest scale which is very robust to the inevitable outliers of the sparse matcher and able to capture large displacements of small objects. Its benefits over other variational methods that also rely on a set of sparse matches are its robustness against very few matches, high levels of noise, and outliers. We validate our proposal using several optical flow variational models. The results consistently outperform the coarse-to-fine approaches and achieve good qualitative and quantitative performance on the standard optical flow benchmarks.

Appendices

Appendix 1: Minimizing the Energy

The numerical minimization algorithm for the general energy (1) is obtained in this paper by decoupling both terms. We linearize the image $I_{t+1}$ near a given optical flow $\mathbf {u}_0=(u_{0,1},u_{0,2})$ and make the following approximation $I_{t+1}(\mathbf {x}+\mathbf {u}(x)) \approx I^{lin}_{t+1}(\mathbf {x}+\mathbf {u}(\mathbf {x}))$, where

$$\begin{aligned} I^{lin}_{t+1}(\mathbf {x}+\mathbf {u}(\mathbf {x}))= & {} I_{t+1}(\mathbf {x}+\mathbf {u}_0(\mathbf {x}))\\&+\,I^x_{t+1}(\mathbf {x}+\mathbf {u}_0(\mathbf {x})) (u_1 - u_{0,1})(\mathbf {x})\\&+\, I^x_{t+1}(\mathbf {x}+\mathbf {u}_0(\mathbf {x})) (u_2 - u_{0,2})(\mathbf {x}), \end{aligned}$$

and $I^x_{t+1}$, $I^y_{t+1}$ denote the partial derivatives of $I_{t+1}$ with respect to x and y, respectively. Let us recall that the two data terms $E_D$ that we have considered in Sect. 4.1 depend on $I_t(\mathbf {x})$ and $I_{t+1}(\mathbf {x}+\mathbf {u}(\mathbf {x}))$; we will denote as $E_{D,lin}$ the same data term but depending on $I_t(\mathbf {x})$ and $I^{lin}_{t+1}(\mathbf {x}+\mathbf {u}(\mathbf {x}))$. In order to decouple the fidelity term $E_{D,lin} (u)$ and the regularization term $E_{R} (\mathbf {u})$ in (1), we introduce an auxiliary variable $\mathbf {v}$ representing the optical flow and we penalize its deviation from $\mathbf {u}$. Thus, the energy to minimize is

$$\begin{aligned} E(\mathbf {u},\mathbf {v}) = E_{D,lin} (\mathbf {v}) + \beta E_{R} (\mathbf {u}) + \frac{1}{2\theta }\int _{{\varOmega }}\Vert \mathbf {u}-\mathbf {v}\Vert ^2, \end{aligned}$$

(6)

depending on the two variables $\mathbf {u}, \mathbf {v}$, where $\theta >0$. The decoupled energy (6) can be minimized by an alternating minimization procedure, alternatively fixing one variable and minimizing with respect to the other one. Section 4.1 presents the different possibilities for the energy.

• 1. For $\mathbf {v}$ fixed, let us consider each of the two different regularization terms, $E^1_R(\mathbf {u})$ and $E^2_R(\mathbf {u})$, presented in Sect. 4.1.

  • 1.1. In the case of $E^1_R(\mathbf {u})$, we reformulate the problem as a min-max problem incorporating the dual variables. Then, our minimization problem can be solved as a saddle-point problem. Following the notation of Osher et al. [22], for $\mathbf {v}=(v_1,v_2)$ fixed, we solve

    $$\begin{aligned}&\int _{\varOmega }\int _{\varOmega }\omega (\mathbf {x},\mathbf {y}) (u_i(\mathbf {x}) - u_i(\mathbf {y})) p(\mathbf {x},\mathbf {y}) {\text {d}}\mathbf {y}{\text {d}}\mathbf {x}\nonumber \\&\quad + \frac{1}{2\theta }\int _{{\varOmega }}\left( u_i -v_i\right) ^{2} {\text {d}}\mathbf {x}, \end{aligned}$$

    (7)

    for $i=1,2$, and p is the dual variable defined on ${\varOmega }\times {\varOmega }$. Let us explain it in detail. First, it is necessary to extend the notion of derivatives to a non-local framework. The non-local derivative can be written as

    $$\begin{aligned} \partial _{y}{u_i(\mathbf {x})} = \frac{u_i(\mathbf {x}) - u_i(\mathbf {y})}{d(\mathbf {x},\mathbf {y})} \end{aligned}$$

    (8)

    where $d(\mathbf {x},\mathbf {y})$ is a positive measure between two points $\mathbf {x},\mathbf {y}$. By taking $d(\mathbf {x},\mathbf {y})$ such that $w(\mathbf {x},\mathbf {y}) = d(\mathbf {x},\mathbf {y})^{-2}$, the non-local gradient $\nabla _{w}u_i(\mathbf {x},\mathbf {y})$ is defined as the vector of all partial derivatives:

    $$\begin{aligned} \nabla _{w}u_i(\mathbf {x},\mathbf {y}) = (u_{i}(\mathbf {x}) - u_{i}(\mathbf {y}))\sqrt{w(\mathbf {x},\mathbf {y})} \quad \mathbf {x},\mathbf {y}\in {\varOmega }.\nonumber \\ \end{aligned}$$

    (9)

    Now, by writing ${\mathbf {p}} := p(\mathbf {x},\mathbf {y})$ for $(\mathbf {x},\mathbf {y})\in {\varOmega }\times {\varOmega }$, the non-local divergence div$_{w}{\mathbf {p}}(\mathbf {x})$ is defined as the adjoint of the non-local gradient:

    $$\begin{aligned} \text {div}_{w}{\mathbf {p}}(\mathbf {x}) = \int _{{\varOmega }} (p(\mathbf {x},\mathbf {y}) - p(\mathbf {y},\mathbf {x}))\sqrt{w(\mathbf {x},\mathbf {y})} {\text {d}}\mathbf {y}.\nonumber \\ \end{aligned}$$

    (10)

Proposition 1

The solution of (7) is given by the following iterative scheme

$$\begin{aligned} p(\mathbf {x},\mathbf {y})^{n+1}= & {} \frac{p(\mathbf {x},\mathbf {y})^{n} + \tau (\overline{u}_{i}^{n}(\mathbf {x}) - \overline{u}_{i}^{n}(\mathbf {y})\sqrt{w(\mathbf {x},\mathbf {y})})}{1 + \tau |\nabla _{w}u_{i}(\mathbf {x},\mathbf {y})|} \end{aligned}$$

(11)

$$\begin{aligned} u_{i}^{n+1}(\mathbf {x})= & {} u_{i}^{n}(\mathbf {x}) - \sigma \left( \frac{(u_{i}^{n}(\mathbf {x}) - v_{i}(\mathbf {x}))}{\theta } - \text {div}_{w} {\mathbf {p}}(\mathbf {x})\right) \nonumber \\ \end{aligned}$$

(12)

$$\begin{aligned} \overline{u}_{i}^{n+1}(\mathbf {x})= & {} 2u_{i}^{n+1}(\mathbf {x}) -u_{i}^{n}(\mathbf {x}) \end{aligned}$$

(13)

where ${u}_{{i}}$ is the primal variable and ${\mathbf {p}}$ is the dual variable.

• 1.2. In the case of $E^2_R(u)$, we use the primal-dual algorithm that Chambolle proposed to minimize the ROF model [13] and which is based on a dual formulation of the TV. Then, our minimization problem can be solved as a saddle-point problem. For $\mathbf {v}$ fixed, we solve

  $$\begin{aligned} \quad \min _\mathbf {u}\max _\mathbf {\xi } \int _{{\varOmega }} \langle D \mathbf {u},\xi \rangle {\text {d}}\mathbf {x}+ \int _{{\varOmega }} \frac{1}{2\theta }\Vert \mathbf {u}-\mathbf {v}\Vert )^{2} {\text {d}}\mathbf {x}\end{aligned}$$

  (14)

  where the dual variables are $\mathbf {\xi } = \left( \begin{array}{cc} \xi _{11} &{} \xi _{12} \\ \xi _{21} &{} \xi _{22} \end{array} \right)$ and satisfy $||\xi ||_{F} \le 1$.

Proposition 2

The solution of (14) is given by the following iterative scheme

$$\begin{aligned} \xi _{i1}^{n+1}= & {} \frac{\xi _{i1}^{n} + \tau \overline{u}_{ix}^{n}}{\max (1,||\xi ||_{2})}, \qquad \xi _{i2}^{n+1} = \frac{\xi _{i2}^{n} + \tau \overline{u}_{iy}^{n}}{\max (1,||\xi ||_{2})}, \end{aligned}$$

(15)

$$\begin{aligned} u_{i}^{n+1}= & {} u_{i}^{n} - \sigma \left( \frac{(u_{i}^{n} - v_{i})}{\theta } - \text {div}\left( \xi _{i1}^{n}, \xi _{i2}^{n}\right) \right) , \end{aligned}$$

(16)

$$\begin{aligned} \overline{u}_{i}^{n+1}= & {} 2u_{i}^{n+1} -u_{i}^{n}, \end{aligned}$$

(17)

where $i=1,2$.

• 2. For $\mathbf {u}$ fixed, let us consider each of the two different data terms, $E^1_D(\mathbf {v})$ and $E^2_D(\mathbf {v})$, presented in Sect. 4.1.

  • 2.1. Case $E^1_D(\mathbf {v})$, Li and Osher [32] present a simple algorithm to find the optimal value of the function $E(x) = \sum \nolimits _{i}^{n} w_i |x-a_i| + F(x)$ when the $w_i$ are nonnegative and F is strictly convex. If F is also differentiable and $F'$ is bijective, it is possible to obtain an explicit formula in terms of the median. For $\mathbf {u}$ fixed, we solve

    $$\begin{aligned} \int _{\varOmega }C(\mathbf {v},\mathbf {x}) {\text {d}}\mathbf {x}+ \frac{1}{2\theta } \int _{{\varOmega }} \Vert \mathbf {u}-\mathbf {v}\Vert ^{2} {\text {d}}\mathbf {x}. \end{aligned}$$

    (18)

    Following the ideas of [54], we solve the discrete version of this problem. Due to the isotropy of the quadratic term, the optimal solution of $C(\mathbf {v},\mathbf {x})$ can be obtained solving a one-dimensional problem. In particular, setting $\mathbf {v}= {\hat{\mathbf {v}}} + \delta \frac{\nabla I(x + v_o)}{|\nabla I(x + v_o)|} + \delta \frac{\nabla ^{+} I(x + v_o)}{|\nabla ^{+} I(x + v_o)|}$ being $\nabla ^{+}I$ an orthogonal vector to the gradient, where $\delta$ is our new variable. Then, we minimize over $\delta$

    $$\begin{aligned} \frac{1}{2\tau }\delta ^2 + \lambda \int _{\varOmega }|\nabla I_{t+1}(\mathbf {x}+ {\hat{\mathbf {v}}}_o)|\left| G({\hat{\mathbf {v}}}) + \delta \right| {\text {d}}\mathbf {y}\end{aligned}$$

    (19)

    where

    $$\begin{aligned} G_{\mathbf {y}}({\hat{\mathbf {v}}}) = \frac{I_t(\mathbf {x}) - I_t(\mathbf {y}) - I_{t+1}(\mathbf {x}+ {\hat{\mathbf {v}}}_o) + I_{t+1}(\mathbf {y}+ {\hat{\mathbf {v}}}_o) + ({\hat{\mathbf {v}}}-{\hat{\mathbf {v}}}_o)^{T}\nabla I_{t+1}(\mathbf {x}+ {\hat{\mathbf {v}}}_o)}{|\nabla I_{t+1}(\mathbf {x}+ {\hat{\mathbf {v}}}_o)|} \end{aligned}$$

Proposition 3

The minimum of (19) with respect to $\delta$ is

$$\begin{aligned} \delta ^{*} = \text {median} \{ b_1,...,b_n,a_0,...a_n \} \end{aligned}$$

(20)

where $b_i = -G_{i}({\hat{\mathbf {v}}})$ and $a_i = (n-2i)\lambda |\nabla I_{t+1}(\mathbf {x}+ \mathbf {v}_o)|$ for all the discrete neighbors i (corresponding to $\mathbf {y}$ above), where n is the number of points in the discrete neighborhood.

• 2.2. Case $E^2_D(\mathbf {v})$. Notice that this term is a particular case of the previous data term. The functional to minimize

  $$\begin{aligned} \int _{W} \lambda |\rho (\mathbf {v})| +\frac{1}{2\theta }\int _{W} \Vert \mathbf {u}-\mathbf {v}\Vert ^{2} {\text {d}}\mathbf {x}, \end{aligned}$$

  (21)

  where $\rho (\mathbf {v}) = I_{t+1}(\mathbf {x},\mathbf {v}_o) + \langle \nabla I_{t+1}(\mathbf {x}+ \mathbf {v}_o),(\mathbf {v}- \mathbf {v}_o) \rangle - I_{t}(\mathbf {x})$, does not depend on spatial derivatives on $\mathbf {v}$. Then, a simple thresholding step gives an explicit solution [62].

Proposition 4

The minimum of (21) with respect to $\mathbf {v}$ is

$$\begin{aligned} \mathbf {v}= \mathbf {u}+ \left\{ \begin{array}{lll} &{}\lambda \theta \nabla I_{t+1} &{}\quad \mathrm {if} \quad \rho (\mathbf {u}) < -\lambda \theta |\nabla I_{t+1}|^{2}\\ &{}- \lambda \theta \nabla I_{t+1} &{}\quad \mathrm {if} \quad \rho (\mathbf {u}) > \; \lambda \theta |\nabla I_{t+1}|^{2}\\ &{}-\rho (\mathbf {u})\frac{\nabla I_{t+1}}{|\nabla I_{t+1}|^{2}} &{}\quad \mathrm { if} \quad |\rho (\mathbf {u})| \le \lambda \theta |\nabla I_{t+1}|^{2} \end{array} \right. \nonumber \\ \end{aligned}$$

(22)

Appendix 2: Implementation Details

Our code is written in C. The numerical scheme to solve the functional E(u), in both steps, is based on the implementation of [41]. Image warpings use bicubic interpolation. The image gradient is computed using centered derivatives. Input images have been normalized between [0,1]. The algorithm parameters are initialized with the same default setting for all the experiments. Both time steps are set to $\tau = \sigma = 0.125$ to ensure convergence. As stopping criterion, the optical flow method uses the infinite norm between two consecutive values of u with a threshold of 0.01. The coupling parameter ${{\theta }}$ is set to 0.3. The smoothness term weight ${\beta }$ is set to 1/40 for the ${\text {TV}}_{{\ell _2}}$-L1 functional and ${\beta } = \frac{N-1}{80}$ for the ${\text {NLTV}}$-${\text {CSAD}}$ one, as suggested by [54], where N is the cardinality of the neighborhood considered in the ${\text {CSAD}}$ term (we use a neighborhood of $7\times 7$ in the data term and then N=49) and we fixed $\sigma _c=2$ and $\sigma _s=2$ for the spatial an color domain of the ${\text {NLTV}}$ term. For the iterated faldoi strategy, we set ${\text {MAX}}{\_}{\text {IT}}$ to 3. The size of the patch ${\mathcal {P}}$ in the local minimization is $11\times 11$. The complexity of our algorithm is ${\mathcal {O}}(n)$, where n is the number of pixels of the image frame. The basic faldoi algorithm takes around 20 seconds for ${\text {TV}}_{\ell _2}$-L1 energy and around 10 minutes for ${\text {NLTV}}$-${\text {CSAD}}$ over an ${\text {MPI}}-{\text {Sintel}}$ image. Notice that our algorithm using ${\text {NLTV}}$-${\text {CSAD}}$ energy is very slow, especially because our implementation does not use parallized code. As ${\text {NLTV}}$-${\text {CSAD}}$ can be easily parallelized, it should take the same time for both functionals using a GPU implementation.