CN101013503A - Method for segmenting abdominal organ in medical image - Google Patents

Abstract

The invention is in the abdominal organs segmentation technology in medical image. Its characteristic is that the method followed with the following steps: through the regular mutual information method, do overall map between the standard images to individual image obtained by CT or MRI scans; separate match of each interested organ using similarity measure; fuzzy connecting segmentation and organ shape amendment; the fuzzy connecting segmentation is divided into: first calculate the gray histogram curve, and treat them as the probability density function of target image; use the point with he greatest probability density region as seed, the histogram curve and seeds as the initial parameters, to get the fuzzy image by fuzzy connecting algorithm, and then separate the target and background pixels by the set threshold. The invention can implement multiple abdominal organs automatic segmentation in objectives image, and can be implemented simply.

Description

A kind of medical image midriff organ segmentation method

Technical field

The invention belongs to Medical Image Processing and application, relate to a kind of medical image midriff organ segmentation method.

Background technology

Cutting apart of abdomen organ has important research meaning and clinical value.At first, from various types of medical images organ being split is to carry out the visual first step.More prior than visual, determine that organ of interest position and zone are significant in the radiotherapy surgical planning.Have only the spatial geometric shape of having grasped organ, just can make radiotherapy scheme accurately and radiotherapy dosage.

People such as Linda G. Shapiro propose to utilize the priori of organ shape dynamically to adjust gray threshold nineteen ninety-five on " Pattern Recognition " magazine to cut apart, discern interested organ one by one according to easy first and difficult later order then.This method is to the single Threshold Segmentation that remains of each organ employing, and the zone that obtains is discontinuous often.People such as John E.Koss proposed to choose image texture features in 1999 on " IEEE Transactions On Medical Imaging " magazine, adopt the Hopfield network to carry out cluster segmentation.Here cut apart the clusters number that the organ that obtains depends on appointment in the Hopfield network.People such as Chien-Cheng Lee 2003 propose to adopt fuzzy diagnosis rule on the context dependent neural net base each zone is discerned and to be merged at " IEEE Transactions On Information Technology In Biomedicine " magazine.The foundation of the fuzzy recognition rule of such cover requires the checking of suitable experience and test of many times.People such as Hyunjin Park 2003 propose then will to be integrated in Bayes's segmentation framework by registration about abdomen organ's anatomy form collection of illustrative plates at " IEEETransactions On Medical Imaging " magazine ^[5]This method needs by hand a plurality of reference mark to be set, and has influenced the robotization of method.

Summary of the invention

The objective of the invention is to provide the organ partitioning algorithm of a kind of abdominal CT of on PC, moving and nuclear-magnetism image, be characterized in carrying out simultaneously the robotization of many organs of image midriff and cut apart, implement simply to need not manual intervention for clinical practice.

The key step of algorithm comprises:

Step (1) in PC, adopts the method for registering based on normalized mutual information, and collection of illustrative plates is made whole registration to the individual images that obtains by CT or nuclear magnetic scanning, to eliminate the whole difference between collection of illustrative plates and target image; Described collection of illustrative plates is meant the image of existing medical images being carried out mark after the anatomy expert is cut apart by hand to each organic region, and it can provide the reference information of human anatomic structure for PC; Described normalized mutual information nMI (S, T*A) represent with following formula:

$\mathrm{nMI}(S,T*A)=\frac{H\left(S\right)+H(T*A)}{H(S,T*A)}$

Wherein, S is a target image,

T*A representative is done image behind the spatial alternation T to collection of illustrative plates A;

H (S) represents the entropy of target image, the entropy of collection of illustrative plates behind H (T*A) representation transformation,

(S T*A) represents the combination entropy between collection of illustrative plates after target image and the conversion to H;

NMI is the abbreviation (normalized mutual information) of the normalized mutual information of English, nMI (S, T*A) normalized mutual information between expression S and the T*A;

Spatial alternation T be the affined transformation that comprises 9 parameters: T=(p, q, r, u, v, w, φ, ω, θ);

Wherein, p, q, r are respectively x behind the spatial alternation, y, the displacement of z direction; R, u, v are respectively x behind the spatial alternation, y, the proportional zoom of z direction; φ, ω, θ are respectively behind the spatial alternation around x, y, the corner of z axle;

If: before the conversion point coordinate be X=(x, y, z),

Then: the coordinate after the conversion is X '=X ', y ', z '),

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ z^{\′}\\ 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& p\\ 0& 1& 0& q\\ 0& 0& 1& r\\ 0& 0& 0& 1\end{array}\right]\×\left[\begin{array}{cccc}\cos \mathrm{\φ}& \sin \mathrm{\φ}& 0& 0\\ -\sin \mathrm{\φ}& \cos \mathrm{\φ}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\×\left[\begin{array}{cccc}0& 0& -\sin \mathrm{\ω}& 0\\ 0& 1& 0& 0\\ \sin \mathrm{\ω}& 0& \cos \mathrm{\ω}& 0\\ 0& 0& 0& 1\end{array}\right]$

$\&times;\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{\&theta;}& \sin \mathrm{\&theta;}& 0\\ 0& -\sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}u& 0& 0& 0\\ 0& v& 0& 0\\ 0& 0& w& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{c}x\\ \mathrm{<figure-callout\; id="1"\; label="y\; z"\; filenames="A20071006306600361.png"\; state="\{\{state\}\}">y</figure-callout>}& \\ z\\ 1\end{array}\right]$

So far, defined target image S and through the normalized mutual information between the collection of illustrative plates A behind the spatial alternation T; Described method for registering based on normalized mutual information is realized successively according to the following steps:

Step (1.1). determine initial spatial alternation T at random ⁰=(T ⁰, q ⁰, r ⁰, u ⁰, v ⁰, w ⁰, φ ⁰, ω ⁰, θ ⁰) in nine parameter values, and act on collection of illustrative plates A with the spatial alternation that obtains thus, obtain the collection of illustrative plates T after the conversion ⁰* A;

Step (1.2). and the associating normalization joint histogram h between the collection of illustrative plates T*A behind calculating target image S and the process spatial alternation T (l, k)

h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}

Wherein, S (x)=l is illustrated in that pixel x corresponding gray is l in the target image,

It is k that T*A (x)=k is illustrated in through x corresponding gray in the collection of illustrative plates T*A pixel of spatial alternation;

Num is the number of pixel x;

Step (1.3). the joint probability distribution p between the collection of illustrative plates T*A behind calculating target image S and the process spatial alternation T _{S, T*A}(l, K):

$p_{S,T*A}(l,k)=\frac{h(l,k)}{\underset{l,k}{\mathrm{\&Sigma;}}h(l,k)}.$

Step (1.4). the marginal probability distribution p of difference computed image S _S(l) and the marginal probability distribution p of T*A _T*A(k):

$p_S\left(l\right)=\underset{k}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k),$

$p_{T*A}\left(k\right)=\underset{l}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k).$

The marginal probability distribution is the notion in the Probability Statistics Theory, for two stochastic variable x, the random vector that y forms (x, y), the probability distribution of its component x is called random vector (x, the marginal distribution about x y); Here be the collection of illustrative plates T*A after image S and the conversion respectively as two stochastic variable l and k, the marginal probability that calculates l and k then respectively distributes;

Step (1.5). calculate the entropy H (S) of target image as follows respectively, the combination entropy H between collection of illustrative plates after the entropy H (T*A) of collection of illustrative plates and target image and the conversion after the conversion (S, T*A);

$H\left(S\right)=-\underset{l}{\mathrm{\&Sigma;}}{p}_S\left(l\right)\log p_S\left(l\right)$

$H(T*A)=-\underset{k}{\mathrm{\&Sigma;}}{p}_{T*A}\left(k\right)\log p_{T*A}\left(k\right)$

$H(S,T*A)=-\underset{k}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k){\log }_{S,T*A}(l,k)$

Step (1.6). the normalized mutual information between the collection of illustrative plates T*A after being calculated as follows target image S and passing through the T conversion:

$\mathrm{nMI}(S,T*A)=\frac{H\left(S\right)+H(T*A)}{H(S,T*A)}$

Step (1.7). (S T*A) brings the Powell optimized Algorithm into spatial alternation T is optimized iteration repeatedly, until till the convergence with normalized mutual information nMI;

The Powell optimized Algorithm is a kind of traditional determinacy optimization method, and is considered to the most effective excellent method that need not to differentiate at present; This basic idea is, for n dimension extreme-value problem, at first asks minimum along n coordinate direction, obtains n conjugate direction afterwards through n time, asks minimum along n conjugate direction then, through repeatedly just trying to achieve minimal value after the iteration;

, utilize the Powell optimized Algorithm that the value of spatial alternation T is optimized here, in the hope of obtaining mutual information nMI (S, minimal value T*A); Spatial alternation T contains nine components, and promptly (p, q, r, u, v, w, φ, ω, θ), nine dimensions in these nine corresponding Powell optimized Algorithm of components difference promptly utilize the Powell optimized Algorithm to ask nine dimension extreme-value problems; When Powell optimized Algorithm search end, promptly after the algorithm convergence, (S, the pairing spatial alternation value of minimal value T*A) is optimal spatial conversion T to the nMI that obtains ^*, T ^*=(p ^*, q ^*, r ^*, u ^*, v ^*, w ^*, φ ^*, ω ^*, θ ^*); T then ^** A is the registration collection of illustrative plates that whole registration will obtain;

