Estimating the moments of a random forcing field of 2D fluid from image sequences using energy minimisation method

In this paper, we consider a version of energy minimisation technique applied to images of a 2D fluid flow. The two Navier–Stokes equations describe the static flow of a 2D fluid in terms of velocity field, (u, v), pressure field, p and forcing field, f. Apart from these two Navier–Stokes equations, we have the incompressibility condition to evaluate the three parameters. While implementing this system, random noise (usually non-Gaussian) creeps into the random force field f̲(x,y)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{f}(x,y)$$\end{document}. We denote this random field by δf̲(x,y)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \underline{f}(x,y)$$\end{document} having zero mean and non-trivial second and third moments. We assume that these two moments are known except for some unknown parameters θ̲\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\theta }$$\end{document} (like mean, variance, co-variance, skewness, etc.) which we wish to estimate. In the proposed technique, we first calculate the approximate shift in the average fluid energy defined as a quadratic function of the velocity field. The energy method then requires that θ̲\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\theta }$$\end{document} should be such that this average increases in the energy due to the random forcing component be minimised. We should, however, note that the standard statistical approach to force field estimation is to calculate the velocity field as a function of the force field and then adopt the statistical moment matching technique. Such an approach assumes spatial ergodicity of the velocity field. This approach to force field estimation is more accurate from the statistical moment matching view point but works only if velocity measurements are made. The former technique of energy minimisation does not require any velocity measurements. Both of these techniques are discussed in this paper and MATLAB simulations presented.


Introduction
Statistical computation of optic flow proposed by Horn and Schunck [1] estimates the velocity field from a sequence of images. The constraint term is derived on the principle of conservation of intensity values over consecutive image patterns. In [1], a method for finding the optical flow pattern is presented which assumes that the apparent velocity of the brightness pattern varies smoothly almost everywhere in the image. Richard et al. [2] presented an approach to measure the fluid flow from some sequence of images. The physics of fluid mechanics is used to propose the algorithm of motion recovery. The paper depicts how extracted information from the physical activity of fluids affects the motion recovery algorithm. Next, the algorithm is applied on synthetic imagery and natural imagery and effectiveness of developed algorithm is discussed. Corpetti et al. [3] proposed an optical flow estimator based on image sequences representing fluid flow phenomenon. This technique is based on the adaptation of the data model and a smoothness term. Further extending his work, Corpetti [4] suggested a data model based on continuity equation for the density of fluid ∂ρ ∂t + ∇ · (ρv) = 0, (1) where ρ is the density of the fluid and v is the 3D velocity vector. Alternatively, through calculations, the continuity equation may be written as dρ dt where dρ dt = ∂ρ ∂t + v · (∇ρ). The equation is similar to 2D optical flow constraint on luminance when the divergence vanishes. Bruhn et al. [5] analysed the smoothing effects in local and global differential methods for optic flow computation. In [6], Ashish et al. presented the technique for modelling complex optical flow patterns using the Navier-Stokes equation for a turbulent fluid. The various terms that appear in the Navier-Stokes equation like controlling force, diffusion terms, nonlinear advection term, etc., have been implemented using discrete step-wise algorithms. In particular, the diffusion term has been implemented in the frequency domain because ∇ 2 becomes a multiplication operator in the spatial frequency domain. Before implementing the advection step, a Gaussian smoothing or the fluid flow has been performed. The diffusion step has been shown to be robust. In [7], Kharbat et al. suggested to replace the individual intensity function of each pixel with a geometric moment of its neighbourhood. Every pixel is made invariant to fluctuations of intensity by doing normalisation. The algorithm was shown to be robust against reflections, shadows and varying illumination. Tianshu et al. [8] established the direct relation between the fluid flow and optical flow for visualisation of typical flow. Ashish et al. [9] proposed robust Hessian-based diffusion, advection and mass conservation for the formulation of stable fluid solver.
The proposed method is applied on synthetic field which is generated with the help of Navier-Stokes equations and it is used to model the complex optical flow using real images.
Feng et al. [10] presented modified Navier-Stokes potential flow method and also suggested models of brightness variation in fluid type motions. In optical flow estimation, the general assumption of constant brightness will not be applicable in case of smoke or cloud type motion. Dirichlet-Neumann operator (DNO) is applied to estimate the velocity potential of such images. Aaron et al. [11] used optical flow method to analyse the fluid transport dynamics. Stream function optical flow (SFOF) regularises variational optimisation problem.
Ashish and Adrian [12] applied anisotropic diffusion and local Hessian-based diffusion of viscous fluids to smooth images. Smoothing has also been achieved by interlacing diffusion with Gaussian kernel integration. In short, the physical laws that drive the Navier-Stokes equation after incorporating appropriate boundary conditions has been applied to image fields and accurate modelling of optical flows has thereby been achieved.
It would be interesting to see the random fluid flow produced by random fluctuation in the driving forces and to statistically characterise the random errors in the simulated image parameters due to the fluid. Our paper takes a step in this direction. We have incorporated the scheme suggested by Ashish [6,12] for control forces having small random noisy components and calculated the statistical correlations in the simulated velocity fields via a Monte-Carlo method. This would enable us to calculate the mean square error between the given optical flow field and the simulated flow. We have assessed the impact of noise in the control force on the mean square error energy.
The rest of the paper is organised as follows. In Sect. 2, the Navier-Stokes and incompressibilty equations are used to define velocity field of fluid. In Sect. 3, the energy due to random forcing component is minimised. Section 4 presents the non-random model for the forcing field. The estimation of parameter vector θ which minimises the average energy shift due to random forcing field δg is presented in Sect. 5. Section 6 includes the simulation results and concluding remarks.

