CN102968779B - Non-subsample Contourlet domain type MRI (Magnetic Resonance Imaging) image enhancing method based on FCM (Fuzzy C-means) cluster - Google Patents

Abstract

The invention relates to a method for enhancing an MRI (Magnetic Resonance Imaging) image through an FCM (Fuzzy C-means) cluster under a non-subsample Contourlet domain, in particular relates to a non-subsample Contourlet domain type MRI image enhancing method based on the FCM cluster. Compared with the prior art, the method provided by the invention has the following advantages: 1, the non-subsample Contourlet transform is carried out, so that the edge, the profile information and the texture detailed information in the MRI image can be effectively captured; and 2, the FCM cluster method is introduced, so that the high-frequency directional sub-band coefficient is divided into noise coefficient, weak edge coefficient and high edge coefficient in a self-adapting way, and the problem that the conventional image enhancing method based on transform domain is carried out according to the threshold selection can be overcome.

Description

Based on the non-downsampling Contourlet territory MRI image enchancing method of FCM cluster

Technical field

The inventive method relates on non-downsampling Contourlet territory, utilize fuzzy C-mean algorithm FuzzyC-Means, FCM cluster, to magnetic resonance imaging Magnetic Resonance Imaging, namely MRI image carries out the method that strengthens, especially a kind of non-downsampling Contourlet territory MRI image enchancing method based on FCM cluster.

Background technology

MRI (magnetic resonance imaging Magnetic Resonance Imaging) is a kind of damage-free type medical imaging, can reconstruct the inner structure of human body, be widely used in medical diagnosis.Details is abundant, MRI image clearly, and doctor can be helped to examine parts such as the main organs of human body, vein, soft tissue and focuses better.Because imaging mechanism has its singularity, the light and shade contrast of MRI image is comparatively large, and with noise, weak minutia is easy to be covered, thus affects the diagnosis of doctor.Strengthen MRI image, outstanding details is beneficial to diagnosis and necessitates.

Existing image enchancing method mainly contains histogram equalization method, unsharp masking method based on spatial domain, and based on the method for transform domain.When strengthening medical image, the image enchancing method based on spatial domain has some shortcomings.Because the light and shade contrast of medical image is comparatively large, from the histogram equalization method that the overall intensity Distribution value of image is considered, obvious water washing effect can be produced, the interpretability of effect diagram picture in enhancing image.Although unsharp masking method effectively can strengthen the minutia in medical image, the method is but very responsive to noise, even if there is the noise of seldom amount to be present in image, also can be enhanced, affects the diagnosis of doctor.Noise, weak edge and strong fringing coefficient are distinguished by choosing suitable threshold value by the image enchancing method based on transform domain, and to noise figure zero setting, strengthen weak fringing coefficient and retain the object that strong fringing coefficient reaches enhancing image.In medical image enhancement, compared to the image enchancing method based on spatial domain, the method based on transform domain has certain advantage.But when lacking the priori of medical image, suitable threshold value becomes and is difficult to select.

Because wavelet transformation has good time domain, frequency-domain analysis characteristic, be widely used in MRI medical image enhancement.Along with going deep into of research, it is found that wavelet analysis effectively cannot capture the geological information of essence in image.The non-downsampling Contourlet conversion proposed for 2005 is a kind of multi-scale geometric analysis method, it compensate for the deficiency of wavelet analysis, accurately can catch edge, profile information and the texture information in image, be suitable for expressing the MRI image having and enrich detailed information.Therefore, be necessary to study a kind of new MRI image enchancing method, that expects that it can effectively utilize non-downsampling Contourlet conversion to detailed information well portrays ability, and can overcome the existing dependence chosen threshold value based on the image enchancing method of transform domain.

Summary of the invention

The object of the invention is to overcome the deficiencies in the prior art, namely when lacking the priori of medical image, the suitable threshold value adopted based on the image enchancing method of transform domain is difficult to the problem selected, propose a kind of non-downsampling Contourlet territory MRI image enchancing method based on FCM cluster, the quality of MRI image can be improved, thus be more conducive to the diagnosis of doctor.

Based on a non-downsampling Contourlet territory MRI image enchancing method for FCM cluster, its special feature is, comprises the steps:

(1) MRI image I (m, n) that size is M × N is inputted, wherein 1≤m≤M, 1≤n≤N, wherein M and N is the natural number being greater than 1, and carries out L layer non-downsampling Contourlet conversion to it, obtains the high frequency direction sub-band coefficients D on each yardstick _l,i(m, n) and low frequency sub-band coefficient, wherein 0≤l≤L-1,1≤i≤k _l, k _lrepresent at yardstick 2 ^-lon the number of high frequency direction subband, D _l,i(m, n) represents at yardstick 2 ^-lon i-th high frequency direction sub-band coefficients, L is 3 ~ 5, L is natural number;

