CN104240203A - Medical ultrasound image denoising method based on wavelet transform and quick bilateral filtering - Google Patents

Abstract

A medical ultrasound image denoising method based on wavelet transform and quick bilateral filtering comprises the following steps of the first step of establishing a medical ultrasound image model, the second step of carrying out wavelet transform on an image obtained in the first step after logarithm transformation to obtain four frequency domains (LL1, LH1, HL1and HH1), continuously carrying out wavelet transform on the low frequency domain LL1 to obtain four frequency domains (LL2, LH2, HL2 and HH2) again and then repeating the step until resolving the maximum layer number, the third step of carrying out quick bilateral filtering on the low frequency portion (LLJ) at the last layer, the fourth step of carrying out threshold value method shrinkage on wavelet coefficients of the high frequency portions (LHj, HLj and HHj) of all layers and the fifth step of carrying out wavelet inverse transformation to obtain the denoised medical ultrasound image. In addition, the J is 1, 2, ..., J. If denoised ultrasound envelope signals are wanted, exponential transformation is carried out on the ultrasound image obtained in the fifth step.

Description

Based on the medical ultrasound image denoising method of wavelet transformation and quick bilateral filtering

Technical field

The present invention is applied to Medical Image Denoising field, designs a kind of denoising method based on wavelet transformation and quick bilateral filtering being applicable to medical ultrasonic image.

Background technology

At medical imaging field, the imaging technique such as ultrasonic imaging, CT, MRI has been applied in medical clinic applications.Because ultrasonic imaging technique has without wound, "dead" infringement, the characteristic such as efficient and convenient, therefore ultrasonic imaging technique is safer relative to other imaging techniques.Especially, in observation pregnant woman body in the clinical practice such as fetus growth situation and diagnosis of abdominal organ lesion, the use of ultrasonic imaging technique is even more important.

Gave a women the statistics suffered from breast cancer according to American Cancer Society in 2013, in the past year, have in American Women 232340 example new infiltrative breast carcinoma cases and 39620 people die from breast cancer.The major technique carrying out detecting for the mammary gland of human body and other organs of belly is ultrasonic imaging technique, namely usually said B ultrasonic image.Therefore medical ultrasonic image quality is improved, for doctor provides more clear muting image to have very important significance.

Due to the restriction of ultrasonic imaging mechanism, the existence of speckle noise has had a strong impact on the quality of ultrasonoscopy, result in ultrasonograph quality poor.The generation of speckle noise is owing to there is a large amount of random scatter phenomenons in the basic resolution element in ultrasonic imaging, and image shows as different fleck relevant in spatial domain, and it will cover the very little characteristics of image of those gray scale difference.For clinician, the Accurate Diagnosis of speckle noise to them causes very large interference, is not particularly that the impact that very abundant doctor causes is larger for experience.Therefore, from the angle of clinical practice, need research to remove the algorithm of speckle noise, provide technical support for doctor makes diagnosis more accurately, reduce the risk of Artificial Diagnosis.

Due to the limitation of hospital resources, the quantity that particularly doctor carries out Artificial Diagnosis patient every day cannot meet the demand of social whole stratum, is namely faced with the few situation of the many doctors of patient.Therefore, the demand of various automatic diagnosis instrument is increasing, the appearance of automatic diagnosis instrument, can save doctor's resource on the one hand, more patient can be facilitated to diagnose on the other hand.Along with developing by leaps and bounds of society economy, people's own health situation but allows of no optimist, so the demand of people to domestic type medical treatment automatic diagnosis instrument is also very large, and such as domestic ultrasound image automatic diagnostic instrument etc.But ultrasonoscopy automatic diagnostic instrument is faced with the not high problem of picture quality equally, and automatic diagnostic instrument needs intellectual analysis ultrasonoscopy being done to the later stage, as feature extraction, rim detection and Images Classification identification etc.Therefore, from the angle of automated diagnostic technology, need research to remove the method for speckle noise, for the later stage Intelligent treatment of image provides technical guarantee, promote the development of automatic diagnostics.

In sum, Research of Medical Ultrasonic Image Denoising method has very important significance:

(1) improve the quality of medical ultrasonic image, improve visual effect;

(2) facilitate doctor to judge for focal area more exactly, reduce the risk of auxiliary diagnosis;

(3) promote the development of ultrasonoscopy automated diagnostic technology, there is immeasurable value.