Velocity field of fluid
In our work, optical flow has been modelled by representing the image intensity field E(t, x) as the fluid energy field and the image velocity field v(t, x) as the fluid velocity field. The rate of change of the energy calculated as the material derivative is given as The velocity field w(x, y) = u(x, y)x + v(x, y)ŷ of a fluid inside the square (x, y) [0, L] × [0, L] satisfies the Navier-Stokes equation [14].
where the density ρ = 1, ν = viscosity and w(x, y) stands for the velocity field in R 2 , (w.∇) is the operator ∂ y , p is the pressure field and f is the external forcing field. Apart from (1), we have the equation of continuity for the incompressible fluid.
The latter ensures the existence of a stream function ψ(x, y) so that or equivalently, By taking the curl of (4) and making use of (6a) and (6b), we arrive at the following equation for the stream function (7) and (8).
this is the same as where g = gẑ = ∇ × f . The derivation of the discretised form of curl operator has been mentioned in the appendix A.

Minimisation of energy function
Now, the standard method of estimating the velocity field w of a flow is [11] by minimising an energy function E of the form where Q(x, y) is a 6 × 6 positive definite matrix and λ(x, y) is a 6 × 1 vector ∀ (x,y) D. When the dynamics is given by (7), then u, v, ∇u, ∇v are expressed in terms of ψ x , ψ y , ψ x x , ψ xy , ψ yy and (9) then becomes a linear quadratic function of ψ. After integration by parts, this can be expressed as where Q 0 : R 2 × R 2 −→ R is a positive definite kernel. This linear quadratic form can be used to minimise the L 2 distance between the simulated fluid velocity field and the true optical flow. Alternatively, it can also be used to minimise a weighted linear combination of the fluid energy and its L 2 distance square from a desired optical flow. In this way, we would be able to simulate a given optical flow while keeping the simulated field energy small. The energy function (9) may equivalently be minimised w.r.t. g(x, y) by solving (7) for ψ. We shall now describe a spatially discretised version of this problem. Divide the square Here, ⊗ denotes the standard Kronecker tensor product. function e n R N has value 1 at (n + 1) th position and 0s at all the other positions for 0 n N − 1. Then, we can write (10) as (approximately since discretisation of the spatial integral has been applied) where Q 0 is an N 2 × N 2 positive definite matrix and is approximately given by and

Modelling of forcing field based on discretisation and optimisation
Now, we make use of the equation of the motion (7) and calculate g using ψ after minimising (10), i.e.
we can express (7) as This gives the value of g In this formalism, we have used a non-random model for the forcing field g and hence also for the velocity and stream function field. Here, the matrices A 1 , A 2 in (16) are calculated for (7) by first replacing ∇ 2 and (∇ 2 ) 2 by N 2 × N 2 matrices: writing n,m and likewise n,m ψ ,y (n , m )e n ⊗ e m = P y ψ, we can express (7) or (8) we have Comparing (16) and (24) gives

Calculating the random shift in the stream function
We assume that the forcing field g is given a small random perturbation δg with statistics depending on an unknown parameter vector θ . We then apply perturbation theory to the equation of random motion (16) and calculate the shift δψ in the stream function caused by the random shift δg in the forcing field g. Let g = g 0 +δg and ψ = ψ 0 + δψ. The unperturbed equation of motion is and the perturbed equation is Neglecting second degree of smaller terms, (26) approximates to The exact perturbed equation (26) can be solved to the next degree of smallness. Let δψ = δψ (0) + δψ (1) This gives the O( 0 ) term, so, thus, E(δψ (0) ) = 0.
Here, E stands for the mathematical expectation with respect to the probability distribution of the random force and perturbation δg or equivalently of δ f (i.e. δg = ∇ × δ f ).
Similarly, for the O( 1 ) term, where Now the average energy of the velocity field is Note: We have removed the perturbation tag because it is now clear that δψ is O( ), δψ (0) is O( ) and δψ (1) is O( 2 ). It should be noted that (34) is a generalisation of the form E(ψ 0 + δψ − ψ (0) ) 2 that represents energy in the error between the simulated flow and the true optical flow. More generally, (34) is also a generalisation of a linear quadratic form that represents a weighted linear combination of the error energy and simulated flow energy. Thus, it accurately represents a problem in which we desire to minimise the error while simultaneously retaining the minimality of the simulated flow energy. In (34), δψ = δψ (0) + δψ (1) , where δψ (0) is of the first order of smallness and δψ (1) is of the second order of smallness. In deriving (34), we retain terms only upto cubic order of smallness i.e. by noting that a product of two δψ (0) s is of second order as is δψ (1) , while a product of a δψ (0) and a δψ (1) is of cubic order, while a product of two δψ (0) s is of fourth order of smallness.

