CN101527037A - Method for denoising stationary wavelet image based on neighborhood windowing - Google Patents

Abstract

The invention belongs to the technical field of image processing of a computer and relates to a method for denoising stationary wavelet image based on neighborhood windowing, which comprises the following steps of: carrying out monoscale transform of a two-dimensional stationary wavelet to a noisy image f (x, y), decomposing the image to be one layer, respectively obtaining four sub-band coefficients of low-frequency coefficient A1, horizontal detail coefficient H1, vertical detail coefficient V1 and diagonal detail coefficient D1; keeping A1 to be constant; adopting a vertical linear filter template, a horizontal linear filter template and a diagonal direction filter template to filter the horizontal detail coefficient H1, the vertical detail coefficient V1 and the diagonal detail coefficient D1 into H1, V1 and D1; carrying out stationary wavelet reconstruction to A1 and high-frequency sub-bands of H1, V1 and D1 after being filtered; and then obtaining the denoised image with *=A1+H1+V1+D1. By adopting the method, the purpose of restoring the original image and improving denoising performance to the image can be achieved better.

Description

Method for denoising stationary wavelet image based on neighborhood window addition

Technical field

The invention belongs to technical field of image processing, be specifically related to a kind of method for denoising stationary wavelet image based on neighborhood window addition.

Background technology

Traditional denoising method mainly is to utilize noise to carry out with signal different characteristic distributions on frequency domain, adopt the low-pass filtering method will mainly be in the noise filtering of high-frequency region, but the details of image also mainly is distributed in high-frequency region simultaneously, though therefore the low-pass filtering method can reach the reduction anti noise, has also destroyed image detail.So keeping image detail as much as possible when reducing picture noise becomes an emphasis problem of studying in the image denoising.Because wavelet transformation has adopted the multiresolution analysis method, so application is wider in image denoising, as the wavelet threshold method that picture noise suppresses that is applicable to of people such as Donoho proposition.But this method is easy to generate vibration to picture signal reconstructed image edge after the detail coefficients behind the orthogonal wavelet transformation is carried out threshold process, causes the edge of image distortion.

Summary of the invention

The present invention be directed to the above-mentioned deficiency of prior art, a kind of method for denoising stationary wavelet image based on neighborhood window addition is provided.This method has been utilized the information redundancy of stationary wavelet conversion, is more conducive to find in the yardstick and the dependence between the wavelet coefficient between yardstick, thereby the parameter variance estimated accuracy that is based upon on the wavelet coefficient neighborhood is greatly improved.And also taken into full account correlativity in the layer of wavelet coefficient on this basis, selected the different window filtering templates that add according to its characteristic for each sub-band coefficients, recovered original image better thereby reach, improvement is to the purpose of the denoising performance of image.

Method for denoising stationary wavelet image based on neighborhood window addition of the present invention comprises the following steps:

Step 1: (x y) carries out the conversion of two-dimentional stationary wavelet single scale, obtains four sub-band coefficients respectively: low frequency coefficient A with noisy image f ₁, level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁

Step 2: with A ₁Remain unchanged, to level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁Adopt veritcal linearity filtering template, the linear filtering template of level and diagonal angle trend pass filtering template to carry out mean filter respectively, after the filtering be

Step 3: with A ₁With filtered high-frequency sub-band

Carry out stationary wavelet reconstruct, can obtain image after the denoising

$\hat{f}=A_1+{\hat{H}}_1+{\hat{V}}_1+{\hat{D}}_1.$

Above-mentioned step 2 can be carried out according to following method:

1) at first respectively according to level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁Frequency characteristic select corresponding mean filter template w respectively _ζζ=H, V, D;

2) carry out windowing respectively for the filtering template;

3) carry out filtering operation: $g(x,y)=\frac{1}{M}\underset{(x,y)\∈S}{\mathrm{\Σ}}f(x,y)w_{\mathrm{\ζ}}w\left(n\right)$ ζ＝H，V，D

Wherein S is so that (x y) is the set of the neighborhood mid point at center, and M is counting in the S, and (x y) is sub-band coefficients after the filtering to g.

In the above-mentioned step 1), for level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁Adopt respectively $\frac{1}{5}\left[\begin{array}{c}1\\ 1\\ 1^{*}\\ 1\\ 1\end{array}\right],\frac{1}{5}\left[\begin{array}{ccccc}1& 1& 1^{*}& 1& 1\end{array}\right],\frac{1}{5}\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1^{*}& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$ As the filtering template.

Method for denoising stationary wavelet image based on neighborhood window addition provided by the invention, made full use of the information redundancy of stationary wavelet conversion, and the characteristics of the interior correlativity of layer of wavelet coefficient, and as adding window template denoising method based on what stationary wavelet decomposed according to providing a kind of, reach higher Y-PSNR, had better image denoising effect.