(2) pass through the high frequency direction sub-band coefficients D on each yardstick _{l, i}(m, n) carries out cluster respectively, is divided into noise, weak edge and strong fringing coefficient;

(3) gain rule formula is utilized, to the high frequency direction sub-band coefficients D obtaining category attribute _l,i(m, n) revises, that is:

In formula, by setup parameter a ₀=0, to noise figure zero setting, by setup parameter a ₂=1.25, appropriateness is carried out to strong fringing coefficient and strengthens, D ' _{l, i}(m, n) is revised high frequency direction sub-band coefficients, A _lfor at yardstick 2 ^-lon the maximum modulus value of high frequency direction sub-band coefficients, a ₁[] is nonlinear gain function, has:

$a_1[D_{l,i}(m,n)/A_l]=\frac{1}{1+e^{(-t\·D_{l,i}(m,n))/A_l}}-\frac{1}{1+e^{(t\·D_{l,i}(m,n))/A_l}};$

In formula, parametric t is for the degree of control MRI image enhaucament, and t value is larger, and MRI image enhaucament degree is larger, and for avoiding MRI image to cross enhancing, arrange t=7, e is natural constant, e=2.71828;

(4) non-downsampling Contourlet inverse transformation is done to the corrected high frequency direction sub-band coefficients obtained in the low frequency sub-band coefficient obtained in step (1) and step (3), MRI image I ' (m, n) that can be enhanced.

Passing through the high frequency direction sub-band coefficients D on each yardstick wherein described in step (2) _{l, i}(m, n) carries out cluster respectively, is divided into noise, weak edge and strong fringing coefficient, specifically carries out according to following process:

1) the high frequency direction sub-band coefficients D on each yardstick is calculated _{l, i}the three-dimensional feature vector Z of (m, n) ^l(m, n), 1≤n≤M, 1≤n≤N, 0≤l≤L-1, definition for characteristics of mean, for standard deviation characteristic, for maximum norm value tag, $Z^l(m,n)={[Z_{\mathrm{mean}}^l(m,n),Z_{\mathrm{std}}^l(m,n),Z_{\max }^l(m,n)]}^T,$ 0≤l≤L-1；

2) by data set { Z ^l(m, n) } two-dimensional position index (m, n) be converted into one dimension index q, namely q=(m-1) N+n, obtains data set { Z ^l(m, n) } adopt the succinct form of one dimension index q have ${\left\{Z_q^l\right\}}_{1\≤q\≤\mathrm{MN}}={\left\{Z^l(m,n)\right\}}_{1\≤m\≤M,1\≤n\≤N};$

3) random initializtion cluster centre and it is right to cluster centre degree of membership relational matrix carry out random initializtion, 1≤q≤MN, { 1,2,3} makes it correspond to high frequency direction sub-band coefficients D to make p ∈ _{l, i}noise belonging to (m, n), weak edge and this three kind of strong edge, namely wherein U _{1, q}u _{2, q}and U _{3, q}correspond respectively to be under the jurisdiction of the probability of noise, weak edge and strong edge classification, and have initialization iteration count variable tem=0;

4) upgrade the value of iteration count variable, have tem=tem+1;

5) utilize with existing degree of membership relational matrix value, upgrade cluster centre value, have:

$V_p^l=\frac{\underset{q=1}{\overset{M\×N}{\mathrm{\Σ}}}{\left(U_{p,q}^l\right)}^2{Z}_q^l}{\underset{q=1}{\overset{M\×N}{\mathrm{\Σ}}}{\left(U_{p,q}^l\right)}^2};$

6) utilize with step 5) calculate new value, upgrade degree of membership relational matrix value, have:

$U_{1,q}^l=\frac{1/d_{1,q}}{\underset{p=1}{\overset{3}{\mathrm{\Σ}}}(1/d_{p,q})},$ $U_{2,q}^l=\frac{1/d_{2,q}}{\underset{p=1}{\overset{3}{\mathrm{\Σ}}}(1/d_{p,q})},$ $U_{3,q}^l=\frac{1/d_{3,q}}{\underset{p=1}{\overset{3}{\mathrm{\Σ}}}(1/d_{p,q})};$

In formula, d _{p, q}for to cluster centre distance, p ∈ { 1,2,3};

