CN101951465A - Probability motion filtration-based real-time video image stabilizing method - Google Patents

Abstract

The invention discloses a probability motion filtration-based real-time video image stabilizing method. The method improves the filtration response speed, reduces the state estimation error in the case of mismatched model and improves the video image stabilization capability of image pick-up equipment in complicated environment, and comprises the following steps of: 1, dividing a current frame into blocks; 2, acquiring the motion vector of each block based on a diamond search method; 3, estimating a global motion parameter in the current frame relative to a previous frame according to the motion vector; 4, performing motion filtration process on the global motion parameter by using a particle filtering method; and 5, performing motion compensation according to the filtration result, wherein the importance sampling in the particle filtering method of the step 4 is achieved by an adaptive Kalman filtering method.

Description

Real-time video image stabilization method based on the probability motion filtering

Technical field

The present invention relates to is a kind of consumer electronics series products such as mobile phone, digital camera, video camera that are widely used in, also can be widely used in simultaneously shooting that aircraft, naval vessel, satellite, guided missile etc. are equipped and the real-time video image stabilization method in the surveillance, specifically be a kind of real-time video image stabilization method based on the probability motion filtering.

Background technology

The video sequence that video camera is taken in moving process, the active movement part that has not only comprised video camera, also introduced simultaneously irregular random motion, the existence of this random motion, can cause the shake of video and fuzzy, had a strong impact on people's visual experience, the purpose of the steady picture of video is exactly to eliminate or weaken this random motion to damaging influence that video pictures caused.The realization of the steady picture of video comprises estimation, motion filtering and three modules of inactive area compensation, and at these three modules, Chinese scholars has proposed a lot of methods.Aspect estimation, have based on piece, based on characteristic point, based on methods such as bit planes; Aspect motion filtering, weighted mean filtering, Kalman filtering (Kalman Filter) etc. are arranged; Aspect the inactive area compensation, have based on the image amplification, based on video reparation methods such as (Video Inpaint).Although proposed so multi-method, the good compromise but these methods still fail to accomplish on speed and effect, such as: aspect estimation, block-based method, speed is fast, but be subjected to illumination effect bigger,, have anti-preferably lighting effect based on the characteristic point method, but it is responsive to local motion, based on the bit plane method, to illumination-insensitive, but speed is slow and memory data output is bigger.

In addition aspect motion filtering, substantially all be to adopt legacy card Kalman Filtering or weighted mean filtering, the response of filtering is bad, and the filtering method that has utilizes future frame information to carry out the forward direction smoothing processing, although filter effect is better, owing to utilize future frame information not handle in real time; Inactive area compensation aspect, a kind of preferably at present method are to utilize video to repair thought to realize, still the fatal shortcoming of this method is the slow and dependence future frame information of speed, can't handle in real time.Therefore,, design one group of real-time and effective estimation, motion filtering and inactive area compensation method and seem extremely important surely as for the device for the video of a robust.

At last, find that by prior art documents Marius Tico is at " IEEEInternational Conference on Image Processing " (pp569-572,2005) deliver " Constraint motion filtering for videostabilization " on (based on the video image stabilization method of constrained motion filtering, the image processing ieee international conference), this article has adopted Kalman's motion filtering method of feature point tracking method and belt restraining restriction, result of the test shows that this method has and surely looks like effect preferably, but the estimation in the literary composition has adopted the characteristic matching based on both full-pixel, the precision height, but speed is slow, with having adopted in the eight-legged essay based on Kalman's motion filtering method of model at the uniform velocity, when motion model mates, filter effect is better, in case but model changes, then filter effect can be subjected to very big influence, and actual camera does not often meet the uniform motion rule in shooting process.

Summary of the invention

Technical problem to be solved by this invention provides a kind of real-time video image stabilization method based on the probability motion filtering, it has improved the response speed of filtering, reduced the state estimation error when model does not match simultaneously, and then improved picture pick-up device the video under the complex environment is surely looked like ability.

In order to solve above technical problem, the invention provides a kind of real-time video image stabilization method based on the probability motion filtering, it comprises the steps: the first step, present frame is carried out piece to be cut apart, second step, obtain the motion vector of each piece based on diamond search (ds), the 3rd step, go out present frame with respect to former frame global motion parameter according to described estimation of motion vectors, the 4th step, kinematic parameter to the overall situation adopts particle filter method to carry out the motion filtering processing, in the 5th step, carries out motion compensation according to the filtering result.

Importance sampling in described the 4th step particle filter method obtains by method for adaptive kalman filtering.

Described second step comprises based on diamond search (ds): step 1, utilize each piece content of integrogram feature description; Step 2 is that matching criterior is carried out the block matching motion estimation to each piece with integrogram characteristic error minimum.

Described the 3rd step comprises: the 3-1 step, whether the decision block motion is unusual, if not unusually then be labeled as active block, the 3-2 step, the mean value of adding up described active block translational movement estimates the global translation amount, in the 3-3 step, estimates rotation and the amount of zoom of each active block with respect to the match block in the former frame according to described global translation amount, in the 3-4 step, add up rotation and amount of zoom that described rotation and amount of zoom estimate the overall situation.