Ornstein-Uhlenbeck process
In turbulent flows, Ornstein-Uhlenbeck (OU) [13] processes are used to model velocity components of fluid particles. The OU process is a stochastic process that satisfies the following stochastic differential equation: where γ is relaxation rate and σ is white noise. The first term in the right-hand side is a drift term that reduces fluctuations in V(t) process. The second term adds stochastic fluctuations in V(t) process and, therefore, counteracts the first term.

Estimating the shift in average error energy for calculating θ
We evaluate the shift in the average energy function E caused by the random shift and calculate the parameter vector θ that minimises this average energy shift. Using the expressions (25) and (27) and the second and third moments E(δg ⊗δg) and E((δg ⊗δg).δg T ) depend on parameter θ (variance) and E −E 0 is minimised w.r.t θ to get the optimal flow.
Remark 2 E 0 stands for the energy in the absence of the random force perturbation.

Results and conclusion
Simulations were done in Matlab R2015 on a PC INTEL CORE I5 to generate the true optical flow V ox and V oy by minimising the distance between the simulated stream function and the true stream function. The true stream function corresponds to optical flow subjected to an energy constraint on the fluid flow. Such a simulated optical flow has small energy and closely resembles to a given optical flow (as shown in Figs. 1, 2). In Fig. 1, the graphs titled V ox , V x , V ox -V x and V oy , V y , V oy -V y are x and y components of true velocity field, simulated velocity field and error velocity field, respectively.

Noise to signal ratio
Let ψ 0 , ψ be the true and simulated flow stream function and likewise (V ox , V oy ) and (V x , V y ) are the true and simulated velocity fields. We present in Table 1 the corresponding noise to signal ratio (NSR)s, here and NSR = ||ψ−ψ 0 || 2 ||ψ 0 || 2 , obtained from true and simulated stream function in x and y direction, respectively. The true and simulated optical flow velocities are of the order 10 5 , while the error between these two flows is of order 10 −2 as shown in Fig. 1. In all cases, we get an NSR < 1 showing the goodness of the energy minimisation technique. As time progresses, the NSR decreases showing that the time-varying transient simulated flow velocity converges close to the true optical flow.

Analysis of random velocity perturbation
Let (V 0x , V 0y ) and (V x , V y ) denote the true velocity field of optical flow and simulated flow respectively. Our method essentially matches V x , V y to V 0x , V 0y subject to the constraint that V x , V y satisfies the Navier-Stokes equation. Figures 1 and 2 justify the proposed statement. The experiment is carried out over the force field. The force field is varied so that we get a velocity which approximates the true optical flow.
The random forcing field, g, is calculated by energy minimisation using the stream function ψ. We then add a small random force δg to g and simulate the correction to the stream function δψ using (34) and (35). The histogram for δg (Fig. 3) shows a nearby normal distribution which justifies our linearised perturbation approximation. The value of parameter θ (mean, variance) obtained from simulation is 11.6 and 9.6171, respectively. The autocorrelation plot for δg (Fig. 4) shows a exponentially decaying function with increasing spatial lag. This δg is a spatial version of the OU(Ornstein-Uhlenbeck) process. This confirms the presence of both damping factors arising from viscosity and fluctuating factors arising from random white noise in δg as in the OU process. Remark The simulation has been done based on the following computational steps: 1. Estimate the true stream function ψ 0 from measurements of the optical flow. 2. Calculate the force field in the Navier Stokes equation for approximating the true stream function ψ 0 by the simulated stream function ψ. To implement this step, we discretised the Navier-Stokes equation into spatial pixels and then optimise in the discrete domain. Optimisation yields the required driving force in the Navier-Stokes equation for optimal modelling. 3. Using perturbation theory outlined in the manuscript, evaluate the effect of a small random perturbation in the optimal force on the change in the modelled velocity field and the corresponding increase in the average modelling error energy.   Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A
The velocity field w after discretisation is a vector in R 2N 2 .
we can write w = C ψ, the matrix C is the discretised version of the curl operation Note that specifically, or equivalently,