CN1924926A - Two-dimensional blur polymer based ultrasonic image division method - Google Patents

Abstract

This invention relates to one two-dimension fuzzy poly spot noise filter and brightness compensation B type of hypersonic cutting method, which comprises the following steps: extending the image brightness information to the update two-dimension poly fuzzy type near the pixel zone based on image brightness information; leading each heter diffusion spot noise filter to the fuzzy poly aim function to provide near information and to strengthening spot robust property; leading two-dimension fuzzy poly aim function with brightness compensation factor based on the noisy on uneven hypothesis and using image aim and background even to strengthen uneven noise robust.

Description

A kind of ultrasonic image division method based on two-dimensional blur polymer

Technical field

The invention belongs to the ultrasonoscopy processing technology field, be specifically related to the ultrasonic image division method of a kind of two-dimensional blur polymer nexus spot noise wave filter and luminance compensation.

Technical background

Image segmentation is meant by certain feature such as brightness or texture, is the process in several zones with image division.In the medical image analysis field, image segmentation is that other is used, as the precondition and the committed step of tissue measurement, anatomical structure analysis and tissue characterization.Although proposed the method for a lot of medical image segmentation, wherein most methods is that signal to noise ratio (S/N ratio) higher relatively nuclear magnetic resonance image or CT image are proposed.Ultrasonic Diagnosis has become one of clinical indispensable diagnostic method because of its advantage such as safe, convenient, painless, inexpensive, but since the existence of much noise in the ultrasonoscopy ultrasonoscopy is cut apart is remained a pendent difficult problem.

In ultrasonic imaging, when the physical dimension of tissue and incident ultrasound waves appearance is near or during less than wavelength, supersonic beam generation scattering.The scatter echo that phase place is different is interfered the generation speckle noise mutually.As the most significant noise in the ultrasonoscopy, the existence of speckle noise has reduced the quality of ultrasonoscopy, covers or has blured in the image diagnosing very important detailed information.Therefore, numerous in recent years scholars have proposed Wiener filtering successively ^[1], adaptive weighted medium filtering ^[2], the denoising of small echo soft-threshold ^[3], anisotropy diffusion ^[4]Deng the method that suppresses speckle noise.Except that speckle noise, the brightness irregularities that is caused by signal dropout in the ultrasonoscopy also is the main cause that signal noise ratio (snr) of image reduces.Brightness irregularities makes same organization internal have the variation of brightness, has also reduced the contrast between different tissues simultaneously, for having made bigger obstacle cutting apart automatically of ultrasonoscopy.Therefore, many brightness irregularities problems that are used for solving medical image were also proposed ^[5,6]Method, but a lot of method remains at nuclear magnetic resonance image and proposes, and for realizing cutting apart automatically of ultrasonoscopy, must suppress the influence to segmentation result of speckle noise and brightness irregularities simultaneously.

Summary of the invention

The objective of the invention is to propose a kind of ultrasonic and brightness irregularities of spot that can suppress simultaneously to cutting apart the ultrasonoscopy automatic division method of influence.

The ultrasonic image division method that the present invention proposes is a kind of based on fuzzy logic and merge that speckle noise suppresses and the ultrasonic image division method of luminance compensation, its concrete steps are as follows: at first, the notion of two-dimentional fuzzy C average (2DFCM) is proposed, when being incorporated into neighborhood information in the cluster process, improve algorithm dirigibility and speed of convergence; Then,, adopt anisotropic diffusion filtering device SRAD to provide neighborhood information, improve the robustness of algorithm speckle noise for 2DFCM at the inhibition of speckle noise; At last, be considered as on the basis of multiplicative noise in the brightness irregularities that will cause,, revise the objective function of 2DFCM according to the characteristics of signal Processing in the ultrasound image acquisition process by signal attenuation, introduce the luminance compensation factor, the two-dimentional fuzzy C-means clustering (2DHFCM) of structure homogeneity.

Below each step is further described in detail.

Related notion: fuzzy C-means clustering (FCM:Fuzzy C-Means) ^[7]

If image pixel value is formed the set X={x of n sample _i| i=1,2 ..., n}, the segmentation problem of image just is converted into the problem that this n sample is divided into c cluster so.If M={m _j| j=1 ..., c} is the center of each cluster, μ _IjBe the membership function of i sample, then with the cluster loss function J of membership function definition to the j class _FCM, can be write as:

$J_{\mathrm{FCM}}=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^b{\left|\right|x_i-m_j\left|\right|}^2---\left(1\right)$

Wherein b＞1 is the constant of a may command cluster result fog-level, and b gets 2 generally speaking.The optimization of image is cut apart, and is exactly to seek cluster centre m by iteration _jWith degree of membership value μ _Ij, make objective function J _FCMGet minimum. $\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}=1$ Restrictive condition under, according to the Lagrange condition extremum method, can obtain objective function J _FCMThe minimizing necessary condition of the condition of getting:

${\mathrm{\μ}}_{\mathrm{ij}}=\frac{{(1/{\left|\right|x_i-m_j\left|\right|}^2)}^{1/(b-1)}}{\underset{k=1}{\overset{c}{\mathrm{\Σ}}}{(1/{\left|\right|x_i-m_k\left|\right|}^2)}^{1/(b-1)}}---\left(2\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b{x}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b}---\left(3\right)$

Traditional F CM algorithm has only utilized the half-tone information of pixel when carrying out image segmentation, so cluster result is comparatively responsive to noise.Adding the neighborhood constraint in the objective function of FCM is the effective means that improves the FCM noiseproof feature ^[8]So-called neighborhood constraint is exactly to calculate current pixel to all kinds of when being subordinate to situation, and the situation that is subordinate to of each pixel that will be adjacent with object pixel is as constraint condition, and it is according to being the high correlation that the natural image neighbor exists.The adding of neighborhood constraint makes cluster result keep the continuity on the image space, has suppressed the influence of noise to cluster result.But the method list of references [8] of structure neighborhood constraint.

1、2DFCM

Suppose that raw image data is expressed as X={x _i| i=1,2 ..., n}, original image is expressed as through the data that obtain behind the wave filter $\tilde{X}=\left\{{\tilde{X}}_i\right|i=\mathrm{1,2},\·\·\·,n\}.$ The segmentation problem of image is exactly will be with two dimensional sample

Be divided into c cluster, each cluster is characterized by a bivector $\stackrel{\‾}{M}=\left\{(m_j,{\tilde{m}}_j)\right|j=1,\·\·\·,c\}.$ Here

For

Cluster centre, j=1 ..., c since the strong correlation between raw image data and its filtering data can use

Classification results influence its respective pixel x _iClassification results, promptly work as x _iTo with m _jWhen having higher degree of membership for the j class at center,

Should to

For the j class at center has similar high degree of membership. Classification to x _iThe influence degree of classification can be by a weighting constant α control.Rule of thumb, α generally gets 6-8.According to above analysis, the objective function of the 2DFCM that the present invention proposes can be written as:

$J_{2\mathrm{DFCM}}=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^b[{(x_i-m_j)}^2+\mathrm{\α}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]---\left(4\right)$

In constraint condition $\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}=1$ Down, finding the solution of objective function can be adopted Lagrange (Lagrange) multiplier method, constructs the objective function of new unconfined condition:

$F=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}}^b[{(x_i-m_j)}^2+\mathrm{\α}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]+\mathrm{\λ}(1-\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\mathrm{\μ}}_{\mathrm{ij}})---\left(5\right)$

λ is the intermediate variable of being introduced by the Lagrange multiplier method.With objective function respectively to μ _Ij, m _jWith

Asking first order derivative and assignment is zero, can obtain the iterative formula of objective function minimalization:

${\mathrm{\μ}}_{\mathrm{ij}}=\frac{{[{(x_i-m_j)}^2+\mathrm{\α}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}{\underset{k=1}{\overset{c}{\mathrm{\Σ}}}{[{(x_k-m_j)}^2+\mathrm{\α}{({\tilde{x}}_k-{\tilde{m}}_j)}^2]}^{-\frac{1}{(b-1)}}}---\left(6\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b{x}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b},{\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b{\tilde{x}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{\mathrm{\μ}}_{\mathrm{ij}}}^b}---\left(7\right)$

The proposition of two-dimensional blur polymer structure has following advantage than traditional F CM: neighborhood information is incorporated in the objective function of 2DFCM, has strengthened the robustness of cluster to noise; On two-dimensional directional, seek cluster centre, potential raising convergence of algorithm speed; Structure is flexible, can adopt different wave filters that the neighborhood constraint information is provided at dissimilar noises.

For suppressing the influence of speckle noise to cluster result, the present invention adopts anisotropic diffusion filtering device SRAD ^[4]2DFCM is provided required neighborhood information.The ultimate principle of anisotropy diffusion is when inhibition edge near zone is level and smooth, strengthens being parallel to the level and smooth of edge direction.In SRAD, filtering is a kind of process of diffusion, and the dynamics of diffusion is determined by following formula:

$c\left(q\right)=\frac{1}{1+[q^2(x,y;t)-q_0^2\left(t\right)]/\left[q_0^2\left(t\right)(1+q_0^2\left(t\right))\right]}---\left(8\right)$

Q (x, y wherein; T) be called the transient change coefficient, it has served as the edge detector in the diffusion process; q ₀(t) be called the noise proportional function, it has served as the diffusion thresholding.The structure of coefficient of diffusion c (q) should be to strengthen smoothly at homogeneous area, and suppresses level and smooth in the direction perpendicular to the edge.Therefore, SRAD effect in actual applications also depends on following two factors: the one, and the validity of rim detection in the noise image.That is: q (x, y; T) can effectively distinguish little transient change coefficient that causes by noise or texture and the big transient change coefficient that causes by the edge; The 2nd, the accuracy that the diffusion threshold value is estimated.That is: q ₀(t) can accurately differentiate and when strengthen or suppress diffusion process.The method that document [4] proposes has solved first problem well, and expands in document [9], directly applies it in the rim detection of ultrasonoscopy.For second problem, still dependence experience is chosen the initial propagations thresholding in the document [4].

For neighborhood constraint information accurately is provided to 2DFCM, the present invention uses the Rayleigh of two-state, Gaussian Mixture distribution that the ultrasonoscopy intensity profile is carried out match, and adopts the EM algorithm to realize the decomposition of mixed distribution; According to the equally distributed zone of speckle noise in the decomposition result predicted picture; By the analysis of homogeneous area statistical property being obtained the diffusion parameter of anisotropy diffusion.The Rayleigh of two-state, Gaussian Mixture distribute as follows:

p(x)＝α ₁p(x|Rayleigh)+α ₂p(x|Gaussian) (9)

α wherein ₁And α ₂Be mixing constant, satisfy α ₁+ α ₂=1, and have:

$p\left(x\right|\mathrm{Rayleigh})=\frac{x}{{\mathrm{\β}}^2}\exp (-\frac{x^2}{{2\mathrm{\β}}^2})---\left(10\right)$

$p\left(x\right|\mathrm{Gaussian})=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp \left[\frac{-{(x-\mathrm{\μ})}^2}{{2\mathrm{\σ}}^2}\right]---\left(11\right)$

But the method list of references [10] that the EM mixed distribution is decomposed.After the EM algorithm convergence, cut apart criterion according to likelihood ^[10], obtain the decision threshold of target area, background area:

$T=\left\{\begin{array}{cc}\sqrt{2[\log ({\mathrm{\&beta;}}^2/{\mathrm{\&sigma;}}^2)+{\mathrm{\&alpha;}}_1/{\mathrm{\&alpha;}}_2]/|1/{\mathrm{\&beta;}}^2-1/{\mathrm{\&sigma;}}^2|}& ({\mathrm{\&alpha;}}_1\&GreaterEqual;{\mathrm{\&alpha;}}_2)\\ \sqrt{2[\log ({\mathrm{\&sigma;}}^2/{\mathrm{\&beta;}}^2)+{\mathrm{\&alpha;}}_2/{\mathrm{\&alpha;}}_1]/|1/{\mathrm{\&sigma;}}^2-1/{\mathrm{\&beta;}}^2|}& ({\mathrm{\&alpha;}}_1<{\mathrm{\&alpha;}}_2)\end{array}\right.---\left(12\right)$

Obtaining behind the decision threshold image segmentation is target area and background area (noise region).Be 30 * 30 fritter then with image division, select the most uniform piece of noise profile, promptly comprise the maximum piece of noise spot, calculate the initial propagations thresholding according to following formula:

$q_0\left(t\right)=\frac{\sqrt{\mathrm{var}\left[z\left(t\right)\right]}}{\stackrel{\&OverBar;}{z\left(t\right)}}---\left(13\right)$

The equally distributed zone of z (t) expression speckle noise wherein, var[z (t)] and z (t) represent variance and the average that this is regional respectively;

2、2DHFCM

The present invention retrains as neighborhood by the diffusion of introducing anisotropy and improves the robustness of cluster to speckle noise in 2DFCM, is the adaptability of enhancement algorithms to actual ultrasonoscopy, also needs to consider to solve the brightness irregularities problem that is caused by signal attenuation.

Brightness irregularities is a common problem that exists in all kinds of medical images, as nuclear magnetic resonance image.Solving the brightness irregularities problem up to now is to be assumed to be a kind of multiplicative noise than successful method ^[5,6], utilize log-compressed to change it into additive noise to simplify follow-up processing procedure then.Suppose X={x ₁, x ₂..., x _n, Y={y ₁, y ₂..., y _nRepresent desirable backward scattering signal data and actual backward scattering signal data respectively, and ideal data refers to not comprise the data of signal attenuation here, and then the relation between ideal data and real data can be expressed as follows:

Y _i＝X _i×B _i，1≤i≤n (14)

Operations such as log-compressed, low-pass filtering, interpolation have been comprised during ultrasound image acquisition, so will observe the relation of ultrasonoscopy grey scale signal and actual backward scattering signal be reduced to the log-compressed relation usually.Log-compressed relation according to gradation of image signal and actual backward scattering signal then can be expressed as follows the relation of ideal image pixel and observed image pixel:

y _i＝x _i+β _i i∈{1，2，...n} (15)

B={ β wherein ₁, β ₂..., β _nBe observation pixel Y={y ₁, y ₂..., y _nLuminance compensation.Because the anisotropic diffusion filtering device has the average retention characteristic to the observation ultrasonoscopy, therefore the luminance compensation to data behind the anisotropic diffusion filtering can adopt the method identical with original image.If through the data after the anisotropy diffusion $\tilde{Y}=\{{\tilde{y}}_1,{\tilde{y}}_2,\&CenterDot;\&CenterDot;\&CenterDot;{\tilde{y}}_N\}$ Has luminance compensation R={ γ ₁, γ ₂..., γ _n, then

Be expressed as follows with the relation of R:

${\tilde{y}}_i={\tilde{x}}_i+{\mathrm{\&gamma;}}_i,\&ForAll;i\&Element;\{\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;n\}---\left(16\right)$

In the objective function of 2DFCM, introduce luminance compensation, be about to the x in the formula (4) _iUse y _i-β _iReplace, With

Replace, just obtain the objective function of 2DHFCM:

$J_{2\mathrm{DHFCM}}=\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]---\left(17\right)$

In constraint condition $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1$ Down, objective function (17) find the solution the method identical that can adopt with 2DFCM, obtain the iterative formula of objective function minimalization:

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{{[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{-\frac{1}{(b-1)}}}{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{[{(y_k-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_k-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}---\left(18\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b(y_i-{\mathrm{\&beta;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b},{\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b({\tilde{y}}_i-{\mathrm{\&gamma;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b}---\left(19\right)$

${\mathrm{\&beta;}}_i=y_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{m}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b},{\mathrm{\&gamma;}}_i={\tilde{y}}_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{\tilde{m}}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b}---\left(20\right)$

Before the 2DHFCM algorithm iteration, the initial parameter that need be predetermined is cluster numbers, initial cluster center and original intensity compensate for estimated.Because 2DHFCM need upgrade three groups of interactional parameters when each parameter iteration, can be caused the wrong iteration of algorithm by error transfer if initial parameter is inaccurate, so the correct convergence of 2DHFCM depends on the accuracy of initial parameter.Here utilize bottom-up lax mark method and Gaussian Mixture Distribution Model to estimate and guarantee the correct convergence of algorithm for 2DHFCM provides initial parameter.Bottom-up lax mark method is to introduce the method for spatial texture information in cluster, this method generally includes two processing procedures: at first predict initial mark probability by a kind of clustering algorithm such as K averaging method, by readjusting classification space continuity is introduced in the final classification results then.Here with 2DFCM image is presorted earlier, then with the original state of its sorting result as 2DHFCM; Then readjust sorting result with 2DHFCM again, reach the purpose of eliminating speckle noise and brightness irregularities simultaneously.Before carrying out 2DFCM, need given cluster numbers and initial cluster center equally, adopt Gaussian Mixture Distribution Model to estimate initial parameter among the present invention.Be about to raw data X and The two-dimensional histogram that constitutes is considered as a plurality of two-dimentional Gaussian distribution and mixes, and the number of mixed distribution is the initial clustering number, and the center position coordinates of mixed distribution is initial cluster center, but implementation procedure list of references [11].

Description of drawings

Fig. 1, to comprising the comparison of speckle noise ultrasonoscopy segmentation result.Wherein Fig. 1 (a) is an ideal image, Fig. 1 (b) is for adding the analog image of speckle noise, Fig. 1 (c) is the filtering result of SRAD, the result that Fig. 1 (d) is cut apart for FCM, the result that Fig. 1 (e) is cut apart for document [8], Fig. 1 (f) is the segmentation result of 2DFCM, and Fig. 1 (g) is the segmentation result of 2DHFCM.

Fig. 2, to comprising the comparison of speckle noise and brightness irregularities ultrasonoscopy segmentation result.Wherein Fig. 2 (a) is for adding the analog image of speckle noise and brightness irregularities, Fig. 2 (b) is the brightness irregularities model, the result that Fig. 2 (c) is cut apart for FCM, the result that Fig. 2 (d) is cut apart for document [8], Fig. 2 (e) is the segmentation result of 2DFCM, Fig. 2 (f) is the segmentation result of 2DHFCM, and Fig. 2 (g) is the luminance compensation figure to original image, and Fig. 2 (h) is the luminance compensation figure to the SRAD filtering image.

Fig. 3, to the comparison of a width of cloth fetal kidney ponding ultrasonoscopy segmentation result.Wherein Fig. 3 (a) is an original image, Fig. 3 (b) is the result of SRAD filtering, and Fig. 3 (c) is the segmentation result of 2DFCM, and Fig. 3 (d) is the segmentation result of 2DHFCM, Fig. 3 (e) is the luminance compensation figure to original image, and Fig. 3 (f) is the luminance compensation figure to the SRAD filtering image.

Fig. 4, to the cluster process explanation of Fig. 3.Fig. 4 (a) is 30 result for cluster numbers, and Fig. 4 (b) is 10 result for cluster numbers, and Fig. 4 (c) is 6 result for cluster numbers, and Fig. 4 (d) is 3 result for cluster numbers, and Fig. 4 (e) is 2 result for cluster numbers.

Fig. 5, to the comparison of a width of cloth fetus arch of aorta ultrasonoscopy segmentation result.Wherein Fig. 5 (a) is an original image, Fig. 5 (b) is the result of SRAD filtering, and Fig. 5 (c) is the segmentation result of 2DFCM, and Fig. 5 (d) is the segmentation result of 2DHFCM, Fig. 5 (e) is the luminance compensation figure to original image, and Fig. 5 (f) is the luminance compensation figure to the SRAD filtering image.

Fig. 6, to the cluster process explanation of Fig. 5.Fig. 6 (a) is 30 result for cluster numbers, and Fig. 6 (b) is 10 result for cluster numbers, and Fig. 6 (c) is 6 result for cluster numbers, and Fig. 6 (d) is 4 result for cluster numbers, and Fig. 6 (e) is 2 result for cluster numbers.

Fig. 7, to the comparison of a width of cloth tumor of breast ultrasonoscopy segmentation result.Wherein Fig. 7 (a) is an original image, and Fig. 7 (b) is the segmentation result of 2DHFCM, and Fig. 7 (c) is the luminance compensation figure to the SRAD filtering image, and Fig. 7 (d) is the result who split image is extracted edge contour.

Embodiment

Below performing step by whole algorithm, further each ingredient in the invention is described

1 pair of original image carries out SRAD filtering

Before filtering, need to adopt the EM algorithm that the probability distribution of image is carried out match, and the thresholding cut apart of computed image target background on this basis, obtain cutting apart figure, the method by windowing finds and comprises the maximum zone of background in the image then, calculates the initial propagations thresholding according to formula (13).Obtain the neighborhood bound data by SRAD again $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}.$

2 adopt Gaussian Mixture Distribution Model to estimate the initial parameter of 2DFCM

At first with according to raw data X={x _i| i=1,2 ..., n} and filtering data $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ Calculate two-dimensional histogram, using variance then is that 0.625 2-d gaussian filters device carries out smoothly histogram, and the main Gaussian distribution number that exists in the smoothed histogram is made as cluster numbers (being made as K), and histogrammic each maximum point is as cluster centre.

3 2DFCM iteration, the objective function of iteration are formula (4), constraint condition $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1.$

Carry out following steps until the change amount of cluster centre before and after twice iteration less than very little constant ε=0.001

Step1) upgrade degree of membership matrix μ with formula (6) _Ij,

Step2) upgrade cluster centre m with formula (7) _j,

Step3) calculate the change amount of upgrading the front and back cluster centre || M _New-M _Old|| ²Here M _OldM before expression is upgraded, M _NewM after expression is upgraded.

4 2DHFCM initialization

With the degree of membership matrix after the 2DFCM convergence and cluster centre initial value as 2DHFCM, and with the initial luminance compensation value of formula (20) calculating.With filtering data $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ Calculate the one dimension histogram, using variance then is that 1.0 one dimension Gaussian filter carries out smoothly histogram, the main Gaussian distribution number that exists in the smoothed histogram is made as the number of categories (being made as c) of expectation.

5 2DHFCM iteration, the objective function of iteration are formula (4), constraint condition $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1.$

The different cluster centre quantity that the execution following steps exist after iteration are less than the number of categories c of expectation.

Step1) upgrade degree of membership matrix μ with formula (18) _Ij

Step2) with formula (19) upgrade cluster centre mj and

Step3) upgrade two luminance compensation vector β with formula (20) _i, γ _i

Step4) if after cluster centre upgraded, the Euclidean distance Xiao Yu  between two cluster centres then merged two cluster centres.Add up the quantity of different cluster centres.

Interpretation is as follows:

To said method of the present invention, carried out emulation testing.The condition of emulation testing is as follows:

Emulation platform is that dominant frequency is the PC of 2.0GHz, and programming tool is Visual c++6.0.Be verification algorithm performance in all cases, adopted a width of cloth only to be subjected to simulation drawing that speckle noise pollutes and existing speckle noise that the simulation drawing of brightness irregularities is arranged again respectively.Produce the simulation software Field_II of the method for simulation drawing for exploitations such as employing Jensen.The validity that simulation drawing is cut apart is weighed with segmentation accuracy, and segmentation accuracy is defined as follows:

$\mathrm{SA}=\frac{N_{\mathrm{CORRECT}}}{N_{\mathrm{TOTAL}}}\&times;100\%---\left(21\right)$

N wherein _CORRECTThe pixel number andN of the correct classification of expression _TOTALAll pixel numbers that presentation video comprises.Fig. 1 (a) is for comprising the ideal image of a circular target, the analog image that polluted by speckle noise of Fig. 1 (b) for producing with Field_II, Fig. 1 (c)-(g) has provided the filtering result of SRAD respectively, and FCM, document [8], the segmentation result of 2DFCM and 2DHFCM.Table 1 has provided FCM respectively, document [8], the SA value of 2DFCM and 2DHFCM segmentation result.Similar by the visible 2DFCM of result with 2DHFCM segmentation result under the situation that only comprises speckle noise, and obviously be better than the segmentation result of FCM and document [8].

Four kinds of methods of table 1 are to the segmentation accuracy contrast of the simulation drawing that only contains speckle noise

	FC M	KFCM_ S2 [8]	FCM2D	HFCM2D
					Fig.2 (c＝2)	92.3 2	96.26	98.88	98.64

The situation that below Fig. 1 (b) is added brightness irregularities is shown in Fig. 2 (a).Fig. 2 (b) is the model of this uneven brightness, and Fig. 2 (c)-(h) has provided FCM respectively, document [8], and the segmentation result of 2DFCM and 2DHFCM, and to original image with to the luminance compensation figure of SRAD filtering image.Table 2 has provided the SA value of four kinds of method segmentation results respectively.Containing algorithm inefficacy under the situation of brightness irregularities by the visible 2DFCM of result, 2DHFCM has obtained higher segmentation accuracy.

Four kinds of methods of table 2 are to not only containing speckle noise but also containing the segmentation accuracy contrast of the simulation drawing of brightness irregularities

	FC M	KFCM_ S2 [8]	FCM2D	HFCM2D
					Fig.2 (c＝2)	66.1 7	69.41	71.79	98.56

Provide algorithm below at actual ultrasonoscopy Application in Segmentation example.Fig. 3 (a) is the fetal kidney ponding image that a width of cloth trans-abdominal ultrasound obtains, and Fig. 3 (b)-(f) has provided the result of SRAD filtering respectively, the segmentation result of 2DFCM and 2DHFCM, and to original image and to the luminance compensation figure of SRAD filtering image.For the cluster process of 2DHFCM is described, Fig. 4 has provided 2DHFCM intermediate result under the different cluster numbers situations in iterative process.Fig. 4 and Fig. 5 have provided has segmentation result than the fetus arch of aorta ultrasonoscopy of complicated shape to cluster process explanation when to a width of cloth.Fig. 6 has provided to piece image the lower galactophore image segmentation result of contrast between second-rate and target background.

By as seen to the segmentation result of simulation ultrasonoscopy and actual ultrasonoscopy, the present invention can suppress in the ultrasonoscopy speckle noise and brightness irregularities simultaneously to the influence of segmentation result, and not needing in cutting apart artificially to adjust empirical parameter, is that a kind of effective ultrasonoscopy is cut apart new method.

List of references

[1]R.C.Gonzalez，R.E.Woods.Digital image processing[M].Beijing：Publishing House ofElectronics Industry，2003.209.

[2]M.Karaman，M.A.Kutay，G.Bozdagi.An adaptive speckle suppression filter for medicalultrasonic imaging[J].IEEE Trans.Medical Image，1995，14(2)：283-292.

[3]L.D.David.De-noising by soft-thresholding[J].IEEE Trans.Information Theory，1995，41(3)：613-627.

[4]YU Yong-jian and T.A.Scott.Speckle reducing anisotropic diffusion[J].IEEE Trans.ImageProcessing，2002，11(11)：1260-1270.

[5]D.L.Pham，J.L.Prince，“Adaptive fuzzy segmentation of magnetic resonance images，”IEEE Trans Med Imag.，vol.18，no.9，pp.737-752，Sept.1999.

[6]G.F.Xiao，M.Brady，J.A.Noble，and Yongyue Zhang，“Segmentation of ultrasound B-modeimages with intensity inhomogeneity correction，”IEEE Trans.Med Imag.，vol.21，no.1，pp.48-57，Jan.2002

[7]J.C.Bezdek，“A convergence theorem for the fuzzy ISODATA clustering algorithm，”IEEETrans.Pattern Anal.Machine Intell.，vol.PAMI-2，pp.1-8，1980.

[8]Chen Song-can，Zhang Dao-qiang，“Robust Image Segmentation Using FCM With SpatialConstraints Based on New Kernel-Induced Distance Measure，”IEEE Trans.Systems，Man.Cybern.，vol.34，no.4，pp.1907-1916，Apr.2004.

[9]YU Yong-jian and T.A.Scott.Edge detection in ultrasound imagery using the instantaneouscoefficient of variation[J].IEEE Trans.Image Processing，2004，13(12)：1640-1655.

[10]M.S.Hassouna，A.A.Farag，S.Hushek and T.Moriarty，“Cerebrovascular segmentationfrom TOF using stochastic models，”Medical Image Analysis.，vol.10，pp.2-18，2006

[11]F.Kurugollu，B.Sankur，A.E.Harmanci，“Color image segmentation using histogrammultithresholding and fusion，”Image and Vision Computing.，vol.19，pp.915-928，2001.

Claims (2)

1, a kind of ultrasonic image division method based on fuzzy clustering is characterized in that concrete steps are: at first, propose the notion of 2DFCM, neighborhood information is incorporated in the cluster process, to improve algorithm dirigibility and speed of convergence; Then,, adopt anisotropic diffusion filtering device SRAD to provide neighborhood information, improve the robustness of algorithm speckle noise for 2DFCM at the inhibition of speckle noise; At last, will be considered as multiplicative noise by the brightness irregularities that signal attenuation causes, according to the characteristics of signal Processing in the ultrasound image acquisition process, revise the objective function of 2DFCM, introduce the luminance compensation factor, the two-dimentional fuzzy C-means clustering of structure homogeneity, 2DHFCM is two-dimentional fuzzy C average here.

2, the ultrasonic image division method based on fuzzy clustering according to claim 1 is characterized in that concrete steps are:

(1) original image is carried out SRAD filtering

Before filtering, need to adopt the EM algorithm that the probability distribution of image is carried out match, and the thresholding cut apart of computed image target background on this basis, obtain cutting apart figure, method by windowing finds and comprises the maximum zone of background in the image then, calculate the initial propagations thresholding according to formula (13), obtain the neighborhood bound data by SRAD again $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\};$

(2) adopt Gaussian Mixture Distribution Model to estimate the initial parameter of 2DFCM

At first according to raw data X={x _i| i=1,2 ..., n} and filtering data $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ Calculate two-dimensional histogram, using variance then is that 0.625 2-d gaussian filters device carries out smoothly histogram, and the main Gaussian distribution number that exists in the smoothed histogram is made as cluster numbers K, and histogrammic each maximum point is as cluster centre;

(3) 2DFCM iteration, the objective function of iteration are formula (4), and constraint condition is $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1,$

Carry out following steps until the change amount of cluster centre before and after twice iteration less than very little constant ε=0.001,

(1) upgrades degree of membership matrix, μ with formula (6) _Ij,

(2) upgrade cluster centre m with formula (7) _j,

(3) calculate the change amount ‖ M that upgrades the front and back cluster centre _New-M _Old‖ ², M here _OldM before expression is upgraded, M _NewM after expression is upgraded;

(4) 2DHFCM initialization, 2DHFCM is the two-dimentional fuzzy C-means clustering of homogeneity here,

With the degree of membership matrix after the 2DFCM convergence and cluster centre initial value as 2DHFCM, and with the initial luminance compensation value of formula (20) calculating; With filtering data $\tilde{X}=\left\{{\tilde{x}}_i\right|i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n\}$ Calculate the one dimension histogram, using variance then is that 1.0 one dimension Gaussian filter carries out smoothly the main Gaussian distribution number that exists in the smoothed histogram being made as the number of categories c of expectation to histogram;

(5) 2DHFCM iteration, the objective function of iteration are formula (17), and constraint condition is $\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}=1,$

Carry out different cluster centre quantity that following steps exist number of categories c less than expectation after iteration,

(1) upgrades degree of membership matrix μ with formula (18) _Ij,

(2) upgrade cluster centre m with formula (19) _j,

(3) upgrade two luminance compensation vectors, β with formula (20) _i, γ _i,

(4) if after cluster centre upgrades, the Euclidean distance Xiao Yu  between two cluster centres then merges two cluster centres the quantity of adding up different cluster centres;

$J_{2\mathrm{DFCM}}=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b[{(x_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{x}}_i-{\tilde{m}}_j)}^2]---\left(4\right)$

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{{[{(x_i-m_j)}^2+{\mathrm{\&alpha;}({\tilde{x}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{[{(x_k-m_j)}^2+\mathrm{\&alpha;}{({\tilde{x}}_k-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}---\left(6\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b{x}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b},$ ${\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b{\tilde{x}}_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b}---\left(7\right)$

$q_0\left(t\right)=\frac{\sqrt{\mathrm{var}\left[z\left(t\right)\right]}}{\stackrel{\&OverBar;}{z\left(t\right)}}---\left(13\right)$

$J_{2\mathrm{DHFCM}}=\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]---\left(17\right)$

${\mathrm{\&mu;}}_{\mathrm{ij}}=\frac{{[{(y_i-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_i-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{[{(y_k-{\mathrm{\&beta;}}_i-m_j)}^2+\mathrm{\&alpha;}{({\tilde{y}}_k-{\mathrm{\&gamma;}}_i-{\tilde{m}}_j)}^2]}^{\frac{1}{(b-1)}}}---\left(18\right)$

$m_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b(y_i-{\mathrm{\&beta;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b},$ ${\tilde{m}}_j=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b({\tilde{y}}_i-{\mathrm{\&gamma;}}_i)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{{\mathrm{\&mu;}}_{\mathrm{ij}}}^b}---\left(19\right)$

${\mathrm{\&beta;}}_i=y_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{m}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b},$ ${\mathrm{\&gamma;}}_i={\tilde{y}}_i-\frac{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b{\tilde{m}}_j}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{\mathrm{ij}}^b}---\left(20\right)$

Here, b=2, α gets 6-8; { β ₁, β ₂..., β _n}=B is to observation pixel Y={y ₁, y ₂..., y _nLuminance compensation, $\tilde{Y}=\{{\tilde{y}}_1,{\tilde{y}}_2,\&CenterDot;\&CenterDot;\&CenterDot;{\tilde{y}}_N\}$ Be the data after the anisotropy diffusion, { γ ₁, γ ₂..., γ _n}=R is

Luminance compensation, the equally distributed zone of z (t) expression speckle noise wherein, var[z (t)] and z (t) represent variance and the average that this is regional respectively.

Images