Determining dense velocity fields for fluid images based on affine motion

Initial patch matching

The procedure of initial patch matching aims to extract a sparse set of patch matches, facilitating affine motion model estimation. Typically, fluid images possess very few stable salient features but a large number of similar textures. Such images also exhibit high spatial and temporal distortions of luminance patterns. Nevertheless, since the spatial relationships among fluid particles will not change significantly during a small temporal step, we adapt the best-buddies similarity (Dekel et al., 2015) to obtain an optimal match for a given patch. Originally, the best-buddies similarity identifies the fraction of best-buddies pairs, where two points are the nearest neighbors of each other, between two sets of points. We applied the best-buddies similarity measuring the degree of matching between two regions, $\mathrm{R}=\{r_i\in {ℝ}^d{\}}_1^M$ and $\mathrm{P}=\{p_j\in {ℝ}^d{\}}_1^N$, where both $r_i$ and $p_j$ are some small patches of approximately $3\times 3$ pixels. Such a similarity is defined as:

(1) $$BBS(\text{R},\text{P})=\frac{1}{min\{M,N\}}{\displaystyle \sum _i^M{\displaystyle \sum _j^{N}S}}(r_i,p_j).$$where M and N are the number of patches contained in regions $\mathrm{R}$ and $\mathrm{P}$, respectively. Function S shows whether or not the two patches, $r_i$ and $p_j$, are best-buddy pairs (Dekel et al., 2015). If $p_j$ is the nearest neighbor of $r_i$ in R and vice versa, the pair of patches $\{r_i,p_j\}$ is regarded as a best-buddies pair. Formally,

