Full Flow: Optical Flow Estimation By Global Optimization over Regular Grids

Introduction

Optical flow is a vital source of information for visual perception. Animals use optical flow to track and control self-motion, to estimate the spatial layout of the environment, and to perceive the shape and motion of objects . In computer vision, optical flow is used for visual odometry, three-dimensional reconstruction, object segmentation and tracking, and recognition.

The classical approach to dense optical flow estimation is to optimize an objective of the form

where $\mathbf{f}$ is the estimated flow field, $E_{\text{data}}$ is a data term that penalizes association of visually dissimilar areas, and $E_{\text{reg}}$ is a regularization term that penalizes incoherent motion . Traditionally, this objective is optimized by iterative local refinement that maintains and updates a single candidate flow . This local refinement does not optimize the objective globally over the full space of flows and is prone to local minima.

Global optimization of the flow objective has generally been considered intractable unless significant restrictions are imposed . Menze et al. achieved impressive results with a discrete optimization approach, but had to heuristically prune the space of flows using descriptor matching.

In this work, we develop a global optimization approach that optimizes the classical flow objective (1) over the full space of mappings between discrete grids. Our work demonstrates that a direct application of global optimization over full regular grids has significant benefits. Since the highly regular structure of the space of mappings is preserved, we can employ optimizations that take advantage of this structure to reduce the computational complexity of the algorithm’s inner loop. The overall approach is simple and does not involve separately-defined descriptor matching modules: simply optimizing the classical flow objective over full grids is sufficient. We show that this minimalistic approach yields state-of-the-art accuracy on both the Sintel and the KITTI 2015 optical flow benchmarks.

Background

The variational approach to optical flow originates with Horn and Schunck . This elegant approach posits a clear global objective (1) and produces a dense flow field connecting the two images. Since the space of flows is so large, the variational objective has traditionally been optimized locally. Starting with a simple initialization, the flow is iteratively updated by gradient-based steps . Through these iterations, a single candidate flow is maintained. While this local refinement approach can be accurate when displacements are small , it does not optimize the objective globally over the full space of flows and is prone to local minima.

Recent methods have used descriptor matching and nearest neighbor search to initialize the continuous refinement . This more sophisticated initialization is known to significantly improve results in the presence of large displacements. However, the descriptor matching module is trained separately, does not optimize a coherent objective over the provided correspondence sets, and can yield globally suboptimal initializations. We show that state-of-the-art accuracy can be achieved by globally optimizing the classical objective (1), with no separately trained or designed descriptors.

A number of approaches to global optimization for optical flow estimation have been proposed. Steinbrücker et al. use an alternating scheme to optimize a quadratic relaxation of the global objective. This formulation relies on the assumption that the regularizer is convex. A number of subsequent approaches use functional lifting to map the problem into a higher-dimensional space, where the optimization reduces to estimating a collection of hypersurfaces . These schemes likewise impose certain assumptions on the model, such as requiring the data term or the regularizer to be convex. In general, these approaches have not been shown to produce state-of-the-art results on modern benchmarks.

Our approach treats objective (1) as a Markov random field and uses discrete optimization techniques. The Markov random field perspective on optical flow estimation dates back to the 80s and discrete optimization techniques have been applied to the problem in different forms since that time . Glocker et al. applied MRF optimization to sets of control points in coarse-to-fine schemes. In contrast, we operate on dense grids with large two-dimensional label spaces. A number of works considered a simplified MRF formulation that decomposes the horizontal and vertical components of the flow . In contrast, we demonstrate the feasibility of operating on much larger models with two-dimensional label spaces. Lempitsky et al. iteratively improved the estimated flow field by generating proposals and integrating them using the QPBO algorithm. In contrast, we optimize over the full space of mappings between discrete grids. Komodakis et al. evaluated MRF optimization on optical flow estimation with small displacements. In contrast, we show that global optimization over full two-dimensional label spaces is tractable and yields state-of-the-art performance on challenging large-displacement problems.

Menze et al. pruned the space of flows using feature descriptors and optimized an MRF on the pruned label space. In contrast, we argue that operating on the full space is both feasible and desirable. First, we avoid heuristic pruning and the reliance on separately-defined feature descriptors that are not motivated by the flow objective itself. Second, pruning destroys the highly regular structure of the space of mappings. We show that optimization over the full space can be significantly accelerated due to the regularity of the space. In particular, the full regular structure enables the use of highly optimized min-convolution algorithms that reduce the complexity of message passing from quadratic to linear .

