CN103247058B - A kind of quick the Computation of Optical Flow based on error Distributed-tier grid - Google Patents

Abstract

The invention belongs to field of machine vision, disclose a kind of quick the Computation of Optical Flow based on error Distributed-tier grid, comprising: step one, input picture, build linear equation Ax=f; Step 2, sets up multi-layer image pyramid; Step 3, carries out front optimization, eliminates high fdrequency component; Step 4, carries out residual error transmission, eliminates low frequency component; Step 5, repeats step 3, four, until residual error is passed to most coarse layer; Step 7, carries out error passback from most coarse layer; Step 8, carries out refined net error correction; Step 9, carries out rear Optimized Iterative, improves stability of solution; Step 10, repeats step 7 ~ nine, until residual error is passed to most sub-layers.Method of the present invention is a kind of effective ways of accelerate equation Optimization Solution, and high frequency error convergence is very fast, and can significantly improve the computing velocity in visual light flow field, compared with the variational method, the speed of convergence of the method can bring up to more than 3.5 times.

Description

A kind of quick the Computation of Optical Flow based on error Distributed-tier grid

Technical field

The invention belongs to field of machine vision, relate to a kind of quick the Computation of Optical Flow based on error Distributed-tier grid.

Background technology

In visual motion analysis theory, when having a relative motion when between camera and scene objects, viewed brightness of image pattern further is referred to as light stream (optical flow).Because the light stream vector field between sequence video image have expressed the change of image, can matching relationship quantitatively between descriptor frame between corresponding image points, and provide image planes side-play amount, therefore for observer provides motion about target and structural information.

The variational method is the main stream approach solving light stream.Its guiding theory is level and smooth in gradient and the meaning of data fidelity builds energy functional, and emphasis obtains fine and close light stream vector field by the Euler-Lagrange equation (partial differential equation) solving optimized energy function.The advantage of variation light stream is, the versatility of model is good, and the degree of accuracy solved is high, and optical flow field is fine and close.Under variation framework, the method that Horn and Schunck proposes can obtain different optical flow computation models according to different assumed condition, its energy functional model comprises level and smooth item and data item, can mathematics expression be:

$\underset{u}{\min }\left\{\underset{\mathrm{\Ω}}{\∫}\mathrm{\α}({|\▿u|}^2+{|\▿v|}^2\mathrm{d\Ω}+\underset{\mathrm{\Ω}}{\∫}{(I_1(x+u\left(x\right))-I_0\left(x\right))}^2\mathrm{d\Ω})\right\}---\left(1\right)$

In formula, I ₀, I ₁represent adjacent two two field pictures of camera acquisition respectively, x=(x, y) ^trepresent certain picture point on image, with u (x)=(u (x), v (x)) ^trepresent the light stream vector of this position, wherein u (x) and v (x) is respectively the horizontal and vertical component of this vector.Symbol ▽ is gradient operator, and Ω represents imaging plane.Integration item is above level and smooth item, and integration item is below data item, and α is the constant coefficient of both adjustments weight.

By above formula gray scale conservation item I ₁(x+u (x))-I ₀x () carries out Taylor expansion, obtain:

I ₁(x+u(x))-I ₀(x)＝I _xu+I _yv+I _t(2)

In formula, x, y represent level and the vertical component of image I respectively, and t represents the time.Lower right corner subscript represents corresponding partial derivative, namely $I_x=\frac{\∂I}{\∂x},I_y=\frac{\∂I}{\∂y},I_t=\frac{\∂I}{\∂t}.$

Bring formula (2) into formula (1), above-mentioned energy functional is converted to:

$\underset{u}{\min }\{\underset{\mathrm{\Ω}}{\∫}({|\▿u|}^2+{|\▿u|}^2)\mathrm{d\Ω}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{(I_{x}u+I_{y}v+I_t)}^2\mathrm{d\Ω}\}$

Utilize variational method, obtaining its corresponding Euler-Lagrange equation is:

$\left\{\begin{array}{c}{I}_x^2\·u+I_x{I}_y\·v-\mathrm{\α\Δu}=-I_x{I}_t\\ I_x{I}_y\·u+I_y^2\·v-\mathrm{\α\Δv}=-I_y{I}_t\end{array}\right.---\left(3\right)$

In formula, u and v is respectively the horizontal and vertical component of light stream vector, and symbol Δ represents Laplace operator.After this system of equations discretize, if with traditional numerical method, as Gauss – Seidel, the methods such as SOR need iteration just can try to achieve comparatively ideal result thousands of times, and therefore the real-time of algorithm is poor.

In conventional numeric iterative algorithm, the convergence process of analytical error can find, high frequency error component is decayed rapidly in an iterative process, and the error component being in low frequency is decayed slowly.Therefore, the same numerical problem of rapid solving can be carried out with the multi-layer net that a series of resolution is different.This is because, refined net utilize less iterations can eliminate high frequency error; The low frequency aberration left over changes high frequency error into relative to coarse grid, will map on thicker grid, can continue on coarse grid, utilize less iterations to eliminate this part low frequency aberration, reach the object of accelerating convergence.Usual setting multi-layer net, transmits layer by layer by error, realizes the Distributed fusion of error.When specific implementation, formula (3) need be converted to system of linear equations.Formula (3) is out of shape a little, following form can be obtained:

β(x)＝f (4)

Wherein:

$\mathrm{\β}=\left(\begin{array}{cc}{I_x}^2& I_x{I}_y\\ I_x{I}_y& {I_y}^2\end{array}\right)+\left(\begin{array}{cc}-\mathrm{\α\Δ}& 0\\ 0& -\mathrm{\α\Δ}\end{array}\right)$

$f=\left(\begin{array}{c}-I_x{I}_t\\ -I_y{I}_t\end{array}\right)$

$x=\left(\begin{array}{c}u\\ v\end{array}\right)$

Visible, not only comprise constant coefficient in function β, also comprise nonlinear Laplace operator the calculating of optical flow field is caused can not directly to apply multi-layer net Algorithm for Solving thus.

Summary of the invention

According to the problems referred to above existed in optical flow computation method, the present invention proposes a kind of quick the Computation of Optical Flow based on error Distributed-tier grid, the speed of convergence of light stream iterative computation can be accelerated, improve the real-time of computer vision system.

Multi-layer net is a kind of effective ways of accelerate equation Optimization Solution.Wherein, equationof structure Ax=f is a primary step.Due in optical flow equation (4), containing nonlinear Laplace operator in function β, therefore function β is converted into not containing the coefficient matrices A of operator, to solve under multi-layer net algorithm frame by the present invention.

Based on a quick the Computation of Optical Flow for error Distributed-tier grid, it is characterized in that comprising the steps:

Step one, input picture, build linear equation Ax=f, method is as follows:

1) calculate spatio-temporal gradient tensor matrix, solving equations discretize will be treated;

Suppose in a certain two field picture plane, (i, j) represents the position of the i-th row jth row pixel, u _i,jand v _i,jrepresent the displacement field in this horizontal x-axis occurred between consecutive frame and longitudinal y-axis respectively.The space partial derivative f of computed image _x, f _yand time partial derivative f _t, the spatio-temporal gradient tensor J matrix of each pixel in image can be obtained:

$J=\left(\begin{array}{ccc}{f}_x{f}_x& f_x{f}_y& f_x{f}_t\\ f_y{f}_x& f_y{f}_y& f_y{f}_t\\ f_t{f}_x& f_t{f}_y& f_t{f}_t\end{array}\right)=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right)---\left(5\right)$

If with [J] _i,jrepresent the value of J matrix at pixel (i, j) place, by h representation space step-length, then Euler-Lagrange equation (formula (3)) can discretely be:

${\left[J_{11}\right]}_{i,j}{u}_{i,j}+{\left[J_{12}\right]}_{i,j}{v}_{i,j}-\frac{\mathrm{\α}}{h^2}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j})=-{\left[J_{13}\right]}_{i,j}---\left(\text{6}\right)$

${\left[J_{12}\right]}_{i,j}{u}_{i,j}+{\left[J_{22}\right]}_{i,j}{v}_{i,j}-\frac{\mathrm{\α}}{h^2}(v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-4v_{i,j})=-{\left[J_{23}\right]}_{i,j}---\left(\text{7}\right)$

2) difference equation (6) is changed into linear constant coefficient system of equations A ₁x=f ₁;

Abbreviation equation (6), obtains the DIFFERENCE EQUATIONS in the D of domain:

$(1+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{i,j})u_{i,j}-\frac{1}{4}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})+\frac{h^2}{4\mathrm{\α}}\left[J_{12}\right]v_{i,j}=-\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{\text{i,j}}$

Always total N+1 is capable to suppose image, and remove zero row and N-th row, the DIFFERENCE EQUATIONS of each pixel of first row has following form:

$\left(\begin{array}{ccccccc}1+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{11}& -\frac{1}{4}& 0& ...& 0& 0& 0\\ -\frac{1}{4}& 1+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{21}& -\frac{1}{4}& ...& 0& 0& 0\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& -\frac{1}{4}& 1+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{(N-2)1}& -\frac{1}{4}\\ 0& 0& 0& ...& 0& -\frac{1}{4}& 1+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{(N-1)1}\end{array}\right)\left(\begin{array}{c}{u}_{11}\\ u_{21}\\ ...\\ u_{(N-1)1}\end{array}\right)$

$-\frac{1}{4}\left(\begin{array}{cccc}1& 0& ...& 0\\ 0& 1& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& 1\end{array}\right)\left(\begin{array}{c}{u}_{12}\\ u_{22}\\ ...\\ u_{(N-1)2}\end{array}\right)+\left(\begin{array}{cccc}\frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_{11}& 0& ...& 0\\ 0& \frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_{21}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& \frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_{(N-1)1}\end{array}\right)\left(\begin{array}{c}{v}_{11}\\ v_{21}\\ ...\\ v_{(N-1)1}\end{array}\right)$

$=\left(\begin{array}{c}-\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{11}+\frac{1}{4}{g}_{01}+\frac{1}{4}{g}_{10}\\ -\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{21}+\frac{1}{4}{g}_{20}\\ ...\\ -\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{(N-1)1}+\frac{1}{4}{g}_{N1}+\frac{1}{4}{g}_{(N-1)0}\end{array}\right)=b_1---\left(8\right)$

In formula, g _ijrepresent the borderline element u of domain D _i, i.e. g _ij=u _ij.

If introduce optical flow components u _ijand v _ijon column vector:

$u_i=\left\{\begin{array}{c}{u}_{1i}\\ u_{2i}\\ ...\\ u_{(N-1)i}\end{array}\right\},v_i=\left\{\begin{array}{c}{v}_{1i}\\ v_{2i}\\ ...\\ v_{(N-1)i}\end{array}\right\}$

And order:

${\left[J_{\mathrm{mn}}\right]}_i=\left(\begin{array}{cccc}{\left[J_{\mathrm{mn}}\right]}_{1i}& 0& ...& 0\\ 0& {\left[J_{\mathrm{mn}}\right]}_{2i}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\left[J_{\mathrm{mn}}\right]}_{(N-1)i}\end{array}\right)---\left(9\right)$

For the J of the J matrix element at pixel (m, n) place _mn, matrix [J in formula (9) _mn] _ibe equivalent to unit matrix and the J on N-1 rank _mnthe product of i-th this column vector of row.Therefore, the DIFFERENCE EQUATIONS (8) of each pixel of first row can turn to:

$(G+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_1)u_1-\frac{1}{4}I{u}_2+\frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_1{v}_1=b_1$

In formula, I is the unit matrix on N-1 rank, b ₁for the column vector of definition in formula (6), G is N-1 rank matrix below:

$G=\left(\begin{array}{ccccccc}1& -\frac{1}{4}& 0& ...& 0& 0& 0\\ -\frac{1}{4}& 1& -\frac{1}{4}& ...& 0& 0& 0\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& -\frac{1}{4}& 1& -\frac{1}{4}\\ 0& 0& 0& ...& 0& -\frac{1}{4}& 1\end{array}\right)$

In like manner, the difference equation of each pixel of secondary series is:

$-\frac{1}{4}I{u}_1+(G+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_2)u_2-\frac{1}{4}I{u}_3+\frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_2{v}_2=b_2$

In formula:

$b_2=\left\{\begin{array}{c}-\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{12}+\frac{1}{4}{g}_{02}\\ -\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{22}\\ ......\\ -\frac{h^2}{4\mathrm{\α}}{\left[J_{13}\right]}_{(N-1)2}+\frac{1}{4}{g}_{N2}\end{array}\right\}$

All transformed by all for image row, then system of equations can be written as A ₁x=f ₁form, in formula:

$A_1=\left(\begin{array}{cccccccccccc}G+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_1& -\frac{1}{4}I& 0& ...& 0& 0& 0& \frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_1& 0& ...& 0& 0\\ -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_2& -\frac{1}{4}I& ...& 0& 0& 0& 0& \frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_2& ...& 0& 0\\ ...& ...& ...& ...& ...& ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{N-2}& -\frac{1}{4}I& 0& 0& ...& \frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_{N-2}& 0\\ 0& 0& 0& ...& 0& -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{\left[J_{11}\right]}_{N-1}& 0& 0& ...& 0& \frac{h^2}{4\mathrm{\α}}{\left[J_{21}\right]}_{N-1}\end{array}\right)$

$x=\left(\begin{array}{c}{u}_1\\ u_2\\ ...\\ u_{N-1}\\ v_1\\ v_2\\ ...\\ v_{N-1}\end{array}\right),f_1=\left(\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_{N-1}\end{array}\right)$

3) difference equation (7) is changed into linear constant coefficient system of equations A ₂x=f ₂;

4) by A ₁x=f ₁with A ₂x=f ₂common formation system of equations:

Ax＝f (10)

Step 2, sets up multi-layer image pyramid;

The pyramidal schematic diagram of multi-layer image as shown in Figure 2.A uppermost tomographic image is former image in different resolution, each tomographic image below represent successively reduce resolution images falls.Because pixel itself is uniform discrete, therefore can multi-layer image be regarded as multi-layer net.

Usually, the size (long and wide) of image is between 300-800.According to the scale factor η (0.5< η <0.95) of setting, former image in different resolution is successively reduced, and size is rounded, obtain multi-layer image pyramid.Usually, number of plies N=4 ~ 5 layer.The initial value of every layer of light stream is set to zero.

Step 3, carries out front optimization, eliminates high fdrequency component;

Current layer i (i=1,2 ..., N-1), take null matrix as light stream initial value to refined net Equation Iterative m time, to eliminate high fdrequency component.If x _ifor initial value, obtain approximate evaluation value be designated as:

${\stackrel{\&OverBar;}{x}}_i={\mathrm{Relax}}^m(x_i,A_i,r_i)$

In formula, Relax ^mrepresent Gauss-Seidel iterative process, m is the iterations of setting.

A _idetermined by formula (10); As i=1, r is determined by formula (10), i.e. r ₁=f.Otherwise r equals the residual error that last layer is passed to this layer.The corresponding residual error of this layer is updated to:

$r_i=A_i{\stackrel{\&OverBar;}{x}}_i-f_i$

Step 4, carries out residual error transmission, eliminates low frequency component;

Front optimizing process is intended to eliminate high frequency error component, and inherited error is mainly low frequency component.Therefore, residual error is passed on coarse grid, to eliminate low frequency component.Be restricted to the surplus r on coarse grid _i+1for:

$r_{i+1}=I_i^{i+1}{r}_i$

In formula, for refined net is to the mapping operator on coarse grid.

Step 5, as i=1 ~ N-1, repeats step 3, four, until residual error is passed to most coarse layer;

Step 6, solves the system of equations on coarse grid:

A _Ne _N＝r _N

In formula, e _nfor the light stream error of Exact Solution in most coarse layer.

Obtained by matrix operation:

e _N＝(A _N) ^-1r _N

Step 7, carries out error passback from most coarse layer (j=N);

Current layer j, passes back to last layer comparatively on refined net j-1, that is: by the amount of error correction of trying to achieve

$e_{j-1}=I_j^{j-1}{e}_j$

In formula, for coarse grid is to the mapping operator on refined net.

Step 8, carries out refined net error correction;

The initial value of calculating is added the error passed back by coarse grid, obtain the solution after refined net correction that is:

${\hat{x}}_{j-1}={\stackrel{\&OverBar;}{x}}_{j-1}+e_{j-1}$

Step 9, carries out rear Optimized Iterative, improves stability of solution;

After refined net corrects, with the solution after correcting for initial value, Optimized Iterative n time after performing, obtains the corrected value after renewal its expression formula is:

$x_{j-1}^{\mathrm{new}}={\mathrm{Relax}}^n({\hat{x}}_{j-1},A_{j-1},f_{j-1})$

Step 10, works as j=N, N-1 ..., when 2, repeat step 7 ~ nine, until residual error is passed to most sub-layers.

Compared with prior art, the present invention has following beneficial effect:

Light stream solving equation launches by row by the present invention, is turned to the form of system of linear equations, makes the complicated calculations process of light stream can accelerate to realize by linear multi-layer trellis algorithm.Secondly, have employed the multi-layer net of different resolution, the error component of different frequency is distributed on different grids, make error become high frequency error relative to this layer of grid.Very fast due to high frequency error convergence, thus the method can significantly improve the computing velocity in visual light flow field, and compared with the variational method, the speed of convergence of the method can bring up to more than 3.5 times.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of method involved in the present invention;

Fig. 2 is multi-layer image pyramid schematic diagram;

Fig. 3 is first group of test pattern RubberWhale of application example of the present invention;

Fig. 4 is first group of test pattern of application example of the present invention.

Embodiment

Below in conjunction with the drawings and specific embodiments, the present invention will be further described.

The microcomputer one implementing to rely on P3800M dominant frequency, 512M internal memory and above configuration of the present invention, as hardware device platform.Running environment is Windows XP operating system and Matlab software platform.

The detailed process of present embodiment is, first according to optical flow estimation to be solved, and structure system of linear equations Ax=f.Then set up 4 ~ 5 tomographic image pyramids, be optimized iteration respectively at every one deck, and iteration residual error is successively mapped on the thicker grid of lower one deck.In most coarse layer, solve the system of equations on most coarse grid, obtain accurate estimation of error.After this, return error in opposite direction, thinner grid successively carries out error correction.Fig. 1 is the process flow diagram of method involved in the present invention, specifically comprises the following steps:

Step one, input picture, builds linear equation Ax=f;

Step 2, sets up multi-layer image pyramid;

Step 3, carries out front optimization, eliminates high fdrequency component;

Step 4, carries out residual error transmission, eliminates low frequency component;

Step 5, as i=1 ~ N-1, repeats step 3, four, until residual error is passed to most coarse layer;

Step 6, solves the system of equations on coarse grid: A _ne _n=r _n

Step 7, carries out error passback from most coarse layer (j=N);

Step 8, carries out refined net error correction;

Step 9, carries out rear Optimized Iterative, improves stability of solution;

Step 10, works as j=N, N-1 ..., when 2, repeat step 7 ~ nine, until residual error is passed to most sub-layers.

Provide an application example of the present invention below.