7) iteration stopping judge, if tem<50 and with renewal amount be greater than 0.001, go to step 4), otherwise stop iteration, obtain optimum degree of membership relational matrix wherein with correspond respectively to be under the jurisdiction of the most fiducial probability of noise, weak edge and strong edge classification;

8) according to optimum degree of membership relational matrix if the most fiducial probability being under the jurisdiction of noise, weak edge or strong edge classification is maximum, then will be classified as such, and on yardstick l, each high frequency direction sub-band coefficients D of (m, n) position _l,i(m, n) is also all considered to belong to this class, by D _l,i(m, n) divides into noise, weak edge or strong fringing coefficient, and n=mod (q, N), m=(q-n)/N+1, mod represents that q is to N remainder.

Wherein step 1) described in each yardstick of calculating on high frequency direction sub-band coefficients D _{l, i}the three-dimensional feature vector Z of (m, n) ^l(m, n), 1≤m≤M, 1≤n≤N, 0≤l≤L-1, carries out according to the following procedure:

1a) the high frequency direction sub-band coefficients on same yardstick is expressed as M × N × k _ldimension matrix D _l, 0≤l≤L-1, D _lthe element of matrix in (m, n) position is k _ldimensional vector, $D_l(m,n)={[D_{l,1}(m,n),D_{l,2}(m,n),...,D_{l,k_l}(m,n)]}^T;$

2a) from D _lextracting three-dimensional feature in (m, n), is characteristics of mean respectively standard deviation characteristic with maximum norm value tag characteristics of mean computing formula be:

$Z_{\mathrm{mean}}^l(m,n)=\frac{1}{k_l}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}\left|D_{l,i}(m,n)\right|;$

Standard deviation characteristic computing formula be:

$Z_{\mathrm{std}}^l(m,n)=\sqrt{\frac{1}{k_l-1}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}{\left(\right|D_{l,i}(m,n)|-\frac{1}{k_l}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}|D_{l,i}(m,n)\left|\right)}^2};$

Maximum norm value tag computing formula be:

$Z_{\max }^l(m,n)=\max {\left\{\right|D_{l,i}(m,n)\left|\right\}}_{1\&le;i\&le;k_l};$

3a) by step 2a) characteristics of mean that calculates standard deviation characteristic with maximum norm value tag be combined as three-dimensional feature vector Z ^l(m, n), namely $Z^l(m.n)={[Z_{\mathrm{mean}}^l(m,n),Z_{\mathrm{std}}^l(m,n),Z_{\max }^l(m,n)]}^T.$

The inventive method compared with prior art has the following advantages:

1. the present invention is owing to adopting non-downsampling Contourlet conversion, effectively can catch the edge in MRI image, profile information and grain details information.

2. high frequency direction sub-band coefficients is divided into noise, weak edge and strong fringing coefficient due to the method introducing FCM cluster by the present invention adaptively, overcomes the existing dependence chosen threshold value based on the image enchancing method of transform domain.

3. simulation result shows, the present invention, relative to traditional MRI image enchancing method, can obtain integral image sharpness high, and details composition is given prominence to, the enhancing result that homogeneous area is comparatively level and smooth.The method can be used for the details such as internal organs, vein, soft tissue and focus in Contrast-enhanced MRI image, is beneficial to the diagnosis of doctor.

Accompanying drawing explanation

Fig. 1 is the implementation procedure schematic diagram of the inventive method;

Fig. 2 uses the inventive method and existing method to the enhancing Comparative result figure of MRI image im1;

Fig. 3 uses the inventive method and existing method to the enhancing Comparative result figure of MRI image im2;

Fig. 4 uses the inventive method and existing method to the enhancing Comparative result figure of MRI image im3.

Embodiment

The inventive method is on non-downsampling Contourlet territory, the method of FCM cluster is adopted adaptively non-downsampling Contourlet coefficient to be divided into noise, weak edge and strong fringing coefficient, and by noise figure zero setting, strong fringing coefficient is remained unchanged substantially, Nonlinear Mapping strengthens weak fringing coefficient, realizes MRI image enhaucament.

Example of the present invention strengthens true MRI image.

Embodiment 1:

With reference to Fig. 1, concrete steps of the present invention are as follows:

Step 1, input MRI image I (m, n), and L layer non-downsampling Contourlet conversion is carried out to it.

Non-downsampling Contourlet conversion is made up of the turriform filtering of non-lower sampling and the trend pass filtering of non-lower sampling.Carry out one deck non-downsampling Contourlet conversion to MRI image I (m, n), its process is:

1) by size be M × N MRI image I (m, n), 1≤m≤M, 1≤n≤N, wherein M and N is the natural number being greater than 1, the turriform bank of filters of input non-lower sampling, obtains low frequency signal and the bandpass signal of MRI image I (m, n) one deck non-downsampling Contourlet conversion;