(2) $$S(r_i,p_j)=\{\begin{array}{cc}1& N(r_i,\mathrm{P})=p_j\mathrm{a}\mathrm{n}\mathrm{d}N(p_j,\mathrm{R})=r_i\\ 0& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$$where $N(r_i,\mathrm{P})=\underset{p\in P}{\arg min}d(r_i,p)$, and $d(r_i,p)$ is the distance between the two patches $r_i$ and $p$, which is measured by some distance metrics. We combine the RGB values, the structure and the location of each patch to measure the distance between two given patches.

Color distance. To measure the color distance between two patches, each patch (of $3\times 3$ pixels) is represented by a $9$-element vector of its RGB values, $\mathcal{C}$. In order to overcome the influence of outliers, we use the 2-Norm of color vector to measure color distance as follows:

(3) $$L_c(r_i,p_j)=||{\mathcal{C}}_{r_i}-{\mathcal{C}}_{p_j}|{|}_2^2.$$

Structure distance. Color distance is incapable of evaluating the texture distance of two patches. As illustrated in Fig. 2, the color distance between patches $a$ and $b$ is same as that between patches $a$ and $c$, regardless of their texture differences. So, we integrate the structures of patches into the distance function.

The RLBP operator (Ojala, Pietikäinen & Mäenpää, 2000) defines a texture descriptor T in a local neighborhood of a monochrome texture image as the joint distribution of the signed differences between the gray value of the center pixel, $g_c$ and those of the circularly symmetric neighborhood, $g_n(n=0,1,\dots ,N-1)$:

(4) $$T=t(s(g_c-g_1),s(g_c-g_2),\dots ,s(g_c-g_{N-1})$$

The RLBP operator is an excellent measurement of the spatial structure of a local image texture. If the RLBP values of two patches are the same, these patches are regarded as having similar structures. For example, given $RLBP(a)=15$, $RLBP(b)=255$ and $RLBP(c)=15$, the structures of patches $a$ and $c$ are regarded as similar.

So, we define the structure distance between two patches as:

(5) $$L_s(r_i,p_j)=\frac{|RLBP(r_i)-RLBP(p_j)|}{MA{X}_{rlab}}$$where $MA{X}_{rlab}$ is the maximum of RLBP number, such as 255.

Location distance. As a small region usually exhibits fewer structure changes within a short period of time, if two small regions are similar, the distance between two matching patches should then be small. We therefore take such a distance into account.

For distance estimation, the location of a patch is represented by that of its center pixel, which is a local coordinate within an image region, normalized to between 0 and 1. The location distance is then measured by:

(6) $$L_l(r_i,p_j)=||{\mathcal{L}}_{r_i}-{\mathcal{L}}_{p_j}|{|}_2^2$$

Finally, the distance function between two patches is defined as:

(7) $$d(r_i,p_j)=L_c(r_i,p_j)+\alpha \cdot L_s(r_i,p_j)+\beta \cdot L_l(r_i,p_j)$$where $\alpha$ and $\beta$ are scalar coefficients.

To proceed, the original image $I_1$ is first divided into regions of $3{\mathcal{R}}_w\times 3{\mathcal{R}}_h$ pixels. Each region is further subdivided into ${\mathcal{R}}_w\times {\mathcal{R}}_h$ distinct patches, where each patch contains $3\times 3$ pixels. Second, for each region of image $I_1$, the optimal matching region from image $I_2$ can be determined by maximizing the best-buddies similarity between the two regions (see Eq. (1)). To improve reliability, a region from $I_2$ can be used to estimate the initial patch matching only if the resulting best-buddies similarity is higher than a given threshold, e.g., 0.8, in our experiments. The procedure for generating initial patch matches is shown in Algorithm 1.

Affine motion model

The expected fine and dense fluid field for a fluid scene can be obtained by estimating a motion model for each patch. Palaniappan et al. (1995) and Zhou et al. (2001) confirmed that the affine motion model is suitable to model small areas of local fluid motion. Therefore, we estimated an affine motion model for each small patch. We use a 2D affine motion model (Andrew, 2003) to fit a local patch motion as follows:

(8) $$q_{{\mathrm{r}}_i^{\mathrm{\prime }}}^k=M_i{q}_{{\mathrm{r}}_i}^k+T_i$$where $({\mathrm{r}}_i,{\mathrm{r}}_i^{\mathrm{\prime }})\in \mathcal{M}$, $q_{{\mathrm{r}}_i}^k$ and $q_{{\mathrm{r}}_i^{\mathrm{\prime }}}^k$ are two matching points in the local patches ${\mathrm{r}}_i$ and ${\mathrm{r}}_i^{\mathrm{\prime }}$, respectively. $M_i$ is an affine transformation matrix for the patch ${\mathrm{r}}_i$, and $T_i$ is a translation vector.

To get $M_i$ and $T_i$, we have to obtain a set of point matches which contains at least six elements for each patch match. If a patch undergoes rotation through fluid motion such as vorticity, it is not suitable for a pair of points at the same location of two matching patches to be considered as a matching pair. Hence, for each patch match $({\mathrm{r}}_i,{\mathrm{r}}_i^{\mathrm{\prime }})\in \mathcal{M}$, both ${\mathrm{r}}_i$ and ${\mathrm{r}}_i^{\mathrm{\prime }}$ are rotated using the rotation invariant binary patterns (Ojala, Pietikäinen & Mäenpää, 2000). Two points at the same location of the two patches are then considered as a matching pair of points. Finally, we obtain a set of point matches for each patch match, denoted by:

(9) $$\mathcal{M}{\mathcal{P}}_{{\mathrm{r}}_j}=\{(q_{{\mathrm{r}}_j}^k,q_{{\mathrm{r}}_j^{\mathrm{\prime }}}^k)\},k=1,2,\dots ,|{\mathrm{r}}_j|$$where $|{\mathrm{r}}_j|$ denotes the number of points in the patch ${\mathrm{r}}_j$.

The transformation matrix $M_j$ and the translation vector $T_j$ can be computed by solving the affine motion model (see Eq. (8)) for a patch ${\mathrm{r}}_j$ with common numerical methods. We use successive over-relaxation (SOR) to solve the system. Finally, we obtain a sparse set of affine motions, denoted as $\mathcal{A}\mathcal{M}=\{(M_i,T_i)\}$.

Interpolation

The few patch matches were transformed into a dense corresponding field by interpolating a sparse set of affine motion models $\mathcal{A}\mathcal{M}=\{(M_i,T_i)\}$.

The distance between two patches is an important factor for computing proper interpolation. We used the Euclidean distance $\mathfrak{D}$ between patches to measure the distance. As the influence of remote patches is either negligible or harmful to the interpolation, we restricted the patches used in the interpolation at a patch $r$ to its $\rho$ neighbors, denoted as $N_{\rho }(r)$, according to distance $\mathfrak{D}$.

If $BBS(\mathrm{R},{\mathrm{R}}^{\mathrm{\prime }})$ is larger, it shows that the results of patch matches inside this region is more reliable. Therefore, we denoted the value $BBS(\mathrm{R},{\mathrm{R}}^{\mathrm{\prime }})$ as the reliability of a match $(r_i,r_i^{\mathrm{\prime }})\in \mathcal{M}$ which is used as a scale element for interpolation.

Combing the distance and reliability elements, an interpolation function is defined as:

(10) $$M_D(r)=\frac{\sum _{r_k\in N_{\rho }(r),(r_k,r_k^{\mathrm{\prime }})\in \mathcal{M}}K(r_k,r)(M_k,T_k)}{\sum _{r_k\in N_{\rho }(r),(r_k,r_k^{\mathrm{\prime }})\in \mathcal{M}}K(r_k,r)}$$where $r_k\in R$, $r_k^{\mathrm{\prime }}\in P$, $K(r_k,r)=exp(-\mathfrak{D}(r_k,r))+\mu BBS(R,P)$ and $\mu$ is a weight coefficient. This allows us to estimate a dense set of affine motion models, denoted as $\overline{\mathcal{A}\mathcal{M}}$ by Eq. (10), where the sparse set of affine motion models $\mathcal{A}\mathcal{M}$ as the inputs.

Fluid flow optimization

After interpolation, each patch has an affine motion model describing its motion. We then obtained dense fluid flows by carrying out affine translation according to Eq. (8) for each patch. Since only local information is used to estimate the dense fluid flows within each small patch independently, motion may change drastically between two adjacent patches or may have some noise. Meanwhile, due to the absence of a global description of fluid motion, there may be a high degree of overfitting involved in each small area. Hence, some global constraints can be applied to improve the initial fluid motion estimation.

A variational energy function is used to optimize optical flow. For instance, in Brox et al. (2004), brightness constancy assumption, gradient constancy assumption and smoothness constraint were involved. Also, minimizing a variational energy function with a full-scale dense initial flows as the input, can generate better results compared to the conventional coarse-to-fine scheme (Brox et al., 2004; Sun, Roth & Black, 2014). This is applicable to a scene with boundary overlapping or thin objects (Chen, Li & Hall, 2016; Revaud et al., 2015). These phenomena often occur in fluid motion scenes, as the diffusion of fluid may generate some thin fluid areas in the surroundings.

Hence, the dense field generated by our method is used as the initialization of a variational energy minimization method to improve motion estimation. Let $\mathbf{x}=(x,y{)}^T$, $\mathbf{w}=(u,v{)}^T$ and $\mathrm{\nabla }=(\mathrm{\partial }x,\mathrm{\partial }y{)}^T$. We use an energy functional expressing as:

(11) $$\begin{gathered} E(\mathbf{w})={\int }_{\mathrm{\Omega }}\{\mathrm{\Psi }(||I_1(\mathbf{x}+\mathbf{w})-I_2(\mathbf{x})|{|}^2)+ \\ \kappa \mathrm{\Psi }(||\mathrm{\nabla }{I}_1(\mathbf{x}+\mathbf{w})-\mathrm{\nabla }{I}_2(\mathbf{x})|{|}^2)\}d\mathbf{x}+ \\ \gamma {\int }_{\mathrm{\Omega }}\mathrm{\Psi }(||\mathrm{\nabla }u|{|}^2+||\mathrm{\nabla }v|{|}^2) \end{gathered}$$where both $\kappa =30$ and $\gamma =80$ are tuning parameters, the three terms on the right hand side of the equation are brightness constancy, gradient constancy, and smoothness constraint, respectively, and $\mathrm{\Psi }(s)$ is a penalty function (Brox et al., 2004). According to the calculus of variations, we compute Euler-Lagrange equations by preforming three fixed-point iterations and the flow updates three times iteratively.

Qualitative analysis

Synthetic image. It is difficult to obtain the ground truth motion field for a real-world fluid phenomenon, therefore, we first synthesized a test image sequence with the given ground truth flow data. To do this, we applied the advection motion and the vortex motion to generate a motion field, which was smoothed by a Gaussian filter. The top, right-hand images in Fig. 3 shows a complicated velocity field that we generated; this image was used as the ground truth for comparison. The figure presents a visual comparison of different methods (second and third rows) against the ground truth velocity field (middle of top row). In the results, the length of each flow vector indicates the magnitude of velocity. To analyze the results, we highlighted part of the flow vectors of the reconstructed smoke motion in a green box, while highlighting part of the reconstructed non-fluid background in a red box. For readability, we also show enlarged views of both boxes. In order to quantitatively evaluate the estimated fluid flow, we choose the average endpoint error (EPE) which computes the Euclidean Distance of the true flow and the estimated one. The proposed method achieves an EPE score of 5.32 surpassing the previous state-of-the-art methods HS, Classic + NL, LODF, EpicFlow, and S2F by $15.98\mathrm{\%}$ (6.17 $\to$ 5.32), $9.59\mathrm{\%}$ (5.83 $\to$ 5.32), $3.0\mathrm{\%}$ (5.48 $\to$ 5.32), $7.52\mathrm{\%}$ (5.72 $\to$ 5.32), $1.31\mathrm{\%}$ (5.39 $\to$ 5.32), respectively.

Generally, all methods can capture both advection and vortex. However, all methods except ours overestimate motion fields for the non-fluid background, generating an improper velocity field and unexpected motions. The Classic+NL method obtains the worst result, generating significantly large velocities for the non-fluid background. In contrast, our method estimates the motion field of the background as having very low (close to zero) velocities. It implies that our method can more accurately estimate fluid flows from images in terms of both fluid phenomena and non-fluid objects.

Smoke. Figure 4 compares fluid motion estimation results. The two adjacent input smoke images are shown in the top row of Fig. 4 as the first two images. Fluid motion estimation results generated from different methods are depicted in rows two to seven of the figure. Specifically, the resultant images of smoke estimation from those methods are presented in the second column, while the velocity field of the highlighted smoke motion region (indicated by the green boxes) is shown in the first column. By visually comparing the results, the reconstructed smoke motion (done by warping interpolation) of our method closely reproduces the expected result (Fig. 4, second image, first row). For results from other methods, for instance, the flow field generated by the HS method is quite noisy. LDOF (Brox, Bregler & Malik, 2009) overestimates flow changes inside the smoke phenomenon, as smoke does not comprise regular structures. More importantly, we produced significantly better flow field estimation for the non-fluid background region, which is important for some applications, e.g., frame detection.

Hurricane. The second fluid phenomenon we tested was a hurricane scene (Fig. 4). The results are depicted in the last two columns. The results are shown as for the smoke motion. Our method captures complicated flow fields extremely well compared to other methods, as highlighted by both the yellow and red boxes. Specifically, besides being capable of capturing complex patterns of flow field, our generated results are less noisy or distorted. Such that our reconstructed image can highly reproduce the ground truth (top right image).

Clouds. Figure 5 depicts another set of experiment results. Here a dense flow field from two input images showing a scene with a cloudy sky is estimated with some rigid moving objects with moving ones, i.e., boats. The first row shows two input images (columns 1, 3) and two of their zoomed partial views comprising the object motions (columns 2, 4 with three small images each). Six methods are involved and their results are shown in 1. Here, the generated flow field (left column) and the reconstructed results (right column), which are obtained by warping the first frame image against the generated flow field, are shown for each method. Both the estimated velocity fields and the interpolation (reconstructed) results show that our method can more accurately capture the subtle cloud motion and the large displacement of boats comparing to other methods. Note that all methods, with the exception of ours and LDOF, overestimate the velocity fields of cloud mainly due to brightness change and a large number of similar textures. However, cloud and boat motions cannot be properly captured by the LDOF method because it lacks the proper segmentations for these features.

Quantitative analysis

The interpolation error (IE) (Baker et al., 2007) is used as the quantitative measure in our experiments because it lacks the ground truth for real world fluid phenomena. For image interpolation, we define IE to be the root-mean-square (RMS) difference between the real second frame ( $I_2(x,y)$) and the interpolated image ( $I_w(x,y)$) by warping the first frame with the flow field result:

(12) $$IE={\left[\frac{1}{N}{\displaystyle \sum _{(x,y)}(I_w(}x,y)-I_2(x,y){)}^2\right]}^{\frac{1}{2}}$$where N is the number of pixels.

Table 1 compares the errors incurred from our method and other five methods for the four fluid phenomena described in our study. Results show that our method incurred less error than other methods.

We further quantitatively measure the errors between the ground truth $I_2$ and the warping interpolation frames by computing image differences. Figure 6 shows the results of all the methods we have tested against the three fluid phenomena (smoke, clouds, and hurricane). Our results demonstrate that warping images produced by our method are more consistent with the ground truths than those produced by other methods.