Ultrasonic imaging principles and methods.Ultrasound wave is the sound wave of frequency higher than 20K hertz, has namely surmounted the upper limit of hearing of the mankind, and it has good directivity and very strong penetration power.Medical ultrasound image technology is exactly by ultrasonic signal as carrier, utilizes the technology such as electronics science, signal transacting to extract the acoustic signals that inside of human body reflects, then do imaging processing [.The realization of ultrasonic diagnostic technique, mainly utilizes the principle of reflection of ultrasonic signal, because ultrasonic signal has good directivity, and the physical characteristicss such as the reflection similar to light, scattering, decay and Doppler effect.By probe, ultrasonic signal is transmitted in body, ultrasonic signal can be propagated in various tissue, and the acoustic reactance of different tissues is discrepant, reflection and scattering is there is in ultrasonic signal on the surface of various tissue, and then by probe reception of echoes signal, strengthen process through some signals, finally obtain concrete ultrasonoscopy.Because the configuration of surface of various tissue is different with to ultrasonic impedance, then binding of pathological knowledge and clinical medicine experience, can the ill locus of patient, character or lesion degree etc. be made and being diagnosed more accurately.Current B ultrasonic diagnostic imaging technology is the technology that domestic use amount is clinically maximum.

The ultimate principle of B ultrasonic imaging is ultrasonic pulse detection of echoes principle.In order to show the power of echoed signal, B ultrasonic imaging adopts intensification modulation mode to show, and the echoed signal namely detected is stronger, and the brightness of B ultrasonic image is larger; The echoed signal detected is more weak, and the brightness of B ultrasonic image is darker.

B ultrasonic instrument primarily of probe, governor circuit, radiating circuit, receiving circuit, signal processing circuit and display etc. 6 part composition.The simple and easy schematic diagram of B ultrasonic basic imaging principle, as shown in Figure 1.

The basic process of B ultrasonic imaging is: governor circuit controls radiating circuit, radiating circuit produces driving pulse, after probe obtains this driving pulse, the Circuits System inner by probe produces ultrasound wave, and ultrasonic signal is transmitted in human body, then through one period of time delay, probe receives the acoustic signals (echoed signal) reflected by human body again, echoed signal through signal transacting such as a series of filtering such as receiving circuit, signal processing circuit, analog to digital conversion, is finally presented on terminal display with the form of two dimensional image again.

Owing to suppressing speckle noise to have very important significance, numerous researcher has dropped into a large amount of energy in this problem.The medical ultrasound image denoising method occurred in recent decades, can simply be divided into 5 types: self-adaptive solution method, anisotropy parameter denoising method, non-local mean denoising method, Noise Elimination from Wavelet Transform method and mixed type denoising method.Although the complexity of adaptive filter method is low by experiment, the often detail section of blurred picture, the inhibition for speckle noise is not very desirable.Anisotropy parameter denoising method, has very strong noise removal capability, but the possibility of result there will be the phenomenon of excess smoothness.Non local denoise algorithm is more satisfactory for the inhibition of speckle noise, but the complexity of this kind of denoising method is higher, not easily meets the requirement of real-time of medical ultrasonic imaging system, is often used for the later stage denoising of medical ultrasonic image.

Summary of the invention

The present invention will overcome the above-mentioned shortcoming of prior art, the feature of the model of nexus spot noise and the processing demands of medical ultrasonic image, proposes a kind of medical ultrasound image denoising method based on wavelet transformation and quick bilateral filtering.

Wavelet transformation has the superiority such as time frequency analysis and multiscale analysis, and it is widely used in image processing field.When processing additive noise problem, the denoising effect of small echo is better, can meet common product demand.But, only utilize the denoising method of wavelet transformation bad to the inhibition of speckle noise in medical ultrasonic image.For quick two-sided filter, it has very strong noise removal capability on the one hand when processing picture noise, can keep image edge details on the other hand.Simultaneously because traditional two-sided filter its efficiency in the process of denoising is lower, working time will be very long, is difficult to use in real-time system.Along with the resolution of image is increasing, this limits the application space of bilateral filtering to a large extent, how to realize quick bilateral filtering, reduces computing time and has great importance.Therefore the present invention utilizes quick two-sided filter to replace traditional two-sided filter, greatly can shorten working time on the basis of not weakening denoising effect.Therefore, the present invention is by the advantage in conjunction with Wavelet Denoising Method and quick bilateral filtering.Concrete thought is as follows, on the basis of traditional Wavelet noise-eliminating method, according to the statistical property of ultrasonoscopy in wavelet field and speckle noise, improves wavelet threshold contraction method, more effectively can remove the speckle noise of HFS.Because the low frequency part of medical ultrasonic image in wavelet field still exists speckle noise, therefore use denoising effect better and the quick two-sided filter that working time is shorter, at the image edge information suppressing can retain while noise in lower frequency region in lower frequency region.

Medical ultrasound image denoising method based on wavelet transformation and quick bilateral filtering of the present invention, comprises the following steps:

Step 1) foundation of medical ultrasonic image model;

The envelope signal of ultrasonic image-forming system collection is divided into two parts, and one is the reflected signal of significant in-vivo tissue, and another part is noise signal; Wherein noise signal can be divided into multiplicative noise and additive noise; Multiplicative noise is relevant with the principle of ultrasonic signal imaging, is mainly derived from random scattered signal; Additive noise thinks system noise; The envelope signal that ultrasonic image-forming system tentatively obtains is f ^pre, universal model is as follows:

f ^pre＝g ^pren ^pre+w ^pre (1)

Here, subscript ^prethe signal that expression system tentatively obtains; Function g ^prerepresent noise-free signal, n ^preand w ^prerepresent multiplicative noise and additive noise respectively, n in formula ^preit is the principal ingredient of noise;

With multiplicative noise n ^precompare, additive noise w ^preproportion is very little, therefore by w ^premodel after ignoring is

f ^pre＝g ^pren ^pre (2)

In order to adapt to the Dynamic Announce scope of ultrasonic image-forming system display screen, log-compressed process is carried out to the envelope signal that ultrasonic image-forming system collects; Formula (2) model be now multiplied will become the model of addition, as follows

log(f ^pre)＝log(g ^pre)+log(n ^pre) (3)

Signal log (the f obtained ^pre) be namely the medical ultrasonic image usually seen;

Step 2) image after the log-transformation that obtains the first step carries out wavelet decomposition, obtains four frequency domain (LL ¹, LH ¹, HL ¹and HH ¹); To lower frequency region LL ¹proceed wavelet decomposition, then obtain four frequency domain (LL ², LH ², HL ²and HH ²); Then this step is repeated, until decompose maximum number of plies J;

Because wavelet transformation is linear transformation, therefore formula (3) model obtains lower surface model after two-dimensional discrete wavelet conversion:

$W_{l,k}^j\left(\log \left(f^{\mathrm{pre}}\right)\right)=W_{l,k}^j\left(\log \left(g^{\mathrm{pre}}\right)\right)+W_{l,k}^j\left(\log \left(n^{\mathrm{pre}}\right)\right)---\left(4\right)$

Wherein with represent the wavelet coefficient containing the wavelet coefficient of noise image, the wavelet coefficient of noise-free picture and speckle noise respectively, wherein subscript j is the Decomposition order of wavelet transformation, and subscript (l, k) is the coordinate in wavelet field; Conveniently represent, formula (4) is reduced to

$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$

For discrete two dimensional image f (n, m), to the step that it carries out 2-d wavelet decomposition be: first one-dimensional discrete wavelet decomposition is carried out to every one-row pixels of image, then one-dimensional discrete wavelet decomposition is carried out to each row of image again, piece image is decomposed into four sub-band signals; The reconstruct of 2-d wavelet just can be able to obtain according to contrary order in addition;

Analyze sub-band component, i.e. corresponding coefficient of wavelet decomposition; LL0 is original signal, and the information of image all concentrates on here; Each wavelet decomposition all can obtain four sub-bands, obtains LL1, LH1, HL1 and HH1 tetra-sub-bands to LL0 after carrying out one-level wavelet decomposition;

LL1 component is the low frequency component obtained after carrying out wavelet decomposition to the columns and rows of original signal LL0, i.e. approximate part after one-level wavelet decomposition, and it contains the maximum low-frequency information of original image;

LH1 is the high fdrequency component in the vertical direction after a wavelet decomposition, and namely it contains the high-frequency informations such as the approximate information on image level direction and the edge in vertical direction;

HL1 is the high fdrequency component in the horizontal direction after a wavelet decomposition, and namely it contains the high-frequency informations such as the approximate information in image vertical direction and the edge in horizontal direction;

HH1 is that namely it contains the high-frequency informations such as the edge in image level and vertical direction to the high fdrequency component on angular direction after a wavelet decomposition;