Two groups of test patterns are selected to verify the fast light flow field algorithm based on error Distributed-tier grid that the present invention proposes.One group is a pair test pattern RubberWhale in the Middlebury java standard library extensively adopted in the world, as shown in Figure 3; Another group is scientific paper " the An improved algorithm for TV-L that the people such as Wedel deliver at them ¹optical flow " middle a pair test pattern adopted, as shown in Figure 4.

During experiment, the fast light flow field algorithm based on error Distributed-tier grid adopting traditional variational algorithm and the present invention to propose respectively calculate often organize image between optical flow field.By the working time of software statistics two kinds of methods, experimental result is as shown in table 1.In table 1, experiment one adopts image shown in Fig. 3, and experiment two adopts image shown in Fig. 4.As shown in Table 1, the operation time that the fast light flow field algorithm based on error Distributed-tier grid compares traditional variational algorithm shortens greatly, and the speed of convergence of two experiments all brings up to more than 3.5 times, significantly improves counting yield.

The contrast of the traditional variational algorithm of table 1 and the method for the invention working time

Claims (1)

1., based on a quick the Computation of Optical Flow for error Distributed-tier grid, it is characterized in that comprising the steps:

Step one, input picture, build linear equation Ax=f, method is as follows:

1) calculate spatio-temporal gradient tensor matrix, solving equations discretize will be treated;

Suppose in a certain two field picture plane, (i, j) represents the position of the i-th row jth row pixel, u _i,jand v _i,jrepresent the displacement field in this horizontal x-axis occurred between consecutive frame and longitudinal y-axis respectively; The space partial derivative f of computed image _x, f _yand time partial derivative f _t, the spatio-temporal gradient tensor J matrix of each pixel in image can be obtained:

$J=\left(\begin{array}{ccc}{f}_x{f}_x& f_x{f}_y& f_x{f}_t\\ f_y{f}_x& f_y{f}_y& f_y{f}_t\\ f_t{f}_x& f_t{f}_y& f_t{f}_t\end{array}\right)=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right)---\left(1\right)$

If with [J] _i,jrepresent the value of J matrix at pixel (i, j) place, by h representation space step-length, the Euler-Lagrange equation variational method obtained by energy functional number is discrete is:

${\left[J_{11}\right]}_{i,j}{u}_{i,j}+{\left[J_{12}\right]}_{i,j}{v}_{i,j}-\frac{\mathrm{\&alpha;}}{h^2}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j})=-{\left[J_{13}\right]}_{i,j}---\left(2\right)$

${\left[J_{12}\right]}_{i,j}{u}_{i,j}+{\left[J_{22}\right]}_{i,j}{v}_{i,j}-\frac{\mathrm{\&alpha;}}{h^2}(v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-4v_{i,j})=-{\left[J_{23}\right]}_{i,j}---\left(3\right)$

In formula, α is for regulating the coefficient of both level and smooth item and data item weight;

2) difference equation (2) is changed into linear constant coefficient system of equations A ₁x=f ₁;

Abbreviation equation (2), obtains the DIFFERENCE EQUATIONS in the D of domain:

$(1+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{i,j})u_{i,j}-\frac{1}{4}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})+\frac{h^2}{4\mathrm{\&alpha;}}\left[J_{12}\right]v_{i,j}=-\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{i,j}$

Always total N+1 is capable to suppose image, and remove zero row and N-th row, the DIFFERENCE EQUATIONS of each pixel of first row has following form:

$\left(\begin{array}{ccccccc}1+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{11}& -\frac{1}{4}& 0& ...& 0& 0& 0\\ -\frac{1}{4}& 1+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{21}& -\frac{1}{4}& ...& 0& 0& 0\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& -\frac{1}{4}& 1+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{(N-2)1}& -\frac{1}{4}\\ 0& 0& 0& ...& 0& -\frac{1}{4}& 1+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{(N-1)1}\end{array}\right)\left(\begin{array}{c}{u}_{11}\\ u_{21}\\ ...\\ u_{(N-1)1}\end{array}\right)$