Step (2). according to the following steps interested each organ is carried out independent registration successively, determine the roughly corresponding region of each interested organ in target image S; Specifically be divided into following a few step:

Step (2.1). determine the interval R of tonal range by following formula for each organ of interest _k

R _k＝[μ _k-λσ _k，μ _k+λσ _k]，

Wherein, μ _k, σ _kBe k the organ of interest of collection of illustrative plates T*A after conversion gray average and the standard variance in the occupied zone in target image S;

λ is a constant, is used for adjusting the uniform gray level scope; Test finds that the span of λ is 1.2～1.3;

Step (2.2). with the final spatial alternation T that obtains in the step (1) ^*Oppose and carry out the initial space conversion T ' of registration for each organ separately, even T '=T ^*=(p ^*, q ^*, r ^*, u ^*, v ^*, w ^*, φ ^*, ω ^*, θ ^*);

Step (2.3). be calculated as follows the similarity measure M of the organic region after the conversion _k(T '):

$M_k\left(T^{\&prime;}\right)=N_k^{\mathrm{in}}\left(T^{\&prime;}\right)-N_k^{\mathrm{out}}\left(T^{\&prime;}\right).$

Wherein, subscript k representative is carried out registration to k organ,

N _k ^InThe organ number of pixels that (T ') representative adopts spatial alternation T ' time organ k to comprise in the corresponding region in the registration collection of illustrative plates,

N _k ^Out(T ') represents the non-organ number of pixels that comprises in this corresponding region;

$\left\{\begin{array}{c}{N}_k^{\mathrm{in}}\left(T^{\&prime;}\right)=N\left\{x\right|x\&Element;C_k\left(T^{\&prime;}\right),S\left(x\right)\&Element;R_k\}\\ N_k^{\mathrm{out}}\left(T^{\&prime;}\right)=N\left\{x\right|x\&Element;C_k\left(T^{\&prime;}\right),S\left(x\right)\&NotElement;R_k\}\end{array}\right.$

Wherein, N{} represents the number of pixels that comprises in the set of computations,

C _kThe zone of (T ') expression organ k correspondence under spatial alternation T,

The gray-scale value of S (x) remarked pixel x in CT or MRI image;

Step (2.4). with the similarity measure M of the organic region after the conversion _k(T ') brings the Powell optimized Algorithm into T ' iterated, till convergence; Detailed process promptly utilizes the Powell optimized Algorithm that the value of spatial alternation T ' is optimized, in the hope of obtaining similarity measure M _kThe minimal value of (T '); Spatial alternation T ' contains nine components, and promptly (p ', q ', r ', u ', v ', w ', φ ', ω ', θ '), nine dimensions in these nine corresponding Powell optimized Algorithm of components difference promptly utilize the Powell optimized Algorithm to ask nine dimension extreme-value problems; When Powell optimized Algorithm search finishes, promptly after the algorithm convergence, the M that obtains _kThe pairing spatial alternation value of the minimal value of (T ') is optimal spatial conversion T ' ^*, T ' ^*=(p ' ^*, q ' ^*, r ' ^*, u ' ^*, v ' ^*, w ' ^*, ω ' ^*, θ ' ^*); C then _k(T ' ^*) organic region that will obtain for independent registration;

Step (3). fuzzy connection is cut apart; Specifically be divided into following a few step:

Step (3.1). establish C _k(T ' ^*) for the fuzzy prime area that k organ cut apart that connects, be calculated as follows the grey level histogram H in this zone _k(l):

H _k(l)＝N{x|x∈C _k，S(x)＝l}，l＝0，1……L

H _k(l) can be regarded as C _k(T ' ^*) gray scale is the number of the pixel x of l in the zone;

Step (3.2). adopt the method for curve fitting that histogram is carried out level and smooth and denoising, obtain the histogram after the match

${\hat{H}}_k\left(l\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_j\exp [-{\left(\frac{l-b_j}{c_j}\right)}^2],$

Wherein, n represents the number of gaussian component, i.e. the number of the organ that will cut apart; For example we cut apart liver and two kinds of organs of spleen, then n=2;

a _j, b _j, c _jRepresent peak value, center and the width of j Gaussian distribution successively;

Step (3.3). the histogram curve after the match can be regarded the grey level probability density function of organ in target image as; Select the initial seed point of at least one pixel of corresponding maximum gray probability as organ;

Step (3.4). investigate the gray probability sum of each seed points adjacent area; Adjacent area be defined as with as the collection of pixels of the seed points of center pixel distance less than the distance of a pixel unit; Described gray probability by have the number of pixels of this gray scale and the ratio of regional interior pixel sum in the research organic region; Have naming a person for a particular job of maximum region gray probability summation and finally be defined as seed points;

Step (3.5). with the histogram that obtains by step (3.2)

The position of the probability density function of pairing organic region intensity profile and final seed points is an initial parameter, calculates the fuzzy strength of joint of each pixel and the final seed points of organ in the entire image with the fuzzy join algorithm of the described usefulness of following formula:

$c_{\mathrm{\&mu;}}(p,q)=\underset{\mathrm{\&rho;}(p,q)\&Element;P(p,q)}{\max }\left[\underset{z\&Element;\mathrm{\&rho;}(p,q)}{\min }\widehat{H_k}\left(l\left(z\right)\right)\right]$

Wherein, p represents final seed points, and q represents in the entire image more arbitrarily,

c _μ(p, q) the fuzzy strength of joint between expression p and the q;

ρ (p, q) any communication path between expression p and the q,

P (p, the q) set of all communication paths between expression p and the q,

Z ∈ ρ (p, q) expression Z be path ρ (p, q) any pixel on,

L (z) is the gray-scale value of pixel Z in target image,

Expression gray-scale value l (z) is at histogram

In pairing probability density value;

This pixel that fuzzy strength of joint characterizes belongs to the degree of membership of organic region, and the degree of membership of each pixel is called fuzzy connection layout picture as the resulting image of gray scale;

Step (3.6). for the fuzzy connection layout picture of each organ is specified an optimal threshold; So-called threshold value is used for exactly and will belongs to the pixel of target organ and the numerical value that background pixel distinguishes, and the regulation gray scale is the pixel that belongs to target organ more than or equal to the pixel of this threshold value, and gray scale is a background pixel less than the pixel of this threshold value; So-called optimal threshold, be exactly thus the zone that assembled of the organ pixel that obtains of threshold value meet real organic region most; Optimal threshold is the important parameter that fuzzy connection layout is looked like to be converted into final segmentation result; Determine that the optimal threshold concrete steps are as follows:

At first set the threshold to 1, progressively reduce threshold value then; Be used for reducing each time the threshold value that obtains and cut apart fuzzy connection layout picture, obtain the cut zone of organ; Observe the change of the shape of cut zone with threshold value, if when a certain subthreshold reduces, shape acutely changes, and then stops to reduce threshold value, will go up a threshold value and be decided to be optimal threshold;

About the measurement way of the alteration of form of cut zone and definite method of optimal threshold, be further explained in detail as follows:

A). the way of the measurement of the alteration of form of cut zone:

Adopt area relative changing value and tight ness rating relative changing value to weighing the alteration of form of cut zone; Tight ness rating (Compactness) is defined as

Area relative changing value and tight ness rating relative changing value are defined as follows: reduce for each threshold value, calculate organ area and tight ness rating before and after reducing respectively; Organ area and tight ness rating before order reduces are original area and former tight ness rating, and organ area after the reduction and tight ness rating are new area and new tight ness rating, and then the relative variation of area and tight ness rating is defined as

B). determine the method for optimal threshold with area relative changing value and tight ness rating relative changing value:

In threshold value reduction process, if when certain once reduces, the relative changing value of area of discovery and tight ness rating thinks then that greater than predefined reference value ρ violent change has taken place shape, thus the outage threshold search, and the threshold value before will reducing is decided to be optimal threshold; Test shows that reference value ρ value 0.15～0.2 can obtain optimal threshold more accurately;

Step (3.7), the optimal threshold that obtains with step (3.6) according to following formula look like to cut apart to fuzzy connection layout, obtain the two-value split image of organ:

If the optimal threshold of k organ is T _k, then the binary segmentation result of this organ is

Wherein, x represents arbitrary pixel,

F _k(x) be the gray-scale value of the fuzzy connection layout picture of k organ at pixel x place,

B _k(x) be the gray-scale value of the bianry image after cutting apart at pixel x place;

Step (4). the organ shape correction, because causing because gray scale is close to blur to connect, two organs being close to mutually in the human abdominal cavity cut apart issuable erroneous segmentation result to solve;

In the medical image, also similar phenomena of the gray scale of two organs of being close to mutually in human abdominal cavity in image often appears; For two such organs, in the medical image on its surface mutually subsides place may not have tangible border, this can make fuzzy join algorithm think that these two organs are to grow together, thereby is divided into same organ mistakenly; Originally research and propose with organ shape and stabilize to according to detecting erroneous segmentation, and adopt range conversion the erroneous segmentation result to be carried out shape corrections in conjunction with the watershed segmentation algorithm; Idiographic flow is as follows:

Step (4.1). for the multilayer medical image, according to the following steps from top to bottom, adopt the method for zone-by-zone analysis to detect the aspect that has erroneous segmentation;

Step (4.1.1).

Set: when anterior layer is the aspect that algorithm is being analyzed, when the sequence number of anterior layer is i

The below adjacent layer is and works as the adjacent below aspect of anterior layer that the sequence number of below adjacent layer is i-1;

O _k ⁱBe the organ cut zone in anterior layer;

O _k ^I-1Be organ cut zone in the adjacent layer of below;

In the anterior layer, making the each point in the organ cut zone is negative value to the distance of zone boundary, for all do not belong to this regional point in the entire image, i.e. all points beyond should the zone, make its distance to the zone boundary be on the occasion of;

Step (4.1.2). to the organ cut zone O in anterior layer i _k ⁱCarry out range conversion, promptly calculate their minimum distance D to the zone boundary to whenever pressing following formula in the entire image _k ⁱ:

Wherein, y has represented in the pixel of zone boundary nearest with x;

B _k ⁱ(x) the final segmentation result B for obtaining in the step (3.7) _k ^*(x) at the pixel value in anterior layer i;

With distance value d _k ⁱ(x) compose to each pixel as gray-scale value, just obtained range image D _k ⁱ

Step (4.1.3). extract regional O in the adjacent layer of below _k ^I-1The border on the each point locations of pixels, then at D _k ⁱIn the D of each point at same position place _k ⁱ(x) value can be thought O _k ^I-1The distance value of borderline pixel; Obtain all O _k ^I-1The distance value of boundary pixel, calculate the standard variance of these distance values, to determine the similarity degree of organ between i layer and i-1 layer segmentation result; Detected segmentation errors if standard variance greater than predefined reference value, is then thought, following one deck is the aspect that has erroneous segmentation; Rule of thumb, this reference value generally is set in 5 pixel distances;

According to the method described above, detect structure at all levels from top to bottom, obtain the aspect that there is erroneous segmentation in all;

Step (4.2). by following formula the bianry image of each aspect of having erroneous segmentation is carried out range conversion;

Set:, make that its distance value is a negative value for all pixels in the organic region; For all do not belong to the zone point in addition in this zone outside the zone in the entire image, make its distance value be negative minimal value, i.e. negative infinite-∞;

Distance map D then _k ⁱ' (x)

Wherein y has represented the frontier point nearest with x; B _k ⁱThe bianry image that promptly has the aspect of erroneous segmentation;

Step (4.3). the range image D that step (4.2) is obtained by the following method _k ⁱ' application watershed algorithm;

Step (4.3.1). entire image is set a gray threshold T; The initial value that makes T is image D _k ⁱ' in except negative minimum value infinite, that is:

$T_0=\min \left(\right\{D_k^i\left(x\right)|D_k^i\left(x\right)\&NotEqual;-\&infin;\}),$

Wherein, the coordinate of any pixel in the x representative image;

Step (4.3.2). with d is that step-length increases threshold value T gradually, promptly

T _n＝T ₀+n·d，

Wherein, T _nThe value of threshold value T when representing the n time increase step-length; N span 1≤n≤K, wherein $K=(\max \left(D_k^i\left(x\right)\right)-T_0)/d,$ Max (D _k ⁱ(x)) presentation video D _k ⁱ' maximal value;

The d value is fixed, and rule of thumb gets d=3;

Step (4.3.3). utilizing increases the threshold value T that step-length obtains the n time _nPress following formula with image D _k ⁱ' be divided into bianry image W _n:

Wherein, the coordinate of any pixel in the x representative image;

W _n(x)=1 expression is submerged in the following part of the water surface, W _n(x)=0 expression is exposed at the above land part of the water surface; Increase step-length for preceding once (promptly the n-1 time) and obtain threshold value T _N-1, press following formula with image D _k ⁱ' be divided into bianry image W _N-1:

Wherein, the coordinate of any pixel in the x representative image;

Step (4.3.4). at T _nPropagation process in, whenever T _nGreater than a local maximum, with regard to key diagram as W _N-1In two connected regions are arranged at W _nIn be merged into a connected domain, then with W _N-1And W _nTwo width of cloth images carry out xor operation, obtain watershed divide image S _n:

Wherein, the coordinate of any pixel in the x representative image,  represents xor operation, S _n(x)=1 remarked pixel x is the position, watershed divide;

Step (4.3.5). to each the n value among n span 1≤n≤K, all can obtain a corresponding watershed divide image S _nIf all pairing watershed divide of n value images are S ₁, S ₂..., Then the result images S of watershed transform is the set of all watershed divide images:

$S={\mathrm{<figure-callout\; id="1"\; label="S"\; filenames="A20071006306600361.png"\; state="\{\{state\}\}">S</figure-callout>}}_1\&CenterDot;S_2\&CenterDot;...\&CenterDot;S_{\max \left(D_k^i\left(x\right)\right)};$

All pixel values are that 1 pixel is with regard to presentation video D among the image S _k ⁱ' in the position at place, watershed divide; By these watershed divide is image D _k ⁱ' be divided into plurality of sub-regions;

Step (4.4). the image D that watershed algorithm is obtained _k ⁱ' in each subregion, calculate its adjacent aspect organ cut zone O respectively with the below _k ^I-1Area overlap ratio:

Be lower than 50% zone and be considered as the erroneous segmentation zone overlapping ratio, be labeled as the background area again, thus with it at range image D _k ⁱ' in remove.

Experiment effect is analyzed:

In order to verify the reliability of originally cutting apart algorithm, experimental selection 5 covers come from the Clinical CT abdomen scanning image of domestic hospital.This 5 sets of data comprises masculinity and femininity, and the object of all ages and classes.Find out that from segmentation result the organ of interest in the test figure has all obtained effectively cutting apart.Collection of illustrative plates has provided roughly distributed areas behind the registration, and through further cutting apart and revising, final segmentation result is compared the two basically identical with manual segmentation result.

Description of drawings

Fig. 1 is an aspect in the VIP-Man collection of illustrative plates of invention employing.

Fig. 2 is an aspect of CT image.

Fig. 3 is the main process figure of computer program.

Fig. 4 be with collection of illustrative plates behind CT integral image registration, the profile stack of extracting organ of interest in the collection of illustrative plates is presented on the CT image.Organ of interest comprises liver, kidney and spleen.

After Fig. 5 is organ registration, organ site contrast in the profile of organ of interest and the CT image in the collection of illustrative plates.

Fig. 6 is the grey level histogram and the curve fitting synoptic diagram of each organ.

Fig. 7 is the fuzzy connection cutting procedure of liver.Be followed successively by from left to right, liver part in the CT image, the fuzzy connection layout picture of liver looks like to carry out liver segmentation results after the Threshold Segmentation to fuzzy connection layout.

Fig. 8 is the watershed algorithm synoptic diagram.

Fig. 9 is for to carry out the shape corrections process to kidney.Be followed successively by from left to right, kidney part in the CT image, the fuzzy connection segmentation result of kidney, kidney is split into the adjacent aspect kidney regional correlation that a plurality of zones and white contours are represented through watershed algorithm, eliminates to overlap segmentation result after the kidney correction that obtains behind the low excessively zone of ratio.

Figure 10 is to the whole segmentation result of CT image.

Embodiment

The present invention proposes a kind of abdomen organ's dividing method, realize that the robotization of the many organs of medical image midriff is cut apart based on the collection of illustrative plates coupling.Invention is adopted by the manual collection of illustrative plates of cutting apart foundation of U.S. visual human CT data process, and the abdomen organ is cut apart.The paper " VIP-man:An image-based whole-body adult male model constructed from color photographs ofthe visible human project for multi-particle Monte Carlo calculations " that can deliver at " Health Physics " magazine in 2000 referring to people such as Xu about the detailed description of this collection of illustrative plates method for building up.The visible Fig. 1 of the example of this collection of illustrative plates.By registration work, what obtain is a collection of illustrative plates that roughly aligns with target image, and target image promptly needs to do the clinical image that organ is cut apart, and its example is seen Fig. 2.The purpose of registration is roughly alignd collection of illustrative plates by deformation and target image exactly, and this alignment concentrates on the alignment of organ of interest, rather than the alignment of image All Ranges, has so just alleviated the pressure and the degree of difficulty of registration work greatly.Thereby allow in registration, to select better simply spatial alternation, only comprise less degree of freedom.Optimizing process has also improved speed accordingly, has strengthened stability.On partitioning algorithm, select the fuzzy partitioning algorithm that connects, it has stronger stability, can adapt to cutting apart of complicated image.An appointment that subject matter is every call parameter that faces is cut apart in fuzzy connection, comprises initial seed points, the gray distribution features of organ of interest etc.Utilize collection of illustrative plates behind the registration, the present invention proposes the robotization designation method of these parameters, avoided artificial or experience commonly used to specify.Find that because the gray scale between adjacent organs is overlapping or overlap phenomenon, the possibility of result that fuzzy connection is cut apart is comprising wrong zone.In order to eliminate these zone errors, improve the accuracy rate of segmentation result, this research based on organ about have the conforming principle of shape in the CT image of two adjacent aspects, the organ shape modification method has been proposed, segmentation result is done further perfect.

The main process figure of the computer program of this method can participate in Fig. 3.

The inventive method realizes by following step:

1. whole registration

The target of whole registration is the whole difference of eliminating between collection of illustrative plates and target image, comprises the difference of image space and body shape (height is fat or thin etc.).Whole registration, searching be a kind of spatial alternation, make collection of illustrative plates through behind this spatial alternation, can roughly align with target image on the whole.Select similarity transformation (comprising) in the whole registration, as objective function, adopt the Powell algorithm to be optimized with normalized mutual information with respect to x, y, z axle 3 translations of axes, rotation and transformations of scale, totally 9 parameters.

Normalized mutual information between collection of illustrative plates and target image (Normalized Mutual Information below is abbreviated as nMI) is defined as

$\mathrm{nMI}(S,T*A)=\frac{H\left(S\right)+H(T*A)}{H(S,T*A)}$

Wherein A represents collection of illustrative plates, and T*A representative is done image behind the spatial alternation T to collection of illustrative plates, and H (S) represents the entropy of target image, the entropy of collection of illustrative plates behind H (T*A) representation transformation, and (S T*A) represents the combination entropy between collection of illustrative plates after target image and the conversion to H.

On spatial alternation, select to comprise the affined transformation of 9 parameters: T=(p, q, r, u, v, w, φ, ω, θ).Wherein, p, q, r are respectively x behind the spatial alternation, y, the displacement of z direction; U, v, w are respectively x behind the spatial alternation, y, the proportional zoom of z direction; φ, ω, θ are respectively rich x behind the spatial alternation, y, the corner of z axle;

Make before the conversion that point coordinate is that (z), the coordinate after the conversion is X '=(x ', y ', z ') to X=, then for x, y

$\left[\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\\ z^{\&prime;}\\ 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& p\\ 0& 1& 0& q\\ 0& 0& 1& r\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}\cos \mathrm{\&phi;}& \sin \mathrm{\&phi;}& 0& 0\\ -\sin \mathrm{\&phi;}& \cos \mathrm{\&phi;}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}0& 0& -\sin \mathrm{\&omega;}& 0\\ 0& 1& 0& 0\\ \sin \mathrm{\&omega;}& 0& \cos \mathrm{\&omega;}& 0\\ 0& 0& 0& 1\end{array}\right]$

$\&times;\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{\&theta;}& \sin \mathrm{\&theta;}& 0\\ 0& -\sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}u& 0& 0& 0\\ 0& v& 0& 0\\ 0& 0& w& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{c}\mathrm{<figure-callout\; id="1"\; label="x\; y\; z"\; filenames="A20071006306600361.png"\; state="\{\{state\}\}">x</figure-callout>}& \\ y\\ z\\ 1\end{array}\right]$