Model

where $\mathcal{N}\subset\Omega^{2}$ is the 4-connected pixel grid. See Figure 1 for illustration. The data term $\rho_{D}(p,f_{p},I_{1},I_{2})$ penalizes flow fields that connect dissimilar pixels $p$ and $(p+f_{p})$. We use truncated normalized cross-correlation :

where $NCC$ is the normalized cross-correlation between two patches, one centered at $p$ in $I_{1}$ and one centered at $(p+f_{p})$ in $I_{2}$, computed in each color channel and averaged. The truncation at zero prevents penalization of negatively correlated patches. If $(p+f_{p})$ is in the buffer zone $\mkern 1.5mu\overline{\mkern-1.5mu\Omega\mkern-1.5mu}\mkern 1.5mu\setminus\Omega$, the data term is set to a constant penalty $\zeta$.

Our optimization approach assumes that the regularization term has the following form:

where $f^{1},f^{2}$ are the two components of vector $f$ and $\rho(\cdot)$ is a penalty function, such as the $L^{1}$ norm or the Charbonnier penalty. Our formulation and the general solution strategy can accomodate non-convex functions $\rho$, such as the Lorentzian and the generalized Charbonnier penalties. The reduction of message passing complexity from quadratic to linear, described in Section 4.2, applies to such functions as well . The highly efficient min-convolution algorithm described in Section 4.3 will assume that the function $\rho$ is convex. Other linear-time algorithms can be used if this assumption doesn’t hold .

Note that the regularization term (4) couples the horizontal and vertical components of the flow. We apply a Laplace weight to attenuate the regularization along color discontinuities:

Optimization

Objective (2) is a discrete Markov random field with a two-dimensional label space . At first glance, global optimization of this model may appear intractable. The number of nodes and the number of labels are both in the tens of thousands. The most sophisticated prior application of discrete optimization to this problem resorted to pruning of the label space to bring the size of the problem under control . We show that the full problem is tractable.

The pixels are then visited in reverse order with analogous message update rules. This completes one forward-backward pass, considered to be one iteration. A number of iterations are performed.

Given updated messages $\mathbf{m}$, we can compute a solution $\mathbf{l}$ greedily . We determine the labels sequentially in a given order. Upon reaching pixel $p$, we choose a label assignment $l_{p}$ that minimizes $\theta_{p}(l_{p})+\sum_{q<p}{\theta_{pq}(l_{p},l_{q})}+\sum_{p<q}{m_{q\rightarrow p}(l_{p})}$, where $p<q$ means that $p$ precedes $q$ in the order.

2 Complexity reduction

A brute-force implementation of a message update requires $O(M^{2})$ operations as there are $M$ elements in each message and computing each element directly requires $O(M)$ operations according to update rule (6). We now show that a message update can be performed using $O(M)$ operations in our model. This builds on the min-convolution acceleration scheme developed by Felzenszwalb and Huttenlocher . A general treatment of the one-dimensional case was presented by Chen and Koltun .

We begin by rewriting the message update rule (6) as

Each $\phi_{pq}(s)$ can be computed using $O(1)$ operations, thus $\phi_{pq}$ can be computed using $O(M)$ operations in total. We now show that, given $\phi_{pq}$, all elements of the message $m_{p\rightarrow q}$ can be computed in $O(M)$ operations as well. Recall that $\theta_{pq}(s,t)=\rho_{S}(t-s)$ and that $\rho_{S}(\cdot)$ has the form given in equation (4). Rearranging terms, we obtain

$T_{pq}$ can be computed using $O(M)$ operations given $\phi_{pq}$. We now show that $\mathcal{D}_{pq}$ can also be computed using $O(M)$ operations in total. Note that $s$ and $t$ are 2D vectors. Abusing notation somewhat, we can rewrite $\mathcal{D}_{pq}(t)$ as a two-dimensional min-convolution: $\displaystyle\mathcal{D}_{pq}(t^{1},t^{2})$ $\displaystyle=\min\limits_{s^{1},s^{2}}$ \displaystyle\Big{(}\phi_{pq}(s^{1},s^{2})+\rho(t^{1}-s^{1})+\rho(t^{2}-s^{2})\Big{)} $\displaystyle=\min\limits_{s^{2}}$ \displaystyle\Big{(}\min_{s^{1}}\big{(}\phi_{pq}(s^{1},s^{2})+\rho(t^{1}-s^{1})\big{)} \displaystyle+\ \rho(t^{2}-s^{2})\Big{)}. This can be decomposed into two sets of $O(\sqrt{M})$ one-dimensional min-convolutions:

