CN104915968A - Navier-Stokes equation-based optical flow velocity estimation method - Google Patents

Abstract

The invention discloses a Navier-Stokes equation-based optical flow velocity estimation method. According to the method of the invention, a novel energy optimization function is constructed, and the Navier Stokes equation is transformed and is led into an energy equation, and the whole energy equation contains the following constraints in four aspects: consistency of an image, the smoothing of a velocity field, the coupling of the velocity field in the term of time and the incompressibility of the velocity field. Optical flow velocity estimation can be performed on a whole sequence simultaneously through utilizing the novel energy optimization equation provided by the invention, while, a large number of variables are needed to be solved, and in order to stably solve the energy optimization function, a variety of optimization methods such as coarse-to-fine multi-resolution optimization, alternating optimization, multiplicator alternating direction optimization and stable point optimization, are adopted.

Description

A kind of optical flow velocity method of estimation based on Navier Stokes equation

Technical field

The invention belongs to computer vision and field of Computer Graphics, specifically a kind of method of flow field velocity in estimated image or three-dimensional data, the method can be used for video tracking, and air motion is estimated and the aspect such as fluid analysis.

Background technology

Along with the develop rapidly of computer technology, people more and more focus on the authenticity of spontaneous phenomenon in computer animation, and in field of Computer Graphics, the simulation as cigarette, fire, water or other fluid is absolutely necessary.One class methods are simulated by physical equation completely, but modulation parameter is very loaded down with trivial details, and another kind of is utilize equipment to gather the fluid data in true environment, then it carried out to analysis and reuses.The present invention is a link of Equations of The Second Kind method, can be used for carrying out velocity field estimation to the image collected or volume data.

The motion of fluid can describe with famous Navier Stokes equation (Navier-Stokes Equation):

$\∂u/\∂t=-(u\·\▿)u-\▿p+f---\left(1\right)$

$\▿\·u=0---\left(2\right)$

Wherein u represents velocity field, and f represents external force, and p represents pressure, and t represents the time.And visible oBject in fluid (as dyestuff, or the particle etc. of cigarette, use I represents) can describe with following formula:

$\∂I/\∂t=-u\·\▿I---\left(3\right)$