$H\left(S\right)=-\underset{l}{\mathrm{\&Sigma;}}{p}_S\left(l\right)\log p_S\left(l\right)$

$H(T*A)=-\underset{k}{\mathrm{\&Sigma;}}{p}_{T*A}\left(k\right)\log p_{T*A}\left(k\right)$

$H(S,T*A)=-\underset{k}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k){\log }_{S,T*A}(l,k)$

Here, p _{S, T*A}(l k) is joint probability distribution between image S and T*A, can represent with normalized joint histogram:

$p_{S,T*A}(l,k)=\frac{h(l,k)}{\underset{l,k}{\mathrm{\&Sigma;}}h(l,k)}.$

p _s(l), p _T*A(k) be that marginal probability between image S and T*A distributes, calculate by following formula:

$p_S\left(l\right)=\underset{k}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k),$

$p_{T*A}\left(k\right)=\underset{l}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k).$

The marginal probability distribution is the notion in the Probability Statistics Theory, for two stochastic variable x, the random vector that y forms (x, y), the probability distribution of its component x is called random vector (x, the marginal distribution about x y); Here be the collection of illustrative plates T*A after image S and the conversion respectively as two stochastic variable l and k, the marginal probability that calculates l and k then respectively distributes;

(l k) is defined as normalization joint histogram h