For each $s^{2}$, $\mathcal{D}_{pq|s^{2}}(t^{1})$ can be computed for all $t_{1}$ by a 1D min-convolution. Then, for each $t^{1}$, $\mathcal{D}_{pq}(t^{1},t^{2})$ can be computed for all $t^{2}$ by a 1D min-convolution. Each min-convolution can be evaluated in $O(\sqrt{M})$ operations, for a total complexity of $O(M)$.

3 Further acceleration

A min-convolution has the following general form:

It is well-known that the min-convolution can be computed using $O(n)$ operations, where $n=\sqrt{M}=2\varsigma+1$ . However, commonly used algorithms require computing intersections of shifted copies of the function $\rho$. While each intersection can be computed in time $O(1)$, this computation can be numerically intensive for some penalty functions. Since this computation is in the inner loop, it can slow the optimization down. We now review an alternative algorithm that can be used to compute the min-convolution without computing intersections. This algorithm relies on the assumption that $\rho$ is convex, which is otherwise not necessary. Related algorithms are reviewed by Eppstein .

The algorithm is based on totally monotone matrix searching . Let $A$ be an $n\!\times\!n$ matrix, such that $A(i,j)=g(j)+\rho(i-j)$. Let $\textup{ind}_{A}(i)$ be the column index of the minimal element in the $i$th row of $A$. The min-convolution $h$ can be defined as $h(i)=A(i,\textup{ind}_{A}(i))$. The challenge is to evaluate $\textup{ind}_{A}$ in time $O(n)$ without explicitly constructing the matrix $A$.

The convexity of $\rho$ implies that $A$ is totally monotone. The totally monotone matrix search algorithm computes $\textup{ind}_{A}$ in $O(n)$ operations by divide-and-conquer. The algorithm first constructs an $\frac{n}{2}\!\times\!n$ submatrix $B$ by removing every other row of $A$. Then $B$ is reduced to an $\frac{n}{2}\!\times\!\frac{n}{2}$ submatrix $C$ by removing columns that do not contain minima of the rows of $B$. The minima $\textup{ind}_{C}$ are computed recursively, after which the missing elements of $\textup{ind}_{A}$ are filled in. As shown by Aggarwal et al. , all steps can be performed in time $O(n)$. Crucially, all steps can be performed without explicitly constructing $A$.

Implementation

To reduce wall-clock time, we implemented a parallelized TRW-S solver. This general-purpose solver along with the rest of our implementation will be made freely available. At each step of TRW-S, a pixel is ready to be processed if all of its predecessors have already been updated during the current iteration. Thus at any time there is a wavefront of pixels that can be processed in parallel. A grid can be swept diagonally. In the first step only one node can be processed, but the size of the wavefront grows rapidly and all nodes on the wavefront can be processed in parallel. This parallelization scheme has previously been explored on special-purpose hardware for stereo matching . We have implemented the scheme on general-purpose processors. Our implementation is evaluated on a workstation with a 6-core Intel i7-4960X CPU. Parallelization with hyper-threading reduces the running time of each iteration of TRW-S from 256 to 39 seconds, a factor of 6.6. Since the size of the wavefront is $\Theta(\varsigma)$ for most of the iteration, increased hardware parallelism is expected to directly translate to reduction in wall-clock time. We also refer the reader to the concurrent work of Shekhovtsov et al. , who developed a parallelized energy minimization scheme that may be applicable to our setting.

Occlusion handling.

Some pixels in $I_{1}$ may not have corresponding points in $I_{2}$. The computed flow field on these occlusion pixels is likely incorrect. We adopt the common tactic of forward-backward consistency checking: compute the forward flow from $I_{1}$ to $I_{2}$ and the backward flow from $I_{2}$ to $I_{1}$, and discard inconsistent matches . Given the forward flow field $\mathbf{f}$ and the backward flow field $\mathbf{f}^{\prime}$, the following criterion is used to determine whether a match is consistent. For each pixel $p$ in $I_{1}$ and its match $(p+f_{p})$ in $I_{2}$, $f_{p}$ is said to be consistent if there is a pair $(q+f^{\prime}_{q},q)\in I_{1}\!\times\!I_{2}$ that is close to the pair $(p,p+f_{p})$. Specifically, for each $f_{p}$, we check if there exists $f^{\prime}_{q}$ for which