The wavelet coefficient without noise cancellation signal after wavelet decomposition obey broad sense laplacian distribution, its probability distribution is as follows

$p_G\left(g\right)=\frac{v}{2\mathrm{s\Γ}(1/v)}\exp (-{\left|\frac{g-u}{s}\right|}^v),s,v>0,---\left(6\right)$

In formula, be gamma function, v is form parameter, and s is scale parameter, and u is location parameter; Work as v=1, during u=0, formula (6) will become laplacian distribution, and it is the particular module of broad sense laplacian distribution;

The wavelet coefficient of speckle noise simultaneously obedience zero-mean gaussian distributes

$p_N\left(n\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_N}\exp (-\frac{n^2}{2{\mathrm{\σ}}_N^2})---\left(7\right)$

σ in formula _nfor the standard deviation of noise in wavelet field;

Step 3) low frequency part (LL to last one deck ^j) carry out quick bilateral filtering process;

Quick two-sided filter is selected to do filtering process to the wavelet coefficient in lower frequency region; Quick two-sided filter is also known as increasing d type two-sided filter, the gray-scale value of pixel on each coordinate is added as three dimensions using the two-dimensional coordinate of image, form the linear convolution of three-dimensional Gaussian kernel function and 3-D view function, be multiplied corresponding on frequency domain, its result carries out Fourier inversion, calculates like this by loaded down with trivial details node-by-node algorithm being converted to Fast Fourier Transform (FFT);

The structure of quick two-sided filter is as follows