h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}，

That is, pixel gray scale in image S is l, and respective pixel is the number of pixels statistics of k in image T*A.

So far, defined target image S and through the mutual information between the collection of illustrative plates A behind the spatial alternation T.

The Powell optimized Algorithm is a kind of traditional determinacy optimization method, and is considered to the most effective excellent method that need not to differentiate at present.This basic idea is, for n dimension extreme-value problem, at first asks minimum along n coordinate direction, obtains n conjugate direction afterwards through n time, asks minimum along n conjugate direction then, through repeatedly just trying to achieve minimal value after the iteration.

, utilize the Powell optimized Algorithm that the value of spatial alternation T is optimized here, in the hope of obtaining mutual information nMI (S, minimal value T*A).Spatial alternation T contains nine components, and promptly (p, q, r, u, v, w, φ, ω, θ), nine dimensions in these nine corresponding Powell optimized Algorithm of components difference promptly utilize the Powell optimized Algorithm to ask nine dimension extreme-value problems.When Powell optimized Algorithm search end, promptly after the algorithm convergence, (S, the pairing spatial alternation value of minimal value T*A) is optimal spatial conversion T to the nMI that obtains ^*, T ^*=(p ^*, q ^*, r ^*, u ^*, v ^*, w ^*, φ ^*, ω ^*, θ ^*); T then ^** A is the registration collection of illustrative plates that whole registration will obtain;

Fig. 4 has showed the effect of collection of illustrative plates behind CT integral image registration.We will be through optimal spatial T ^*The profile of the organ of interest in the collection of illustrative plates after the conversion extracts and superposes and is presented on the CT image, to show effect.

2. organ registration

The basic thought of organ registration is further to eliminate the position difference between unsolved corresponding organ in the whole registration.Basic skills is to the registration on individual images respectively of each organ of interest in the collection of illustrative plates, also will guarantee not occurrence positions conflict between each organ simultaneously.

Independent registration to each organ is still selected similarity transformation, adopts the Powell optimized Algorithm.Brand-new objective function of organ registration definition, all there is a corresponding region in each organ of interest in the collection of illustrative plates in target image, organ registration will be made great efforts to make and be comprised organ pixel as much as possible and the least possible non-organ pixel in this corresponding region.Owing to can't determine also whether each pixel belongs to certain organ here, therefore, adopt a kind of supposition standard.If total n organ of interest then to k organ of interest wherein, is speculated as between its gray area,

R _k＝[μ _k-λσ _k，μ _k+λσ _k]，

Wherein, μ _k, σ _kBe k organ of interest of collection of illustrative plates gray average and the standard variance in the occupied zone behind registration.

Wherein λ is a constant, is used for adjusting the uniform gray level scope.Test finds that the span of λ is 1.2～1.3.As the pixel gray scale at interval R _kIn, think that then it belongs to organ k pixel; Otherwise, think that it belongs to non-organ pixel.The target of registration is by optimization algorithm, and each organ is determined a new similarity transformation respectively; Here the optimization algorithm of Cai Yonging is the Powell algorithm.New conversion makes that the organ pixel is many as far as possible in the organ corresponding region, and non-organ pixel is the least possible.The similarity measure M of definition organ registration _k(T) as follows:

$M_k\left(T\right)=N_k^{\mathrm{in}}\left(T\right)-N_k^{\mathrm{out}}\left(T\right).$

Wherein, subscript k representative is carried out registration to k organ, N _k ^InThe organ number of pixels that comprises in the corresponding region of organ k when (T) spatial alternation T is adopted in representative, N _k ^Out(T) represent the non-organ number of pixels that comprises in the corresponding region.