$\begin{array}{c}-\frac{1}{4}\left(\begin{array}{cccc}1& 0& ...& 0\\ 0& 1& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& 1\end{array}\right)\left(\begin{array}{c}{u}_{12}\\ u_{22}\\ ...\\ u_{(N-1)2}\end{array}\right)+\left(\begin{array}{cccc}\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_{11}& 0& ...& 0\\ 0& \frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_{21}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& \frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_{(N-1)1}\end{array}\right)\left(\begin{array}{c}{v}_{11}\\ v_{21}\\ ...\\ v_{(N-1)1}\end{array}\right)\\ =\begin{array}{}\end{array}\left(\begin{array}{c}-\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{11}+\frac{1}{4}{g}_{01}+\frac{1}{4}{g}_{10}\\ -\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{21}+\frac{1}{4}{g}_{20}\\ ...\\ -\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{(N-1)1}+\frac{1}{4}{g}_{N1}+\frac{1}{4}{g}_{(N-1)0}\end{array}\right)=b_1\end{array}---\left(4\right)$

In formula, g _ijrepresent the borderline element u of domain D _i, i.e. g _ij=u _ij;

Introduce optical flow components u _ijand v _ijon column vector:

$u_i=\left\{\begin{array}{c}{u}_{1i}\\ u_{2i}\\ ...\\ u_{(N-1)i}\end{array}\right\},v_i=\left\{\begin{array}{c}{v}_{1i}\\ v_{2i}\\ ...\\ v_{(N-1)i}\end{array}\right\}$

And order:

${\left[J_{\mathrm{mn}}\right]}_i=\left(\begin{array}{cccc}{\left[J_{\mathrm{mn}}\right]}_{1i}& 0& ...& 0\\ 0& {\left[J_{\mathrm{mn}}\right]}_{2i}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ...& {\left[J_{\mathrm{mn}}\right]}_{(N-1)i}\end{array}\right)---\left(5\right)$

For the J of the J matrix element at pixel (m, n) place _mn, matrix [J in formula (5) _mn] _ibe equivalent to unit matrix and the J on N-1 rank _mnthe product of i-th this column vector of row; Therefore, the DIFFERENCE EQUATIONS (4) of each pixel of first row can turn to:

$(G+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_1)u_1-\frac{1}{4}{\mathrm{Iu}}_2+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_1{v}_1=b_1$

In formula, I is N-1 rank unit matrix, b ₁for the column vector of definition in formula (2), G is N-1 rank matrix below:

$G=\left(\begin{array}{ccccccc}1& -\frac{1}{4}& 0& ...& 0& 0& 0\\ -\frac{1}{4}& 1& -\frac{1}{4}& ...& 0& 0& 0\\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& -\frac{1}{4}& 1& -\frac{1}{4}\\ 0& 0& 0& ...& 0& -\frac{1}{4}& 1\end{array}\right)$

The difference equation of each pixel of secondary series is:

$-\frac{1}{4}{\mathrm{Iu}}_1+(G+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_2)u_2-\frac{1}{4}{\mathrm{Iu}}_3+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_2{v}_2=b_2$

In formula:

$b_2=\left\{\begin{array}{c}-\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{12}+\frac{1}{4}{g}_{02}\\ -\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{22}\\ ......\\ -\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{13}\right]}_{(N-1)2}+\frac{1}{4}{g}_{N2}\end{array}\right\}$

All transformed by all for image row, then system of equations can be written as A ₁x=f ₁form, in formula:

$A_1=\left(\begin{array}{cccccccccccc}G+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_1& -\frac{1}{4}I& 0& ...& 0& 0& 0& \frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_1& 0& ...& 0& 0\\ -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_2& -\frac{1}{4}I& ...& 0& 0& 0& 0& \frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_2& ...& 0& 0\\ ...& ...& ...& ...& ...& ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{N-2}& -\frac{1}{4}I& 0& 0& ...& \frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_{N-2}& 0\\ 0& 0& 0& ...& 0& -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{11}\right]}_{N-1}& 0& 0& ...& 0& \frac{h^2}{4\mathrm{\&alpha;}}{\left[J_{21}\right]}_{N-1}\end{array}\right)$