$\mathrm{BI}(x,y)=\frac{\mathrm{IY}(x,y)}{\mathrm{EY}(x,y)}=\frac{\mathrm{interp}(G\⊗\mathrm{IX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}{\mathrm{interp}(G\⊗\mathrm{EX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}---\left(8\right)$

In formula, (x, y) is the coordinate of input image pixels point; $\mathrm{IX}(x,y,z)=\left\{\begin{array}{c}z,z=I(x,y)\\ 0,z\≠I(x,y)\end{array}\right.$ The three-dimensional image matrix that representing input images I obtains after increasing dimension, $\mathrm{EX}(x,y,z)=\left\{\begin{array}{c}1,z=I(x,y)\\ 0,z\≠I(x,y)\end{array}\right.$ Represent three-dimensional weight matrix; Interp is interpolating function; G is the spatial neighbor degree factor G after linearization _swith gray scale similarity factor G _rproduct, i.e. gaussian kernel function; represent the linear convolution of matrix; s _sand s _rrepresent spatial domain sampling rate and gray scale territory sampling rate respectively;

Quick two-sided filter is characterized by and increases dimension matrix and increase the linear convolution tieing up kernel function, respectively three-dimensional Gaussian filtering is carried out to three-dimensional matrice IX and EX, be reduced to two-dimensional matrix IY, EY after again two filter result being carried out linear interpolation, IY is divided by EY pointwise and obtains restored image BI; Its process simplification is: after first interpolation, point removes, and ensure that the numerical value that can not to exist in EY matrix element close to zero, can Recovery image well; In three dimensions, the numerical solution of sampling, convolutional calculation, interpolation etc. achieves the acceleration of several order of magnitude;

Step 4) HFS (LH to every one deck ^j, HL ^jand HH ^j, j=1,2 ..., J) wavelet coefficient carry out threshold method shrink process;

In Wavelet noise-eliminating method, the selection of threshold function table directly can have influence on final image denoising result; When Threshold selection is less, the noise figure that a part is greater than this threshold value can be taken as useful signal and remain, and the image after this just causes denoising still exists much noise; When Threshold selection is larger, the useful information that a lot of coefficient is very little can be used as noise and zero setting, this becomes very level and smooth by making the image after denoising, loses a lot of detailed information; Therefore select appropriate wavelet threshold function extremely important; The people such as Donoho devise a general wavelet shrinkage threshold function table, namely here namely M is the overall number of wavelet coefficient in corresponding wavelet field; But this threshold function table is performed poor in the denoising of medical ultrasonic image, in order to obtain better denoising effect, universal threshold function is done following improvement by the present invention

$T_j=a_j{\mathrm{\σ}}_N\sqrt{2\log M}---\left(9\right)$

J in formula (=1,2 ..., J) and be the Decomposition order at wavelet coefficient place, J is the maximum decomposition level number of wavelet transformation, in the present invention, the auto-adaptive parameter a of j layer _jelect 2 as ^j-j+1;

In Wavelet noise-eliminating method, a first selected given threshold value, then shrinks wavelet coefficient according to certain rule, just completes the denoising to wavelet coefficient; An i.e. given threshold value, the coefficient that all absolute values are less than this threshold value is taken as noise, then does zero setting process to it; Wavelet coefficient absolute value being greater than to threshold value reduces by certain method, then obtains the new value after reducing;

Without the wavelet coefficient of noise cancellation signal obey broad sense laplacian distribution, the speckle noise part in wavelet field gaussian distributed; Select v=1, u=0, then formula (6) becomes laplacian distribution

$p_G\left(g\right)=\frac{1}{2s}\exp (-\frac{\left|g\right|}{s}),s>0---\left(10\right)$

In order to obtain the Signal estimation value in wavelet field, use the method that Bayesian MAP is estimated; In the computation process of posterior probability, use Bayesian formula as follows

$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f\right|g)\·p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N(f-g)\·p_G\left(g\right)\end{array}---\left(11\right)$

Bring formula (7), formula (10) into above formula (11), obtain

$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}\·\frac{1}{2\sqrt{2\mathrm{\π}}s{\mathrm{\σ}}_N}\×\\ \exp \left(\frac{2{\mathrm{\σ}}_N^2\left|g\right|-s{(f-g)}^2}{2s{\mathrm{\σ}}_N^2}\right)\end{array}---\left(12\right)$

In order to obtain maximum a posteriori probability, by ln (p _g|F(g|f)) g is asked to the equation zero setting of first order derivative, finally obtain

$\hat{g}=\mathrm{sign}\left(f\right)\·\max \left(\right|f|-\frac{{\mathrm{\σ}}_N^2}{s},0)---\left(13\right)$

for the estimation of g, and suppose f and without noise cancellation signal g jack per line; So just obtain new contraction method

$\hat{g}=\left\{\begin{array}{cc}0,& f<=T_j\\ \mathrm{sign}\left(f\right)\&CenterDot;\max \left(\right|f|-\frac{{\mathrm{\&sigma;}}_N^2}{s},0),& f>T_j\end{array}\right.---\left(14\right)$

Only have scale parameter s to be unknown in formula, can be determined by following formula

$s={\left[0.5({\mathrm{\&sigma;}}_{F,j}^2-{\mathrm{\&sigma;}}_N^2)\right]}^2---\left(15\right)$

Wherein σ _f,jfor noise image wavelet coefficient is in the standard deviation of j layer;

By observing the curve map of wavelet shrinkage function, can find out that the wavelet shrinkage function improved shows on curve image herein more level and smooth, especially when wavelet coefficient is greater than in the interval range of wavelet threshold;

Step 5) do wavelet inverse transformation process, obtain the medical ultrasonic image after denoising; If obtain the ultrasonic envelope signal after denoising, exponential transform is done to the ultrasonoscopy that the 5th step obtains.

Advantage of the present invention is:

One aspect of the present invention is improved general threshold value contraction method, has very strong noise removal capability to the tiny noise of high-frequency domain; On the other hand, because the present invention has carried out quick bilateral filtering process to low frequency part, therefore for the speckle noise (being present in low frequency part) that particle is larger, there is very strong rejection ability equally.Due to the time of the present invention consumption very major part be due to introduce the wave filter for the treatment of low frequency part noise, and traditional two-sided filter is long for working time in the process of denoising, so the introducing of quick two-sided filter can greatly improve operational efficiency while guarantee denoising effect in the present invention.Simultaneously for the feature of medical ultrasonic image, the method for this combination well can not only suppress speckle noise, can also retain the detail section at focus edge in image etc. simultaneously, doctor can better be helped to carry out illness analysis.

Accompanying drawing explanation

Fig. 1 B ultrasonic Image-forming instrument basic principle schematic

Fig. 2 image secondary of the present invention wavelet decomposition schematic diagram

Fig. 3 is soft/schematic diagram of hard threshold method

The wavelet shrinkage method curve that Fig. 4 the present invention improves

Fig. 5 overall flow figure of the present invention

The secondary decomposing schematic representation of Fig. 6 medical ultrasonic image of the present invention

The denoising effect of Fig. 7 emulating image of the present invention

The denoising result of Fig. 8 clinical ultrasound image of the present invention

Embodiment

With reference to accompanying drawing 2-8:

For making the object, technical solutions and advantages of the present invention more clear, below just technical scheme of the present invention is further described.

Medical ultrasound image denoising method based on wavelet transformation and quick bilateral filtering of the present invention, as shown in Figure 5, comprises the following steps:

Step 1) set up the model of medical ultrasonic image.

If think that the factor that ultrasonic image-forming system can affect acoustic power to those makes appropriate dynamic compensation, then the envelope signal of ultrasonic image-forming system collection is made up of two parts, and one is the reflected signal of significant in-vivo tissue, and another part is noise signal.Wherein noise signal can be divided into multiplicative noise and additive noise.Multiplicative noise is relevant with the principle of ultrasonic signal imaging, is mainly derived from random scattered signal.Additive noise thinks system noise, as the noise etc. of sensor.The envelope signal that ultrasonic image-forming system tentatively obtains is f ^pre, its universal model is as follows

f ^pre＝g ^pren ^pre+w ^pre (1)

Here, subscript ^prethe signal that expression system tentatively obtains.Function g ^prerepresent noise-free signal, n ^preand w ^prerepresent multiplicative noise and additive noise respectively, n in formula ^preit is the principal ingredient of noise.

With multiplicative noise n ^precompare, additive noise w ^preproportion is very little, therefore by w ^premodel after ignoring is

f ^pre＝g ^pren ^pre (2)

In order to adapt to the Dynamic Announce scope of ultrasonic image-forming system display screen, log-compressed process is carried out to the envelope signal that ultrasonic image-forming system collects.Formula (2) model be now multiplied will become the model of addition, as follows

log(f ^pre)＝log(g ^pre)+log(n ^pre) (3)

Now, the signal log (f obtained ^pre) be namely the medical ultrasonic image usually seen.If the target of process is common medical ultrasonic image, namely did log-transformation, then omitted this step.Speckle noise log (n after log-transformation ^pre) can be similar to be expressed as Gaussian noise.

Step 2) wavelet decomposition is carried out to the medical ultrasonic image after log-transformation.

Because wavelet transformation is linear transformation, therefore formula (3) model obtains lower surface model after two-dimensional discrete wavelet conversion:

$W_{l,k}^j\left(\log \left(f^{\mathrm{pre}}\right)\right)=W_{l,k}^j\left(\log \left(g^{\mathrm{pre}}\right)\right)+W_{l,k}^j\left(\log \left(n^{\mathrm{pre}}\right)\right)---\left(4\right)$

Wherein with represent the wavelet coefficient containing the wavelet coefficient of noise image, the wavelet coefficient of noise-free picture and speckle noise respectively.Wherein subscript j is the Decomposition order of wavelet transformation, and subscript (l, k) is the coordinate in wavelet field.Conveniently represent, formula (4) is reduced to

$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$

For discrete two dimensional image f (n, m), to the step that it carries out 2-d wavelet decomposition be: first one-dimensional discrete wavelet decomposition is carried out to every one-row pixels of image, then one-dimensional discrete wavelet decomposition is carried out to each row of image again, so just piece image is decomposed into four sub-band signals.The schematic diagram of secondary wavelet decomposition is carried out as shown in Figure 2 to image.

In like manner, the schematic diagram that 2-d wavelet is reconstructed and Fig. 2 closely similar, namely process according to contrary order, can obtain, therefore repeat no more here.Next by some the sub-band components in Fig. 2, namely corresponding coefficient of wavelet decomposition does simple analysis.LL0 is original signal, and the information of image all concentrates on here.Each wavelet decomposition all can obtain four sub-bands, obtains LL1, LH1, HL1 and HH1 tetra-sub-bands to LL0 after carrying out one-level wavelet decomposition.

LL1 component is the low frequency component obtained after carrying out wavelet decomposition to the columns and rows of original signal LL0, i.e. approximate part after one-level wavelet decomposition, it contains the maximum low-frequency information of original image.

LH1 is the high fdrequency component in the vertical direction after a wavelet decomposition, and namely it contains the high-frequency informations such as the approximate information on image level direction and the edge in vertical direction.

HL1 is the high fdrequency component in the horizontal direction after a wavelet decomposition, and namely it contains the high-frequency informations such as the approximate information in image vertical direction and the edge in horizontal direction.

HH1 is that namely it contains the high-frequency informations such as the edge in image level and vertical direction to the high fdrequency component on angular direction after a wavelet decomposition.

Fig. 6 is the secondary decomposing schematic representation of medical ultrasonic image, can find out that the speckle noise of medical ultrasonic image is mainly distributed in HFS.Low frequency part remains the bulk information of image, but still can there is a fraction of noise.

It is considered herein that the wavelet coefficient without noise cancellation signal after wavelet decomposition obey broad sense laplacian distribution, its probability distribution is as follows

$p_G\left(g\right)=\frac{v}{2\mathrm{s\&Gamma;}(1/v)}\exp (-{\left|\frac{g-u}{s}\right|}^v),s,v>0,---\left(6\right)$

In formula, be gamma function, v is form parameter, and s is scale parameter, and u is location parameter; Work as v=1, during u=0, formula (6) will become laplacian distribution, and it is the particular module of broad sense laplacian distribution;

The wavelet coefficient of speckle noise simultaneously obedience zero-mean gaussian distributes

$p_N\left(n\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_N}\exp (-\frac{n^2}{2{\mathrm{\&sigma;}}_N^2})---\left(7\right)$

σ in formula _nfor the standard deviation of noise in wavelet field.

Step 3) quick bilateral filtering process is carried out to the low frequency part of last one deck.