Description of drawings

Fig. 1 the present invention is based on the overview flow chart of the method for denoising stationary wavelet image of neighborhood window addition.

Fig. 2 stationary wavelet conversion synoptic diagram.Fig. 2 (a) stationary wavelet filtering; Fig. 2 (b) is the stationary wavelet interpolation process.

Fig. 3 denoising master drawing of the present invention.Fig. 3 (a) is the former figure of denoising master drawing; Fig. 3 (b) adds the image of making an uproar for the former figure of master drawing; Fig. 3 (c)～(g) is an image after the use denoising method denoising of the present invention, and wherein, (c) window area is 5, adds rectangular window; (d) window area is 5, adds Hanning window; (e) window area is 5, adds quarter window; (f) window area is 5, adds the hamming window; (g) window area is 5, adds the Brackman window.

Embodiment

Below by drawings and Examples the present invention is further described.

1. image transformation

(x y) carries out the stationary wavelet conversion, obtains the wavelet coefficient matrix, decomposes 1 layer, and wavelet basis is the sym8 small echo will to add the image f that makes an uproar.

The stationary wavelet conversion all produces the wavelet coefficient of similar number at each yardstick, digital picture f (its decomposition formula is (synoptic diagram is as shown in Figure 2) for x, two-dimentional stationary wavelet conversion y):

$A_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}=\underset{n1}{\mathrm{\&Sigma;}}\underset{n2}{\mathrm{\&Sigma;}}{h_0}^{\&UpArrow;2j}(n_1-2k_1){h_0}^{\&UpArrow;2j}(n_2-2k_2)A_{j-1,n_1,{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}_2}$

$H_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}=\underset{n1}{\mathrm{\&Sigma;}}\underset{n2}{\mathrm{\&Sigma;}}{h_0}^{\&UpArrow;2j}(n_1-2k_1){g_0}^{\&UpArrow;2j}(n_2-2k_2)A_{j-1,n_1,{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}_2}$

$V_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}=\underset{n1}{\mathrm{\&Sigma;}}\underset{n2}{\mathrm{\&Sigma;}}{g_0}^{\&UpArrow;2j}(n_1-2k_1){h_0}^{\&UpArrow;2j}(n_2-2k_2)A_{j-1,n_1,{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}_2}$

$D_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}=\underset{n1}{\mathrm{\&Sigma;}}\underset{n2}{\mathrm{\&Sigma;}}{g_0}^{\&UpArrow;2j}(n_1-2k_1){g_0}^{\&UpArrow;2j}(n_2-2k_2)A_{j-1,n_1{,\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}_2}$

Wherein j is a decomposition scale, { h _kAnd { g _kBe respectively low pass and Hi-pass filter, h ₀ ^{↑ 2j}And g ₀ ^{↑ 2j}Be illustrated in h ₀, g ₀Insert 2j-1 individual zero between 2.J-1 yardstick tomographic image A _J-1Result after one deck wavelet decomposition is: low frequency coefficient

The level detail coefficient

The vertical detail coefficient

With the diagonal detail coefficient

Corresponding restructing algorithm is:

$A_{j-1,n_1,{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}_2}=\frac{1}{4}\underset{i=0}{\overset{3}{\mathrm{\&Sigma;}}}(\underset{k_1}{\mathrm{\&Sigma;}}\underset{k_2}{\mathrm{\&Sigma;}}{h}_1(n_1-2k_1-i)h_1(n_2-2k_2-i)A_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}$

$+\underset{k_1}{\mathrm{\&Sigma;}}\underset{k_2}{\mathrm{\&Sigma;}}{h}_1(n_1-2k_1-i)g_1(n_2-2k_2-i)H_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}$

$+\underset{k_1}{\mathrm{\&Sigma;}}\underset{k_2}{\mathrm{\&Sigma;}}{g}_1(n_1-2k_1-i)h_1(n_2-2k_2-i)V_{j,k_1,{\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">k</figure-callout>}}_2}$

$+\underset{k_1}{\mathrm{\&Sigma;}}\underset{k_2}{\mathrm{\&Sigma;}}{g}_1(n_1-2k_1-i)g_1(n_2-2k_2-i)D_{j,k_1,k_2})$

Here (x y) carries out the stationary wavelet conversion, decomposes 1 layer, obtains low frequency coefficient A respectively to image f ₁With high frequency detail coefficients H ₁, V ₁And D ₁

2. neighborhood window addition filtering

Wavelet transformation can recursively use low pass and Hi-pass filter to realize by the low frequency coefficient to same subband, means that wavelet coefficient is correlated with in a small neighbourhood, is called the interior correlativity of layer of wavelet coefficient.In the neighborhood of the bigger wavelet coefficient of value, may have one group of bigger wavelet coefficient.

To the wavelet coefficient individual processing in each subband, treatment step is as follows:

1) with A ₁Remain unchanged, respectively according to H ₁, V ₁And D ₁Frequency characteristic select corresponding mean filter template, after the filtering be

H wherein ₁The low-frequency information of image signal level direction and the high-frequency information of vertical direction have been comprised, and Gaussian noise is bigger at high frequency region noise energy proportion, carry out filtering so selected veritcal linearity filtering template, suc as formula (2), so both eliminate the noise signal of vertical direction, kept image edge information again largely; V ₁Then comprised the high-frequency information of image signal level direction and the low-frequency information of vertical direction, therefore selected the linear filtering template of level for use, suc as formula (4); Comprised high-frequency information, therefore adopted diagonal angle trend pass filtering template, as the formula (6) the angular direction.

$\frac{1}{3}\left[\begin{array}{c}1\\ 1^{*}\\ 1\end{array}\right]$ (formula 1) $\frac{1}{5}\left[\begin{array}{c}1\\ 1\\ 1^{*}\\ 1\\ 1\end{array}\right]$ (formula 2)

$\frac{1}{3}\left[\begin{array}{ccc}1& 1^{*}& 1\end{array}\right]$ (formula 3) $\frac{1}{5}\left[\begin{array}{ccccc}1& 1& 1^{*}& 1& 1\end{array}\right]$ (formula 4)

$\frac{1}{3}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1^{*}& 0\\ 0& 0& 1\end{array}\right]$ (formula 5) $\frac{1}{5}\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1^{*}& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$ (formula 6);

