CN103606165B - Phase locked medical image processing method Network Based - Google Patents

Abstract

The invention discloses a kind of phase locked medical image processing method Network Based, mainly solve the problem that existing medical image processing method diagnosis discrimination is not high and local detail process is bad.Implementation step is: (1) reads in primitive medicine image, extracts the feature of pixel, carries out watershed segmentation；(2) similarity between segmentation rear region is calculated；(3) with region as node, similarity, as network weight, sets up network；(4) randomly generating the initial phase of network, being iterated phase place, until stablizing, extracting the node phase after stablizing；(5) by the phase place after stablizing, node is carried out just classification then computing node slope, according to slope, node is finely divided class；(6) redistribute by the node of misclassification, and classification sequence number is tagged on pixel；(7) different gray values, the image after being split are composed to inhomogeneous pixel.The present invention compared with the conventional method, improves the diagnosis discrimination of medical image.

Description

Phase locked medical image processing method Network Based

Technical field

The invention belongs to image processing field, relate to the process of medical image, can be used for monitoring of diseases distribution, research pathogeny and disease auxiliary diagnosis.

Background technology

Medical image includes CT, positron emission chromatography imaging technique PET, single photon emission tomographic SPECT, NMR (Nuclear Magnetic Resonance)-imaging MR and the image that other medical imaging device obtains.Complexity and multiformity due to medical image, and the property difference of the image-forming principle of medical image and tissue itself, and the formation of image is subject to the impact of such as noise, field offset effect, local volume effect and histokinesis etc., medical image is compared with normal image, has the features such as fuzzy, inhomogeneities.It addition, human dissection organizational structure and complex-shaped, the difference between Different Individual is big, and the problem also making medical image processing is increasingly complex.

In recent years, along with medical imaging is in clinical medical successful Application, the emerging technology of some medical image processing, such as fuzzy mathematics, mathematical morphology, digital topology, artificial intelligence etc., make remarkable progress, new Medical Imaging Technology emerges in an endless stream.Wherein, calculating speed is fast, segmentation is accurate and is accurately positioned the characteristics such as image border because having for the watershed algorithm based on mathematical morphology that Digabel et al. proposes, and receives significant attention.The ultimate principle of watershed algorithm is the topological landforms that image is regarded as in geodesy, and gray level image is considered as landform appearance, the grey scale pixel value height above sea level to putting.Each local minizing point, and its influence area is referred to as reception basin, the border of reception basin is then referred to as watershed.Watershed algorithm is by finding watershed that image is split.But, this method has the drawback that there is over-segmentation phenomenon.

Summary of the invention

Present invention aims to the problem that medical image is fuzzy, uneven, individual variation is big, it is proposed to a kind of phase locked medical image processing method Network Based, to solve the over-segmentation problem that existing method exists, improve the diagnosis discrimination of medical image.

The technical thought of the present invention is: with dividing ridge method, the image after feature extraction is split, and the zonule after segmentation is corresponded to the node in network, sets up network model；By agitator formula, node is carried out phase place iteration, complete the merging of overdivided region, improve the effect of image segmentation.Implementation step includes as follows:

(1) extract the pixel feature of primitive medicine image, and the pixel feature extracted is carried out watershed segmentation, obtain N number of region, N >=1000；

(2) similarity between calculating each two region, and each region and the similarity of himself, obtain N × (N+1)/2 Similarity value；

(3) with N number of region as each node in network, with N × (N+1)/2 Similarity value as the weights in network；

(4) in [0,2 π], randomly generate N number of number, as the initial phase of each node in network, adopt agitator formula that the N number of node in network is carried out phase place iteration, when all node N meet formula θ_i(k)-θ_i(k-3) during < 0.01, finishing iteration, θ_iK () represents the phase place of node i kth time iteration, iteration total degree is designated as m；

(5) [0,2 π] is divided into n₁Individual subinterval, statistics phase place falls the node number in each subinterval, with n₁Individual subinterval sequence number is as abscissa, to drop on subinterval interior nodes number as vertical coordinate, draws PHASE DISTRIBUTION curve, and statistics is between adjacent trough all nodes on curve, is classified as a class, and N number of node is divided into L altogether₁Class, L₁Numerical value identical with trough number, the classification sequence number of N number of node is stored in the vectorial F that length is N, if the arbitrary node i in N number of node belongs to jth class, then the i-th value f in F_i=j, j=1 ..., L₁；