Generally based on the denoising method of small echo, the wavelet coefficient namely retaining lower frequency region (LL) is constant, only does threshold process to the wavelet coefficient of high-frequency domain (LH, HL, HH).But, perform poor when the method is applied to medical ultrasound image denoising.Through many experiments, find that the wavelet coefficient in lower frequency region still has a lot of speckle noise, in order to the speckle noise in more effectively filtering lower frequency region, the present invention selects quick two-sided filter to do filtering process to the wavelet coefficient in lower frequency region.

The structure of traditional two-sided filter is as follows

$\begin{array}{c}h\left(x\right)=k^{-1}\left(x\right)\&times;\\ \underset{\mathrm{\&xi;}\&Element;\mathrm{\&Omega;}\left(x\right)}{\&Integral;}f\left(\mathrm{\&xi;}\right)\&CenterDot;c(\mathrm{\&xi;},x)\&CenterDot;s(f\left(\mathrm{\&xi;}\right),f\left(x\right))\mathrm{d\&xi;}\end{array}---\left(16\right)$

In formula for normalized factor, Ω (x) is the window area centered by pixel x.

Two-sided filter is combined into by the kernel of two wave filters, i.e. the kernel c (ξ, x) of region filters and the kernel s (f (ξ), f (x)) of codomain wave filter.In formula, c (ξ, x) represents in Ω (x) scope, the distance function of surrounding pixel point ξ and center pixel x.S (f (ξ), f (x)) represents in Ω (x) scope, the similar function between the pixel value f (ξ) of surrounding pixel point and pixel value f (x) of central pixel point.The end value of pixel value f (x) of pixel after bilateral filtering centered by h (x).

Although traditional two-sided filter also can ensure good denoising effect, but its working time is long, be unfavorable for the application in real-time system, so quick two-sided filter is introduced the effect reaching medical ultrasonic image filtration in conjunction with Wavelet Denoising Method by the present invention, and can preserving edge information.

Quick two-sided filter is also known as increasing d type two-sided filter, its thought adds the gray-scale value of pixel on each coordinate as three dimensions using the two-dimensional coordinate of image, form the linear convolution of three-dimensional Gaussian kernel function and 3-D view function, be multiplied corresponding on frequency domain, its result carries out Fourier inversion, calculating loaded down with trivial details node-by-node algorithm being converted to Fast Fourier Transform (FFT) like this, effectively improving counting yield.

The structure of quick two-sided filter is as follows