$\left\{\begin{array}{c}{N}_k^{\mathrm{in}}\left(T\right)=N\left\{x\right|x\&Element;C_k\left(T\right),S\left(x\right)\&Element;R_k\}\\ N_k^{\mathrm{out}}\left(T\right)=N\left\{x\right|x\&Element;C_k\left(T\right),S\left(x\right)\&NotElement;R_k\}\end{array}\right.$

Wherein, N{} represents the number of pixels that comprises in the set of computations, C _k(T) prime area of expression organ k correspondence under spatial alternation T, the gray-scale value of S (x) remarked pixel in CT or MRI image.

The concrete optimizing process of Powell optimized Algorithm is: with the similarity measure M of the organic region after the conversion _k(T ') brings the Powell optimized Algorithm into, utilizes the Powell optimized Algorithm that the value of spatial alternation T ' is optimized, in the hope of obtaining similarity measure M _kThe minimal value of (T '); Spatial alternation T ' contains nine components, and promptly (p ', q ', r ', u ', v ', w ', φ ', ω ', θ '), nine dimensions in these nine corresponding Powell optimized Algorithm of components difference promptly utilize the Powell optimized Algorithm to ask nine dimension extreme-value problems; When Powell optimized Algorithm search finishes, promptly after the algorithm convergence, the M that obtains _kThe pairing spatial alternation value of the minimal value of (T ') is optimal spatial conversion T ' ^*, T ' ^*=(p ' ^*, q ' ^*, r ' ^*, u ' ^*, v ' ^*, w ' ^*, φ ' ^*, ω ' ^*, θ ' ^*); C then _k(T ' ^*) organic region that will obtain for independent registration; Fig. 5 has showed through after the organ registration, organ site contrast in the profile of the organ of interest in the collection of illustrative plates and the CT image.

3. fuzzy the connection cut apart

Fuzzy join algorithm is proposed by J.Udupa, and it obtains the degree of membership that each pixel belongs to this organ by the fuzzy strength of joint of each pixel and organ seed points in the definition image.The degree of membership of each pixel is promptly blured the connection layout picture as image gray.Be applied to organ and cut apart in order to blur method of attachment, at first need to specify the gray distribution features of each organ of interest and seed points position (seed points of organ be meant fuzzy join algorithm cut apart the initial pixel point when beginning to cut apart).After these parameters were determined, fuzzy join algorithm can be calculated the fuzzy connection layout picture of each organ of interest respectively; Fuzzy connection layout similarly is a gray level image, and wherein each gray values of pixel points has been represented fuzzy the be connected degree of this point with seed points.After fuzzy connection layout picture generates, also need further to specify a threshold value, be translated into two-value split image, as final segmentation result to organ.

Collection of illustrative plates behind the registration has provided approximate location and the scope of organ in target image, is called the prime area of organ.The grey level histogram of pixel in the statistics prime area.With k organ is example, and the computing formula of grey level histogram is,

H _k(l)＝N{x?|?x∈C _k，S(x)＝l}，?l＝0，1……L，

Promptly gray scale is the number of the pixel of l in the prime area of organ k.

The organ intensity histogram has reflected that grey scale pixel value and pixel belong to the relation between the possibility size of organ.Because what comprise in the prime area all is not the organ pixel, also has non-organ pixel.Therefore adopt the method for curve fitting that histogram is carried out level and smooth and denoising.Curve fitting is considered as the stack of several Gaussian functions with histogram, adopts minimum variance to carry out match.

The organ fitting formula is

${\hat{H}}_k\left(l\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_j\exp [-{\left(\frac{l-b_j}{c_j}\right)}^2],$

Wherein, n represents the number of gaussian component, i.e. the number of the organ that will cut apart; For example we cut apart liver and two kinds of organs of spleen, then n=2; a _j, b _j, c _jRepresent peak value, center and the width of j Gaussian distribution successively.Between gaussian component number and fitting precision, need carry out balance, obtain high as far as possible fitting precision with the least possible gaussian component.Here, adopt the residual value of adjusting square (Adjusted R-square is called for short AR) to seek optimum gaussian component number.

$\mathrm{AR}=1-\frac{\mathrm{SSE}*L}{\mathrm{SST}*(L-n*3)},$

Wherein

$\mathrm{SSE}=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{(H_k\left(l\right)-{\hat{H}}_k\left(l\right))}^2$

$\mathrm{SST}=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{(H_k\left(l\right)-{\stackrel{\&OverBar;}{H}}_k)}^2,$

${\stackrel{\&OverBar;}{H}}_k^i=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}{H}_k\left(l\right)/(L+1).$

Concrete Gauss curve fitting process is:

(1) makes n=0;

(2) make n=n+1;

(3) adopt minimum SSE method, match formula

(4) calculate AR;

(5) when AR＜ _ARAnd execution in step (2) is returned in n＜5;

In above-mentioned steps, T _ARBe according to the experience preset threshold, be used to control balance fitting precision and gaussian component number.Test shows, T _AR=0.95～0.97 is a more rational threshold setting.

Fig. 6 has showed the grey level histogram and the curve fitting synoptic diagram of each organ.Histogram curve after the match can be regarded the grey level probability density function of organ in target image as.Select the initial seed point of at least one pixel of corresponding maximum gray probability as organ.The number of pixels of corresponding maximum gray probability may further be investigated the gray probability sum of pixel adjacent area again greater than 1.Adjacent area is commonly defined as and the collection of pixels of center pixel distance less than a pixel unit.Described area grayscale probability by have the number of pixels of this gray scale and the ratio of regional interior pixel sum in the research organic region; Have naming a person for a particular job of maximum region gray probability summation and finally be defined as seed points.Stipulate that simultaneously seed points must be positioned at the inside, prime area of organ.

Determine the gray distribution features and the initial seed point of organ, can calculate each organ fuzzy connection layout picture separately.With the histogram that obtains above The position of the probability density function of pairing organic region intensity profile and final seed points is an initial parameter, calculates the fuzzy strength of joint of each pixel and the final seed points of organ in the entire image with the fuzzy join algorithm of the described usefulness of following formula:

$c_{\mathrm{\&mu;}}(p,q)=\underset{\mathrm{\&rho;}(p,q)\&Element;P(p,q)}{\max }\left[\underset{z\&Element;\mathrm{\&rho;}(p,q)}{\min }{\hat{H}}_k\left(l\left(z\right)\right)\right]$

Wherein p represents final seed points, and q represents in the entire image more arbitrarily, c _μ(p, q) the fuzzy strength of joint between expression p and the q; ρ (p, q) any communication path between expression p and the q, P (p, the q) set of all communication paths between expression p and the q, z ∈ ρ (p, q) expression Z be path ρ (l (z) is the gray-scale value of pixel Z in target image for p, q) any pixel on,

Expression gray-scale value l (z) is at histogram In pairing probability density value;

This pixel that fuzzy strength of joint characterizes belongs to the degree of membership of organic region, and the degree of membership of each pixel is called fuzzy connection layout picture as the resulting image of gray scale.And organ really will be split, also need to specify a fuzzy strength of joint threshold value, fuzzy connection layout is looked like to carry out two-value cut apart.So-called threshold value is used for a numerical value that object pixel and background pixel are distinguished exactly; In fuzzy connection layout picture, we think that gray-scale value is the organ pixel more than or equal to the pixel of threshold value, and gray-scale value is a background pixel less than threshold pixels.So-called optimal threshold, be exactly thus the zone that assembled of the organ pixel that obtains of threshold value meet real organic region most; Optimal threshold is the important parameter that fuzzy connection layout is looked like to be converted into final segmentation result.

In order to find best fuzzy connection segmentation threshold, following searching method is proposed.At first set the threshold to 1, progressively reduce threshold value then; Be used for reducing each time the threshold value that obtains and cut apart fuzzy connection layout picture, obtain the cut zone of organ; Dividing method is as follows:

If the threshold value of k organ is T _k, then the Threshold Segmentation result of this organ is

Wherein x represents arbitrary pixel, F _k(x) be the gray-scale value of the fuzzy connection layout picture of k organ at pixel x place.

For threshold value reduction each time, observe the alteration of form of segmentation result with threshold value.If when a certain subthreshold reduced, shape acutely changed, then outage threshold search will be gone up a threshold value and will be decided to be optimal threshold.Here two parameters have been used in the measurement of alteration of form, i.e. area relative changing value and tight ness rating relative changing value.Tight ness rating (Compactness) is defined as

Reduce for each threshold value, calculate organ area and tight ness rating before and after reducing respectively; Organ area and tight ness rating before order reduces are called original area and former tight ness rating, and organ area after the reduction and tight ness rating are called new area and new tight ness rating, and then the relative variation of area and tight ness rating is defined as

If when a certain subthreshold reduced, the relative changing value of area of discovery and tight ness rating thought then that greater than predefined reference value ρ violent change has taken place shape, thus the outage threshold search, and the threshold value before will reducing is decided to be optimal threshold.Test shows that reference value ρ value 0.15～0.2 can obtain optimal threshold more accurately.

Next, will blur connection layout with top definite optimal threshold and look like to be divided into bianry image, as final output.Concrete grammar is as follows: the optimal threshold of establishing k organ is T _k ^*, then the binary segmentation result of this organ is

Wherein x represents arbitrary pixel, F _k(x) be the gray-scale value of the fuzzy connection layout picture of k organ, B at pixel x place _k ^*(x) be the gray-scale value of the bianry image after cutting apart at pixel x place.

Fig. 7 is that example has been showed fuzzy connection cutting procedure with the liver.Be followed successively by from left to right, liver part in the CT image, the fuzzy connection layout picture of liver looks like to carry out liver segmentation results after the Threshold Segmentation to fuzzy connection layout.

4. organ shape correction

In the medical image, also similar phenomena of the gray scale of two organs of being close to mutually in human abdominal cavity in image often appears; For two such organs, in the medical image on its surface mutually subsides place may not have tangible border, this can make fuzzy join algorithm think that these two organs are to grow together, thereby is divided into same organ mistakenly; Originally research and propose with organ shape and stabilize to according to detecting erroneous segmentation, and adopt range conversion the erroneous segmentation result to be carried out shape corrections in conjunction with the watershed segmentation algorithm; Idiographic flow is as follows:

1). for the multilayer medical image, adopt the method for zone-by-zone analysis to detect the aspect that has erroneous segmentation;

Detection order is from the top aspect to below aspect; Move for each aspect, the aspect that claims algorithm analyzing is to work as anterior layer, and adjacent below aspect is called the below adjacent layer with working as anterior layer; If working as the sequence number of anterior layer is i, then the sequence number of below adjacent layer is i-1; In anterior layer, establish O _k ⁱRepresent the organ cut zone; To O _k ⁱCarry out range conversion, promptly in the zoning each discrete point to the distance of zone boundary; For having a few in the zone, make that its distance value is a negative value; For having a few outside the zone (scope comprises that all do not belong to this regional point in the entire image), make its distance value on the occasion of, obtaining range image like this is D _k ⁱ,

The discrete pixel at any place in the x representative image wherein, y have been represented point nearest with x on the zone boundary; B _k ⁱ(x) be the final segmentation result B of fuzzy connection _k ^*(x) at the the pixel value when anterior layer.

Below in the adjacent layer, suppose O _k ^I-1Represent the organ cut zone; Extract O _k ^I-1The position of boundary pixel, then at D _k ⁱIn the D at same position place _k ⁱ(x) value can be thought O _k ^I-1Borderline pixel and O _k ⁱThe distance value of borderline respective pixel; When obtaining all O _k ^I-1This distance value of boundary pixel after, calculate the standard variance of these distance values, this standard variance has reflected the similarity degree of organ between i layer and i-1 layer segmentation result.If organ just expands uniformly from one deck to another layer or shrinks, then the standard variance value should be less than the predefine reference value.Therefore, can adopt the size of organ frontier distance standard variance between adjacent aspect to weigh wrong unusual outburst area whether occurred.If greater than reference value, then thinking, standard variance detected segmentation errors.According to the test experience, this reference value generally is set in 5 pixel distances.

According to the method described above, detect structure at all levels from top to bottom, obtain the aspect that there is erroneous segmentation in all.

2). the bianry image to the aspect that has erroneous segmentation carries out range conversion;

The range conversion here and top different; For having a few in the zone, make that its distance value is a negative value; For having a few outside the zone (scope comprises that all do not belong to this regional point in the entire image), make the negative minimal value of its distance value; Obtain range image D like this _k ⁱ' be:

The discrete pixel at any place in the x representative image wherein, y have been represented point nearest with x on the zone boundary; B _k ⁱThe bianry image that promptly has the aspect of erroneous segmentation;

3). to image D _k ⁱ' application watershed algorithm.