The probability particle filter method based on adaptive Kalman filter (Particle Filter) of the present invention's design, by regulating the observation noise variance, improve the response speed of Kalman filtering, come the state estimation error of correction card Kalman Filtering when motion model does not match by the probability particle filter method simultaneously.The present invention has simultaneously adopted the piece method for estimating based on the integrogram feature, and this piece method for estimating, at first current frame image is divided into some, as 8x8 or 16x16, every is mated in the particular search scope in reference frame, matching criterior between piece and the piece, conventional method adopts absolute error, mean square error etc., these criterions are easy to calculate, but it is responsive to picture noise and illumination variation, the present invention has adopted the integrogram characteristic error as matching criterior, this feature has fast operation, and have picture noise and the insensitive advantage of illumination variation, strengthened the anti-illumination of picture pick-up device and the ability of picture noise, improved the robustness of global motion parameter Estimation.In addition, after obtaining the piece kinematic parameter, need estimate the global motion parameter of present frame based on these kinematic parameters with respect to former frame, in order to estimate the global motion parameter exactly, need judge the movement tendency of each piece, and need reject the unusual piece of those movement tendencies, generally speaking, the movement tendency of several piece is consistent mostly, but the piece that has meeting is unusual owing to movement tendency appears in factors such as object local motion or interference, the present invention is based on most piece movement tendency agreement principles, think that the motion vector of most several piece obeys a Gaussian Profile, and near pairing of the motion vector that is positioned at the Gaussian Profile average is only active block, and all the other are the abnormal motion piece.Then, after judging active block, need utilize the kinematic parameter of these pieces to estimate the global motion parameter, classical way is to adopt least square method to find the solution an overdetermined equation group, promptly

Wherein Be the global motion parameter, X is the coordinate of piece in the former frame, X ' is the coordinate of corresponding blocks in the present frame, when finding the solution this equation, need ask matrix inversion, the dimension of this matrix is generally all higher, therefore computing is complicated, and might matrix to occur irreversible owing to the computer numerical problem, causes obtaining the global motion parameter, and utilization of the present invention estimates the translation that each piece of present frame is produced with respect to the former frame piece, rotation and amount of zoom, then, adopt the method for average statistical to obtain the global motion parameter, this method realizes simple, avoided complicated higher dimensional matrix inversion operation, and estimated accuracy is higher.

Description of drawings

Below in conjunction with the drawings and specific embodiments the present invention is described in further detail.

Fig. 1 (a)-1 (h) is a processing method The general frame of the present invention.

Wherein Fig. 1 (a) is integrogram feature description figure; Fig. 1 (b) is the integrogram feature description figure of moving mass; Fig. 1 (c) is a piece estimation flow chart; Fig. 1 (d) is abnormal motion piece and lost motion piece mark figure; Fig. 1 (e) is the anglec of rotation tolerance of single active block and the estimation schematic diagram of amount of zoom; Fig. 1 (f) is the global motion parameter Estimation flow chart based on active block; Fig. 1 (g) is motion filtering implementation procedure figure; Fig. 1 (h) is the The general frame of the steady picture of video.

Fig. 2 (a)-2 (b) is two kinds of motion filtering effect comparison diagrams among the present invention.

Wherein Fig. 2 (a) is Kalman's motion filtering method design sketch of belt restraining restriction, and Fig. 2 (b) is the probability particle filter method design sketch based on adaptive Kalman filter.

Fig. 3 (a)-3 (f) surely looks like design sketch for video among the present invention.

Wherein Fig. 3 (a) is an original video sequence, Fig. 3 (b) is a steady picture back sequence, Fig. 3 (c) is the motion filtering design sketch of x direction displacement, Fig. 3 (d) is the motion filtering design sketch of y direction displacement, Fig. 3 (e) is the motion filtering design sketch of anglec of rotation tolerance, and Fig. 3 (f) is the motion filtering design sketch of amount of zoom.

Embodiment

The The general frame of a kind of real-time video image stabilization method based on the probability motion filtering that Fig. 1 (a)-1 (h) proposes for the present invention.The concrete implementation detail of each several part is as follows:

Fig. 1 (a) is integrogram feature description figure

1. integrogram feature (corresponding English be Integral Image Feature)

Shown in Fig. 1 (a), so-called integrogram feature, in the drawings in the white rectangle zone in grey scale pixel value sum and the black rectangle zone difference of grey scale pixel value sum be the integrogram feature, adopted 5 kinds of integrogram features among the figure.The benefit of this feature is embodied in: the influence that the pixel value sum in (1) zone can effectively resist picture noise, (2) in the adjacent rectangle zone pixel and the effective customer service illumination effect of difference, so this feature all has better anti-jamming capability to picture noise and illumination.In order to calculate these features fast, need to convert original image to integral image, conversion formula is:

$\mathrm{ii}(x,y)=\underset{x^{\′}\≤x,y^{\′}\≤y}{\mathrm{\Σ}}i(x^{\′},y^{\′})---\left(1\right)$