2) by the directional filter banks of the bandpass signal of MRI image I (m, n) input non-lower sampling, the high frequency direction sub-band coefficients of MRI image I (m, n) one deck non-downsampling Contourlet conversion is obtained;

3) by MRI image I (m, the low frequency signal of non-downsampling Contourlet conversion n) is as new input picture, repeat above-mentioned steps 1) and 2), obtain the high frequency direction sub-band coefficients D that MRI image I (m, n) carries out on each yardstick of L layer non-downsampling Contourlet conversion _l,i(m, n) and low frequency sub-band coefficient, 0≤l≤L-1,1≤i≤k _l, k _lrepresent at yardstick 2 ^-lon the number of high frequency direction subband, D _{l, i}(m, n) represents at yardstick 2 ^-lon i-th high frequency direction sub-band coefficients, L is 3 ~ 5, and L is natural number.

Step 2, by the high frequency direction sub-band coefficients D on each yardstick _l,i(m, n) carries out cluster respectively, is divided into noise, weak edge and strong fringing coefficient.Its process is:

1) the high frequency direction sub-band coefficients D on each yardstick is calculated _l,ithree-dimensional feature vector Z ' (m, n) of (m, n), 1≤m≤M, 1≤n≤N, 0≤l≤L-1, carries out according to the following procedure:

1a) the high frequency direction sub-band coefficients on same yardstick is expressed as M × N × k _ldimension matrix D _l, 0≤l≤L-1, D _lthe element of matrix in (m, n) position is k _ldimensional vector, $D_l(m,n)={[D_{l,1}(m,n),D_{l,2}(m,n),...,D_{l,k_l}(m,n)]}^T;$

2a) from D _lextracting three-dimensional feature in (m, n), is characteristics of mean respectively standard deviation characteristic with maximum norm value tag characteristics of mean computing formula be:

$Z_{\mathrm{mean}}^l(m,n)=\frac{1}{k_l}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}\left|D_{l,i}(m,n)\right|---\left(1\right)$

Standard deviation characteristic computing formula be:

$Z_{\mathrm{std}}^l(m,n)=\sqrt{\frac{1}{k_l-1}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}{\left(\right|D_{l,i}(m,n)|-\frac{1}{k_l}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}|D_{l,i}(m,n)\left|\right)}^2}---\left(2\right)$

Maximum norm value tag computing formula be:

$Z_{\max }^l(m,n)=\max {\left\{\right|D_{l,i}(m,n)\left|\right\}}_{1\&le;i\&le;k_l}---\left(3\right)$

3a) by step 2a) characteristics of mean that calculates standard deviation characteristic with maximum norm value tag be combined as three-dimensional feature vector Z ^l(m, n), namely $Z^l(m,n)={[Z_{\mathrm{mean}}^l(m,n),Z_{\mathrm{std}}^l(m,n),Z_{\max }^l(m,n)]}^T.$

2) by data set { Z ^l(m, n) } two-dimensional position index (m, n) be converted into one dimension index q, namely q=(m-1) N+n, obtains data set { Z ^l(m, n) } adopt the succinct form of one dimension index q have ${\left\{Z_q^l\right\}}_{1\&le;q\&le;\mathrm{MN}}={\left\{Z^l(m,n)\right\}}_{1\&le;m\&le;M,1\&le;n\&le;N};$

3) random initializtion cluster centre and it is right to cluster centre degree of membership relational matrix carry out random initializtion, 1≤q≤MN, { 1,2,3} makes it correspond to high frequency direction sub-band coefficients D to make p ∈ _{l, i}noise belonging to (m, n), weak edge and this three kind of strong edge, namely wherein U _{1, q}u _{2, q}and U _{3, q}correspond respectively to be under the jurisdiction of the probability of noise, weak edge and strong edge classification;

And have initialization iteration count variable tem=0;

4) upgrade the value of iteration count variable, have tem=tem+1;

5) utilize with existing degree of membership relational matrix value, upgrade cluster centre value, have:

$V_p^l=\frac{\underset{q=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{\left(U_{p,q}^l\right)}^2{Z}_q^l}{\underset{q=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{\left(U_{p,q}^l\right)}^2}---\left(4\right)$

6) utilize with step 5) calculate new value, upgrade degree of membership relational matrix value, have:

$U_{1,q}^l=\frac{1/d_{1,q}}{\underset{p=1}{\overset{3}{\mathrm{\&Sigma;}}}(1/d_{p,q})},$ $U_{2,q}^l=\frac{1/d_{2,q}}{\underset{p=1}{\overset{3}{\mathrm{\&Sigma;}}}(1/d_{p,q})},$ $U_{3,q}^l=\frac{1/d_{3,q}}{\underset{p=1}{\overset{3}{\mathrm{\&Sigma;}}}(1/d_{p,q})}---\left(5\right)$

In formula, d _{p, q}for to cluster centre distance, p ∈ { 1,2,3};

7) iteration stopping judge, if tem<50 and with renewal amount be greater than 0.001, go to step 4), otherwise stop iteration, obtain optimum degree of membership relational matrix wherein with correspond respectively to be under the jurisdiction of the most fiducial probability of noise, weak edge and strong edge classification;

8) according to optimum degree of membership relational matrix if the most fiducial probability being under the jurisdiction of noise, weak edge or strong edge classification is maximum, then will be classified as such, and on yardstick l, each high frequency direction sub-band coefficients D of (m, n) position _{l, i}(m, n) is also all considered to belong to this class, by D _{l, i}(m, n) divides into noise, weak edge or strong fringing coefficient, and n=mod (q, N), m=(q-n)/N+1, mod represents that q is to N remainder.

Step 3, utilizes gain rule formula, to the high frequency direction sub-band coefficients D obtaining category attribute _l,i(m, n) revises.Its process is:

In formula, by setup parameter a ₀=0, to noise figure zero setting, by setup parameter a ₂=1.25, appropriateness is carried out to strong fringing coefficient and strengthens, D ' _{l, i}(m, n) is revised high frequency direction sub-band coefficients, A _lfor at yardstick 2 ^-lon the maximum modulus value of high frequency direction sub-band coefficients, a ₁[] is nonlinear gain function, has:

$a_1[D_{l,i}(m,n)/A_l]=\frac{1}{1+e^{(-t\&CenterDot;D_{l,i}(m,n))/A_l}}-\frac{1}{1+e^{(t\&CenterDot;D_{l,i}(m,n))/A_l}}---\left(7\right)$

In formula, parametric t is for the degree of control MRI image enhaucament, and t value is larger, and MRI image enhaucament degree is larger, and for avoiding MRI image to cross enhancing, arrange t=7, e is natural constant, e=2.71828.

Step 4, do non-downsampling Contourlet inverse transformation to the corrected high frequency direction sub-band coefficients obtained in the low frequency sub-band coefficient of the L layer non-downsampling Contourlet conversion obtained in step 1 and step 3, its process is:

1) successively to high frequency direction sub-band coefficients D ' _l,i(m, n), 0≤l≤L-1,1≤i≤k _l, do the reconstruct of non-lower sampling directional filter banks, obtain MRI and strengthen image I ' (m, n) L, L-1 ..., the bandpass signal of 1 layer of non-downsampling Contourlet conversion;

2) image I ' (m is strengthened to low frequency sub-band coefficient and MRI, n) bandpass signal of L layer does the reconstruct of non-lower sampling turriform bank of filters, obtain the low-pass signal that MRI strengthens image I ' (m, n) L-1 layer non-downsampling Contourlet conversion;

3) image I ' (m is strengthened to MRI, n) low-pass signal of Γ layer non-downsampling Contourlet conversion and MRI strengthen image I ' (m, n) bandpass signal of Γ layer non-downsampling Contourlet conversion does the reconstruct of non-lower sampling turriform bank of filters, obtain MRI and strengthen image I ' (m, n) low-pass signal of Γ-1 layer of non-downsampling Contourlet conversion, makes Γ=L-1 successively, L-2, ..., 1;

The MRI finally obtaining Accurate Reconstruction strengthens image I ' (m, n), and namely MRI strengthens the low-pass signal of image I ' (m, n) the 0th layer of non-downsampling Contourlet conversion.

Below by way of the validity of Simulation experiments validate the inventive method.

Simulated conditions: enhancing emulation experiment is carried out to the true MRI image of 256 × 256 sizes containing little noise.

Emulation content: have selected histogram equalization method HE, unsharp masking method UM, wavelet field Enhancement Method WT and method NSCT of the present invention and contrast.

In experiment, four layers of decomposition are carried out to MRI image.Wavelet transformation adopts ' 9-7 ' wavelet basis.Contourlet transformation and NSCT conversion all adopt ' 9-7 ' Pyramid transform and ' c-d ' directional filter banks, and by thin yardstick to thick yardstick, the decomposition number of high frequency direction subband is respectively 8,8,4 and 4.