Introduce the watershed divide technology by the method for rising.Watershed divide (Watershed) algorithm is regarded the gray-scale value of pixel in the image as on the topomap sea level elevation, the corresponding mountain peak of the part that then gray scale is high in the image, the corresponding mountain valley of the part that gray scale is low.Utilize Fig. 8 (a) to discuss, for easy, the one dimension sectional view of the image that only draws.Pixel in the horizontal ordinate representative image, the gray-scale value of ordinate represent pixel.Water level raises gradually to suppose to have water to gush out also from the lowest point.If the water water level of gushing out from two the lowest point has exceeded mountain peak therebetween, these water will converge.If stop converging of water, then need on the mountain peak, build dam, and the height of dam increases with the rising of water level.This process is not all had by water logging along with whole mountain peaks and finishes.All dams of building are divided into a lot of zones to the view picture topomap, and these dams have just constituted the watershed divide in this sheet soil.Each zone of being come out by watershed segmentation is corresponding to the sub regions in the image, so watershed algorithm is with range image D _k ⁱ' be divided into a plurality of subregions.With watershed algorithm with range image D _k ⁱ' be divided into a plurality of subregions detailed process as follows: still with Fig. 8 (a), three mountain peaks that occur among the figure are represented from image D _k ⁱ' in a part of taking out arbitrarily, the pixel value that gray scale is high is represented on the mountain peak, the pixel value that gray scale is low is represented in the mountain valley.If gray threshold T represents the height of the water surface; The initial value that makes T is image D _k ⁱ' in except that negative minimum value infinite, that is:

$T_0=\min$ $\left(\right\{D_k^i\left(x\right)|D_k^i\left(x\right)\&NotEqual;-\&infin;\}),$

Wherein, the coordinate of any pixel in the x representative image.

With d is that step-length increases threshold value T gradually, promptly

T _n＝T ₀+n·d，

Wherein, T _nThe value of threshold value T when representing the n time increase step-length;

N span 1≤n≤K; $K=(\max \left(D_k^i\left(x\right)\right)-T_0)/d,$ Max (D _k ⁱ(x)) presentation video D _k ⁱ' maximal value;

The d value is fixed, and rule of thumb gets d=3.

Utilizing increases the threshold value T that step-length obtains the n time _n(shown in Fig. 8 (c)) can be by following formula with image D _k ⁱ' be divided into bianry image W _n:

Wherein, the coordinate of any pixel in the x representative image.W _n(x)=1 expression is submerged in the following part of the water surface, W _n(x)=0 expression is exposed at the above land part of the water surface.

Increase step-length for preceding once (promptly the n-1 time) and obtain threshold value T _N-1(shown in Fig. 8 (b)) can be by following formula with image D _k ⁱ' be divided into bianry image W _N-1:

Wherein, the coordinate of any pixel in the x representative image.

At T _nPropagation process in, whenever T _nGreater than a local maximum, with regard to key diagram as W _N-1Two interior connected regions are (as the connected domain C among Fig. 8 (b) _jAnd C _J+1) at W _nIn be merged into a connected domain (as the C among Fig. 8 (c) _{J, j+1}), then with W _N-1And W _nTwo width of cloth images carry out xor operation, obtain watershed divide image S _n:

Wherein, the coordinate of any pixel in the x representative image,  represents xor operation, S _n(x)=1 remarked pixel x is the position, watershed divide, shown in Fig. 8 (d).

Because n span 1≤n≤K, each n value of scope all can obtain a corresponding watershed divide image S hereto _nThe then all pairing watershed divide of n value images are S ₁, S ₂...,

The result images S of watershed transform is the set of all watershed divide images:

$S={\mathrm{<figure-callout\; id="1"\; label="S"\; filenames="A20071006306600361.png"\; state="\{\{state\}\}">S</figure-callout>}}_1\&CenterDot;S_2\&CenterDot;...\&CenterDot;S_{\max \left(D_k^i\left(x\right)\right)}.$

Then all pixel values are that 1 pixel is with regard to presentation video D among the image S _k ⁱ' in the position at place, watershed divide; By these watershed divide is image D _k ⁱ' be divided into plurality of sub-regions.

4). to each subregion that watershed algorithm obtains, calculate its adjacent aspect organ cut zone O respectively with the below _k ^I-1Area overlap ratio; Area overlaps ratio and is defined as:

Based on the principle of organ shape unanimity in adjacent aspect, think that overlapping ratio is lower than 50% zone for should be the erroneous segmentation zone, should be labeled as the background area again,, thus with it at range image D _k ⁱ' in remove.

Fig. 9 is that example has illustrated the shape corrections process with the kidney.Be followed successively by from left to right, kidney part in the CT image, the fuzzy connection segmentation result of kidney, kidney is split into the below adjacent aspect kidney regional correlation that a plurality of zones and white contours are represented through watershed algorithm, eliminates to overlap segmentation result after the kidney correction that obtains behind the low excessively zone of ratio.

So far, just all finished the belly medical image cut apart the visible Figure 10 of final segmentation result.

Claims (1)

1. the method that medical image midriff organ is cut apart specifically is characterised in that, this method realizes in PC successively according to the following steps:

The key step of algorithm comprises:

Step (1) in PC, adopts the method for registering based on normalized mutual information, and collection of illustrative plates is made whole registration to the individual images that obtains by CT or nuclear magnetic scanning, to eliminate the whole difference between collection of illustrative plates and target image; Described collection of illustrative plates is meant the image of existing medical images being carried out mark after the anatomy expert is cut apart by hand to each organic region, and it can provide the reference information of human anatomic structure for PC; Described normalized mutual information nMI (S, T*A) represent with following formula:

$\mathrm{nMI}(S,T*A)=\frac{H\left(S\right)+H(T*A)}{H(S,T*A)}$

Wherein, S is a target image,

T*A representative is done image behind the spatial alternation T to collection of illustrative plates A;

H (S) represents the entropy of target image, the entropy of collection of illustrative plates behind H (T*A) representation transformation,

(S T*A) represents the combination entropy between collection of illustrative plates after target image and the conversion to H;

NMI is the abbreviation (normalized mutual information) of the normalized mutual information of English, nMI (S, T*A) normalized mutual information between expression S and the T*A;

Spatial alternation T be the affined transformation that comprises 9 parameters: T=(p, q, r, u, v, w, φ, ω, θ);

Wherein, p, q, r are respectively x behind the spatial alternation, y, the displacement of z direction; R, u, v are respectively x behind the spatial alternation, y, the proportional zoom of z direction; φ, ω, θ are respectively behind the spatial alternation around x, y, the corner of z axle;

If: before the conversion point coordinate be X=(x, y, z),

Then: the coordinate after the conversion is X '=(x ', y ', z '),

$\left[\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\\ z^{\&prime;}\\ 1\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& p\\ 0& 1& 0& q\\ 0& 0& 1& r\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}\cos \mathrm{\&phi;}& \sin \mathrm{\&phi;}& 0& 0\\ -\sin \mathrm{\&phi;}& \cos \mathrm{\&phi;}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}0& 0& -\sin \mathrm{\&omega;}& 0\\ 0& 1& 0& 0\\ \sin \mathrm{\&omega;}& 0& \cos \mathrm{\&omega;}& 0\\ 0& 0& 0& 1\end{array}\right]$

$\&times;\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{\&theta;}& \sin \mathrm{\&theta;}& 0\\ 0& -\sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{cccc}u& 0& 0& 0\\ 0& v& 0& 0\\ 0& 0& w& 0\\ 0& 0& 0& 1\end{array}\right]\&times;\left[\begin{array}{c}x\\ y\\ z\\ 1\end{array}\right]$