Wherein, x and y are the pixel coordinate in the integral image, and x ' and y ' are the pixel coordinate in the original image, ii is an integral image, i is an original image, and integral image has been arranged, and the computational short cut of integrogram feature is the simple plus-minus to rectangle summit locational integrogram pixel value so.

Fig. 1 (b) is the integrogram feature description figure of moving mass

1. piece integrogram feature and piece matching criterior

The content of piece can be described by one group of size single integrogram feature different with the position, piece for a 16x16, integrogram feature in it has thousands of, if come the content of description block with whole integrogram features, the precision height, but real-time is poor, just can describe out the content of a piece well by about 100 the integrogram features of test discovery.Behind good each piece of integrogram feature description, then the matching error criterion of two pieces can design as follows:

$\mathrm{ISAD}({\mathrm{\Δ}}_x,{\mathrm{\Δ}}_y)=\underset{k=0}{\overset{M}{\mathrm{\Σ}}}|F_{B(i,j)}\left(k\right)-F_{B(i+{\mathrm{\Δ}}_x,j+{\mathrm{\Δ}}_y)}\left(k\right)|---\left(2\right)$

Wherein, k is the call number of integrogram feature, and M is an integrogram feature sum, Δ _xAnd Δ _yBe displacement, (i is (i, piece j), B (i+ Δ for the present frame center j) to B _x, the j+ Δ _y) for the former frame center be (i+ Δ _x, the j+ Δ _y) piece, F _{B (i, j)}(k) and

(k) be respectively piece B (i, j) and B (i+ Δ _x, the j+ Δ _y) last k integrogram feature, if this error is minimum, then show this two piece couplings in search procedure.

Fig. 1 (c) is a piece estimation flow chart

1. after the piece estimation had been determined the matching error criterion, searching algorithm searched out corresponding blocks with present frame block-matching error minimum according to this criterion in the former frame image.In order to ensure the accuracy of coupling, in search procedure, select the piece of feature rich to go to search for as far as possible, and, directly carry out the invalid block mark the not abundant piece of those features, do not participate in motion match.The feature rich degree can be calculated by following formula:

$T=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}{F}_B\left(k\right)---\left(3\right)$

Wherein, F _B(k) be k integrogram feature on the piece B, M is an integrogram feature sum.Generally speaking, the texture information in the piece is relevant with T value size, and the T value is big more, and this piece texture information is obvious more, and feature is abundant more, otherwise, otherwise.If T is less than a minimum threshold, then the feature of this piece cannot not be defined as abundantly, and Direct Mark is an invalid block.Equally, the piece bigger to those matching errors in estimation also carries out the invalid block mark, because error is bigger than normal, illustrates that coupling may be inaccurate.The part piece is carried out invalid mark, can improve the precision of global motion parameter Estimation.On searching method, the present invention has adopted traditional diamond search (ds), and this method speed is fast, search precision is high.

Fig. 1 (d) is abnormal motion piece and lost motion piece mark figure

1. the abnormal motion piece is rejected

The rejecting abnormalities moving mass is the same with the purpose of mark invalid block, all be in order to estimate the global motion parameter more accurately, generally speaking, the movement tendency of several piece all is consistent mostly, but the piece that has can be because the local motion or the interference of object, it is unusual movement tendency to occur, must find out these pieces and also be rejected.In the present invention, adopted a kind of statistical method to realize rejecting to the abnormal motion piece, at first, add up the average and the variance of the motion vector of all pieces, then with this average and gauss of distribution function of variance structure, if the motion vector distribution of certain piece then is labeled as the abnormal motion piece with this piece beyond 3 times of variances of this Gaussian Profile average.The piece that indicates the canescence line among Fig. 1 (d) is the dyskinesia piece, and the piece of no lines mark is the bigger piece of error in the abundant or matching process of feature, and the piece that indicates the black line is an active block.The judgement performing step of abnormal motion piece is as follows:

Step 1: the average of statistics block kinematic parameter and variance

The average of piece kinematic parameter and variance computing formula are as follows:

$U_x=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{V}_x\left(l\right)---\left(4\right)$

$S_x^2=\frac{1}{L-1}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{[V_x\left(l\right)-U_x]}^2---\left(5\right)$

$U_y=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{V}_y\left(l\right)---\left(6\right)$

$S_y^2=\frac{1}{L-1}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{[V_y\left(l\right)-U_y]}^2---\left(7\right)$

Wherein, l is the index value of piece, and L is the total quantity of piece, V _x(l) be the moving displacement of l piece on the x direction, V _y(l) be moving displacement on the y direction, U _xBe the moving displacement average of all pieces on the x direction,

Be the moving displacement variance on the x direction, U _yBe the moving displacement average of all pieces on the y direction,

Be the moving displacement variance on the y direction.

Step 2: the dyskinesia of checking each piece

For k piece, if its x and y direction displacement V _x(k) and V _y(k) satisfy following any one formula, then this piece is the abnormal motion piece, needs to reject, and does not participate in average and variance calculating next time, and formula is as follows:

$\left\{\begin{array}{c}{[V_x\left(k\right)-U_x]}^2>T*S_x^2\\ {[V_y\left(k\right)-U_y]}^2>T*S_y^2\end{array}\right.---\left(8\right)$

Wherein, T is a threshold value, generally gets 2.5-4.

Step 3: circulation step 1 and step 2, when iteration satisfies a maximum iteration time, or the current average that counts compares variation when very little with previous average, stops iteration.

Fig. 1 (e) is the anglec of rotation tolerance of single active block and the estimation schematic diagram of amount of zoom

1. the anglec of rotation of single active block tolerance and amount of zoom are estimated

After judging active block, need to estimate the global motion parameter based on the kinematic parameter of these pieces.Present frame, then can be described by the piece motion process among Fig. 1 (e) if be reacted on the piece with respect to the global motion process of former frame, that is: the piece of former frame is from O ₁Point moves to the O of present frame ₂Point then, rotates to O ₃Point, last, zoom to O ₄The point.Before the anglec of rotation tolerance of estimating single active block and amount of zoom, need count the global translation amount earlier, the global translation amount can be calculated by following formula:

$G_x=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{V}_x\left(m\right)---\left(9\right)$

$G_y=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{V}_y\left(m\right)---\left(10\right)$

Wherein, G _xAnd G _yBe the global translation amount of present frame with respect to former frame, M is the sum of active block, and m is the index of active block, V _x(m) and V _y(m) be respectively the moving displacement amount of m active block on x direction and y direction.

Piece is at O ₁G has been passed through in the position _xAnd G _yAfter the translation, arrive O ₂Point rotates to O then ₃Point, anglec of rotation θ (m) can calculate by following formula:

$\mathrm{\&theta;}\left(m\right)=\&angle;{\mathrm{<figure-callout\; id="2"\; label="O"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">O</figure-callout>}}_2{\mathrm{OO}}_4$

$=\&angle;{\mathrm{BOO}}_4-\&angle;{\mathrm{<figure-callout\; id="2"\; label="BOO"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">BOO</figure-callout>}}_2$

$=a\tan \left(\frac{O_{4x}}{O_{4y}}\right)-a\tan \left(\frac{O_{2x}}{O_{2y}}\right)---\left(11\right)$

$=a\tan \left(\frac{O_{4x}}{O_{4y}}\right)-a\tan \left(\frac{O_{1x}+G_x}{O_{1y}+G_y}\right)$

Wherein, O _4yAnd O _4xFor the present frame center is (O _4y, O _4x) piece, O _1yAnd O _1xBe the best matching blocks center that in former frame, searches, G _xAnd G _yBe the global translation amount of estimating.

In like manner, can derive the amount of zoom of each active block:

$S\left(m\right)=\frac{r_{\mathrm{oo}4}}{r_{\mathrm{oo}3}}=\frac{r_{\mathrm{oo}4}}{{\mathrm{<figure-callout\; id="2"\; label="r\; oo"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">r</figure-callout>}}_{{\mathrm{oo}}_{}}}$

$=\frac{\sqrt{O_{4\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">y</figure-callout>}}^2+O_{4x}^2}}{\sqrt{O_{2\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">y</figure-callout>}}^2+O_{2x}^2}}---\left(12\right)$

$=\frac{\sqrt{O_{4\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">y</figure-callout>}}^2+O_{4x}^2}}{\sqrt{{(O_{1y}+G_y)}^2+{(O_{1x}+G_x)}^2}}$

Fig. 1 (f) estimates flow chart based on the global parameter of active block

1. global motion parameter Estimation

After the translational movement that obtains each active block, anglec of rotation tolerance and amount of zoom, then the global motion parameter also can be obtained by this tittle statistics, wherein, the global translation amount is estimated by formula (9) and formula (10), and overall anglec of rotation tolerance and amount of zoom can be estimated by following formula:

$G_{\mathrm{\&theta;}}=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{\&theta;}\left(m\right)---\left(13\right)$

$G_S=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}S\left(m\right)---\left(14\right)$

Wherein, G _θBe overall rotation amount, G _SBe overall amount of zoom, θ (m) is the anglec of rotation tolerance of m active block, and S (m) is the amount of zoom of m active block.

Fig. 1 (g) is motion filtering implementation procedure figure.

So-called motion filtering is exactly the filtering jittering component, makes processed video become steady.Because the method for amplifying based on image has been adopted in the compensation of the inactive area among the present invention, therefore need carry out amplitude limit to the output of filtering algorithm, promptly compensation rate can not be greater than certain threshold value.The present invention has used for reference the Kalman filtering thought of belt restraining in the Marius Tico paper, but the method in this paper exists filter response bad, and motion model is not when matching, state estimation error problem bigger than normal, at these two problems, the present invention has designed a kind of probability particle filter method based on adaptive Kalman filter, regulate the fast-response that the observation noise variance improves Kalman filtering, state estimation error problem bigger than normal when coming the customer service motion model not match by the probability particle filter algorithm simultaneously by self adaptation.

Particle filter is a kind of probability filtering algorithm based on Monte Carlo method, and it can effectively solve the state estimation problem under non-linear, non-Gauss's situation.Particle filter algorithm distributes by the posteriority that the heavy particle of one group of cum rights approaches state, and when the number of particle was tending towards infinity, then these particles can approach the true posteriority distribution of state in theory.Particle filter algorithm generally comprises importance sampling, resampling and three steps of state estimation.In the importance sampling step, at first, t-1 particle constantly

Pass through motion model

Be delivered to next constantly; Then, at given measured value Z _t, the particle that each is transmitted

Be endowed weights

At last, these weights of normalization

Resampling the step particle Will be by resampling, the selecteed number of times of each particle is directly proportional with its weights size, and the selecteed possibility of the particle that weights are big more is big more.After the posteriority of the state of estimating distributed, the state of target can utilize lowest mean square root error estimation (MMSE) to estimate.Probability particle filter method performing step based on adaptive Kalman filter is as follows:

Step 1: importance sampling

In the importance sampling step, by the sample states of method for adaptive kalman filtering acquisition probability particle filter, the specific implementation process is as follows:

$X_{t|t-1}^{\left(i\right)}={\mathrm{AX}}_{t-1}^{\left(i\right)},i=1,...,N---\left(15\right)$

$P_{t|t-1}^{\left(i\right)}={\mathrm{AP}}_{t-1}^{\left(i\right)}{A}^T+{\mathrm{\&sigma;}}_e^2{\mathrm{BB}}^T---\left(16\right)$

$G=P_{t|t-1}^{\left(i\right)}C{(C^T{P}_{t|t-1}^{\left(i\right)}C+{\mathrm{\&sigma;}}_{t,u}^2)}^{-1}---\left(17\right)$

$P_t^{\left(i\right)}=(I-{\mathrm{GC}}^T)P_{t|t-1}^{\left(i\right)}---\left(18\right)$

$X_t^{\left(i\right)}=X_{t|t-1}^{\left(i\right)}+G[Z_t-C^T{X}_{t|t-1}^{\left(i\right)}]---\left(19\right)$

${\tilde{X}}_t^{\left(i\right)}=X_t^{\left(i\right)}+N(0,P_t^{\left(i\right)})---\left(20\right)$

Wherein, i is the index value of particle, and N is a total number of particles, B=[1 1] ^T, C=[1 0] ^T,

Be i particle t-1 state constantly,

Be i particle in t predicted state constantly, Be i particle t-1 state covariance matrix constantly,

Be that i particle predicted battle array in t state covariance constantly, G is the gain battle array of Kalman filtering, Be that i particle estimated battle array in t state covariance constantly,

Be the estimated state of i particle in t moment Kalman filtering,

Be the state after i particle resamples, Z _tBe measured value,

Be the state-noise variance,

Be t observation noise variance constantly, this value is to utilize prior information to set a fixed value in traditional kalman filter method, but be difficult to obtain this prior information in the real system, and the size that should be worth directly influences the response speed of Kalman filtering, among the present invention, designed the method that a kind of self adaptation is regulated the observation noise variance.

When shake is very serious, be reflected on the measured value bigger for its variance, at this moment, if also carry out Kalman filtering with fixing observation noise variance, the fast-response of filtering is affected, in order to improve the fast-response of filtering, the present invention has designed the strategy that following self adaptation is regulated the observation noise variance.At first, the measured value in statistics a period of time window changes severe degree, and the description formula of severe degree is as follows:

$M_z=\frac{1}{L}\underset{l=t-L-1}{\overset{t}{\mathrm{\&Sigma;}}}{Z}_l---\left(21\right)$

$S_z=\frac{1}{L}\underset{l=t-L-1}{\overset{L}{\mathrm{\&Sigma;}}}|Z_l-M_z|---\left(22\right)$

Wherein, L is a time window length, Z _lBe l measured value constantly, M _zBe the measured value average in the time window, S _zFor the measured value in the time window changes severe degree, work as S _zBig more, then in this time window, it is strong that measured value changes Shaoxing opera, otherwise, otherwise.Secondly, change severe degree according to measured value, self adaptation is regulated the observation noise variance, regulates formula and is:

$\left\{\begin{array}{cc}{\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2=\mathrm{\&alpha;}{\mathrm{\&sigma;}}_{t-1,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2,& \mathrm{if}{S}_z>S_{\max }\\ {\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2={\mathrm{\&sigma;}}_{t-1,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2,& \mathrm{if}{S}_{\min }\&le;S_z\&le;S_{\max }\\ {\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2=\mathrm{\&beta;}{\mathrm{\&sigma;}}_{t-1,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2,& \mathrm{if}{S}_z<S_{\min }\end{array}\right.---\left(23\right)$

Wherein,

Be t observation noise variance constantly,

Be t-1 observation noise variance constantly; S _zFor the measured value in the time window changes severe degree, α and β are adjustment factor, and α generally gets 0.7-0.9, and β generally gets 1.1-1.3; S _MaxBe the max-thresholds of measured value variation severe degree, S _MinFor measured value changes the severe degree minimum threshold.That is to say and work as S _zGreater than S _MaxThe time, reduce t observation noise variance constantly; Work as S _zLess than S _MinThe time, increase t observation noise variance constantly.The purpose of formula (23) is exactly to change when violent in measured value, increases the gain of Kalman filtering by turning the observation noise variance down, thereby improves the response speed of state estimation.

Excessive or too small for the observation noise variance that prevents to adjust, adopt following formula to retrain, make the self adaptation adjustment between a minimum and maximum value of observation noise variance, it is as follows to regulate the formula design:

$\left\{\begin{array}{cc}{\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2={\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="max"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">max</figure-callout>}}^2,& \mathrm{if}{\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2>{\mathrm{\&sigma;}}_{\max }^2\\ {\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2={\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="min"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">min</figure-callout>}}^2,& \mathrm{if}{\mathrm{\&sigma;}}_{t,\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="HSA00000133084200052.png,HSA00000133084200053.png"\; state="\{\{state\}\}">u</figure-callout>}}^2<{\mathrm{\&sigma;}}_{\min }^2\end{array}\right.---\left(24\right)$

That is to say Adjustment process in worthwhile during greater than maximum set threshold value, get maximum set threshold value, when less than minimum setting threshold, get minimum setting threshold.Maximum set threshold value generally gets 40 ²-50 ²Between value, minimum setting threshold generally gets 5 ²-10 ²Between value.

Step 2: calculate the particle weights

$w_t^{\left(i\right)}=\frac{1}{\sqrt{2{\mathrm{\&pi;\&sigma;}}_{t,u}}}\exp [-\frac{{({\tilde{X}}_t^{\left(i\right)}-Z_t)}^2}{2{\mathrm{\&sigma;}}_{t,u}^2}]---\left(25\right)$

Step 3: normalization particle weights

${\tilde{w}}_t^{\left(i\right)}=w_t^{\left(i\right)}/\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_t^{\left(i\right)}---\left(26\right)$

Step 4: resampling particle

${\{X_t^{\left(i\right)},1/N\}}_{i=1}^N=\mathrm{Resample}\left[{\{{\tilde{X}}_t^{\left(i\right)},{\tilde{w}}_t^{\left(i\right)}\}}_{i=1}^N\right]---\left(27\right)$

Step 5: state estimation

When state estimation, if the difference of selecteed particle state and actual observed value is too big, when promptly motion compensation quantity is too big, need retrain processing to compensation rate, the constraint process is achieved as follows:

$X_t^{\left(i\right)}=\left\{\begin{array}{c}{X}_t^{\left(i\right)}+\mathrm{sign}(Z_n-C^T{X}_t^{\left(i\right)})\left[\right|Z_t-C^T{X}_t^{\left(i\right)}|-D]P_t^{\left(i\right)}C{\left[C^T{P}_t^{\left(i\right)}C\right]}^{-1},\mathrm{if}|Z_t-C^T{X}_t^{\left(i\right)}|>D\\ X_t^{\left(i\right)},\mathrm{otherwise}\end{array}\right.---\left(28\right)$

Wherein, Z _tBe t measured value constantly, C=[1 0],

Be that i particle estimated battle array in t state covariance constantly, D is the max-thresholds of compensation rate.

Step 6: state output

After each particle state carried out amplitude limiting processing, final output state was obtained by following formula:

${\hat{X}}_t=\frac{1}{N}{X}_t^{\left(i\right)}---\left(29\right)$

Wherein,

Be t estimated state constantly,

Be i particle state after resampling and the amplitude limiting processing.

Fig. 2 (a)-2 (b) switches between rapid movement and static two states back and forth for the effect comparison diagram of two kinds of motion filtering methods among the present invention, the dbjective state among this figure.Fig. 2 (a) is Kalman's motion filtering method design sketch of belt restraining restriction, Fig. 2 (b) is the probability particle filter method design sketch based on adaptive Kalman filter, dotted line 2 is filtered design sketch among two figure, solid line 1 is the design sketch before the filtering, the longitudinal axis is represented the displacement (Horizontal Position) of horizontal x direction, and transverse axis is represented frame number (Frame).From figure, can obviously find out, no matter on fast-response still is the state estimation of model when not matching, all be better than Kalman's motion filtering method of belt restraining restriction based on the probability particle filter method of adaptive Kalman filter.

Fig. 3 (a)-3 (f) surely looks like design sketch for video among the present invention.Wherein, Fig. 3 (a) is an original video sequence, Fig. 3 (b) is a steady picture back sequence, Fig. 3 (c) is the motion filtering design sketch of x direction displacement, Fig. 3 (d) is the motion filtering design sketch of y direction displacement, Fig. 3 (e) is the motion filtering design sketch of anglec of rotation tolerance, Fig. 3 (f) is the motion filtering design sketch of amount of zoom, Fig. 3 (c) wherein, (d), (e) and the dotted line (f) 2 be filtered design sketch, solid line 1 is the design sketch before the filtering, the longitudinal axis of Fig. 3 (c) is represented horizontal x direction displacement (Horizontal Position), transverse axis is represented frame number (Frame), the longitudinal axis of Fig. 3 (d) is represented the displacement (Vertical Position) of vertical y direction, transverse axis is represented frame number (Frame), and Fig. 3 (e) longitudinal axis is represented anglec of rotation tolerance (Rotation Angle), and transverse axis is represented frame number (Frame), Fig. 3 (f) longitudinal axis is represented amount of zoom (Scale), and transverse axis is represented frame number (Frame).From these as a result the figure as can be seen, the state estimation of the filtering method among the present invention when filter response and model do not match all has better effects.

Claims (15)

1. real-time video image stabilization method based on the probability motion filtering, it is characterized in that the first step is carried out piece to present frame and cut apart, second step, obtain the motion vector of each piece based on diamond search (ds), in the 3rd step, go out the global motion parameter of present frame with respect to former frame according to described estimation of motion vectors, the 4th step, kinematic parameter to the overall situation adopts particle filter method to carry out the motion filtering processing, in the 5th step, carries out motion compensation according to the filtering result.

2. the real-time video image stabilization method based on the probability motion filtering according to claim 1 is characterized in that, the importance sampling in described the 4th step particle filter method obtains by method for adaptive kalman filtering.

3. the real-time video image stabilization method based on the probability motion filtering according to claim 2, it is characterized in that, described method for adaptive kalman filtering is to come self adaptation to regulate observation noise variance in the Kalman filtering by the variation severe degree of judging the global motion parameter, wherein regulates the observation noise variance and adopts following formula:

$\left\{\begin{array}{cc}{\mathrm{\&sigma;}}_{t,u}^2={\mathrm{\&alpha;\&sigma;}}_{t-1,u}^2,& \mathrm{if}{S}_z>S_{\max }\\ {\mathrm{\&sigma;}}_{t,u}^2={\mathrm{\&sigma;}}_{t-1,u}^2,& \mathrm{if}{S}_{\min }\&le;S_z\&le;S_{\max }\\ {\mathrm{\&sigma;}}_{t,u}^2={\mathrm{\&beta;\&sigma;}}_{t-1,u}^2,& \mathrm{if}{S}_z<S_{\min }\end{array}\right.---\left(1\right)$

Wherein, Be t observation noise variance constantly,

Be t-1 observation noise variance constantly; S _zFor the measured value in the time window changes severe degree, work as S _zBig more, then in this time window, it is strong that measured value changes Shaoxing opera, otherwise, otherwise; α and β are adjustment factor; S _MaxBe the max-thresholds of measured value variation severe degree, S _MinFor measured value changes the severe degree minimum threshold.

4. the real-time video image stabilization method based on the probability motion filtering according to claim 3 is characterized in that, and is described

Adjustment process in worthwhile during greater than maximum set threshold value, get maximum set threshold value, when less than minimum setting threshold, get minimum setting threshold.

5. according to claim 3 or 4 described real-time video image stabilization methods, it is characterized in that described maximum set threshold value generally gets 40 based on the probability motion filtering ²-50 ²Between value, minimum setting threshold generally gets 5 ²-10 ²Between value.

6. according to claim 3 or 4 described real-time video image stabilization methods based on the probability motion filtering, it is characterized in that described α gets 0.7-0.9, β gets 1.1-1.3.

7. according to the described real-time video image stabilization method of claim 1, it is characterized in that described second step comprises based on diamond search (ds) based on the probability motion filtering: step 1, utilize each piece content of integrogram feature description; Step 2 is that matching criterior is carried out the block matching motion estimation to each piece with integrogram characteristic error minimum.

8. according to the described real-time video image stabilization method of claim 7, it is characterized in that, before described step 2, also comprise described feature rich degree is judged, calculate by following formula based on the probability motion filtering:

$T=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{F}_B\left(k\right)---\left(2\right)$

Wherein, F _B(k) be k integrogram feature on the piece B, M is an integrogram feature sum, the T representation feature enriches degree, the T value is big more, feature is abundant more, if T less than a minimum threshold, then the feature of this piece cannot not be defined as abundantly, Direct Mark is an invalid block, this piece is not carried out the block matching motion of described step 2 and estimates.

9. according to claim 7 or 8 described real-time video image stabilization methods based on the probability motion filtering, it is characterized in that described step 2 also comprises the block-matching error determining step, this step is that the matching error of former and later two pieces is calculated according to following formula:

$\mathrm{ISAD}({\mathrm{\&Delta;}}_x,{\mathrm{\&Delta;}}_y)=\underset{k=0}{\overset{M}{\mathrm{\&Sigma;}}}|F_{B(i,j)}\left(k\right)-F_{B(i+{\mathrm{\&Delta;}}_x,j+{\mathrm{\&Delta;}}_y)}\left(k\right)|---\left(3\right)$

K is the call number of integrogram feature, and M is an integrogram feature sum, Δ _xAnd Δ _yBe displacement, (i is (i, piece j), B (i+ Δ for the present frame center j) to B _x, the j+ Δ _y) for the former frame center be (i+ Δ _x, the j+ Δ _y) piece, F _{B (i, j)}(k) and

Be respectively piece B (i, j) and B (i+ Δ _x, the j+ Δ _y) last k integrogram feature, ISAD (Δ _x, Δ _y) be the matching error value, if this matching error value then is labeled as invalid block greater than a maximum set value, this piece does not participate in the global motion parameter Estimation in described the 3rd step.

10. according to the described real-time video image stabilization method of claim 1 based on the probability motion filtering, it is characterized in that, described the 3rd step comprises: the 3-1 step, whether the decision block motion is unusual, if not unusually then be labeled as active block, the 3-2 step, the mean value of adding up described active block translational movement estimates the global translation amount, the 3-3 step, estimate rotation and the amount of zoom of each active block according to described global translation amount with respect to the match block in the former frame, in the 3-4 step, add up rotation and amount of zoom that described rotation and amount of zoom estimate the overall situation.

11. the real-time video image stabilization method based on the probability motion filtering according to claim 10 is characterized in that, whether decision block comprised the steps: unusually during described 3-1 went on foot

Step 1: statistics block is at x, the average of the kinematic parameter on the y direction and variance, and concrete computing formula is as follows:

$U_x=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{V}_x\left(l\right)---\left(4\right)$

$S_x^2=\frac{1}{L-1}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{[V_x\left(l\right)-U_x]}^2---\left(5\right)$

$U_y=\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{V}_y\left(l\right)---\left(6\right)$

$S_y^2=\frac{1}{L-1}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}{[V_y\left(l\right)-U_y]}^2---\left(7\right)$

Wherein, l is the index value of piece, and L is the total quantity of piece, V _x(l) be the moving displacement of l piece on the x direction, V _y(l) be moving displacement on the y direction, U _xBe the moving displacement average of all pieces on the x direction,

Be the moving displacement variance on the x direction, U _yBe the moving displacement average of all pieces on the y direction,

Be the moving displacement variance on the y direction;

Step 2: check the dyskinesia of each piece, for k piece, if its x and y direction displacement V _x(k) and V _y(k) satisfy following any one formula, then this piece is the abnormal motion piece, otherwise is labeled as active block, and formula is as follows:

$\left\{\begin{array}{c}{[V_x\left(k\right)-U_x]}^2>T*S_x^2\\ {[V_y\left(k\right)-U_y]}^2>T*S_y^2\end{array}\right.---\left(8\right)$

Wherein, T is a threshold value, generally gets 2.5-4.As this piece is the abnormal motion piece, does not then participate in the statistics of next average and variance;

Step 3: circulation step 1 and step 2, when iteration satisfies a maximum iteration time, or the current average that counts compares variation when very little with previous average, stops iteration.

12. the real-time video image stabilization method based on the probability motion filtering according to claim 10, it is characterized in that, the mean value of the described active block translational movement of described 3-2 step statistics estimates the global translation amount, and described global translation amount is included in the global translation amount G on the x direction _xWith the global translation amount G on the y direction _y, concrete computing formula is as follows:

$G_x=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{V}_x\left(m\right)---\left(9\right)$

$G_y=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{V}_y\left(m\right)---\left(10\right)$

Wherein, M is the sum of active block, and m is the index of active block, V _x(m) and V _y(m) be respectively the moving displacement of m active block on x direction and y direction.

13., it is characterized in that the anglec of rotation θ (m) of m the active block in described 3-3 step calculates by following formula according to claim 10 or 12 described real-time video image stabilization methods based on the probability motion filtering:

$\mathrm{\&theta;}\left(m\right)=a\tan \left(\frac{O_{4x}}{O_{4y}}\right)-a\tan \left(\frac{O_{1x}+G_x}{O_{1y}+G_y}\right)---\left(11\right)$

Wherein, O _4yAnd O _4xFor the present frame center is (O _4y, O _4x) piece, O _1yAnd O _1xBe the best matching blocks center that in former frame, searches, G _xAnd G _yFor estimate described on x direction and y direction the global translation amount, m is from 1 to M, M is the sum of active block.

14., it is characterized in that the amount of zoom S (m) of m the active block in described 3-3 step calculates by following formula according to claim 10 or 12 described real-time video image stabilization methods based on the probability motion filtering:

$S\left(m\right)=\frac{\sqrt{O_{4y}^2+O_{4x}^2}}{\sqrt{{(O_{1y}+G_y)}^2+{(O_{1x}+G_x)}^2}}---\left(12\right)$

Wherein, O _4yAnd O _4xFor the present frame center is (O _4y, O _4x) piece, O _1yAnd O _1xBe the best matching blocks center that in former frame, searches, G _xAnd G _yBe the global translation amount of estimating.

15. the real-time video image stabilization method based on the probability motion filtering according to claim 10 is characterized in that, described overall situation rotation of described 3-4 step statistics and amount of zoom can calculate by following formula:

$G_{\mathrm{\&theta;}}=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{\&theta;}\left(m\right)---\left(13\right)$

$G_S=\frac{1}{M}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}S\left(m\right)---\left(14\right)$

G _θBe overall rotation amount, G _SBe overall amount of zoom, θ (m) is the anglec of rotation of m active block, and S (m) is the amount of zoom of m active block.

Images