$x=\left(\begin{array}{c}{u}_1\\ u_2\\ ...\\ u_{N-1}\\ v_1\\ v_2\\ ...\\ v_{N-1}\end{array}\right),f_1=\left(\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_{N-1}\end{array}\right)$

3) difference equation (3) is changed into linear constant coefficient system of equations A ₂x=f ₂;

4) by A ₁x=f ₁with A ₂x=f ₂common formation system of equations:

Ax＝f (6)

Step 2, set up multi-layer image pyramid, method is as follows:

The uppermost tomographic image of multi-layer image pyramid is former image in different resolution, each tomographic image below represent successively reduce resolution images falls; Because pixel itself is uniform discrete, therefore can multi-layer image be regarded as multi-layer net;

Usually, the length of image and wide between 300 ~ 800; According to the scale factor η of setting, 0.5< η <0.95, successively reduces former image in different resolution, and size is rounded, and obtains multi-layer image pyramid; Usually, number of plies N=4 ~ 5 layer, the initial value of every layer of light stream is set to zero;

Step 3, carries out front optimization, and eliminate high fdrequency component, method is as follows:

At current layer i, i=1,2 ..., N-1 take null matrix as light stream initial value to refined net Equation Iterative m time, to eliminate high fdrequency component; If x _ifor initial value, obtain approximate evaluation value be designated as:

${\stackrel{\&OverBar;}{x}}_i={\mathrm{Relax}}^m(x_i,A_i,r_i)$

In formula, Relax ^mrepresent Gauss-Seidel iterative process, m is the iterations of setting;

A _idetermined by formula (6); As i=1, r is determined by formula (6), i.e. r ₁=f; Otherwise r equals the residual error that last layer is passed to this layer; The corresponding residual error of this layer is updated to:

$r_i=A_i{\stackrel{\&OverBar;}{x}}_i-f_i$

Step 4, carries out residual error transmission, and eliminate low frequency component, method is as follows:

Front optimizing process is intended to eliminate high frequency error component, and inherited error is mainly low frequency component, is therefore passed on coarse grid by residual error, to eliminate low frequency component; Be restricted to the surplus r on coarse grid _i+1for:

$r_{i+1}=I_i^{i+1}{r}_i$

In formula, for refined net is to the mapping operator on coarse grid;

Step 5, as i=1 ~ N-1, repeats step 3, four, until residual error is passed to most coarse layer;

Step 6, solves the system of equations on coarse grid:

A _Ne _N＝r _N

In formula, e _nfor the light stream error of Exact Solution in most coarse layer;

Obtained by matrix operation:

e _N＝(A _N) ^-1r _N

Step 7, from most coarse layer, carry out error passback, method is as follows:

Current layer is j, the amount of error correction of trying to achieve is passed back to last layer comparatively on refined net j-1, that is:

$e_{j-1}=I_j^{j-1}{e}_j$

In formula, for coarse grid is to the mapping operator on refined net;

Step 8, carry out refined net error correction, method is as follows:

The initial value of calculating is added the error passed back by coarse grid, obtain the solution after refined net correction that is:

${\hat{x}}_{j-1}={\stackrel{\&OverBar;}{x}}_{j-1}+e_{j-1}$

Step 9, carries out rear Optimized Iterative, and improve stability of solution, method is as follows:

After refined net corrects, with the solution after correcting for initial value, Optimized Iterative n time after performing, obtains the corrected value after renewal its expression formula is:

$x_{j-1}^{\mathrm{new}}={\mathrm{Relax}}^n({\hat{x}}_{j-1},A_{j-1},f_{j-1})$

Step 10, works as j=N, N-1 ..., when 2, repeat step 7 ~ nine, until residual error is passed to most sub-layers.