So far, defined target image S and through the normalized mutual information between the collection of illustrative plates A behind the spatial alternation T; Described method for registering based on normalized mutual information is realized successively according to the following steps:

Step (1.1). determine initial spatial alternation T at random ⁰=(p ⁰, q ⁰, r ⁰, u ⁰, v ⁰, w ⁰, φ ⁰, ω ⁰, θ ⁰) in nine parameter values, and act on collection of illustrative plates A with the spatial alternation that obtains thus, obtain the collection of illustrative plates T after the conversion ⁰* A;

Step (1.2). the associating normalization joint histogram h between the collection of illustrative plates T*A behind calculating target image S and the process spatial alternation T (l, k),

h(l，k)＝Num{x|S(x)＝l，T*A(x)＝k}

Wherein, S (x)=l is illustrated in that pixel x corresponding gray is l in the target image,

It is k that T*A (x)=k is illustrated in through x corresponding gray in the collection of illustrative plates T*A pixel of spatial alternation;

Num is the number of pixel x;

Step (1.3). the joint probability distribution p between the collection of illustrative plates T*A behind calculating target image S and the process spatial alternation T _{S, T*A}(l, k):

$p_{S,T*A}(l,k)=\frac{h(l,k)}{\underset{l,k}{\mathrm{\&Sigma;}}h(l,k)}.$

Step (1.4). the marginal probability distribution p of difference computed image S _S(l) and the marginal probability distribution p of T*A _T*A(k):

$P_S\left(l\right)=\underset{k}{\mathrm{\&Sigma;}}{P}_{S,T*A}(l,k),$

$p_{T*A}\left(k\right)=\underset{l}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k).$

The marginal probability distribution is the notion in the Probability Statistics Theory, for two stochastic variable x, the random vector that y forms (x, y), the probability distribution of its component x is called random vector (x, the marginal distribution about x y); Here be the collection of illustrative plates T*A after image S and the conversion respectively as two stochastic variable l and k, the marginal probability that calculates l and k then respectively distributes;

Step (1.5). calculate the entropy H (S) of target image as follows respectively, the combination entropy H between collection of illustrative plates after the entropy H (T*A) of collection of illustrative plates and target image and the conversion after the conversion (S, T*A);

$H\left(S\right)=-\underset{l}{\mathrm{\&Sigma;}}{p}_S\left(l\right)\log p_S\left(l\right)$

$H(T*A)=-\underset{k}{\mathrm{\&Sigma;}}{p}_{T*A}\left(k\right)\log p_{T*A}\left(k\right)$

$H(S,T*A)=-\underset{k}{\mathrm{\&Sigma;}}{p}_{S,T*A}(l,k)\mathrm{lo}{g}_{S,T*A}(l,k)$

Step (1.6). the normalized mutual information between the collection of illustrative plates T*A after being calculated as follows target image S and passing through the T conversion:

$\mathrm{nMI}(S,T*A)=\frac{H\left(S\right)+H(T*A)}{H(S,T*A)}$

Step (1.7). (S T*A) brings the Powell optimized Algorithm into spatial alternation T is optimized iteration repeatedly, until till the convergence with normalized mutual information nMI;

The Powell optimized Algorithm is a kind of traditional determinacy optimization method, and is considered to the most effective excellent method that need not to differentiate at present; This basic idea is, for n dimension extreme-value problem, at first asks minimum along n coordinate direction, obtains n conjugate direction afterwards through n time, asks minimum along n conjugate direction then, through repeatedly just trying to achieve minimal value after the iteration;

, utilize the Powell optimized Algorithm that the value of spatial alternation T is optimized here, in the hope of obtaining mutual information nMI (S, minimal value T*A); Spatial alternation T contains nine components, and promptly (p, q, r, u, v, w, φ, ω, θ), nine dimensions in these nine corresponding Powell optimized Algorithm of components difference promptly utilize the Powell optimized Algorithm to ask nine dimension extreme-value problems; When Powell optimized Algorithm search end, promptly after the algorithm convergence, (S, the pairing spatial alternation value of minimal value T*A) is optimal spatial conversion T to the nMI that obtains ^*, T ^*=(p ^*, q ^*, r ^*, u ^*, v ^*, w ^*, φ ^*, ω ^*, θ ^*); T then ^** A is the registration collection of illustrative plates that whole registration will obtain;

Step (2). according to the following steps interested each organ is carried out independent registration successively, determine the roughly corresponding region of each interested organ in target image S; Specifically be divided into following a few step:

Step (2.1). determine the interval R of tonal range by following formula for each organ of interest _k

R _k＝[μ _k-λσ _k，μ _k+λσ _k]，

Wherein, μ _k, σ _kBe k the organ of interest of collection of illustrative plates T*A after conversion gray average and the standard variance in the occupied zone in target image S;

λ is a constant, is used for adjusting the uniform gray level scope; Test finds that the span of λ is 1.2～1.3;

Step (2.2). with the final spatial alternation T that obtains in the step (1) ^*Oppose and carry out the initial space conversion T ' of registration for each organ separately, even T '=T ^*=(p ^*, q ^*, r ^*, u ^*, v ^*, w ^*, φ ^*, ω ^*, θ ^*);

