CN105160674A - Improved quick bidimensional empirical mode decomposition method - Google Patents

Abstract

An improved quick bidimensional empirical mode decomposition method comprises: determining to-be-screened data; detecting a local limit point of data ri-1(x,y); obtaining a mean envelope mmean(x,y) from hi(x,y) by using a mean filter; after obtaining the mean envelope, calculating a new signal hi'(x,y); calculating whether the generated new signal meets a screening stop criterion SD or not, and if SD is smaller than a given threshold value xi, determining that hi'(x,y) is a BIMF (Bidimensional Intrinsic Mode Function) of an ith layer; and setting hi(x,y) to be equal to hi'(x,y). Compared with an existing bidimensional empirical mode decomposition method, the improved quick bidimensional empirical mode decomposition method has the characteristics that the mean envelope is obtained by applying one fast convolution based mean filter, so that compared with a conventional bidimensional empirical mode decomposition method for obtaining a mean envelope based on an RBF (Radial Basis Function), the time complexity is greatly lowered, and the calculation rate of an algorithm is increased under the circumstance that a similar experimental result is obtained; and the window size of the mean filter is calculated by applying a Manhattan distance, so that the whole algorithm is adaptively completed without the need for manual intervention.

Description

A kind of fast two-dimensional empirical mode decomposition method of improvement

Technical field

The present invention is a kind of fast two-dimensional empirical mode decomposition method of improvement, belongs to the novel Time-Frequency Analysis Method in digital image processing field.

Background technology

Along with the development of computer technology and the continuous expansion in each application thereof, digital image processing techniques have also been obtained swift and violent development.Empirical mode decomposition algorithm (EMD) is a kind of novel Time-Frequency Analysis Method being applicable to analyze non-linear/non-stationary signal, the method is considered to the important breakthrough linearly and to be steadily assumed to be traditional Time-Frequency Analysis Method such as basic Fourier analysis and wavelet transformation, it is the time scale feature based on self, completely by data-driven, can adaptive be limited solid-state mode function (IMF) and a residual error trend term sum by " sieve " process by the signal decomposition of complexity.Image can regard a 2D signal as, and two-dimensional empirical mode decomposition algorithm (BEMD) is the popularization of EMD in 2D signal process, becomes the study hotspot in image processing field in the last few years.

Traditional two-dimensional empirical mode decomposition algorithm, in the process of screening, first needs to ask for Local Extremum, then asks for up/down envelope by local extremum, then ask for average envelope, finally judges whether average envelope meets screening stopping criterion.Common two-dimensional empirical mode decomposition algorithm, all by one group of radial basis function (RadialBasisFunction after searching local extremum, be called for short RBF) obtain envelope, but this method needs to use method of interpolation to be similar to obtains envelope, therefore, the computing time of algorithm is significantly increased.In addition, the time using the method to obtain envelope approximately account for 75% in whole screening process, can not meet the real-time of computation requirement.Therefore, be necessary to improve this algorithm, two-dimensional solid-state modular function (BIMF) can be generated more efficiently, and there is better availability, the real-time in practical application can be met.

Summary of the invention

The object of the invention is, when in two-dimensional empirical mode decomposition algorithm application and image procossing, improves algorithm counting yield and makes it meet real-time in practical application; Overcome in existing method shortcomings such as there is envelope formation speed is slow, border pollution, propose a kind of fast two-dimensional empirical mode decomposition method of improvement.

To achieve these goals, the present invention takes following technical scheme:

One, the data will carrying out sieving are determined.The data of carrying out sieving are the trial image of standard; As Lena image, cameraman image, Elaine image etc.If pending original image is f (m, n) (wherein m is the row-coordinate of image, n image column coordinate), make i=1, the data so will carrying out sieving are r ₀(x, y)=f (m, n) (wherein x is the row-coordinate of image, y image column coordinate).

Two, data r is detected _i-1the Local Extremum of (x, y).Detect data r _i-1in (x, y), x and y is respectively row-coordinate and the row coordinate of data to be tested; The detection of extreme point comprises the detection to Local modulus maxima detection and local minizing point in data.For above-mentioned detection data r _i-1(x, y) adopts 2-D pattern to represent, so wherein a _mnrepresent r _i-1the pixel of middle m capable n row, adopts neighborhood method to detect extreme point (generally adopting the window of 3 × 3), by the Local modulus maxima that detects stored in 2-D array MAXMA, by the local minizing point that detects stored in 2-D array MINMA.So make r _i-1middle any point is a _xy(x, y), with a _xyadjacent pixel is $\&ForAll;a_{uv}=\left\{a_{uv}|x-1<u<x+1,y-1<v<y+1\right\},$ If a _xy> a _uv, then a _xy∈ MAXMA; If a _xy< a _uv, then a _xy∈ MINMA; If above condition does not all meet, so this point does not belong to Local Extremum.

