CN103247058A - Fast optical flow field calculation method based on error-distributed multilayer grid - Google Patents

Abstract

The invention belongs to the field of machine vision, and discloses a fast optical flow field calculation method based on an error-distributed multilayer grid. The method comprises the following steps: step one, inputting an image, and constructing a linear equation, that is Ax is equal to f; step two, establishing a multilayer image pyramid; step three, performing prior optimization, and eliminating high-frequency components; step four, performing residual error transmission, and eliminating low-frequency components; step five, repeating steps three and four until residual errors are transmitted to the thickest layer; step six, solving equation sets on thick girds; step seven, performing error passing back from the thickest layer; step eight, performing error correction on thin grids; step nine, performing subsequent optimization iteration, and improving stability of solutions; and step ten, repeating steps seven, eight and nine until the residual errors are transmitted to the thinnest layer. The method provided by the invention is an effective method for accelerating optimization solution of equations, can quick converge high-frequency errors, and can remarkably improve computing speed of a visual optical flow field; and compared with a variational method, the convergence rate of the method can be improved above 3.5 times.

Description

A kind of quick optical flow field computing method based on the distributed multi-layer net of error

Technical field

The invention belongs to field of machine vision, relate to a kind of quick optical flow field computing method based on the distributed multi-layer net of error.

Background technology

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

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

$\underset{u}{\min }\{\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\Ω}\}---\left(1\right)$

In the formula, I ₀, I ₁Adjacent two two field pictures of representing camera acquisition respectively, and x=(x, y) ^TCertain picture point on the presentation video is with u (x)=(u (x), v (x)) ^TThe light stream vector of representing this position, wherein u (x) and v (x) are respectively level and the vertical component of this vector.Symbol ▽ is gradient operator, and Ω represents imaging plane.The integration item of front is level and smooth, and the integration item of back is data item, and α is for regulating the constant coefficient of the two weight.

With following formula gray scale conservation item I ₁(x+u (x))-I ₀(x) carry out Taylor expansion, obtain:

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

In the formula, x, y be level and the vertical component of presentation video I respectively, and t represents the time.Lower right corner subscript is represented 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+{|\▿v|}^2)\mathrm{d\Ω}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}{(I_{x}u+I_{y}v+I_t)}^2\mathrm{d\Ω}\}$

Utilize variational method, obtain its corresponding Euler-Lagrange equation and be:

$\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 the formula, u and v are respectively level and the vertical component of light stream vector, and the symbol Δ is represented Laplace operator.After this system of equations discretize, if with traditional numerical solution method, as Gauss – Seidel, methods such as SOR need iteration just can try to achieve comparatively ideal result thousands of times, so the real-time of algorithm is relatively poor.

In traditional numerical value iterative algorithm, the convergence process of analytical error can find, the high frequency error component is decay rapidly in iterative process, and the error component that is in low frequency is decayed slowly.Therefore, the different multi-layer net of available a series of resolution comes the same numerical problem of rapid solving.This is because utilize less iterations can eliminate high frequency error at refined net; The low frequency aberration of leaving over changes high frequency error into respect to coarse grid, will map on the thicker grid, can continue to utilize less iterations to eliminate this part low frequency aberration at coarse grid, has reached the purpose of accelerating convergence.Usually set multi-layer net, error is transmitted layer by layer, realize the distributed estimation of error.When specific implementation, need formula (3) is converted to system of linear equations.Formula (3) is out of shape a little, can obtains following form:

β(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)$

As seen, not only comprise constant coefficient among the function β, also comprise nonlinear Laplace operator

Causing the calculating of optical flow field can not directly use the multi-layer net algorithm thus finds the solution.

Summary of the invention

According to the problems referred to above that exist in the optical flow computation method, the present invention proposes a kind of quick optical flow field computing method based on the distributed multi-layer net of error, can accelerate the light stream convergence speed of iteration calculation, improve the real-time of computer vision system.

Multi-layer net is a kind of effective ways that equation optimization is found the solution that accelerate.Wherein, equationof structure Ax=f is a primary step.Because in optical flow equation (4), contain nonlinear Laplace operator among the function β, so the present invention is converted into the coefficient matrices A that does not contain operator with function β, in order under the multi-layer net algorithm frame, find the solution.

A kind of quick optical flow field computing method based on the distributed multi-layer net of error is characterized in that comprising the steps:

Step 1, input picture makes up linear equation Ax=f, and method is as follows:

(1) calculates space-time gradient tensor matrix, will wait the group discretize of solving an equation;

Suppose on a certain two field picture plane (i, the j) position of the capable j row of expression i pixel, u _{I, j}And v _{I, j}Represent this displacement field on the horizontal x axle that takes place between consecutive frame and vertical y axle respectively.The space partial derivative f of computed image _x, f _yAnd time partial derivative f _t, can obtain the space-time gradient tensor J matrix of each pixel in the image:

$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, j}Expression J matrix pixel (with h representation space step-length, then can disperse is Euler-Lagrange equation (formula (3)) for i, the value of j) locating:

${\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(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(7\right)$

(2) difference equation (6) is changed into linear equation with constant coefficient group A ₁X=f ₁

Abbreviation equation (6) obtains finding the solution the DIFFERENCE EQUATIONS in the regional D:

$(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]}_{i,j}$

Always total N+1 is capable to suppose image, and it is capable to remove the 0th row and N, and the DIFFERENCE EQUATIONS of first each pixel of 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{\α}}{{[J}_{21}]}_{(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{\α}}{{[J}_{13}]}_{11}+\frac{1}{4}{g}_{01}+\frac{1}{4}{g}_{10}\\ -\frac{h^2}{4\mathrm{\α}}{{[J}_{13}]}_{21}+\frac{1}{4}{g}_{20}\\ \·\·\·\\ \frac{h^2}{4\mathrm{\α}}{{[J}_{13}]}_{(N-1)1}+\frac{1}{4}{g}_{N1}+\frac{1}{4}{g}_{(N-1)0}\end{array}\right)=b_1---\left(8\right)$

In the formula, g _IjThe borderline element u of regional D is found the solution in expression _i, i.e. g _Ij=u _Ij

If introduce light stream component u _IjAnd v _IjOn column vector:

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

And order:

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

For pixel (m, the J of the J matrix element of n) locating _Mn, matrix [J in the formula (9) _Mn] _iThe unit matrix and the J that are equivalent to the N-1 rank _MnI is listed as the product of this column vector.Therefore, the DIFFERENCE EQUATIONS (8) of first each pixel of row can turn to:

$(G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_1)u_1-\frac{1}{4}{\mathrm{Iu}}_2+\frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_1{v}_1=b_1$

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

$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}{\mathrm{Iu}}_1+(G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_2)u_2-\frac{1}{4}{\mathrm{Iu}}_3+\frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_2{v}_2=b_2$

In the formula:

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

All row of image are all transformed, and then system of equations can be written as A ₁X=f ₁Form, in the formula:

$A_1=\left(\begin{array}{cccccccccccc}G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_1& \frac{1}{4}I& 0& \·\·\·& 0& 0& 0& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_1& 0& \·\·\·& 0& 0\\ \frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_2& \frac{1}{4}I& \·\·\·& 0& 0& 0& 0& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_2& \·\·\·& 0& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& \·\·\·& \frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_{N-2}& \frac{1}{4}I& 0& 0& \·\·\·& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_{N-2}& 0\\ 0& 0& 0& \·\·\·& 0& \frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_{N-1}& 0& 0& \·\·\·& 0& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_{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 equation with constant coefficient group A ₂X=f ₂

(4) by A ₁X=f ₁With A ₂X=f ₂Common formation system of equations:

Ax=f (10)

Step 2 is set up the multi-layer image pyramid;

The pyramidal synoptic diagram of multi-layer image as shown in Figure 2.A uppermost tomographic image is former image in different resolution, and each following tomographic image is represented successively to dwindle falls the resolution image.Because pixel itself is uniform discrete, therefore can regard multi-layer image as multi-layer net.

Usually, size of images (long and wide) is between 300-800.According to the scale factor η (0.5＜η＜0.95) that sets former image in different resolution is successively dwindled, and size is rounded, obtain the 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 is carried out preceding optimization, eliminates high fdrequency component;

When anterior layer i (i=1,2 ..., N-1), be the light stream initial value with the null matrix to refined net equation iteration m time, to eliminate high fdrequency component.If x _iBe initial value, obtain the approximate evaluation value

Be designated as:

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

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

A _iDetermined by formula (10); When i=1, r is definite 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{\‾}{x}}_i-f_i$

Step 4 is carried out the residual error transmission, eliminates low frequency component;

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

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

In the formula,

Be the mapping operator of refined net to the coarse grid.

Step 5, when i=1～N-1, repeating step three, four is passed to the thickest layer until residual error;

Step 6, find the solution the system of equations on the coarse grid:

A _Ne _N＝r _N

In the formula, e _NBe the light stream error of accurately finding the solution on the thickest layer.

Obtain by matrix operation:

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

Step 7 begins to carry out the error passback from the thickest layer (j=N);

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

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

In the formula,

Be the mapping operator of coarse grid to the refined net.

Step 8 is carried out refined net error correction;

The initial value that calculates is added the error of being passed back by coarse grid, obtain the solution after refined net is proofreaied and correct

That is:

${\hat{x}}_{j-1}={\hat{x}}_{j-1}+e_{j-1}$

Step 9 is carried out the back and is optimized iteration, improves stability of solution;

After refined net is proofreaied and correct, be initial value with the solution after proofreading and correct, carry out the back and optimize iteration n time, the corrected value after obtaining upgrading

Its expression formula is:

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

Step 10 is worked as j=N, N-1 ..., 2 o'clock, repeating step seven～nine was passed to the thinnest layer until residual error.

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

The present invention launches the light stream solving equation by row, with its form that turns to system of linear equations, make the complicated calculations process of light stream to accelerate to realize by linear multi-layer net algorithm.Secondly, adopt the multi-layer net of different resolution, the error component of different frequency has been distributed on the different grids, made error become high frequency error with respect to this layer grid.Because high frequency error convergence is very fast, thereby this method can significantly improve the computing velocity in visual light flow field, compares with the variational method, and the speed of convergence of this method can be brought up to more than 3.5 times.

Description of drawings

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

Fig. 2 is multi-layer image pyramid synoptic 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

The present invention will be further described below in conjunction with the drawings and specific embodiments.

One of the microcomputer of enforcement support P3800M dominant frequency of the present invention, 512M internal memory and above configuration is as the hardware device platform.Running environment is Windows XP operating system and Matlab software platform.

The detailed process of present embodiment is at first according to light stream model to be found the solution, to construct system of linear equations Ax=f.Set up 4～5 tomographic image pyramids then, be optimized iteration respectively at each layer, and the iteration residual error successively is mapped to down on the thicker grid of one deck.At the thickest layer, find the solution the system of equations on the coarse grid, obtain accurate estimation of error.After this, return error in opposite direction, on thinner grid, successively carry out error correction.Fig. 1 is the process flow diagram of method involved in the present invention, specifically may further comprise the steps:

Step 1, input picture makes up linear equation Ax=f;

Step 2 is set up the multi-layer image pyramid;

Step 3 is carried out preceding optimization, eliminates high fdrequency component;

Step 4 is carried out the residual error transmission, eliminates low frequency component;

Step 5, as i=1～N-1, repeating step three, four is passed to the thickest layer until residual error;

Step 6 is found the solution the system of equations on the coarse grid: A _Ne _N=r _N

Step 7 begins to carry out the error passback from the thickest layer (j=N);

Step 8 is carried out refined net error correction;

Step 9 is carried out the back and is optimized iteration, improves stability of solution;

Step 10 is worked as j=N, N-1 ..., 2 o'clock, repeating step seven～nine was passed to the thinnest layer until residual error.

Provide an application example of the present invention below.

Select two groups of test patterns that the quick optical flow field algorithm based on the distributed multi-layer net of error that the present invention proposes is verified.One group is a pair of test pattern RubberWhale in the Middlebury java standard library that extensively adopts in the world, as shown in Figure 3; Another group is scientific paper " the An improved algorithm for TV-L that people such as Wedel deliver at them ¹Optical flow " the middle a pair of test pattern that adopts, as shown in Figure 4.

During experiment, the quick optical flow field algorithm based on the distributed multi-layer net of error that adopts traditional variational algorithm and the present invention to propose respectively calculate every group of image between optical flow field.By the working time of two kinds of methods of software statistics, experimental result is as shown in table 1.In the table 1, image shown in Figure 3 is adopted in experiment one, and image shown in Figure 4 is adopted in experiment two.As shown in Table 1, shorten the operation time of comparing traditional variational algorithm based on the quick optical flow field algorithm of the distributed multi-layer net of error greatly, and the speed of convergence of two experiments is all brought up to more than 3.5 times, has significantly improved counting yield.

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

Claims (1)

1. the quick optical flow field computing method based on the distributed multi-layer net of error is characterized in that comprising the steps:

Step 1, input picture makes up linear equation Ax=f, and method is as follows:

(1) calculates space-time gradient tensor matrix, will wait the group discretize of solving an equation;

Suppose on a certain two field picture plane (i, the j) position of the capable j row of expression i pixel, u _{I, j}And v _{I, j}Represent this displacement field on the horizontal x axle that takes place between consecutive frame and vertical y axle respectively; The space partial derivative f of computed image _x, f _yAnd time partial derivative f _t, can obtain the space-time gradient tensor J matrix of each pixel in the image:

$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, j}Expression J matrix pixel (with h representation space step-length, the Euler-Lagrange equation that the variational method is obtained by the energy functional number is discrete to be for i, the value of j) locating:

${\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(2\right)$

${\left[J_{12}\right]}_{i,j}{u}_{i,j}+{\left[J_{22}\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_{23}\right]}_{i,j}---\left(3\right)$

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

(2) difference equation (2) is changed into linear equation with constant coefficient group A ₁X=f ₁

Abbreviation equation (2) obtains finding the solution the DIFFERENCE EQUATIONS in the regional D:

$(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]}_{i,j}$