(6) with m generation of node and m-3 for the slope of all nodes of phase calculation, N number of slope value is obtained, to L₁Each class of apoplexy due to endogenous wind carries out following operation: [-0.01,0.01] is divided into n₂Individual subinterval, statistics slope value falls the node number in each subinterval, with n₂Individual subinterval sequence number, as abscissa, to drop on subinterval node number as vertical coordinate, draws L₁Bar slope distribution curve, adds up being between adjacent trough all nodes on each curve, is classified as a class, and all nodes are divided into L altogether₂Class, L₂With L₁The trough sum of bar slope distribution curve is identical, updates the classification sequence number of each node in vector F；

(7) detection L₂The nodes of each class of individual apoplexy due to endogenous wind, if the nodes of certain class is less than threshold value T=N/200, then redistributes such all nodes, the similarity of the arbitrary node i of this apoplexy due to endogenous wind Yu class j is designated as v_ij, calculate the L of i and all classes₂Individual similarity, to L₂Individual similarity is ranked up, maximizing, and node i is given class that maximum is corresponding again, updates the classification sequence number of vector F interior joint i；

(8) by the classification sequence number of vector F storage, namely the classification sequence number belonging to the region of pixel, is tagged on each pixel；

(9) compose different gray values, the image after being split to inhomogeneous pixel, and export.

The present invention compared with prior art has the advantage that

1, the present invention adopts agitator formula that the region after watershed algorithm segmentation is carried out phase place iteration, and operand is little, calculates speed fast；

2, the region after segmentation is merged by the present invention by extracting phase place and slope, compensate for the over-segmentation defect of watershed algorithm, improves the diagnosis discrimination of medical image；

3, the present invention fully takes into account the complexity of medical image, adopts the method for sequencing of similarity that the region sorting out mistake is reclassified, it is possible to focal area is accurately positioned.

Simulation result shows, mammogram figure adopts the medical image processing method of the present invention, improves the discrimination of image, contribute to the Accurate Diagnosis of medical image.

Accompanying drawing explanation

Fig. 1 is the flow chart of the present invention；

Fig. 2 is the original mammogram figure that the present invention emulates employing；

Fig. 3 is the mammogram figure after Fig. 2 is carried out watershed segmentation by the present invention；

Fig. 3 is carried out the mammogram figure after region merging technique by the present invention and existing two kinds of methods by Fig. 4.

Detailed description of the invention

With reference to Fig. 1, the present invention phase locked medical image processing method Network Based, comprise the steps:

Step 1, extracts the feature of primitive medicine image pixel point and pixel feature is carried out watershed segmentation.

(1a) inputting primitive medicine image, this example selects original mammogram figure as shown in Figure 2, is sized to 1024 × 1024；

(1b) feature based on gray level co-occurrence matrixes is extracted, for any pixel i, extract contrast, concordance, energy [0,45,90,135] 12 dimensional features on four direction, then extract 10 dimensional features based on un-decimated wavelet transform conversion, above-mentioned 12 dimensional features and this two feature of 10 dimensional features are merged composition 22 dimensional vectors, as the feature of ith pixel point；

(1c) all pixels in image are repeated step (1b), obtain the feature of all pixels of primitive medicine image；

(1d) the pixel feature extracted is carried out watershed segmentation, obtain N number of region, N >=1000, the mammogram figure after segmentation, as shown in Figure 3.

Step 2, the similarity between zoning.

(2a) calculate the feature in all regions, obtain the feature x in N number of region₁,...x_N, wherein, the feature x of ith zone_iThe meansigma methods of all pixel features equal in ith zone；

(2b) with the feature x in N number of region₁,...x_N, calculate the similarity between each two region, and each region and himself similarity, the computing formula of similarity is:

$w_{\mathrm{ij}}=e^{-\frac{{(x_i-x_j)}^2}{{2\mathrm{\&sigma;}}^2}},i=1,...N,j=1,...N,$

Wherein, w_ijRepresent the similarity between ith zone and jth region, x_iIt is the feature of ith zone, x_jBeing the feature in jth region, σ is scale parameter, and span is (0,1).

Step 3, with N number of region as each node in network, with N × (N+1)/2 Similarity value as the weights in network, sets up network.

Step 4, adopts oscillator phase synchronization formula to carry out phase place iteration.

(4a) arithmetic mean of instantaneous value of all weights N × (N+1)/2 in computing networkGenerator matrix E, the value e of the i-th row jth row in this matrix E_ijFor:

$\left\{\begin{array}{c}{e}_{\mathrm{ij}}=1,w_{\mathrm{ij}}\&GreaterEqual;\stackrel{\&OverBar;}{w}\\ e_{\mathrm{ij}}=0,w_{\mathrm{ij}}<\stackrel{\&OverBar;}{w}\end{array}\right.,i=1,...N,j=1,...N;$

(4b) in [0,2 π], produce N number of random number, as the initial phase of N number of node, be designated as θ (0), θ (0)={ θ₁(0),...θ_i(0),...,θ_N(0) }, wherein, θ_i(0) it is the initial phase of i-th node；

(4c) the initial phase θ according to i-th node_i(0), obtaining oscillator phase synchronization formula is:

${\mathrm{\&theta;}}_i(k+1)={\mathrm{\&theta;}}_i\left(k\right)+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Gamma;}}{{\mathrm{\&Gamma;}}_{\max }}\sin ({\mathrm{\&theta;}}_j\left(k\right)-{\mathrm{\&theta;}}_i\left(k\right)),i=1,...N,$

Wherein, N is the nodes of network, and k is iterations, θ_iK () is the phase place of node i kth time iteration, Γ is stiffness of coupling, works as e_ijWhen=1, Γ=3, work as e_ijWhen=0, Γ=-1,

(4d) from k=0, the initial phase θ (0) of N number of node is substituted in above-mentioned oscillator phase synchronization formula, N number of node is carried out phase place iteration；

(4e) as iterations k > 3, it is judged that whether node reaches stable, if all N number of nodes all meet formula θ_i(k)-θ_i(k-3) < 0.01, then node phase reaches stable, and finishing iteration, iteration total degree is designated as m；

(4f) when iterations k≤3, or also not all node all meets formula θ_i(k)-θ_i(k-3) < 0.01, then the kth of N number of node for phase theta₁(k),...θ_NK () substitutes in above-mentioned oscillator phase synchronization iterative formula, calculate kth+1 generation phase place of N number of node, and add 1 by iterations k, returns step (4e).

Step 5, classifies to node according to phase place.

[0,2 π] is divided into n₁Individual subinterval, statistics phase place falls the node number in each subinterval, with n₁Individual subinterval sequence number, as abscissa, to drop on subinterval interior nodes number as vertical coordinate, draws PHASE DISTRIBUTION curve, and statistics is in all nodes between adjacent trough on curve, is classified as a class, and N number of node is divided into L altogether₁Class, L₁Numerical value identical with trough number, the classification sequence number of N number of node is stored in the vectorial F that length is N, if the arbitrary node i in N number of node belongs to jth class, then the i-th value f in F_i=j, j=1 ... L₁。

Step 6, classifies to node according to slope.

(6a) by slope S for all node N of phase calculation of m generation of node and m-3:

S={s₁,...s_i,...s_N,

Wherein, the slope s of i-th node_iFor:

$s_i=\frac{{\mathrm{\&theta;}}_i\left(m\right)-{\mathrm{\&theta;}}_i(m-3)}{3},i=1,...N;$

(6b) to L₁Each class of apoplexy due to endogenous wind carries out following operation: [-0.01,0.01] is divided into n₂Individual subinterval, statistics slope value falls the node number in each subinterval, with n₂Individual subinterval sequence number, as abscissa, to drop on subinterval node number as vertical coordinate, draws L₁Bar slope distribution curve, adds up being between adjacent trough all nodes on each curve, is classified as a class, and all nodes are divided into L altogether₂Class, L₂With L₁The trough sum of bar slope distribution curve is identical, updates the classification sequence number of each node in vector F.

Step 7, redistributes the node of possible misclassification.

(7a) L is detected₂The nodes of individual apoplexy due to endogenous wind any type, if such nodes is less than threshold value T=N/200, then calculates the L of this apoplexy due to endogenous wind i-th node and all classes₂Individual similarity v_ij:

$v_{\mathrm{ij}}=\underset{n=1,f_n=j}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{\mathrm{in}},j=1,...L_2;$

(7b) to L₂Individual similarity is ranked up, maximizing, i-th node is given class that maximum is corresponding again, and updates the node classification sequence number in vector F；

(7c) to all L₂Individual class repeats step (7a)-(7b), is reclassified by all nodes class interior joint less than threshold value T.

Step 8, by the classification sequence number of vector F storage, namely the classification sequence number belonging to the region of pixel, is tagged on each pixel.

Step 9, composes different gray values, the image after being split to inhomogeneous pixel, and exports.

(9a) the gray value h of the i-th class is calculated_i:

$h_i=255\&times;\frac{L_2-i+1}{L_2},i=1,...L_2,$

By h_iIt is assigned to all pixels of the i-th class；

(9b) make i=i+1, repeat step (9a), until i=L₂, now L₂All pixels all assignment of individual class are complete, and the gray level image formed by all pixels after assignment is as final segmentation image；

(9c) final segmentation image is exported.

The effect of the present invention can be further illustrated by the following emulation experiment to mammogram figure:

1, simulated conditions

The emulation of the present invention is at windowsXP, SPI, CPUInterCore2Duo, fundamental frequency 2.33Ghz, and software platform is that Matlab2007b runs.The original mammogram that emulation is selected derives from common data sets MIAS, obtains 330 original mammograms altogether, and wherein Fig. 2 is a width canceration mammogram therein.

2, emulation content and result

Emulation 1, carries out watershed segmentation by the present invention to canceration mammogram shown in Fig. 2, obtains the mammogram after Fig. 3 segmentation.

Emulation 2, by the present invention and existing Kmeas method and existing Fcm method, the image after Fig. 3 is split carries out region classification, and result is Fig. 4 such as.Wherein Fig. 4 a is the result adopting the network phase synchronous method of the present invention that Fig. 3 carries out region classification；Fig. 4 b is the result adopting existing Kmeas method that Fig. 3 carries out region classification；Fig. 4 c is the result adopting existing Fcm method that Fig. 3 carries out region classification.

From Fig. 4 a it can be seen that the present invention can detect the lesion region of mammary gland exactly, and compared with Fig. 4 b and Fig. 4 c, the lesion region that the present invention detects is without assorted point, and region consistency is higher, and in the local detail process to tumor region, advantage is comparatively obvious.

Above-mentioned simulation result shows, the present invention can improve the order of accuarcy of image segmentation compared to traditional method, it is thus achieved that better segmentation result, is a kind of effective medical image processing method.

Above example is carried out under premised on technical solution of the present invention, gives detailed embodiment and concrete operating process, but protection scope of the present invention is not limited to above-described embodiment.

Claims (5)

1. a phase locked medical image processing method Network Based, comprises the steps:

(1) extract the pixel feature of primitive medicine image, and the pixel feature extracted is carried out watershed segmentation, obtain N number of region, N >=1000；

(2) similarity between calculating each two region, and each region and the similarity of himself, obtain N × (N+1)/2 Similarity value；

(3) with N number of region as each node in network, with N × (N+1)/2 Similarity value as the weights in network；

(4) in [0,2 π], randomly generate N number of number, as the initial phase of each node in network, adopt agitator formula that the N number of node in network is carried out phase place iteration:

(4a) the arithmetic mean of instantaneous value w of all weights N × (N+1)/2 in computing network, generator matrix E, the value e of the i-th row jth row in this matrix E_ijFor:

$\left\{\begin{array}{c}{e}_{ij}=1,w_{ij}\&GreaterEqual;\stackrel{\&OverBar;}{w}\\ e_{ij}=0,w_{ij}<\stackrel{\&OverBar;}{w}\end{array}\right.,i=1,...N,j=1,...N;$

(4b) in [0,2 π], produce N number of random number, as the initial phase of N number of node, be designated as θ (0), θ (0)={ θ₁(0),...θ_i(0),...,θ_N(0) }, wherein, θ_i(0) it is the initial phase of i-th node；

(4c) the initial phase θ according to i-th node_i(0), obtaining oscillator phase synchronization iterative formula is:

${\mathrm{\&theta;}}_i(k+1)={\mathrm{\&theta;}}_i\left(k\right)+\underset{j=1}{\overset{N}{\&Sigma;}}\frac{\mathrm{\&Gamma;}}{{\mathrm{\&Gamma;}}_{max}}sin\left({\mathrm{\&theta;}}_j\right(k)-{\mathrm{\&theta;}}_i(k\left)\right),i=1,...N,$

Wherein, N is the nodes of network, and k is iterations, θ_iK () is the phase place of the kth time iteration of i-th node, Γ is stiffness of coupling, works as e_ijWhen=1, Γ=3, work as e_ijWhen=0, Γ=-1,

(4d) from k=0, the initial phase θ (0) of N number of node is substituted in above-mentioned oscillator phase synchronization iterative formula, N number of node is carried out phase place iteration；

When all node N meet formula θ_i(k)-θ_i(k-3) < when 0.01, finishing iteration, θ_iK () represents the phase place of node i kth time iteration, iteration total degree is designated as m；

(5) [0,2 π] is divided into n₁Individual subinterval, statistics phase place falls the node number in each subinterval, with n₁Individual subinterval sequence number is as abscissa, to drop on subinterval interior nodes number as vertical coordinate, draws PHASE DISTRIBUTION curve, and statistics is between adjacent trough all nodes on curve, is classified as a class, and N number of node is divided into L altogether₁Class, L₁Numerical value identical with trough number, the classification sequence number of N number of node is stored in the vectorial F that length is N, if the arbitrary node i in N number of node belongs to jth class, then the i-th value f in F_i=j, j=1 ..., L₁；

(6) with m generation of node and m-3 for the slope of all nodes of phase calculation, N number of slope value is obtained, to L₁Each class of apoplexy due to endogenous wind carries out following operation: [-0.01,0.01] is divided into n₂Individual subinterval, statistics slope value falls the node number in each subinterval, with n₂Individual subinterval sequence number, as abscissa, to drop on subinterval node number as vertical coordinate, draws L₁Bar slope distribution curve, adds up being between adjacent trough all nodes on each curve, is classified as a class, and all nodes are divided into L altogether₂Class, L₂With L₁The trough sum of bar slope distribution curve is identical, updates the classification sequence number of each node in vector F；

(7) detection L₂The nodes of each class of individual apoplexy due to endogenous wind, if the nodes of certain class is less than threshold value T=N/200, then redistributes such all nodes, the similarity of the arbitrary node i of this apoplexy due to endogenous wind Yu class j is designated as v_ij, calculate the L of i and all classes₂Individual similarity, to L₂Individual similarity is ranked up, maximizing, and node i is given class that maximum is corresponding again, updates the classification sequence number of vector F interior joint i；

(8) by the classification sequence number of vector F storage, namely the classification sequence number belonging to the region of pixel, is tagged on each pixel；

(9) compose different gray values, the image after being split to inhomogeneous pixel, and export.

2. the Network Based phase locked medical image processing method according to claims 1, wherein the pixel feature extracting primitive medicine image described in step (1), carry out according to following steps:

(2a) extract based on the feature of gray level co-occurrence matrixes, for any pixel i, extract contrast, concordance, the energy 12 dimensional features on [0 °, 45 °, 90 °, 135 °] four direction；

(2b) for any pixel i, 10 dimensional features based on un-decimated wavelet transform conversion are extracted；

(2c) above-mentioned 12 dimensional features and this two feature of 10 dimensional features are merged composition 22 dimensional vectors, as the feature of ith pixel point；

(2d) all pixels in image are repeated step (2a)-(2c), obtain the feature of all pixels of primitive medicine image.

3. the Network Based phase locked medical image processing method according to claims 1, the wherein similarity between the zoning described in step (2), carries out as follows:

(3a) calculate the feature in all regions, obtain the feature x in N number of region₁,…x_N, wherein, the feature x of ith zone_iThe meansigma methods of all pixel features equal in ith zone；

(3b) with the feature x in N number of region₁,…x_N, calculate the similarity between each two region, and each region and himself similarity, the computing formula of similarity is:

$w_{ij}=e^{-\frac{{(x_i-x_j)}^2}{2{\mathrm{\&sigma;}}^2}},i=1,...N,j=1,...N,$

Wherein, w_ijRepresent the similarity between ith zone and jth region, x_iIt is the feature of ith zone, x_jBeing the feature in jth region, σ is scale parameter, and span is (0,1).

4. the Network Based phase locked medical image processing method according to claims 1, the wherein slope calculating all nodes described in step (6), carry out according to equation below:

$s_i=\frac{{\mathrm{\&theta;}}_i\left(m\right)-{\mathrm{\&theta;}}_i(m-3)}{3},i=1,...N,$

Wherein, s_iIt it is the slope of i-th node.

5. the Network Based phase locked medical image processing method according to claims 1, the wherein L of the computing node i described in step (7) and all classes₂Individual similarity, carries out according to equation below:

$v_{ij}=\underset{n=1,f_n=j}{\overset{N}{\&Sigma;}}{w}_{in},j=1,...L_2,$

Wherein, v_ijIt it is the similarity of node i and class j.