$\mathrm{BI}(x,y)=\frac{\mathrm{IY}(x,y)}{\mathrm{EY}(x,y)}=\frac{\mathrm{interp}(G\&CircleTimes;\mathrm{IX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}{\mathrm{interp}(G\&CircleTimes;\mathrm{EX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}---\left(8\right)$

In formula, (x, y) is the coordinate of input image pixels point; $\mathrm{IX}(x,y,z)=\left\{\begin{array}{c}z,z=I(x,y)\\ 0,z\&NotEqual;I(x,y)\end{array}\right.$ The three-dimensional image matrix that representing input images I obtains after increasing dimension, $\mathrm{EX}(x,y,z)=\left\{\begin{array}{c}1,z=I(x,y)\\ 0,z\&NotEqual;I(x,y)\end{array}\right.$ Represent three-dimensional weight matrix; Interp is interpolating function; G is the spatial neighbor degree factor G after linearization _swith gray scale similarity factor G _rproduct, i.e. gaussian kernel function; represent the linear convolution of matrix; s _sand s _rrepresent spatial domain sampling rate and gray scale territory sampling rate respectively;

As can be seen from formula (8), quick two-sided filter is characterized by and increases dimension matrix and increase the linear convolution tieing up kernel function, respectively three-dimensional Gaussian filtering is carried out to three-dimensional matrice IX and EX, be reduced to two-dimensional matrix IY, EY after again two filter result being carried out linear interpolation, IY is divided by EY pointwise and obtains restored image BI.Its process simplification is: after first interpolation, point removes, and ensure that the numerical value that can not to exist in EY matrix element close to zero, can Recovery image well.In three dimensions, the numerical solution of sampling, convolutional calculation, interpolation etc. achieves the acceleration of several order of magnitude.

Step 4) threshold value shrink process is carried out to the wavelet coefficient of the HFS of every one deck.

In Wavelet noise-eliminating method, the selection of threshold function table directly can have influence on final image denoising result.When Threshold selection is less, the noise figure that a part is greater than this threshold value can be taken as useful signal and remain, and the image after this just causes denoising still exists much noise; When Threshold selection is larger, the useful information that a lot of coefficient is very little can be used as noise and zero setting, this becomes very level and smooth by making the image after denoising, loses a lot of detailed information.Therefore select appropriate wavelet threshold function extremely important.The people such as Donoho devise a general wavelet shrinkage threshold function table, namely here namely M is the overall number of wavelet coefficient in corresponding wavelet field.But this threshold function table is performed poor in the denoising of medical ultrasonic image, in order to obtain better denoising effect, universal threshold function is done following improvement by the present invention

$T_j=a_j{\mathrm{\&sigma;}}_N\sqrt{2\log M}---\left(9\right)$

J in formula (=1,2 ..., J) and be the Decomposition order at wavelet coefficient place, J is the maximum decomposition level number of wavelet transformation, in the present invention, the auto-adaptive parameter a of j layer _jelect 2 as ^j-j+1.

In Wavelet noise-eliminating method, a first selected given threshold value, then shrinks wavelet coefficient according to certain rule, just completes the denoising to wavelet coefficient.An i.e. given threshold value, the coefficient that all absolute values are less than this threshold value is taken as noise, then does zero setting process to it; Wavelet coefficient absolute value being greater than to threshold value reduces by certain method, then obtains the new value after reducing.

Classical wavelet shrinkage method has Soft thresholding and hard threshold method as Fig. 3, but in Soft thresholding, larger wavelet coefficient is always reduced by threshold value, therefore shrink after signal mathematical expectation with shrink before different, so process after image relative smooth some.The shortcoming of hard threshold method be wavelet coefficient near null value territory by unexpected zero setting, result in the uncontinuity of wavelet data, and this makes the variance of signal larger, these conversion are larger for the details impact in image.But in actual applications, when particularly noise level is very high, the image after hard threshold method process can produce concussion around point of discontinuity, the denoising effect of effect diagram picture.

Threshold value contraction method due to classics can not meet the requirement to medical ultrasound image denoising, so the present invention improves contraction method.

Without the wavelet coefficient of noise cancellation signal obey broad sense laplacian distribution, the speckle noise part in wavelet field gaussian distributed.In order to simplify calculating, select v=1 herein, then formula (6) becomes laplacian distribution

$p_G\left(g\right)=\frac{1}{2s}\exp (-\frac{\left|g\right|}{s}),s>0---\left(10\right)$

In order to obtain the Signal estimation value in wavelet field, use the method that Bayesian MAP is estimated.In the computation process of posterior probability, use Bayesian formula as follows

$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f\right|g)\&CenterDot;p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N(f-g)\&CenterDot;p_G\left(g\right)\end{array}---\left(11\right)$

Bring formula (7), formula (10) into above formula (11), obtain

$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}\&CenterDot;\frac{1}{2\sqrt{2\mathrm{\&pi;}}s{\mathrm{\&sigma;}}_N}\&times;\\ \exp \left(\frac{2{\mathrm{\&sigma;}}_N^2\left|g\right|-s{(f-g)}^2}{2s{\mathrm{\&sigma;}}_N^2}\right)\end{array}---\left(12\right)$

In order to obtain maximum a posteriori probability, by ln (p _g|F(g|f)) g is asked to the equation zero setting of first order derivative, finally obtain

$\hat{g}=\mathrm{sign}\left(f\right)\&CenterDot;\max \left(\right|f|-\frac{{\mathrm{\&sigma;}}_N^2}{s},0)---\left(13\right)$

for the estimation of g, and suppose f and without noise cancellation signal g jack per line.So just obtain new contraction method

$\hat{g}=\left\{\begin{array}{cc}0,& f<=T_j\\ \mathrm{sign}\left(f\right)\&CenterDot;\max \left(\right|f|-\frac{{\mathrm{\&sigma;}}_N^2}{s},0),& f>T_j\end{array}\right.---\left(14\right)$

Only have scale parameter s to be unknown in formula, can be determined by following formula

$s={\left[0.5({\mathrm{\&sigma;}}_{F,j}^2-{\mathrm{\&sigma;}}_N^2)\right]}^2---\left(15\right)$

Wherein σ _f,jfor noise image wavelet coefficient is in the standard deviation of j layer;

By observing the curve map Fig. 4 of the wavelet shrinkage function that the present invention improves, can find out that the wavelet shrinkage function improved shows on curve image herein more level and smooth, especially when wavelet coefficient is greater than in the interval range of wavelet threshold.

Step 5) do wavelet inverse transformation process, obtain the medical ultrasonic image after denoising.

The wavelet coefficient after denoising just can be obtained through threshold value shrink process and bilateral filtering process, in order to obtain the ultrasonoscopy after denoising, need to carry out wavelet inverse transformation to wavelet coefficient, thus the image that can obtain being beneficial to after the denoising that doctor analyzes, also demonstrate the present invention by experiment and really can meet requirement for medical ultrasound image denoising.

Experiment results.

In order to evaluate denoising method in this paper objectively, using Y-PSNR (PSNR), structural similarity (SSIM), FoM (Pratt ' s Figure of Merit) and working time as image quality evaluation standard.The computing formula of Y-PSNR is as follows

$\mathrm{PSNR}(X,\hat{X})=101g\left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(17\right)$

In formula, for the estimated value of signal X, MSE is obtained by formulae discovery below

$\mathrm{MSE}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(X_{i,j}-{\hat{X}}_{i,j})}^2---\left(18\right)$

Here M, N represent length and the width of 2D signal X respectively.

Structural similarity can quantize two width images difference structurally, and formula is defined as follows

$\mathrm{SSIM}(X,\hat{X})=\frac{(2{\mathrm{\&mu;}}_X{\mathrm{\&mu;}}_{\hat{X}}+c_1)({2\mathrm{\&sigma;}}_{X,\hat{X}}+c_2)}{({\mathrm{\&mu;}}_X^2+{\mathrm{\&mu;}}_{\hat{X}}^2+c_1)({\mathrm{\&sigma;}}_X^2+{\mathrm{\&sigma;}}_{\hat{X}}^2+c_2)}---\left(19\right)$

In formula, μ _x, with average and the variance of reference picture and estimated image respectively. be X and covariance, c ₁and c ₂for constant.Work as c ₁and c ₂when being all chosen as positive number, the span of SSIM is [0 1], and wherein 1 is best result, represents that the structure of two width figure is identical.

FoM can compare the rim detection quality of denoising image objectively, and formula is defined as follows

$\mathrm{FoM}(X,\hat{X})=\frac{1}{\max (N_X,N_{\hat{X}})}\underset{i=1}{\overset{N_X}{\mathrm{\&Sigma;}}}\frac{1}{1+{\mathrm{\&alpha;d}}_i^2}---\left(20\right)$

In formula, N _xwith represent desirable with the actual edge pixel number detected respectively.α is constant (usually getting α=1/9), d _ibe expressed as the distance of the i-th edge pixel point to nearest ideal edge pixel.The span of FoM is [0 1], and wherein 1 is best result, is expressed as the image border detected consistent with desirable image border.Here Canny detection algorithm (standard deviation value σ=3 of Gaussian filter) is used during detection of edge pixels.

Certainly in order to better represent advantage of the present invention, next mainly by contrast experiment, quantize to compare each picture appraisal standard to make objects and advantages of the present invention more clear.

In order to the effect of quantitative predication denoising of the present invention, first carry out the experiment of emulating image.Obtain as the data in Fig. 7 and table 1.

The Performance comparision of table 1 denoising method

Image for emulating obtains (speckle noise of medical ultrasonic image meets Gaussian distribution) by adding Gaussian noise to noise-free picture.The data quantitatively drawn by emulation experiment can be found out, the effect in speckle noise is being suppressed to be not fine iff utilizing the denoising method of wavelet transformation, certain speckle noise is still there is in low frequency part, go to do bilateral filtering process to low frequency part so introduce bilateral filtering, certain denoising effect is significantly improved, but can find that working time is longer, be unfavorable for the application in real-time system.So the present invention improves this, introduce quick two-sided filter and replace traditional two-sided filter, by relatively finding, this method shortens the time of operation greatly while not weakening denoising effect.

Clinical ultrasound image is tested, selection be the ultrasonoscopy Fig. 8 of mammary glands in women tissue containing focus.

Owing to there is not muting ultrasonoscopy in reality, therefore, the quality index such as PSNR here cannot use effectively, so just introduce another index, and non-reference picture quality index (NIQE).Obtain table 2

The Performance comparision of table 2 denoising method

By relatively finding that the present invention is being applied in medical ultrasonic image, denoising effect is not significantly weakened, and obtained by emulation experiment the efficiency that the present invention greatly can improve operation, therefore can be applied greatly in the middle of Real-time System.The deficiency that before solving, in method, working time is long.

Content described in this instructions embodiment is only enumerating the way of realization of inventive concept; protection scope of the present invention should not be regarded as being only limitted to the concrete form that embodiment is stated, protection scope of the present invention also and conceive the equivalent technologies means that can expect according to the present invention in those skilled in the art.

Claims (1)

1., based on the medical ultrasound image denoising method of wavelet transformation and quick bilateral filtering, comprise the following steps:

Step 1) foundation of medical ultrasonic image model;

The envelope signal of ultrasonic image-forming system collection is divided into two parts, and one is the reflected signal of significant in-vivo tissue, and another part is noise signal; Wherein noise signal can be divided into multiplicative noise and additive noise; Multiplicative noise is relevant with the principle of ultrasonic signal imaging, is mainly derived from random scattered signal; Additive noise thinks system noise; The envelope signal that ultrasonic image-forming system tentatively obtains is f ^pre, universal model is as follows:

f ^pre＝g ^pren ^pre+w ^pre (1)

Here, subscript ^prethe signal that expression system tentatively obtains; Function g ^prerepresent noise-free signal, n ^preand w ^prerepresent multiplicative noise and additive noise respectively, n in formula ^preit is the principal ingredient of noise;

With multiplicative noise n ^precompare, additive noise w ^preproportion is very little, therefore by w ^premodel after ignoring is

f ^pre＝g ^pren ^pre (2)

In order to adapt to the Dynamic Announce scope of ultrasonic image-forming system display screen, log-compressed process is carried out to the envelope signal that ultrasonic image-forming system collects; Formula (2) model be now multiplied will become the model of addition, as follows

log(f ^pre)＝log(g ^pre)+log(n ^pre) (3)

Signal log (the f obtained ^pre) be namely the medical ultrasonic image usually seen;

Step 2) image after the log-transformation that obtains the first step carries out wavelet decomposition, obtains four frequency domain (LL ¹, LH ¹, HL ¹and HH ¹); To lower frequency region LL ¹proceed wavelet decomposition, then obtain four frequency domain (LL ², LH ², HL ²and HH ²); Then this step is repeated, until decompose maximum number of plies J;

Because wavelet transformation is linear transformation, therefore formula (3) model obtains lower surface model after two-dimensional discrete wavelet conversion:

$W_{l,k}^j\left(\log \left(f^{\mathrm{pre}}\right)\right)=W_{l,k}^j\left(\log \left(g^{\mathrm{pre}}\right)\right)+W_{l,k}^j\left(\log \left(n^{\mathrm{pre}}\right)\right)---\left(4\right)$

Wherein with represent the wavelet coefficient containing the wavelet coefficient of noise image, the wavelet coefficient of noise-free picture and speckle noise respectively; Wherein subscript j is the Decomposition order of wavelet transformation, and subscript (l, k) is the coordinate in wavelet field; Conveniently represent, formula (4) is reduced to

$F_{l,k}^j=G_{l,k}^j+N_{l,k}^j---\left(5\right)$

For discrete two dimensional image f (n, m), to the step that it carries out 2-d wavelet decomposition be: first one-dimensional discrete wavelet decomposition is carried out to every one-row pixels of image, then one-dimensional discrete wavelet decomposition is carried out to each row of image again, piece image is decomposed into four sub-band signals; The reconstruct of 2-d wavelet just can be able to obtain according to contrary order in addition;

Analyze sub-band component, i.e. corresponding coefficient of wavelet decomposition; LL0 is original signal, and the information of image all concentrates on here; Each wavelet decomposition all can obtain four sub-bands, obtains LL1, LH1, HL1 and HH1 tetra-sub-bands to LL0 after carrying out one-level wavelet decomposition;

LL1 component is the low frequency component obtained after carrying out wavelet decomposition to the columns and rows of original signal LL0, i.e. approximate part after one-level wavelet decomposition, and it contains the maximum low-frequency information of original image;

LH1 is the high fdrequency component in the vertical direction after a wavelet decomposition, and namely it contains the high-frequency informations such as the approximate information on image level direction and the edge in vertical direction;

HL1 is the high fdrequency component in the horizontal direction after a wavelet decomposition, and namely it contains the high-frequency informations such as the approximate information in image vertical direction and the edge in horizontal direction;

HH1 is that namely it contains the high-frequency informations such as the edge in image level and vertical direction to the high fdrequency component on angular direction after a wavelet decomposition;

The wavelet coefficient without noise cancellation signal after wavelet decomposition obey broad sense laplacian distribution, its probability distribution is as follows

$p_G\left(g\right)=\frac{v}{2\mathrm{s\&Gamma;}(1/v)}\exp (-{\left|\frac{g-u}{s}\right|}^v),s,v>0,---\left(6\right)$

In formula, be gamma function, v is form parameter, and s is scale parameter, and u is location parameter; Work as v=1, during u=0, formula (6) will become laplacian distribution, and it is the particular module of broad sense laplacian distribution;

The wavelet coefficient of speckle noise simultaneously obedience zero-mean gaussian distributes

$p_N\left(n\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_N}\exp (-\frac{n^2}{2{\mathrm{\&sigma;}}_N^2})---\left(7\right)$

σ in formula _nfor the standard deviation of noise in wavelet field;

Step 3) low frequency part (LL to last one deck ^j) carry out quick bilateral filtering process;