This test can be performed by finding nearest neighbors across two point sets: $\{(p,p+f_{p})\}_{p\in\Omega}$ and $\{(q+f^{\prime}_{q},q)\}_{q\in\Omega}$.

Postprocessing.

We optimize the model described in Section 3, remove inconsistent matches as described in the previous paragraph, and then interpolate the results to get subpixel-resolution flow. We use the EpicFlow interpolation scheme , which has become a common postprocessing step in recent state-of-the-art pipelines . Since an interpolation step is necessary to obtain subpixel-accurate flow, the discrete optimization need not operate at the highest resolution. We found that optimizing the presented model on $1/3$-resolution images still yields state-of-the-art performance. We attribute this both to the power of the presented global optimization approach and to the effectiveness of the EpicFlow interpolation scheme. The postprocessing is illustrated in Figure 2.

Experiments

The presented approach is implemented in Matlab, with a C++ wrapper for the parallelized TRW-S solver. Our Matlab code is less than 50 lines long, not including the general-purpose solver.

The experiments are performed on two challenging optical flow datasets, MPI Sintel and KITTI 2015 . We use a workstation with a 6-core Intel i7-4960X 3.6GHz CPU and 64GB of RAM. The computation of the 1.3 billion values in the unary cost volume takes 51 seconds. We use $3\!\times\!3$ patches for NCC in $1/3$-resolution images. Performing 3 iterations of TRW-S using the general optimization framework described in Section 4 takes about 2 minutes on either Sintel or KITTI images, downsampled by a factor of 3, with any penalty function $\rho$. When the penalty function is the $L^{1}$ norm, we can accelerate the optimization further with the $L^{1}$ distance transform , which reduces the running time to about 30 seconds for 3 iterations of TRW-S. EpicFlow interpolation takes 3 seconds. For each dataset, we train the parameters on $5\%$ of the training set by grid search. We use the same parameters for the ‘final’ and ‘clean’ sequences in the Sintel dataset.

In experiments reported in this section, we use the $L^{1}$ norm for regularization ($\rho(x)=|x|$) and no truncation ($\tau=\infty$). This decision is motivated by the controlled experiments reported in Section 6.2.

MPI Sintel is a dataset for large-displacement optical flow . There are two types of sequences in the dataset, clean and final. The clean sequences exhibit a variety of illumination and shading effects including specular reflectance and soft shadows. The final sequences additionally have motion blur, depth of field, and atmospheric effects.

The experimental results are provided in Table 1. We use the $10$ metrics reported by Bailer et al. , including all, noc, occ, d0-10, and s40+ for both clean and final test sequences. All the errors are measured as endpoint error (EPE), which is the Euclidean distance between the estimated flow vector and the ground truth. Since some error metrics are extremely close for different methods, and because the average EPE is sensitive to outliers (the top methods generally have errors of $20$ to $40$ pixels on a number of challenging sequences), we highlight every method that achieves within $1\%$ of the lowest reported error as one of the top methods according to that error metric.

Our approach outperforms EpicFlow , TF+OFM , NNF-Local , PH-Flow , and Classic+NL on almost all metrics. Our approach ranks 2nd on the key EPE-all metric for both final and clean sequences. Compared to EpicFlow, our approach reduces EPE-all by $6.2\%$ on the final sequences and by $12.5\%$ on the clean sequences.

KITTI 2015.

KITTI Optical Flow 2015 is an optical flow dataset that comprises outdoor images of dynamic scenes captured from a car . The car is equipped with a LiDAR sensor and color cameras. Ground-truth flow is obtained by rigid alignment of the static environment and by fitting CAD models to moving objects. Ground-truth correspondences are sparse. The dataset contains $200$ training sequences and $200$ test sequences. A flow vector is considered an outlier if its endpoint error is 3 pixels or higher. Table 2 lists the most accurate methods on this dataset, along with the classical Horn-Schunck algorithm for reference. Note that SOF was developed concurrently with our work and uses substantially more information at training time, at the cost of generality.