Three, in order to make that algorithm does not need manual intervention and self-adaptation completes, therefore use MAXMA and MINMA obtained above to determine the size w of mean filter window.Use l ₁normed linear space manhatton distance under norm calculates MAXMA=(q ₁, q ₂, q ₃..., q _n) ^tin distance Distanse between each element _maxma, wherein q _nrepresent the n-th maximum point, n is positive integer; Then identical method is used to calculate MINMA=(p ₁, p ₂, p ₃..., p _n) ^t, in distance Distanse between each element _minma, wherein p _nrepresent the n-th minimum point, n is positive integer;

Distanse _maxma(q _i,q _j)＝|q _i-q _j|(1)

Distanse _minma(p _i,p _j)＝|p _i-p _j|(2)

Wherein { i|1 < i < n}, { j|1 < j < n}, i ≠ j.Then by Distanse _maxmaand Distanse _minmastored in array d _{max_l1}with in; Finally by calculating

d ₁＝min{min{d _{max_l1}},min{d _{min_l1}}}(3)

d ₂＝max{max{d _{max_l1}},max{d _{min_l1}}}(4)

Obtain the size w of two filtering windows of candidate, and by d ₁d ₂be defined as Type-1 and Type-2 respectively.

Four, the 2-D data of local extremum will have been detected as the object being about to carry out sieving

h _i(x,y)＝r _i-1(5)

Five, use mean filter from h _iaverage envelope m is obtained in (x, y) _mean(x, y).Mean filter is actually the process of a convolution, therefore uses the fast convolution algorithm based on svd to improve the counting yield of wave filter.The mask θ of the m × n used in Arbitrary Matrix M, first uses the decomposition based on singular value to obtain result

$\mathrm{SVD}\left(\mathrm{\&theta;}\right)=U\&times;\mathrm{\&Sigma;}\&times;V^{*}---\left(6\right)$

Wherein U is a unit matrix, and Σ is a diagonal matrix, V ^*for the conjugate matrices of unit matrix V, then obtained the order of mask θ by diagonal matrix Σ result after obtaining convolution the most at last

Wherein i-th row of U (i :) representing matrix U, V (:, i) i-th row of representing matrix V, * represents and is carrying out convolution operation.Final acquisition average envelope m _mean(x, y)=R.

Six, after acquisition average envelope, new signal h is calculated _i' (x, y).

h _i'(x,y)＝h _i(x,y)-m _mean(x,y)(8)

Seven, calculate the new signal generated and whether meet screening stopping criterion SD, if SD is less than given threshold xi (such as 0.05), then have h _i' (x, y) be i-th layer of BIMF; Otherwise, make h _i(x, y)=h _i' (x, y), jump to the 5th step;

$SD=\frac{{\mathrm{\&Sigma;}}_X^W{\mathrm{\&Sigma;}}_Y^H{|h_i^{\&prime;}(x,y)-h_i(x,y)|}^2}{{\mathrm{\&Sigma;}}_X^W{\mathrm{\&Sigma;}}_Y^H{\left|h_i(x,y)\right|}^2}---\left(9\right)$

Eight, residual error surplus r is judged _i(x, y)=r _i-1(x, y)-h _i' whether (x, y) be monotonic quantity, if so, then whole decomposable process terminates; If not, then make i=i+1 and jump to second step.

Compared with existing two-dimensional empirical mode decomposition method, tool of the present invention has the following advantages: this method is applied a mean filter based on fast convolution and obtained average envelope, average envelope is obtained based on RBF in conventional two-dimensional empirical mode decomposition method, greatly reduce time complexity, obtaining the computation rate that improve algorithm in similar experiment result situation, and application manhatton distance carrys out the window size of computation of mean values wave filter, such that whole algorithm does not need manual intervention and self-adaptation completes.

Accompanying drawing explanation

Fig. 1 fast two-dimensional empirical mode decomposition of the present invention process flow diagram;

Fig. 2 the present invention decomposes the 256 × 256Lena image used by example;

Fig. 3 the present invention decomposes the 256 × 256Elaina image used by example;

Fig. 4 detects extreme point result to one-dimensional data;

Fig. 5 detects Local Extremum result to 2-D data (Fig. 2); Wherein (a) represents Local modulus maxima testing result, and (b) represents local minizing point's testing result;

Fig. 6 uses the Contrast on effect of fast convolution and traditional convolution respectively to Fig. 2; A () uses traditional 2-D convolution of 27 × 27 windows; B () uses traditional 2-D convolution of 57 × 57 windows; C () uses the quick 2-D convolution of 27 × 27 windows; D () uses the quick 2-D convolution of 57 × 57 windows.