In experiment, WT and NSCT all adopts four layers of decomposition to MRI image.The wavelet basis that WT selects is ' 9-7 ' wavelet basis.NSCT adopts classical ' 9-7 ' Pyramid transform and ' c-d ' directional filter banks, and by thin yardstick to thick yardstick, the decomposition number of high frequency direction subband is respectively 8,8,4 and 4.

The evaluation index of the enhancing image that the present invention adopts is:

(1) information entropy.Namely the average information comprised in image.The definition of information entropy is:

$H=-\underset{j=0}{\overset{J-1}{\mathrm{\&Sigma;}}}{p}_j{\log }_{p_j}---\left(8\right)$

In formula, H represents the entropy of image, and J represents the number of greyscale levels that image is total, p _jrepresent that gray-scale value is the pixel count N of j _jthe pixel count N total with image _gratio, that is: p _j=N _j/ N _g.Information entropy reflect strengthen quantity of information contained by image number.The entropy strengthening image is larger, shows that the quantity of information strengthened in image is more.

(2) average.The i.e. gray-scale value average of image pixel.Average reflection be the mean flow rate of image.After image enhaucament, average should increase.The definition of average is:

$\mathrm{mean}=\frac{1}{\mathrm{MN}}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}F(m,n)---\left(9\right)$

In formula, mean represents the gray-scale value average of image pixel, and the size of image F (m, n) is M × N, 1≤m≤M, 1≤n≤N.

(3) standard deviation.Namely the gray value standard of image pixel is poor.What standard deviation reflected is the degree of scatter that image intensity value departs from average.Detailed information in image is abundanter, and the standard deviation of image is also larger.The definition of standard deviation is:

$\mathrm{std}=\sqrt{\frac{1}{\mathrm{MN}}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{[F(m,n)-\mathrm{mean}]}^2}---\left(10\right)$

In formula, std represents that the gray value standard of image pixel is poor.

Simulation result:

(1) according to several image metric indexs of described emulation content simulation as table 1.

Table 1 Comparison of experiment results

From the experimental data of table 1, the average of the enhancing result adopting HE method to obtain and standard deviation are far away higher than additive method, this is that water washing effect in the enhancing result owing to adopting HE method to obtain and noise are stronger, as for MRI image im1, the average of MRI image is 48.81, adopts the be enhanced average of image of HE method to be 133.38, is respectively 50.03 far away higher than UM method, the average of enhancing image that adopts WT method and method NSCT of the present invention to obtain, 50.22,51.14.Although UM method is comparatively responsive to noise, due to original image, to contain noise little, from evaluation index, a small amount of noise does not make a big impact to the enhancing result adopting UM method to obtain, as for MRI image im2, adopt the be enhanced standard deviation of image of UM method to be 32.26, be respectively 34.91 compared to the standard deviation of the enhancing image adopting WT method and method NSCT of the present invention to obtain, 36.97, not greatly different.Method NSCT of the present invention compares with other Enhancement Method in information entropy has advantage, as MRI image iml, adopts HE method, UM method, the information entropy of the enhancing image that WT method and method NSCT of the present invention obtain is respectively 5.09,6.07,6.17,6.19; And the inventive method obtains average and variance is respectively 51.14 and 66.10, all exceed UM method and WT method.

(2) the present invention and existing method is used to the enhancing Comparative result figure of MRI image iml as Fig. 2.Wherein the upper figure of Fig. 2 is MRI image iml; The left figure of Fig. 2 is the enhancing image adopting HE method to obtain; The right figure of Fig. 2 is the enhancing image adopting UM method to obtain; Fig. 2 lower-left figure is the enhancing image adopting WT method to obtain; Fig. 2 bottom-right graph is the enhancing image adopting method NSCT of the present invention to obtain.

(3) the present invention and existing method is used to the enhancing Comparative result figure of MRI image im2 as Fig. 3.Wherein the upper figure of Fig. 3 is MRI image im2; The left figure of Fig. 3 is the enhancing image adopting HE method to obtain; The right figure of Fig. 3 is the enhancing image adopting UM method to obtain; Fig. 3 lower-left figure is the enhancing image adopting WT method to obtain; Fig. 3 bottom-right graph is the enhancing image adopting method NSCT of the present invention to obtain.

(4) the present invention and existing method is used to the enhancing Comparative result figure of MRI image im3 as Fig. 4.Wherein the upper figure of Fig. 4 is MRI image im3; The left figure of Fig. 4 is the enhancing image adopting HE method to obtain; The right figure of Fig. 4 is the enhancing image adopting UM method to obtain; Fig. 4 lower-left figure is the enhancing image adopting WT method to obtain; Fig. 4 bottom-right graph is the enhancing image adopting method NSCT of the present invention to obtain.