Quick two-sided filter is selected to do filtering process to the wavelet coefficient in lower frequency region; Quick two-sided filter is also known as increasing d type two-sided filter, the gray-scale value of pixel on each coordinate is added as three dimensions using the two-dimensional coordinate of image, form the linear convolution of three-dimensional Gaussian kernel function and 3-D view function, be multiplied corresponding on frequency domain, its result carries out Fourier inversion, calculates like this by loaded down with trivial details node-by-node algorithm being converted to Fast Fourier Transform (FFT);

The structure of quick two-sided filter is as follows

$\mathrm{BI}(x,y)=\frac{\mathrm{IY}(x,y)}{\mathrm{EY}(x,y)}=\frac{\mathrm{interp}(G\&CircleTimes;\mathrm{IX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}{\mathrm{interp}(G\&CircleTimes;\mathrm{EX},\frac{x}{s_s},\frac{y}{s_s},\frac{I(x,y)}{s_r})}---\left(8\right)$

In formula, (x, y) is the coordinate of input image pixels point; $\mathrm{IX}(x,y,z)=\left\{\begin{array}{c}z,z=I(x,y)\\ 0,z\&NotEqual;I(x,y)\end{array}\right.$ The three-dimensional image matrix that representing input images I obtains after increasing dimension, $\mathrm{EX}(x,y,z)=\left\{\begin{array}{c}1,z=I(x,y)\\ 0,z\&NotEqual;I(x,y)\end{array}\right.$ Represent three-dimensional weight matrix; Interp is interpolating function; G is the spatial neighbor degree factor G after linearization _swith gray scale similarity factor G _rproduct, i.e. gaussian kernel function; represent the linear convolution of matrix; s _sand s _rrepresent spatial domain sampling rate and gray scale territory sampling rate respectively;

Quick two-sided filter is characterized by and increases dimension matrix and increase the linear convolution tieing up kernel function, respectively three-dimensional Gaussian filtering is carried out to three-dimensional matrice IX and EX, be reduced to two-dimensional matrix IY, EY after again two filter result being carried out linear interpolation, IY is divided by EY pointwise and obtains restored image BI; Its process simplification is: after first interpolation, point removes, and ensure that the numerical value that can not to exist in EY matrix element close to zero, can Recovery image well; In three dimensions, the numerical solution of sampling, convolutional calculation, interpolation etc. achieves the acceleration of several order of magnitude;

Step 4) HFS (LH to every one deck ^j, HL ^jand HH ^j, j=1,2 ..., J) wavelet coefficient carry out threshold method shrink process;

In Wavelet noise-eliminating method, the selection of threshold function table directly can have influence on final image denoising result; When Threshold selection is less, the noise figure that a part is greater than this threshold value can be taken as useful signal and remain, and the image after this just causes denoising still exists much noise; When Threshold selection is larger, the useful information that a lot of coefficient is very little can be used as noise and zero setting, this becomes very level and smooth by making the image after denoising, loses a lot of detailed information; Therefore select appropriate wavelet threshold function extremely important; The people such as Donoho devise a general wavelet shrinkage threshold function table, namely here namely M is the overall number of wavelet coefficient in corresponding wavelet field; Then, this threshold function table is performed poor in the denoising of medical ultrasonic image, and in order to obtain better denoising effect, universal threshold function is done following improvement by the present invention

$T_j=a_j{\mathrm{\&sigma;}}_N\sqrt{2\log M}---\left(9\right)$

J in formula (=1,2 ..., J) and be the Decomposition order at wavelet coefficient place, J is the maximum decomposition level number of wavelet transformation, in the present invention, the auto-adaptive parameter a of j layer _jelect 2 as ^j-j+1;

In Wavelet noise-eliminating method, a first selected given threshold value, then shrinks wavelet coefficient according to certain rule, just completes the denoising to wavelet coefficient; An i.e. given threshold value, the coefficient that all absolute values are less than this threshold value is taken as noise, then does zero setting process to it; Wavelet coefficient absolute value being greater than to threshold value reduces by certain method, then obtains the new value after reducing;

Without the wavelet coefficient of noise cancellation signal obey broad sense laplacian distribution, the speckle noise part in wavelet field gaussian distributed; Select v=1, u=0, then formula (6) becomes laplacian distribution

$p_G\left(g\right)=\frac{1}{2s}\exp (-\frac{\left|g\right|}{s}),s>0---\left(10\right)$

In order to obtain the Signal estimation value in wavelet field, use the method that Bayesian MAP is estimated; In the computation process of posterior probability, use Bayesian formula as follows

$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}{p}_{F|G}\left(f\right|g)\&CenterDot;p_G\left(g\right)=\\ \frac{1}{p_F\left(f\right)}{p}_N(f-g)\&CenterDot;p_G\left(g\right)\end{array}---\left(11\right)$

Bring formula (7), formula (10) into above formula (11), obtain

$\begin{array}{c}{p}_{G|F}\left(g\right|f)=\frac{1}{p_F\left(f\right)}\&CenterDot;\frac{1}{2\sqrt{2\mathrm{\&pi;}}s{\mathrm{\&sigma;}}_N}\&times;\\ \exp \left(\frac{2{\mathrm{\&sigma;}}_N^2\left|g\right|-s{(f-g)}^2}{2s{\mathrm{\&sigma;}}_N^2}\right)\end{array}---\left(12\right)$

In order to obtain maximum a posteriori probability, by ln (p _g|F(g|f)) g is asked to the equation zero setting of first order derivative, finally obtain

$\hat{g}=\mathrm{sign}\left(f\right)\&CenterDot;\max \left(\right|f|-\frac{{\mathrm{\&sigma;}}_N^2}{s},0)---\left(13\right)$

for the estimation of g, and suppose f and without noise cancellation signal g jack per line; So just obtain new contraction method

$\hat{g}=\left\{\begin{array}{cc}0,& f<=T_j\\ \mathrm{sign}\left(f\right)\&CenterDot;\max \left(\right|f|-\frac{{\mathrm{\&sigma;}}_N^2}{s},0),& f>T_j\end{array}\right.---\left(14\right)$

Only have yardstick s to be unknown in formula, can be determined by following formula

$s={\left[0.5({\mathrm{\&sigma;}}_{F,j}^2-{\mathrm{\&sigma;}}_N^2)\right]}^2---\left(15\right)$

Wherein σ _f,jfor noise image wavelet coefficient is in the standard deviation of j layer;

By observing the curve map of wavelet shrinkage function, can find out that the wavelet shrinkage function improved shows on curve image herein more level and smooth, especially when wavelet coefficient is greater than in the interval range of wavelet threshold;

Step 5) do wavelet inverse transformation process, obtain the medical ultrasonic image after denoising; If obtain the ultrasonic envelope signal after denoising, exponential transform is done to the ultrasonoscopy that the 5th step obtains.