2) carry out windowing for the filtering template;

For improving the filtering performance of template, the filtering template is carried out windowing, select rectangular window, Hanning window, quarter window, hamming window and Blackman window here.The characteristic of various window function w (n) is as follows:

◆ rectangular window (Rectangular Window)

The time domain form can be expressed as:

Frequency domain characteristic is: $W_R\left(e^{\mathrm{j\&omega;}}\right)=e^{-j\left(\frac{N-1}{2}\right)\mathrm{\&omega;}}\frac{\sin \left(\frac{\mathrm{\&omega;N}}{2}\right)}{\sin \left(\frac{\mathrm{\&omega;}}{2}\right)}$

◆ Hanning window

The time domain form can be expressed as:

$\stackrel{\&CenterDot;}{w}\left(k\right)=0.5(1-\cos \left(2\mathrm{\&pi;}\frac{k}{n+1}\right))$ k＝1，2，…，N

Frequency domain characteristic is:

$W\left(\mathrm{\&omega;}\right)=\{0.5W_R\left(\mathrm{\&omega;}\right)+0.25[W_R(\mathrm{\&omega;}-\frac{2\mathrm{\&pi;}}{N-1})+W_R(\mathrm{\&omega;}+\frac{2\mathrm{\&pi;}}{N-1})\left]\right\}{e}^{-\mathrm{j\&omega;}\left(\frac{N-1}{2}\right)}$

Wherein, W _R(ω) be the amplitude-frequency characteristic function of rectangular window function.

◆ quarter window

Quarter window is the simplest frequency spectrum function W (e ^{J ω}) be non-negative a kind of window function.The time domain form of quarter window function can be expressed as:

When n is odd number

$w\left(k\right)=\left\{\begin{array}{cc}\frac{2k}{n+1},& 1\&le;k\&le;\frac{n+1}{2}\\ \frac{2(n-k+1)}{n+1},& \frac{n+1}{2}\&le;k\&le;n\end{array}\right.$

When n is even number

$w\left(k\right)=\left\{\begin{array}{cc}\frac{2k-1}{n},& 1\&le;k\&le;\frac{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}{2}\\ \frac{2(n-k+1)}{n},& \frac{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="A20091006845700092.png"\; state="\{\{state\}\}">n</figure-callout>}}{2}\&le;k\&le;n\end{array}\right.$

Frequency domain characteristic is:

$W_R\left(e^{\mathrm{j\&omega;}}\right)=e^{-\mathrm{j\&omega;}\left(\frac{N-1}{2}\right)}\frac{2}{N-1}{\left(\frac{\sin \left(\frac{\mathrm{\&omega;}(N-1)}{4}\right)}{\sin \left(\frac{\mathrm{\&omega;}}{2}\right)}\right)}^2$

◆ the hamming window function

The time domain form can be expressed as

$w\left(k\right)=0.54-0.46\cos \left(2\mathrm{\&pi;}\frac{k}{N-1}\right)$ k＝1，2，…，N

Frequency domain characteristic is

$W\left(\mathrm{\&omega;}\right)=0.54W_R\left(\mathrm{\&omega;}\right)+0.23[W_R(\mathrm{\&omega;}-\frac{2\mathrm{\&pi;}}{N-1})+W_R(\mathrm{\&omega;}+\frac{2\mathrm{\&pi;}}{N-1})]$

Wherein, W _R(ω) be the amplitude-frequency characteristic function of rectangular window function.

◆ the Blacknam window function

The time domain form can be expressed as

$w\left(k\right)=0.42-0.5\cos \left(2\mathrm{\&pi;}\frac{k-1}{N-1}\right)+0.08\cos \left(4\mathrm{\&pi;}\frac{k-1}{N-1}\right)$ k＝1，2，…，N

Frequency domain characteristic is

$W\left(\mathrm{\&omega;}\right)=0.42W_R\left(\mathrm{\&omega;}\right)+0.25[W_R(\mathrm{\&omega;}-\frac{2\mathrm{\&pi;}}{N-1})+W_R(\mathrm{\&omega;}+\frac{2\mathrm{\&pi;}}{N-1})]+0.04[W_R(\mathrm{\&omega;}-\frac{4\mathrm{\&pi;}}{N-1})+W_R(\mathrm{\&omega;}+\frac{4\mathrm{\&pi;}}{N-1})]$

Wherein, W _R(ω) be the amplitude-frequency characteristic function of rectangular window function.

4) carry out filtering operation.Use corresponding filtering template that each subband is carried out filtering operation.