Step (2.3). be calculated as follows the similarity measure M of the organic region after the conversion _k(T '):

$M_k\left(T^{\&prime;}\right)=N_k^{\mathrm{in}}\left(T^{\&prime;}\right)-N_k^{\mathrm{out}}\left(T^{\&prime;}\right).$

Wherein, subscript k representative is carried out registration to k organ,

N _k ^InThe organ number of pixels that (T ') representative adopts spatial alternation T ' time organ k to comprise in the corresponding region in the registration collection of illustrative plates,

N _k ^Out(T ') represents the non-organ number of pixels that comprises in this corresponding region;

$\left\{\begin{array}{c}{N}_k^{\mathrm{in}}\left(T^{\&prime;}\right)=N\left\{x\right|x\&Element;C_k\left(T^{\&prime;}\right),S\left(x\right)\&Element;R_k\}\\ N_k^{\mathrm{out}}\left(T^{\&prime;}\right)=N\left\{x\right|x\&Element;C_k\left(T^{\&prime;}\right),S\left(x\right)\&NotElement;R_k\}\end{array}\right.$

Wherein, N{} represents the number of pixels that comprises in the set of computations,

C _kThe zone of (T ') expression organ k correspondence under spatial alternation T,

The gray-scale value of S (x) remarked pixel x in CT or MRI image;

Step (2.4). with the similarity measure M of the organic region after the conversion _k(T ') brings the Powell optimized Algorithm into T ' iterated, till convergence; Detailed process promptly utilizes the Powell optimized Algorithm that the value of spatial alternation T ' is optimized, in the hope of obtaining similarity measure M _kThe minimal value of (T '); Spatial alternation T ' contains nine components, and promptly (p ', q ', r ', u ', v ', w ', φ ', ω ', θ '), nine dimensions in these nine corresponding Powell optimized Algorithm of components difference promptly utilize the Powell optimized Algorithm to ask nine dimension extreme-value problems; When Powell optimized Algorithm search finishes, promptly after the algorithm convergence, the M that obtains _kThe pairing spatial alternation value of the minimal value of (T ') is optimal spatial conversion T ^*', T ' ^*=(p ' ^*, q ' ^*, r ' ^*, u ' ^*, v ' ^*, w ' ^*, φ ' ^*, ω ' ^*, θ ' ^*); C then _k(T ' ^*) organic region that will obtain for independent registration;

Step (3). fuzzy connection is cut apart; Specifically be divided into following a few step:

Step (3.1). establish C _k(T ' ^*) for the fuzzy prime area that k organ cut apart that connects, be calculated as follows the grey level histogram H in this zone _k(l):

H _k(l)＝N{x|x∈C _k，S(x)＝l}，l＝0，1……L

H _k(l) can be regarded as C _k(T ' ^*) gray scale is the number of the pixel x of l in the zone;

Step (3.2). adopt the method for curve fitting that histogram is carried out level and smooth and denoising, obtain the histogram after the match

${\hat{H}}_k\left(l\right)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_j\exp [-{\left(\frac{l-b_j}{c_j}\right)}^2],$

Wherein, n represents the number of gaussian component, i.e. the number of the organ that will cut apart; For example we cut apart liver and two kinds of organs of spleen, then n=2;

a _j, b _j, c _jRepresent peak value, center and the width of j Gaussian distribution successively;

Step (3.3). the histogram curve after the match can be regarded the grey level probability density function of organ in target image as; Select the initial seed point of at least one pixel of corresponding maximum gray probability as organ;

Step (3.4). investigate the gray probability sum of each seed points adjacent area; Adjacent area be defined as with as the collection of pixels of the seed points of center pixel distance less than the distance of a pixel unit; Described gray probability by have the number of pixels of this gray scale and the ratio of regional interior pixel sum in the research organic region; Have naming a person for a particular job of maximum region gray probability summation and finally be defined as seed points;

Step (3.5). with the histogram that obtains by step (3.2) The position of the probability density function of pairing organic region intensity profile and final seed points is an initial parameter, calculates the fuzzy strength of joint of each pixel and the final seed points of organ in the entire image with the fuzzy join algorithm of the described usefulness of following formula:

$c_{\mathrm{\&mu;}}(p,q)=\underset{\mathrm{\&rho;}(p,q)\&Element;p(p,q)}{\max }\left[\underset{z\&Element;\mathrm{\&rho;}(p,q)}{\min }{\hat{H}}_k\left(l\left(z\right)\right)\right]$

Wherein, p represents final seed points, and q represents in the entire image more arbitrarily,

c _μ(p, q) the fuzzy strength of joint between expression p and the q;

ρ (p, q) any communication path between expression p and the q,

P (p, the q) set of all communication paths between expression p and the q,

Z ∈ ρ (p, q) expression Z be path ρ (p, q) any pixel on,

L (z) is the gray-scale value of pixel Z in target image,

Expression gray-scale value l (z) is at histogram

In pairing probability density value;

This pixel that fuzzy strength of joint characterizes belongs to the degree of membership of organic region, and the degree of membership of each pixel is called fuzzy connection layout picture as the resulting image of gray scale;

Step (3.6). for the fuzzy connection layout picture of each organ is specified an optimal threshold; So-called threshold value is used for exactly and will belongs to the pixel of target organ and the numerical value that background pixel distinguishes, and the regulation gray scale is the pixel that belongs to target organ more than or equal to the pixel of this threshold value, and gray scale is a background pixel less than the pixel of this threshold value; So-called optimal threshold, be exactly thus the zone that assembled of the organ pixel that obtains of threshold value meet real organic region most; Optimal threshold is the important parameter that fuzzy connection layout is looked like to be converted into final segmentation result; Determine that the optimal threshold concrete steps are as follows:

At first set the threshold to 1, progressively reduce threshold value then; Be used for reducing each time the threshold value that obtains and cut apart fuzzy connection layout picture, obtain the cut zone of organ; Observe the change of the shape of cut zone with threshold value, if when a certain subthreshold reduces, shape acutely changes, and then stops to reduce threshold value, will go up a threshold value and be decided to be optimal threshold;

About the measurement way of the alteration of form of cut zone and definite method of optimal threshold, be further explained in detail as follows:

A). the way of the measurement of the alteration of form of cut zone:

Adopt area relative changing value and tight ness rating relative changing value to weighing the alteration of form of cut zone; Tight ness rating (Compactness) is defined as

Area relative changing value and tight ness rating relative changing value are defined as follows: reduce for each threshold value, calculate organ area and tight ness rating before and after reducing respectively; Organ area and tight ness rating before order reduces are original area and former tight ness rating, and organ area after the reduction and tight ness rating are new area and new tight ness rating, and then the relative variation of area and tight ness rating is defined as

B). determine the method for optimal threshold with area relative changing value and tight ness rating relative changing value:

In threshold value reduction process, if when certain once reduces, the relative changing value of area of discovery and tight ness rating thinks then that greater than predefined reference value ρ violent change has taken place shape, thus the outage threshold search, and the threshold value before will reducing is decided to be optimal threshold; Test shows that reference value ρ value 0.15～0.2 can obtain optimal threshold more accurately;

Step (3.7), the optimal threshold that obtains with step (3.6) according to following formula look like to cut apart to fuzzy connection layout, obtain the two-value split image of organ:

If the optimal threshold of k organ is T _k, then the binary segmentation result of this organ is

Wherein, x represents arbitrary pixel,

F _k(x) be the gray-scale value of the fuzzy connection layout picture of k organ at pixel x place,

B _k(x) be the gray-scale value of the bianry image after cutting apart at pixel x place;

Step (4). the organ shape correction, because causing because gray scale is close to blur to connect, two organs being close to mutually in the human abdominal cavity cut apart issuable erroneous segmentation result to solve;

In the medical image, also similar phenomena of the gray scale of two organs of being close to mutually in human abdominal cavity in image often appears; For two such organs, in the medical image on its surface mutually subsides place may not have tangible border, this can make fuzzy join algorithm think that these two organs are to grow together, thereby is divided into same organ mistakenly; Originally research and propose with organ shape and stabilize to according to detecting erroneous segmentation, and adopt range conversion the erroneous segmentation result to be carried out shape corrections in conjunction with the watershed segmentation algorithm; Idiographic flow is as follows:

Step (4.1). for the multilayer medical image, according to the following steps from top to bottom, adopt the method for zone-by-zone analysis to detect the aspect that has erroneous segmentation;

Step (4.1.1).

Set: when anterior layer is the aspect that algorithm is being analyzed, when the sequence number of anterior layer is i

The below adjacent layer is and works as the adjacent below aspect of anterior layer that the sequence number of below adjacent layer is i-1;

O _k ⁱBe the organ cut zone in anterior layer;

O _k ^I-1Be organ cut zone in the adjacent layer of below;

In the anterior layer, making the each point in the organ cut zone is negative value to the distance of zone boundary, for all do not belong to this regional point in the entire image, i.e. all points beyond should the zone, make its distance to the zone boundary be on the occasion of;

Step (4.1.2). to the organ cut zone O in anterior layer i _k ⁱCarry out range conversion, promptly calculate their minimum distance D to the zone boundary to whenever pressing following formula in the entire image _k ⁱ:

Wherein, y has represented in the pixel of zone boundary nearest with x;

B ⁱ _k(x) the final segmentation result B for obtaining in the step (3.7) _k ^*(x) at the pixel value in anterior layer i;

With distance value D _k ⁱ(x) compose to each pixel as gray-scale value, just obtained range image D _k ⁱ

Step (4.1.3). extract regional O in the adjacent layer of below _k ^I-1The border on the each point locations of pixels, then at D _k ⁱIn the D of each point at same position place _k ⁱ(x) value can be thought O _k ^I-1The distance value of borderline pixel; Obtain all O _k ^I-1The distance value of boundary pixel, calculate the standard variance of these distance values, to determine the similarity degree of organ between i layer and i-1 layer segmentation result; Detected segmentation errors if standard variance greater than predefined reference value, is then thought, following one deck is the aspect that has erroneous segmentation; Rule of thumb, this reference value generally is set in 5 pixel distances;

According to the method described above, detect structure at all levels from top to bottom, obtain the aspect that there is erroneous segmentation in all;

Step (4.2). by following formula the bianry image of each aspect of having erroneous segmentation is carried out range conversion;

Set:, make that its distance value is a negative value for all pixels in the organic region; For all do not belong to the zone point in addition in this zone outside the zone in the entire image, make its distance value be negative minimal value, i.e. negative infinite-∞;

Range image D then _k ⁱ' (x)

Wherein y has represented the frontier point nearest with x; B _k ⁱThe bianry image that promptly has the aspect of erroneous segmentation;

Step (4.3). the range image D that step (4.2) is obtained by the following method _k ⁱ' application watershed algorithm;

Step (4.3.1). entire image is set a gray threshold T; The initial value that makes T is image D _k ⁱ' in except negative minimum value infinite, that is:

$T_0=\min \left(\right\{D_k^i\left(x\right)|D_k^i\left(x\right)\&NotEqual;-\&infin;\}),$

Wherein, the coordinate of any pixel in the x representative image;

Step (4.3.2). with d is that step-length increases threshold value T gradually, promptly

T _n＝T ₀+n·d，

Wherein, T _nThe value of threshold value T when representing the n time increase step-length; N span 1≤n≤K, wherein $K=(\max \left(D_k^i\left(x\right)\right)-T_0)/d$ , max (D _k ⁱ(x)) presentation video D _k ⁱ' maximal value;

The d value is fixed, and rule of thumb gets d=3;

Step (4.3.3). utilizing increases the threshold value T that step-length obtains the n time _nPress following formula with image D _k ⁱ' be divided into bianry image W _n:

Wherein, the coordinate of any pixel in the x representative image;

W _n(x)=1 expression is submerged in the following part of the water surface, W _n(x)=0 expression is exposed at the above land part of the water surface; Increase step-length for preceding once (promptly the n-1 time) and obtain threshold value T _N-1, press following formula with image D _k ⁱ' be divided into bianry image W _N-1:

Wherein, the coordinate of any pixel in the x representative image;

Step (4.3.4). at T _nPropagation process in, whenever T _nGreater than a local maximum, with regard to key diagram as W _N-1In two connected regions are arranged at W _nIn be merged into a connected domain, then with W _N-1And W _nTwo width of cloth images carry out xor operation, obtain watershed divide image S _n:

Wherein, the coordinate of any pixel in the x representative image,  represents xor operation, S _n(x)=1 remarked pixel x is the position, watershed divide;

Step (4.3.5). to each the n value among n span 1≤n≤K, all can obtain a corresponding watershed divide image S _nIf all pairing watershed divide of n value images are S ₁, S ₂..., S _{Max (Dki (x))}, then the result images S of watershed transform is the set of all watershed divide images:

$S=S_1\&CenterDot;S_2\&CenterDot;...\&CenterDot;S_{\max \left(D_k^i\left(x\right)\right);}$

All pixel values are that 1 pixel is with regard to presentation video D among the image S _k ⁱ' in the position at place, watershed divide; By these watershed divide is image D _k ⁱ' be divided into plurality of sub-regions;

Step (4.4). the image D that watershed algorithm is obtained _k ⁱ' in each subregion, calculate its adjacent aspect organ cut zone O respectively with the below _k ^I-1Area overlap ratio:

Be lower than 50% zone and be considered as the erroneous segmentation zone overlapping ratio, be labeled as the background area again, thus with it at range image D _k ⁱ' in remove.

Images