Fig. 7 uses the time of fast convolution and traditional convolution to contrast to Fig. 2 respectively;

Analysis result after Fig. 8 decomposes Fig. 3;

Output rusults after Fig. 9 decomposes Fig. 3.

Embodiment

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

Step one: first, determines the data will carrying out sieving.The data of generally carrying out sieving are the trial image of standard, as Lena image, and cameraman image, Elaine image etc.

Step 2: secondly, detects data r _i-1the Local Extremum of (x, y).The detection of extreme point comprises the detection to Local modulus maxima detection and local minizing point in data.For 2-D data neighborhood method is adopted to detect extreme point (generally adopting the window of 3 × 3), by the Local modulus maxima that detects stored in 2-D array MAXMA, by the local minizing point that detects stored in 2-D array MINMA.So make r ₀middle any point is a _xy(x, y), with a _xyadjacent pixel is $\&ForAll;a_{uv}=\left\{a_{uv}|x-1<u<x+1,y-1<v<y+1\right\},$ If a _xy> a _uv, then a _xy∈ MAXMA; If a _xy< a _uv, then a _xy∈ MINMA; If above condition does not all meet, so this point does not belong to Local Extremum.In order to make that algorithm does not need manual intervention and self-adaptation completes, therefore use MAXMA and MINMA obtained above to determine the size w of mean filter window.Use l ₁normed linear space manhatton distance under norm calculates MAXMA=(q ₁, q ₂, q ₃..., q _n) ^tin distance Distanse between each element _maxma; Then identical method is used to calculate MINMA=(p ₁, p ₂, p ₃..., p _n) ^tin distance Distanse between each element _minma;

Distanse _maxma(q _i,q _j)＝|q _i-q _j|(1)

Distanse _minma(p _i,p _j)＝|p _i-p _j|(2)

Wherein { i|1 < i < n}, { j|1 < j < n}, i ≠ j.Then by Distanse _maxmaand Distanse _minmastored in array with in; Finally by calculating

d ₁＝min{min{d _{max_l1}},min{d _{min_l1}}}(3)

d ₂＝max{max{d _{max_l1}},max{d _{min_l1}}}(4)

Obtain the size w of two filtering windows of candidate, and by d ₁d ₂be defined as Type-1 and Type-2 respectively.In 2-D data after having detected local extremum, as the object being about to carry out sieving

h _i(x,y)＝r _i-1(5)

Step 3: next, uses mean filter from h _iaverage envelope m is obtained in (x, y) _mean(x, y).Mean filter is actually the process of a convolution, therefore uses the fast convolution algorithm based on svd to improve the counting yield of wave filter.The mask θ of the m × n used in Arbitrary Matrix M, first uses the decomposition based on singular value to obtain result

$\mathrm{SVD}\left(\mathrm{\&theta;}\right)=U\&times;\mathrm{\&Sigma;}\&times;V^{*}---\left(6\right)$

Wherein U is a unit matrix, and Σ is a diagonal matrix, V ^*for the conjugate matrices of unit matrix V, then obtain the order of mask θ by diagonal matrix Σ result after obtaining convolution the most at last

Wherein i-th row of U (i :) representing matrix U, V (:, i) i-th row of representing matrix V, * represents and is carrying out convolution operation.Final acquisition average envelope m _mean(x, y)=R, after acquisition average envelope, calculates new signal h _i' (x, y).

h _i'(x,y)＝h _i(x,y)-m _mean(x,y)(8)

Step 4: last, calculates the new signal generated and whether meets screening stopping criterion SD, if SD is less than given threshold xi (such as 0.05), then have h _i' (x, y) be i-th layer of BIMF; Otherwise, make h _i(x, y)=h _i' (x, y), jump to step 3;

$SD=\frac{{\mathrm{\&Sigma;}}_X^W{\mathrm{\&Sigma;}}_Y^H{|h_i^{\&prime;}(x,y)-h_i(x,y)|}^2}{{\mathrm{\&Sigma;}}_X^W{\mathrm{\&Sigma;}}_Y^H{\left|h_i(x,y)\right|}^2}---\left(9\right)$

Judge residual error surplus r _i(x, y)=r _i-1(x, y)-h _i' whether (x, y) be monotonic quantity, if so, then whole decomposable process terminates; If not, then make i=i+1 and jump to step 2.

As shown in figs 1-9, the present invention has quoted from respective image and has analyzed, and exports corresponding analysis result simultaneously; This method is applied a mean filter based on fast convolution and is obtained average envelope, average envelope is obtained based on RBF in conventional two-dimensional empirical mode decomposition method, greatly reduce time complexity, obtaining the computation rate that improve algorithm in similar experiment result situation, and application manhatton distance carrys out the window size of computation of mean values wave filter, such that whole algorithm does not need manual intervention and self-adaptation completes.