Always total N+1 is capable to suppose image, and it is capable to remove the 0th row and N, and the DIFFERENCE EQUATIONS of first each pixel of 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{\α}}{{[J}_{21}]}_{(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{\α}}{{[J}_{13}]}_{11}+\frac{1}{4}{g}_{01}+\frac{1}{4}{g}_{10}\\ -\frac{h^2}{4\mathrm{\α}}{{[J}_{13}]}_{21}+\frac{1}{4}{g}_{20}\\ \·\·\·\\ \frac{h^2}{4\mathrm{\α}}{{[J}_{13}]}_{(N-1)1}+\frac{1}{4}{g}_{N1}+\frac{1}{4}{g}_{(N-1)0}\end{array}\right)=b_1---\left(4\right)$

In the formula, gi _jThe borderline element ui of regional D, i.e. g are found the solution in expression _Ij=u _Ij

Introduce light stream component u _IjAnd v _IjOn column vector:

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

And order:

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

For pixel (m, the J of the J matrix element of n) locating _Mn, matrix [J in the formula (5) _Mn] _iThe unit matrix and the J that are equivalent to the N-1 rank _MnI is listed as the product of this column vector; Therefore, the DIFFERENCE EQUATIONS (4) of first each pixel of row can turn to:

$(G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_1)u_1-\frac{1}{4}{\mathrm{Iu}}_2+\frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_1{v}_1=b_1$

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

$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{\α}}{{[J}_{11}]}_2)u_2-\frac{1}{4}{\mathrm{Iu}}_3+\frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_2{v}_2=b_2$

In the formula:

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

All row of image are all transformed, and then system of equations can be written as A ₁X=f ₁Form, in the formula:

$A_1=\left(\begin{array}{cccccccccccc}G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_1& \frac{1}{4}I& 0& \·\·\·& 0& 0& 0& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_1& 0& \·\·\·& 0& 0\\ -\frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_2& -\frac{1}{4}I& \·\·\·& 0& 0& 0& 0& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_2& \·\·\·& 0& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& \·\·\·& \frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_{N-2}& \frac{1}{4}I& 0& 0& \·\·\·& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_{N-2}& 0\\ 0& 0& 0& \·\·\·& 0& \frac{1}{4}I& G+\frac{h^2}{4\mathrm{\α}}{{[J}_{11}]}_{N-1}& 0& 0& \·\·\·& 0& \frac{h^2}{4\mathrm{\α}}{{[J}_{21}]}_{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 equation with constant coefficient group A ₂X=f ₂

(4) by A ₁X=f ₁With A ₂X=f ₂Common formation system of equations:

Ax=f (6)

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

The uppermost tomographic image of multi-layer image pyramid is former image in different resolution, and each following tomographic image is represented successively to dwindle falls the resolution image; Because pixel itself is uniform discrete, therefore can regard multi-layer image as multi-layer net;

Usually, the length of image and wide between 300～800; According to the scale factor η that sets, 0.5＜η＜0.95 successively dwindles former image in different resolution, and size is rounded, and obtains the 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 is carried out preceding optimization, eliminates high fdrequency component, and method is as follows:

Working as anterior layer i, i=1,2 ..., N-1 is the light stream initial value with the null matrix to refined net equation iteration m time, to eliminate high fdrequency component; If x _iBe initial value, obtain the approximate evaluation value

Be designated as:

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

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

A _iDetermined by formula (6); When i=1, r is definite 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{\‾}{x}}_i-f_i$

Step 4 is carried out the residual error transmission, eliminates low frequency component, and method is as follows:

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

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

In the formula, I _i ^I+1Be the mapping operator of refined net to the coarse grid;

Step 5, when i=1～N-1, repeating step three, four is passed to the thickest layer until residual error;

Step 6, find the solution the system of equations on the coarse grid:

A _Ne _N＝r _N

In the formula, e _NBe the light stream error of accurately finding the solution on the thickest layer;

Obtain by matrix operation:

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

Step 7 begins to carry out the error passback from the thickest layer, and method is as follows:

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

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

In the formula, Be the mapping operator of coarse grid to the refined net;

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

The initial value that calculates is added the error of being passed back by coarse grid, obtain the solution after refined net is proofreaied and correct

That is:

${\hat{x}}_{j-1}={\stackrel{\‾}{x}}_{j-1}+e_{j-1}$

Step 9 is carried out the back and is optimized iteration, improves stability of solution, and method is as follows:

After refined net is proofreaied and correct, be initial value with the solution after proofreading and correct, carry out the back and optimize iteration n time, the corrected value after obtaining upgrading

Its expression formula is:

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

Step 10 is worked as j=N, N-1 ..., 2 o'clock, repeating step seven～nine was passed to the thinnest layer until residual error.

Images