Early stage optical flow analysis method only considers formula (3), the i.e. consistance coupling of image, this method (and afterwards on this basis improve one's methods) roughly can obtain velocity field, but this velocity field does not utilize any constraint of Navier Stokes equation, accurate not enough for fluid.Formula (2) was incorporated in the middle of the model of optical flow velocity estimation by Gregson afterwards, made the flow field obtained meet incompressible characteristic.This method is come, for estimating speed field, can obtain good result for use two consecutive frames, but the motion of fluid is have very strong coupling in time, and this coupling is described by formula (1).

Summary of the invention

For overcoming above-mentioned shortcoming, the object of the present invention is to provide a kind of optical flow velocity field method of estimation of global optimum, the method can carry out optical flow velocity estimation to whole sequence simultaneously, meet the characteristics of motion of fluid, and whole sequence also meets the characteristics of motion of fluid between adjacent two frames.

In order to achieve the above object, the present invention proposes a kind of optical flow velocity method of estimation based on Navier Stokes equation.First proposed a brand-new energy optimizing model based on Na Wei-Stokes, then employ very elaborate optimisation strategy and carry out velocity field and solve.Implementation method is: first, each frame structure multi-resolution pyramid of convection cell sequence, and the pyramidal number of plies is specified by user, but the speed that standard is lowest resolution image does not exceed a length in pixels; Then, each level of optimization from coarse to fine; In each level, optimize each frame speed field according to vertical sequence alternate, the iterations of alternative optimization is specified by user; When each frame speed field of alternative optimization, energy model corresponding for this frame speed place is split into three parts: image consistency and velocity field level and smooth, the time coupling of velocity field, the incompressibility of velocity field, with multiplier alternating direction optimization method to this three parts alternative optimization, iterations is specified by user.

Accompanying drawing explanation

Fig. 1 illustrates the main flow figure of the optical flow velocity method of estimation that the present invention is based on Navier Stokes equation;

Fig. 2 illustrates the algorithm of multiplier alternating direction optimization method (ADMM) of the present invention.

Embodiment

As shown in Figure 1, the optical flow velocity method of estimation based on Navier Stokes equation of the present invention adopts following steps:

(1) to each frame structure multi-resolution pyramid of image sequence.Take into account the accuracy of arithmetic speed and result, we used the zoom factor of 0.5 to construct pyramid, and the pyramidal number of plies is specified according to actual conditions by user, the standard of specifying allows the length of optical flow velocity field in the image of lowest resolution not exceed a pixel.

(2) successively secondaryly from lowest resolution to highest resolution to be optimized.Optimized first lowest resolution, starts in computing most, and the initial value of all velocity fields is set to 0, and comes for each frame speed field calculates initial value by the following formula of optimization:

E(u _i)＝E _data(u _i)+α ²E _sm(u _i)

E _datafor image consistency energy, concrete formula is:

$E_{\mathrm{data}}\left(u_i\right)={\∫}_{\mathrm{\Ω}}{\left|\right|I_i\left(x\right)-I_{i+1}(x+u_i)\left|\right|}_2^2\mathrm{d\Ω,}$

Locus is x, and integral domain is Ω, I _iit is the i-th two field picture.

E _smfor velocity field smoothed energy, concrete formula is α is weight coefficient, and wherein, N is sequence length, i.e. number of image frames, u _ibe the i-th frame speed field, alpha, gamma is weight coefficient.

Although away from this initial value also differs from optimum solution, it enough supports follow-up optimizing process.

(3) each frame speed field is alternately optimized.When optimization i-th frame speed field, the velocity field of other frames is all considered as known quantity, is formed specially for the energy theorem of the i-th frame speed field thus:

$E\left(u_i\right)=E_{\mathrm{data}}\left(u_i\right)+{\mathrm{\&alpha;}}^2{E}_{\mathrm{sm}}\left(u_i\right)+\left\{\begin{array}{c}{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{E}_{\mathrm{tem}}(u_i,u_{i+1}),i=1,\\ {\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{E}_{\mathrm{tem}}(u_{i-1},u_i),i=N-1,\\ {\mathrm{\&alpha;}}^2\mathrm{\&gamma;}(E_{\mathrm{tem}}(u_{i-1},u_i)+E_{\mathrm{tem}}(u_i,u_{i+1})),1<i<N-1,\end{array}\right.---\left(4\right)$

$\begin{array}{cc}s.t.& \&dtri;\&CenterDot;u_i=0\end{array}$

(4) formula (4) is out of shape, is divided into three parts and is optimized with multiplier alternating direction optimization method.First be that the constraints conversion of formula (4) is become unconfined formula, divergence constraint transfers energy term to then three parts are respectively:

F ₁(u _i)＝E _data(u _i)+α ²E _sm(u _i),

$F_2\left(u_i\right)=\left\{\begin{array}{c}{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{E}_{\mathrm{tem}}(u_i,u_{i+1}),i=1,\\ {\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{E}_{\mathrm{tem}}(u_{i-1},u_i),i=N-1,\\ {\mathrm{\&alpha;}}^2\mathrm{\&gamma;}(E_{\mathrm{tem}}(u_{i-1},u_i)+E_{\mathrm{tem}}(u_i,u_{i+1})),1<i<N-1,\end{array}\right.$

F ₃(u _i)＝ηE _div(u _i).

Wherein η is coefficient.Multiplier alternating direction optimization method can operate (proximal operator) and represents with near-end, see the algorithm in Fig. 2.

Near-end handling function is defined as ${\mathrm{prox}}_{\mathrm{\&lambda;F}}\left(z\right)={\arg \min }_{x}F\left(x\right)+\frac{\mathrm{\&lambda;}}{2}{\left|\right|z-x\left|\right|}_2^2$

A the optimization method of () Part I is ${\mathrm{prox}}_{{\mathrm{\&lambda;F}}_1}\left(b\right)={\arg \min }_{u_i}{F}_1\left(u_i\right)+\frac{\mathrm{\&lambda;}}{2}{\left|\right|b-u_i\left|\right|}_2^2\text{,}$ Here the x3-y1 in b corresponding diagram 2 in algorithm is corresponding after formula launches;

${\min }_{u_i}{\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|I_i\left(x\right)-I_{i+1}(x+u_i)\left|\right|}_2^2\mathrm{d\&Omega;}+\frac{\mathrm{\&lambda;}}{2}{\left|\right|b-u_{\text{i}}\left|\right|}_2^2$

First obtain Euler-Lagrange equation by the variational method, then use stable point alternative manner and first order Taylor to launch to solve.

B the optimization method of () Part II is ${\mathrm{prox}}_{{\mathrm{\&lambda;F}}_2}\left(c\right)={\arg \min }_{u_i}{F}_2\left(u_i\right)+\frac{\mathrm{\&lambda;}}{2}{\left|\right|c-u_i\left|\right|}_2^2\text{,}$ Here the x3-y2 in c respective figure 2 in algorithm is corresponding after formula launches;

${\min }_{u_i}{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|P(u_{i+1}(x+u_i)-f_i\left(x\right))-u_i\left(x\right)\left|\right|}_2^2\mathrm{d\&Omega;}+{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|P(u_{i-1}(x-u_i)+f_i\left(x\right))-u_i\left(x\right)\left|\right|}_2^2\mathrm{d\&Omega;}+\frac{\mathrm{\&lambda;}}{2}{\left|\right|x-c\left|\right|}_2^2$

Adopt the method for stable point iteration, during kth+1 iteration, its Euler-Lagrange equation is:

During i=1

$2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}(u_i^{k+1}\left(x\right)-P(u_{i+1}(x+u_i^k)-f_i\left(x\right)))+\mathrm{\&lambda;}(u_i^{k+1}-c)=0$

During i=N-1

$2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}(u_i^{k+1}\left(x\right)-P(u_{i-1}(x-u_i^k)+f_i\left(x\right)))+\mathrm{\&lambda;}(u_i^{k+1}-c)=0$

During 1<i<N-1

$\begin{array}{c}2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}(u_i^{k+1}\left(x\right)-P(u_{i+1}(x+u_i^k)-f_i\left(x\right)))+2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}(u_i^{k+1}\left(x\right)-P(u_{i-1}(x-u_i^k)+f_i\left(x\right)))\\ +\mathrm{\&lambda;}(u_i^{k+1}-c)=0\end{array}$

Computation sequence is ${\mathrm{warpU}}^{+}=u_{i+1}(x+u_i^k),\mathrm{warp}{U}^{-}=u_{i-1}(x-u_i^k),{\hat{u}}^{+}=\mathrm{warp}{U}^{+}-f_i,$ ${\hat{u}}^{-}=\mathrm{warp}{U}^{-}-f_i,u^{+}=P\left({\hat{u}}^{+}\right),u^{-}=P\left({\hat{u}}^{-}\right),$ Last optimum solution is:

During i=1, $u_i^{k+1}=\frac{2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{u}^{+}+\mathrm{\&lambda;c}}{2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}+\mathrm{\&lambda;}}$

During i=N-1, $u_i^{k+1}=\frac{2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{u}^{-}+\mathrm{\&lambda;c}}{2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}+\mathrm{\&lambda;}}$

During 1<i<N-1, $u_i^{k+1}=\frac{2{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}{\left(}^{u^{+}+u^{-}}+\mathrm{\&lambda;c}}{4{\mathrm{\&alpha;}}^2\mathrm{\&gamma;}+\mathrm{\&lambda;}}.$

C the optimization method of () Part III is ${\mathrm{prox}}_{{\mathrm{\&lambda;F}}_3}\left(d\right)={\arg \min }_{u_i}{F}_3\left(u_i\right)+\frac{\mathrm{\&lambda;}}{2}{\left|\right|d-u_i\left|\right|}_2^2\text{,}$ Here in d corresponding diagram 2 in algorithm when coefficient η is tending towards infinity (being now the target that formula (4) is optimized), its optimum solution is equivalent to u _i=P (d).

Claims (2)

1., based on an optical flow velocity method of estimation for Navier Stokes equation, it is characterized in that comprising the following steps:

(1) be each two field picture structure multi-resolution pyramid, optimize the optical flow velocity field of each level from coarse to finely;

(2) each frame speed field of whole sequence alternately being optimized, when optimizing a certain frame speed field, regarding the velocity field of other frames as constant;

(3) the current velocity field that will optimize, the energy theorem of its correspondence can turn three parts into, uses multiplier alternating direction optimization method (ADMM) to be optimized; Part I energy is image consistency energy and velocity field smoothed energy, and Part II energy is velocity field coupling energy in time, and Part III energy is the incompressibility energy of velocity field;

(4) stable point alternative manner is used to optimize Part I energy;

(5) stable point alternative manner is used to optimize Part II energy;

(6) Poisson equation is used to carry out velocity projections to optimize Part III energy.

2., as claimed in claim 1 based on the optical flow velocity method of estimation of Navier Stokes equation, it is characterized in that: adopt following light stream energy optimization model, be specially:

$E(u_1,u_2,...,u_{N-1})=\underset{i=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{E}_{\mathrm{data}}\left(u_i\right)+{\mathrm{\&alpha;}}^2(\underset{i}{\overset{N-1}{\mathrm{\&Sigma;}}}{E}_{\mathrm{sm}}\left(u_i\right)+\mathrm{\&gamma;}\underset{i=1}{\overset{N-2}{\mathrm{\&Sigma;}}}{E}_{\mathrm{tem}}(u_i,u_{i+1}))---\left(1\right)$

$s.t.\&ForAll;i\&Element;[0,N-1],\&dtri;\&CenterDot;u_i=0$

Wherein, N is sequence length, i.e. number of image frames, u _ibe the i-th frame speed field, alpha, gamma is weight coefficient;

E _datafor image consistency energy, concrete formula is:

$E_{\mathrm{data}}\left(u_i\right)={\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|I_i\left(x\right)-I_{i+1}(x+u_i)\left|\right|}_2^2\mathrm{d\&Omega;},$

I _ibe the i-th two field picture, locus is x, and integral domain is Ω,

E _smfor velocity field smoothed energy, concrete formula is

E _temfor time coupling energy, concrete formula is,

$E_{\mathrm{tem}}(u_i,u_{i+1})={\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|P(u_{i+1}(x+u_i)-f_i\left(x\right))-u_i\left(x\right)\left|\right|}_2^2\mathrm{d\&Omega;}---\left(2\right)$

$E_{\mathrm{tem}}(u_{i-1},u_i)={\&Integral;}_{\mathrm{\&Omega;}}{\left|\right|P(u_{i-1}(x-u_i)+f_i\left(x\right))-u_i\left(x\right)\left|\right|}_2^2\mathrm{d\&Omega;}---\left(3\right)$

Speed is thrown into incompressible subspace by P () expression, uses mode and the u of Poisson projection _out=P (u _in) be calculated as u _i+1(x+u _i), u _i-1(x-u _i) represent former frame, a rear frame speed field advection respectively after velocity field, f _irepresent the external force of the i-th frame.