Visible with reference to Fig. 2, Fig. 3 and Fig. 4, adopt the enhancing image that HE method obtains: the left figure of left figure and Fig. 4 of Fig. 2 left figure, Fig. 3 creates obvious water washing effect, reduces image readability; Adopt the enhancing image that UM method obtains: the right figure of right figure and Fig. 4 of Fig. 2 right figure, Fig. 3, while enhancing image detail feature, also provide enhanced the noise in image; Adopt the enhancing image that WT method obtains: figure, Fig. 3 lower-left, Fig. 2 lower-left Tu Hetu4 lower-left figure and the enhancing image adopting NSCT method of the present invention to obtain: Fig. 2 bottom-right graph, Fig. 3 bottom-right graph and Fig. 4 bottom-right graph details are clear, comparatively level and smooth at homogeneous area, and compared to MRI image, the overall sharpness strengthening image obtains effective raising.The lower middle portion in abdominal cavity in fig. 2, the right part in abdominal cavity in Fig. 3 and Fig. 4, by contrasting with the enhancing result adopting additive method to obtain, be not difficult to find, the minutia of the enhancing result adopting NSCT method of the present invention to obtain is the abundantest, highlight in MRI image the weak edge details being difficult to observe, and the homogeneous area strengthening image is comparatively level and smooth, the overall sharpness of image is higher.

The present invention is compared to existing traditional images Enhancement Method, no matter from the evaluation of objective parameter, or all there is superiority from the visual quality of image, integral image sharpness can be obtained high, details composition is given prominence to, the enhancing result that homogeneous area is comparatively level and smooth.This enhancing result is conducive to the diagnosis of doctor, is a kind of effective and feasible MRI algorithm for image enhancement.

Claims (2)

1., based on a non-downsampling Contourlet territory MRI image enchancing method for fuzzy C-means clustering, it is characterized in that, comprise the steps:

(1) MRI image I (m, n) that size is M × N is inputted, wherein 1≤m≤M, 1≤n≤N, wherein M and N is the natural number being greater than 1, and carries out L layer non-downsampling Contourlet conversion to it, obtains the high frequency direction sub-band coefficients D on each yardstick _l,i(m, n) and low frequency sub-band coefficient, wherein 0≤l≤L-1,1≤i≤k _l, k _lrepresent at yardstick 2 ^-lon the number of high frequency direction subband, D _l,i(m, n) represents at yardstick 2 ^-lon i-th high frequency direction sub-band coefficients, L is 3 ~ 5, L is natural number;

(2) pass through the high frequency direction sub-band coefficients D on each yardstick _l,i(m, n) carries out cluster respectively, is divided into noise, weak edge and strong fringing coefficient;

Specifically carry out according to following process:

1) the high frequency direction sub-band coefficients D on each yardstick is calculated _l,ithe three-dimensional feature vector Z of (m, n) ^l(m, n), 1≤m≤M, 1≤n≤N, 0≤l≤L-1, definition for characteristics of mean, for standard deviation characteristic, for maximum norm value tag, $Z^l(m,n)={[Z_{\mathrm{mean}}^l(m,n),Z_{\mathrm{std}}^l(m,n),Z_{\max }^l(m,n)]}^T,0\&le;l\&le;L-1;$

2) by data set { Z ^l(m, n) } two-dimensional position index (m, n) be converted into one dimension index q, namely q=(m-1) N+n, obtains data set { Z ^l(m, n) } adopt the succinct form of one dimension index q have ${\left\{Z_q^l\right\}}_{1\&le;q\&le;\mathrm{MN}}={\left\{Z^l(m,n)\right\}}_{1\&le;m\&le;M,1\&le;n\&le;N};$

3) random initializtion cluster centre and it is right to cluster centre degree of membership relational matrix carry out random initializtion, 1≤q≤MN, { 1,2,3} makes it correspond to high frequency direction sub-band coefficients D to make p ∈ _l,inoise belonging to (m, n), weak edge and this three kind of strong edge, namely wherein U _{1, q}u _{2, q}and U _{3, q}correspond respectively to be under the jurisdiction of the probability of noise, weak edge and strong edge classification, and have initialization iteration count variable tem=0;

4) upgrade the value of iteration count variable, have tem=tem+1;

5) utilize with existing degree of membership relational matrix value, upgrade cluster centre value, have:

$V_p^l=\frac{\underset{q=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{\left(U_{p,q}^l\right)}^2{Z}_q^l}{\underset{q=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{\left(U_{p,q}^l\right)}^2};$

6) utilize with step 5) calculate new value, upgrade degree of membership relational matrix value, have:

$U_{1,q}^l=\frac{1/d_{1,q}}{\underset{p=1}{\overset{3}{\mathrm{\&Sigma;}}}(1/d_{p,q})},U_{2,q}^l=\frac{1/d_{2,q}}{\underset{p=1}{\overset{3}{\mathrm{\&Sigma;}}}(1/d_{p,q})},U_{3,q}^l=\frac{1/d_{3,q}}{\underset{p=1}{\overset{3}{\mathrm{\&Sigma;}}}(1/d_{p,q})};$

In formula, d _p,qfor to cluster centre distance,

7) iteration stopping judge, if tem < 50 and with renewal amount be greater than 0.001, go to step 4), otherwise stop iteration, obtain optimum degree of membership relational matrix wherein with correspond respectively to be under the jurisdiction of the most fiducial probability of noise, weak edge and strong edge classification;

8) according to optimum degree of membership relational matrix if the most fiducial probability being under the jurisdiction of noise, weak edge or strong edge classification is maximum, then will be classified as such, and on yardstick l, each high frequency direction sub-band coefficients D of (m, n) position _l,i(m, n) is also all considered to belong to this class, by D _l,i(m, n) divides into noise, weak edge or strong fringing coefficient, and n=mod (q, N), m=(q-n)/N+1, mod represents that q is to N remainder;

(3) gain rule formula is utilized, to the high frequency direction sub-band coefficients D obtaining category attribute _l,i(m, n) revises, that is:

In formula, by setup parameter a ₀=0, to noise figure zero setting, by setup parameter a ₂=1.25, appropriateness is carried out to strong fringing coefficient and strengthens, D ' _l,i(m, n) is revised high frequency direction sub-band coefficients, A _lfor at yardstick 2 ^-lon the maximum modulus value of high frequency direction sub-band coefficients, a ₁[] is nonlinear gain function, has:

$a_1[D_{l,i}(m,n)/A_l]=\frac{1}{1+e^{(-t\&CenterDot;D_{l,i}(m,n))/A_l}}-\frac{1}{1+e^{(t\&CenterDot;D_{l,i}(m,n))/A_l}};$

In formula, parametric t is for the degree of control MRI image enhaucament, and t value is larger, and MRI image enhaucament degree is larger, and for avoiding MRI image to cross enhancing, arrange t=7, e is natural constant, e=2.71828;

(4) non-downsampling Contourlet inverse transformation is done to the corrected high frequency direction sub-band coefficients obtained in the low frequency sub-band coefficient obtained in step (1) and step (3), MRI image I'(m, the n that can be enhanced).

2. a kind of non-downsampling Contourlet territory MRI image enchancing method based on fuzzy C-means clustering as claimed in claim 1, is characterized in that:

Wherein step 1) described in each yardstick of calculating on high frequency direction sub-band coefficients D _l,ithe three-dimensional feature vector Z of (m, n) ^l(m, n), 1≤m≤M, 1≤n≤N, 0≤l≤L-1, carries out according to the following procedure:

1a) the high frequency direction sub-band coefficients on same yardstick is expressed as M × N × k _ldimension matrix D _l, 0≤l≤L-1, D _lthe element of matrix in (m, n) position is k _ldimensional vector, $D_l(m,n)={[D_{l,1}(m,n),D_{l,2}(m,n),...,D_{l,k_l}(m,n)]}^T;$

2a) from D _lextracting three-dimensional feature in (m, n), is characteristics of mean respectively standard deviation characteristic with maximum norm value tag characteristics of mean computing formula be:

$Z_{\mathrm{mean}}^l(m,n)=\frac{1}{k_l}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}\left|D_{l,i}(m,n)\right|;$

Standard deviation characteristic computing formula be:

$Z_{\mathrm{std}}^l(m,n)=\sqrt{\frac{1}{k_l-1}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}{\left(\right|D_{l,i}(m,n)|-\frac{1}{k_l}\underset{i=1}{\overset{k_l}{\mathrm{\&Sigma;}}}|D_{l,i}(m,n)\left|\right)}^2};$

Maximum norm value tag computing formula be:

$Z_{\max }^l(m,n)=\max {\left\{\right|D_{l,i}(m,n)\left|\right\}}_{1\&le;i\&le;k_t};$

3a) by step 2a) characteristics of mean that calculates standard deviation characteristic with maximum norm value tag be combined as three-dimensional feature vector Z ^l(m, n), namely $Z^l(m,n)={[Z_{\mathrm{mean}}^l(m,n),Z_{\mathrm{std}}^l(m,n),Z_{\max }^l(m,n)]}^T.$