$g(x,y)=\frac{1}{M}\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}f(x,y)w_{\mathrm{\&zeta;}}w\left(n\right)$ ζ＝H，V，D

Wherein S is so that (x y) is the set of the neighborhood mid point at center, and M is counting in the S, and (x y) is sub-band coefficients after the filtering to g.

3. image reconstruction

With A ₁With

Carry out stationary wavelet reconstruct, can obtain image after the denoising

$\hat{f}=A_1+{\hat{H}}_1+{\hat{V}}_1+{\hat{D}}_1.$

4. experimental result

In order to verify the validity of denoising method of the present invention, (shown in Fig. 3 (a)) tests to concrete picture.Adopt the sym8 small echo to carry out Flame Image Process in the experiment, to image in addition standard variance be 15 noise, image is decomposed 1 layer with stationary wavelet, and uses the filtering template of different windowings to carry out mean filter adding the image of making an uproar.As the criterion of anti-acoustic capability quality, experimental result is as shown in table 1 with PSNR (PeakSignal to Noise Ratio).

Use the PSNR/db of various window area filtering to compare under the different windowing situations of table 1

The data that provide from table 1 use the method for denoising stationary wavelet image based on neighborhood window addition that provides among the present invention can obtain higher Y-PSNR as can be seen.Simultaneously from Fig. 3 (c)～(g) handle the back image also as can be seen this method obtained denoising effect preferably.

Claims (3)

1. the method for denoising stationary wavelet image based on neighborhood window addition comprises the following steps:

Step 1: (x y) carries out the conversion of two-dimentional stationary wavelet single scale, decomposes 1 layer, obtains four sub-band coefficients respectively: low frequency coefficient A with noisy image f ₁, level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁

Step 2: with A ₁Remain unchanged, to level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁Adopt veritcal linearity filtering template, the linear filtering template of level and diagonal angle trend pass filtering template to carry out mean filter respectively, after the filtering be

Step 3: with A ₁With filtered high-frequency sub-band

Carry out stationary wavelet reconstruct, can obtain image after the denoising $\hat{f}=A_1+{\hat{H}}_1+{\hat{V}}_1+{\hat{D}}_1.$

2. the method for denoising stationary wavelet image based on neighborhood window addition according to claim 1, step 2 wherein can be carried out according to following method:

1) at first respectively according to level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁Frequency characteristic select corresponding mean filter template w respectively _ζζ=H, V, D;

2) carry out windowing respectively for the filtering template;

3) carry out filtering operation: $g(x,y)=\frac{1}{M}\underset{(x,y)\&Element;S}{\mathrm{\&Sigma;}}f(x,y)w_{\mathrm{\&zeta;}}w\left(n\right)$ ζ＝H，V，D

Wherein S is so that (x y) is the set of the neighborhood mid point at center, and M is counting in the S, and (x y) is sub-band coefficients after the filtering to g.

3. the method for denoising stationary wavelet image based on neighborhood window addition according to claim 2 is in the step 1) wherein, for level detail coefficient H ₁, vertical detail coefficient V ₁With diagonal detail coefficient D ₁Adopt respectively $\frac{1}{5}\left[\begin{array}{c}1\\ 1\\ 1*\\ 1\\ 1\end{array}\right],$ $\frac{1}{5}\left[\begin{array}{ccccc}1& 1& 1*& 1& 1\end{array}\right],$ $\frac{1}{5}\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1*& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$ As the filtering template.

Images