Claims (1)

1. the fast two-dimensional empirical mode decomposition method improved, is characterized in that: the implementation process of the method is as follows,

1) data will carrying out sieving are determined; The data of carrying out sieving are the trial image of standard; If pending original image is f (m, n), wherein m is the row-coordinate of image, and n image column coordinate, makes i=1, and the data so will carrying out sieving are r ₀(x, y)=f (m, n), wherein x is the row-coordinate of image, y image column coordinate;

2) data r is detected _i-1the Local Extremum of (x, y); Detect data r _i-1in (x, y), x and y is respectively row-coordinate and the row coordinate of data to be tested; The detection of extreme point comprises the detection to Local modulus maxima detection and local minizing point in data; For above-mentioned detection data r _i-1(x, y) adopts 2-D pattern to represent, so wherein a _mnrepresent r _i-1the pixel of middle m capable n row, adopts neighborhood method to detect extreme point, by the Local modulus maxima that detects stored in 2-D array MAXMA, by the local minizing point that detects stored in 2-D array MINMA; So make r _i-1middle any point is a _xy(x, y), with a _xyadjacent pixel is $\&ForAll;a_{uv}=\left\{a_{uv}|x-1<u<x+1,y-1<v<y+1\right\},$ If a _xy> a _uv, then a _xy∈ MAXMA; If a _xy< a _uv, then a _xy∈ MINMA; If above condition does not all meet, so this point does not belong to Local Extremum;

3) in order to make that algorithm does not need manual intervention and self-adaptation completes, therefore use MAXMA and MINMA obtained above to determine the size w of mean filter window; Use l ₁normed linear space manhatton distance under norm calculates MAXMA=(q ₁, q ₂, q ₃..., q _n) ^tin distance Distanse between each element _maxma, wherein q _nrepresent the n-th maximum point, n is positive integer; Then identical method is used to calculate MINMA=(p ₁, p ₂, p ₃..., p _n) ^t, in distance Distanse between each element _minma, wherein p _nrepresent the n-th minimum point, n is positive integer;

Distanse _maxma(q _i,q _j)＝|q _i-q _j|(1)

Distanse _minma(p _i,p _j)＝|p _i-p _j|(2)

Wherein { i|1 < i < n}, { j|1 < j < n}, i ≠ j; Then by Distanse _maxmaand Distanse _minmastored in array with in; Finally by calculating

d ₁＝min{min{d _{max_l1}},min{d _{min_l1}}}(3)

d ₂＝max{max{d _{max_l1}},max{d _{min_l1}}}(4)

Obtain the size w of two filtering windows of candidate, and by d ₁d ₂be defined as Type-1 and Type-2 respectively;

4) the 2-D data of local extremum will have been detected as the object being about to carry out sieving

h _i(x,y)＝r _i-1(5)

5) use mean filter from h _iaverage envelope m is obtained in (x, y) _mean(x, y); Mean filter is actually the process of a convolution, therefore uses the fast convolution algorithm based on svd to improve the counting yield of wave filter; The mask θ of the m × n used in Arbitrary Matrix M, first uses the decomposition based on singular value to obtain result

SVD(θ)＝U×Σ×V ^*(6)

Wherein U is a unit matrix, and Σ is a diagonal matrix, V ^*for the conjugate matrices of unit matrix V, then obtained the order of mask θ by diagonal matrix Σ result after obtaining convolution the most at last

Wherein i-th row of U (i :) representing matrix U, V (:, i) i-th row of representing matrix V, * represents and is carrying out convolution operation; Final acquisition average envelope m _mean(x, y)=R;

6) after acquisition average envelope, new signal h ' is calculated _i(x, y);

h′ _i(x,y)＝h _i(x,y)-m _mean(x,y)(8)

7) calculate the new signal generated and whether meet screening stopping criterion SD, if SD is less than given threshold xi, then have h ' _i(x, y) is i-th layer of BIMF; Otherwise, make h _i(x, y)=h ' _i(x, y), jumps to step 5);

$SD=\frac{{\mathrm{\&Sigma;}}_X^W{\mathrm{\&Sigma;}}_Y^H{|h_i^{\&prime;}(x,y)-h_i(x,y)|}^2}{{\mathrm{\&Sigma;}}_X^W{\mathrm{\&Sigma;}}_Y^H{\left|h_i(x,y)\right|}^2}---\left(9\right)$

8) residual error surplus r is judged _i(x, y)=r _i-1(x, y)-h ' _iwhether (x, y) is monotonic quantity, and if so, then whole decomposable process terminates; If not, then make i=i+1 and jump to step 2).