Qualitative results.

In Figure 4, we compare our visual results to EpicFlow and DiscreteFlow on three scenes from MPI Sintel and three scenes from KITTI 2015. On MPI Sintel, our approach performs well on regions with large displacements (we rank first on s40+ in Table 1). This is also reflected in the visual results. See the flapping wings in scene 1 and the flying butterfly in scene 2. In scene 3, all three methods fail but our approach and DiscreteFlow recover more of the flow field than EpicFlow. On KITTI 2015, our approach is visually similar to DiscreteFlow in most street scenes (for example, scene 1). In some cases, our approach is visually more accurate (scene 2), but not on others (outliers on the white line in scene 3). Both our approach and DiscreteFlow are visually more accurate than EpicFlow.

2 Controlled experiments

The generality of the presented optimization framework enables a comprehensive evaluation of different data terms and regularization terms. We perform such an evaluation using $5\%$ of the MPI Sintel training set (final pass). In all conditions, we optimize variants of the model presented in Section 3 using the method presented in Sections 4 and 5. We evaluate two data terms: the patch-based truncated NCC term given in equation (3) and the classical pixelwise Horn-Schunck data term given by the squared Euclidean distance in RGB color space. We also evaluate three penalty functions for the regularization term (equation 4): $L^{1}$ ($\rho(x)=|x|$), squared $L^{2}$ ($\rho(x)=x^{2}$), and Charbonnier ($\rho(x)=\sqrt{x^{2}+\varepsilon^{2}}$, where $\varepsilon=5$). For each penalty function, we evaluate a truncated regularizer ($\tau$ is determined by grid search) and a non-truncated convex form ($\tau=\infty$). All free parameters are determined by grid search.

The results are shown in Table 3, which provides the average EPE over the images used for the evaluation for each combination of the three factors (data term, penalty function, truncation). The results suggest that the data term is of primary importance: the patch-based truncated NCC term is much more effective than the pixelwise Horn-Schunck data term, irrespective of the regularizer. Note that the non-truncated HS+$L^{2}$ condition corresponds to global optimization of the classical Horn-Schunck model. The results for the non-truncated NCC+$L^{2}$ condition indicate that by replacing the pixelwise Horn-Schunck data term with patch-based truncated NCC, retaining the classical non-truncated quadratic regularizer, and globally optimizing the objective we come within $10\%$ of the accuracy achieved by our top-performing variant. The key factors are global optimization and a patch-based data term.

Figure 3 provides a more detailed visualization of the results. For each condition, the figure shows the average endpoint error for each tested image, sorted by magnitude. This figure shows models with truncation in the regularizer. The plots indicate that for most tested images the accuracy achieved by global optimization is high in all conditions, irrespective of the tested factors. The different conditions, specifically the two data terms, are distinguished by their robustness when accuracy is low. The patch-based NCC data term limits the error on challenging images much more effectively than the pixelwise HS data term.

To summarize, the presented optimization approach was designed to support global optimization with very general data and regularization terms. The generality of the presented framework enabled a controlled evaluation of global optimization with different data terms and regularizers. The results indicate that within a global optimization framework the detailed form of the regularizer is less important than the data term, the classical quadratic regularizer yields competitive performance, and the highest accuracy is achieved using the $L^{1}$ penalty.

Conclusion

We have shown that optimizing a classical Horn-Schunck-type objective globally over full regular grids is sufficient to initialize continuous interpolation and obtain state-of-the-art accuracy on challenging modern optical flow benchmarks. In particular, this demonstrates that state-of-the-art accuracy on large-displacement optical flow estimation can be achieved without externally defined descriptors. The flow objective itself is sufficiently powerful to produce accurate mappings even in the presence of large displacements. We have shown that optimizing the objective globally over the full space of mappings between regular grids is feasible and that the regular structure of the space enables significant optimizations.

Our Matlab implementation is less than 50 lines long, excluding the general-purpose TRW-S solver. We hope that the simplicity of our approach will support further advances. More advanced data terms can easily be integrated into our global optimization framework and are likely to yield even more accurate results. In addition, we believe that there is scope for further improvement in continuous interpolation accuracy, building on the impressive performance of the interpolation scheme of Revaud et al. . The output of the presented global optimization approach can serve as a canonical initialization for benchmarking such continuous interpolation schemes.

References