CN1234091C - Medical image segmentation method based on horizontal collection and watershed method - Google Patents

Abstract

The present invention relates to a medical image segmentation method based on level set and a Watershed method. The present invention comprises the steps that anisotropic diffusion filtration is used for removing noise; images are excessively divided by the Watershed method; a stack data structure is built, which is used for positioning a grid point in a narrow band, which has the minimal time T; the images are finally divided by a Fast Marching method. In the present invention, medical image segmentation is carried out by the Watershed method and the improved Fast Marching method; the segmentation speed can be greatly improved, and simultaneously, the present invention has wide adaptability no matter CT images or MR images. The present invention has important application value for the field of computer-aided diagnosis and treatment, etc.

Description

Medical image cutting method based on horizontal collection and watershed method

Technical field

The present invention relates to Flame Image Process, particularly carry out the method for medical image segmentation in conjunction with Watershed (watershed divide) method and level set (Level set).

Background technology

In research and application to image, people are often only interested in some part in the image.Image segmentation utilizes exactly that Partial Feature extracts some interesting areas in the image in the image information, interesting areas is commonly referred to target or prospect (remainder is called background), could carry out feature extraction and measurement to target on this basis, thereby make more high-rise analysis and understanding become possibility.Image segmentation can with feature comprise gradation of image, color, texture, partial statistics characteristic or spectrum signature etc.

Image segmentation is by the committed step of Flame Image Process to graphical analysis, also is the subject matter in the computer vision field Level Visual, and it is again a classic problem simultaneously.Because the importance and the difficulty of this problem, image segmentation problem is attracting domestic and international researchist to make great efforts for it always for many years, has proposed thousands of various types of partitioning algorithms so far.But to the Partial Feature cut zone that can only utilize owing to image segmentation in the image information, therefore the whole bag of tricks must have limitation and specific aim, can only select the whole bag of tricks at the demand of various practical application area, up to the present also not have a general dividing method.

Level Set Method mainly progressively grows up from research fields such as interface propagation, and it is to handle the sealing moving interface effective computational tool of how much change in topology in the evolutionary process in time.Osher and Sethian at first propose the level set of the moving interface of the time that relies on and describe.The basic thought of Level Set Method is the level set that is expressed as two-dimentional toroidal function that the plane closed curve is implicit, and the point set that promptly has the same functions value is by the implicit motion of finding the solution curve of the evolution of level set function curved surface.Although it is complicated that this conversion problem that makes becomes in form, bring lot of advantages finding the solution of problem, its biggest advantage is that the change in topology of curve can access very natural processing, and can obtain unique weak solution that satisfies entropy condition.

Following fundamental equation is satisfied in the evolution of level set function:

Φ _t+F|Φ|＝0 (1)

Wherein Φ is a level set function, its zero level set representations objective contour curve, that is:

Γ(t)＝{x|Φ(x，t)＝0} (2)

| Φ | the gradient norm of expression level set function; F is the velocity function on the surface normal direction, the motion of control curve, and general F comprises two: with image-related item (as gradient information) and the item (as the curvature of curve) relevant with the geometric configuration of curve.

Suppose that Γ is the initial position of a lineoid, F is the speed on normal direction.In the LevelSet algorithm, we look Γ is a higher-dimension function phi (x, y, zero level collection z).Then, according to the chain type rule, can produce EVOLUTION EQUATION suc as formula the lineoid of the motion of (1).

In the level set algorithm prevailing a kind of be exactly so-called Fast Marching (advancing fast) method.Consider the special circumstances of a kind of interface motion now, that is exactly movement velocity F＞0 at interface.Suppose T be the interface through a specified point (x, time y), T just satisfied following equation like this:

|T|F＝1 (3)

This formula understands that briefly the gradient of time of arrival and the movement velocity at interface are inversely proportional to.In a broad sense, have two kinds of methods can be used for being similar to the time dependent position, interface of motion: a kind of is to solve by difference quotient in iteration and the digital approximation formula (1); Another kind is the solution that makes up formula (3) middle time of arrival of T.And Fast Marching method depends on a kind of method in back.

Formula (3) is a kind of form of famous Eikonal equation, and Sethian at first proposes, and obtain T time of arrival in the formula (3), is equivalent to and finds the solution following quadratic equation.

${\left[\begin{array}{c}\max {(D_{\mathrm{ij}}^{-x}T,0)}^2+\min {(D_{\mathrm{ij}}^{+x}T,0)}^2+\\ \max {(D_{\mathrm{ij}}^{-y}T,0)}^2+\min {(D_{\mathrm{ij}}^{+y}T,0)}^2+\end{array}\right]}^{1/2}=1/F_{i,j}---\left(4\right)$

Here D ^-And D ⁺Be respectively backward difference and forward-difference operator.

Below we provide the concrete steps of the Fast Marching method that Sethian at first proposes:

Initialization:

1) moving point: moving point is exactly the fixing point of time T in all net points.In our algorithm, the seed points of user's appointment just, (x y)=0, sees Fig. 2 to time T, is example with the two-dimensional case;

2) arrowband point: all points in the arrowband are called the arrowband point.In our algorithm, the point of the 4-adjacency of all seed points just, time T (x, y)=(x y), sees Fig. 2 to 1/F, is example with the two-dimensional case;

3) point of distance: except moving point and arrowband point, all other net point is a point of distance, and T (x, y)=TIME_MAX, see Fig. 2, with the two-dimensional case example;

Red grid among Fig. 2 is represented moving point, and green grid is represented the arrowband point, and the grid of white is represented point of distance.The beginning of Fig. 2 (a) expression Fast Marching method, the process of Fig. 2 (b) and Fig. 2 (c) expression Fast Marching method.

Circulation:

1) begin circulation: setting up an office, (i j) is the point that has minimum time T in the arrowband;

2) (i is moving point j), and deletes from the arrowband identification point;

3) identification point (i, 4-abutment points j): as fruit dot (i, abutment points j) is a moving point, does not then change the time; If its abutment points is the arrowband point, then upgrade the time of abutment points according to formula (4); If its abutment points is a point of distance, then identifying this abutment points is the arrowband point, simultaneously the time of upgrading this abutment points according to formula (4);

4) if surpass specified threshold value the time of arrival of a bit, end loop then, otherwise jump to 1);

We readily appreciate that in the worst case, Fast Marching method need travel through all pixels of piece image from the concrete steps of Fast Marching method.R.Malladi has at first defined following velocity function:

$F(x,y)=e^{-\mathrm{\α}|\▿G_{\mathrm{\σ}}*I(x,y)|}---\left(5\right)$

Here (x y) is original image to I, and is a gradient operator, G _σ ^{* I (x, y)}It is level and smooth that expression is carried out Gauss to original image, and 0＜α＜1st, weight coefficient.

Summary of the invention

The purpose of this invention is to provide a kind of algorithm of medical image segmentation fast and accurately, can cut apart the image of mode such as CT and MR.

For achieving the above object, the medical image cutting method based on horizontal collection and watershed method comprises step:

Anisotropic diffusion filtering is used to remove noise;

Adopt the Watershed method, image is carried out over-segmentation;

The foundation of heap data structure is used for locating the net point that the arrowband has minimum time T;

Adopt Fast Marching method, image is carried out last cutting apart.

The present invention utilizes Watershed method and improved Fast Marching method to carry out medical image segmentation, not only can improve the speed of cutting apart widely, also has adaptability widely simultaneously, no matter is CT or MR image.There is important use to be worth in fields such as computer-aided diagnosis and treatments.

Description of drawings

Fig. 1 utilizes Watershed method and improvement Fast marching method to carry out the structured flowchart of image segmentation;

Fig. 2 is activity, arrowband and point of distance synoptic diagram;

Fig. 3 is that Fast Marching method is crossed over the image boundary synoptic diagram;

Fig. 4 is an adjacent area synoptic diagram after the Watershed conversion;

Fig. 5 is a knee joint CT image segmentation;

Fig. 6 is a brain MR image: gyrus is cut apart;

Fig. 7 is a brain MR image: white matter of brain is cut apart;

Fig. 8 is a brain MR image: white matter of brain is cut apart.

Embodiment

Core concept of the present invention is Watershed algorithm and improved Fast Marching method to be combined carry out image segmentation.Owing to introduced the Watershed algorithm image carried out over-segmentation, this method only need be calculated the time of arrival that seed points arrives the zone boundary at its place, for other pixel of intra-zone, then need not the time of arrival of seed points to calculate, the arithmetic speed of algorithm just is much improved like this.And we have redefined the velocity function of Fast Marching method according to the similarity of the statistical property between the zone.

Describe medical image segmentation algorithm of the present invention in detail below in conjunction with accompanying drawing.As a kind of concrete implementation, structured flowchart is seen Fig. 1.Mainly comprise four steps: anisotropic diffusion filtering, the Watershed method is carried out over-segmentation (number of regions of cutting apart surpasses the practical object number that comprises in the image), the foundation of heap data structure, improved Fast Marching method is carried out last cutting apart.Below it is made introductions all round.

Step 1: anisotropic diffusion filtering

Medical image contains random noise usually, and noise might cause the rim detection mistake, and for reducing the influence of random noise to partitioning algorithm, we can carry out filtering to image earlier.But, also can smoothly fall the edge usually when utilizing low-pass filtering operator smoothed image noise, will cause boundary information destroyed like this, thereby bring difficulty for the aftertreatment of image.For avoiding the generation of this situation, in our algorithm, will adopt anisotropic diffusion filtering to remove noise, this method can also keep marginal position and strengthen the edge when removing noise.

Anisotropic diffusion filtering is proposed by Perona and Malik the earliest, the anisotropic diffusion equation below they have proposed on traditional isotropy diffusion equation basis, and use it for the edge of image enhancing.

$\frac{\∂I(x,y,t)}{\∂t}=\mathrm{div}[g\left(\right||\▿I)\▿I]---\left(6\right)$

Wherein | I| presentation video gradient magnitude, coefficient of diffusion are function gs () relevant with gradient.G () is a monotonic decreasing function that is used for controlling level and smooth degree, satisfies g (0)=1, and when x → ∞ g (x) → 0.This means that coefficient of diffusion reduces with the increase of image gradient, can guarantee like this to spread with fast speed at intra-zone ( I is little), in then no longer diffusion of marginal point ( I is big), thereby play the effect that the edge strengthens, so g () is also referred to as the edge and stops function.

Step 2: the Watershed method is carried out over-segmentation

Watershed (watershed divide) algorithm is a kind of method based on mathematical morphology, is proposed by Digabel and Lantuejou the earliest.The Watershed algorithm also can be categorized as the dividing method based on the zone.The thought source of Watershed method is in geography: it regards the gradient magnitude image as a width of cloth topomap, and the corresponding sea level elevation of gradient magnitude, the zone of different Grad is just corresponding to basin between mountain peak and mountain valley in the image.Imagination is made a call to a hole in the position of each local minizing point, then topomap is immersed gradually in the lake, and the basin of global minimum point is water inlet earlier.Water level raises gradually and covers the basin, when the water in adjacent two basins is about to merge, at this moment builds the dam interception between two basins.This process is basin, many mountain valleys with image division, and the watershed divide is exactly dykes and dams of separating these basins.The basic step of Watershed algorithm is as follows:

1) nonlinear diffusion filtering smoothed image;

2) calculate the smoothly gradient map of back image;

3) the maximal value MAX of compute gradient figure;

4) select suitable water logging threshold value DT=MAX*1/n;

5) to each pixel, if current pixel and its 8-neighbours' value all less than DT, then merges current pixel and its 8-neighbours

6) merged zonule (as: counting of comprising of zone＜10).

Watershed algorithm will produce the over-segmentation figure (number of regions of cutting apart surpasses the practical object number that comprises in the image) of original image, and the number of regions of over-segmentation depends on the selection of flooding elevation of water.

Step 3: the foundation of heap data structure

The crucial part of Fast Marching method is to locate how fast the net point that has minimum time T in the arrowband.Here we use following heap data structure:

struct

{

Double time; Seed points is to the time of arrival of net point

Int x; The directions X coordinate of net point

Int y; The Y direction coordinate of net point

}HeapNode；

We store the time of arrival T of seed points to other net point, preserve the position (coordinate of the coordinate of directions X and Y direction) of this net point in network simultaneously.Be some operations below to the heap data structure:

Heap_init: initialization heap data structure, seed points as moving point, is calculated the time of arrival that seed points arrives marginal point, and by depositing in proper order from small to large in the heap data structure;

Heap_push: point of distance is converted into the arrowband point, and by depositing in proper order from small to large in the heap data structure;

Heap_pop: deletion has the net point of minimum time from heap, and changes it into moving point;

Heap_update: the time of arrival of upgrading seed points net point in the arrowband.

Step 4: improved Fast Marching method is carried out last cutting apart

From formula (5) as can be known, (x y) only relies on boundary information in the image to velocity function F, just the gradient information of image and do not utilize the global information of image-region fully.Therefore, for obscurity boundary in the image or discrete shape, just be difficult to obtain desirable segmentation effect.As shown in Figure 3:

Fig. 3 (a) expression knee image original graph, the green arrow indication be respectively that seed points is selected and the place of obscurity boundary; Fig. 3 (b) is illustrated in the fuzzy local segmentation result of image boundary and has crossed over the border, and this is closely bound up with the selection of the definition of velocity function and parameter alpha, and the key of Fast Marching method is choosing of velocity function.Because its velocity function has been chosen the gradient information of image, the boundary information of image just, and do not utilize the global information of image-region fully must cause at gradation of image more approachingly, and the local segmentation result of obscurity boundary makes a mistake.

We know that in general, image is to be made of Pork-pieces zonule, and every zonule all has identical character, and is more approaching as gray-scale value, and texture structure is similar etc.How effectively to utilize the information of these adjacent zonules most important to last segmentation result.For this reason, we introduce the Watershed conversion, earlier image over-segmentation (number of regions of cutting apart surpasses the practical object number that comprises in the image) are become many such zonules.

Image is carried out the benefit that over-segmentation has three aspects.The one, for Fast Marching method, suppose that seed points (moving point) is positioned at the inside of a zonule, we only need to calculate the time of arrival that seed points arrives its boundary pixel point, other pixel of inside, zonule hereto, because it has identical character (in our algorithm, it is more approaching to refer to gray values of pixel points), so seed points need not to calculate to time of arrival of these pixels, the speed of Fast Marching method improves greatly like this; The 2nd, because the border of these zonules might be potential correct segmentation result after the over-segmentation, so arranging the border one of last segmentation result, we are positioned on the border in the zone that obtains by over-segmentation, thereby reduced the blindness of cutting apart, improved the accuracy of cutting apart; The 3rd, the similarity of the partial statistics characteristic of these adjacent zonules has good reference value for the velocity function of Fast Marching method.

The zone that note Watershed transfer pair image carries out after the over-segmentation is R ₁, R ₂, R ₃..., R _n, n is the number of regions after the over-segmentation.Each region R _iStatistical nature be designated as: { SF _{I, j}| 1≤i≤n, 1≤j≤m}, wherein m represents the number of range statistics feature.SF _{I, 1}Be the gray average of pixel in this zone, SF _{I, 2}It is the gray variance of pixel in this zone.The number of the feature of record can decide according to the characteristic of image, and tangible texture is for example arranged in the image, and can concentrate in provincial characteristics increases the feature of distinguishing texture.In our experiment, only use SF _{I, 1}And SF _{I, 2}Two features.

That shown in Figure 4 is latter two adjacent areas of Watershed conversion R ₁And R ₂, seed points (moving point) s ₁Be located at region R ₁In.What red line was represented among the figure is region R ₁And R ₂Common border, p ₁, p ₂, p ₃... be borderline pixel, and q ₁, q ₂, q ₃... be region R ₁Other borderline pixel.

As seed points s ₁When arriving the zone boundary, we can calculate region R ₁And R ₂Statistical nature SF _{1, j}And SF _{2, j}, and be calculated as follows out similarity between them:

${\mathrm{Sm}}_{i,j}=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}|S{F}_{i,k}-{\mathrm{SF}}_{j,k}|---\left(7\right)$

Here Sm _{I, j}The expression region R _iAnd region R _jBetween similarity, m represents region R _iAnd region R _jThe statistical nature number, SF _{I, m}The expression region R _iStatistical nature, SF _{J, m}The expression region R _jStatistical nature.Obtain region R _iAnd region R _jBetween similarity, we just can defined range R _iAnd region R _jThe velocity function of pixel on the common border:

$F_{i,j}=e^{-\mathrm{\β}\left|{\mathrm{Sm}}_{i,j}\right|}---\left(8\right)$

Here be the positive number of 0≤β≤1.

In sum, the step of the Fast Marching method after the improvement is as follows:

1. initialization:

1) moving point: the interior pixels point in the zone at seed points of identifying user appointment (can have a plurality of) and seed points place is a moving point, and time T (x, y)=0;

2) arrowband point: the boundary pixel point that identifies all seed points regions is the arrowband point.Time T (x, y)=1/F (x, y);

3) point of distance: except moving point and arrowband point, all other net point is a point of distance, and T (x, y)=TIME_MAX;

2. circulation:

1) begin circulation: setting up an office, (i is the point that has minimum time T in the arrowband j), and the zone that might as well establish its place is R _i, R _j

2) zoning R _iAnd R _jSimilarity, (i j) is moving point to identification point, and deletes from the arrowband;

3) identification point (i, all frontier points of region j): (i, the frontier point of region j) they are moving point, then do not change the time as fruit dot; If the frontier point of its region is the arrowband point, then upgrade the time of frontier point according to formula (4); If the frontier point of its region is a point of distance, then identifying this frontier point is the arrowband point, simultaneously the time of upgrading this frontier point according to formula (4);

4) if surpass specified threshold value the time of arrival of a bit, end loop then, otherwise jump to 1);

Operation result

We have realized algorithm described in the invention with C Plus Plus, and have done experiment on several different data sets.All experiments all are at a PIII 800, the 128MB internal memory, and operating system is to finish on the PC of Windows 2000.

Fig. 5,6,7,8 is several groups of examples.Fig. 5 (a) is original CT knee joint image, Fig. 5 (b) utilizes the transition of Watershed method to cut apart figure, Fig. 5 (c) be utilize improved FastMarching method cut apart figure at last, think ratio with Fig. 3, the Fast Marching method after the improvement can obtain satisfied segmentation result in the place of obscurity boundary.Fig. 6 (a) is a brain MR original image, and Fig. 6 (b) is Watershed over-segmentation figure, and Fig. 6 (c) is last gyrus segmentation result.Fig. 7 (a) is a brain MR original image, and Fig. 7 (b) is Watershed over-segmentation figure, and Fig. 7 (c) is that last white matter of brain is cut apart figure.Fig. 8 is as Fig. 7.

These several groups of examples have fully proved the validity of our algorithms.

Claims (3)

1. medical image cutting method based on horizontal collection and watershed Watershed method comprises step:

Original image is carried out anisotropic diffusion filtering, be used to remove the noise of image, safeguard image edge information simultaneously;

Adopt watershed divide Watershed method, image is carried out over-segmentation;

The foundation of heap data structure, nodes records in the pile structure be T time of arrival, described time of arrival, T was the time of arrival that seed points arrives other net point, and T (x, y)=1/F (x, y), wherein (x, y) expression net point coordinate, F (x, y) be velocity function, with the net point that has minimum time T in the heap data structural orientation arrowband;

Adopt improved fast moving Fast Marching method, image is carried out last cutting apart, its step is as follows;

1. initialization:

1) moving point: the interior pixels point in the zone at the seed points of identifying user appointment and seed points place is a moving point;

2) arrowband point: the boundary pixel point that identifies all seed points regions is the arrowband point;

3) point of distance: except moving point and arrowband point, all other net point is a point of distance;

2. circulation:

1) begin circulation: setting up an office, (i j) is the point that has minimum time T in the arrowband;

2) region R _iAnd R _jBe point (i, j) region, zoning R _iAnd R _jSimilarity, (i j) is moving point to identification point, and deletes from the arrowband;

3) (i, j) all frontier points of region: (i, the frontier point of region j) are moving point to identification point, then do not change the time as fruit dot; If the frontier point of its region is the arrowband point, upgrade the time of frontier point; If the frontier point of its region is a point of distance, then identifying this frontier point is the arrowband point, the time of upgrading this frontier point;

4) if seed points surpasses specified threshold value to certain any time of arrival, end loop then, otherwise restart circulation;

2. by the described method of claim 1, it is characterized in that the time of described renewal frontier point is undertaken by following formula:

${[\max {(D_{\mathrm{ij}}^{-x}T,0)}^2+\min {(D_{\mathrm{ij}}^{+x}T,0)}^2+\max {(D_{\mathrm{ij}}^{-y}T,0)}^2+\min {(D_{\mathrm{ij}}^{+y}T,0)}^2]}^{1/2}=1/F_{\mathrm{ij}}$

Wherein, F _i, j be point (i, initial velocity function j), according to $F_{i,j}=e^{-\mathrm{\β}\left[{\mathrm{Sm}}_{I,J}\right|}$ Calculate, β is a constant, Sm _{I, j}Be region R _iAnd R _jSimilarity, its calculating is carried out according to following formula:

${\mathrm{Sm}}_{i,j}=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}|S{F}_{i,k}-{\mathrm{SF}}_{j,k}|$

Here D ^-And D ⁺Be respectively backward difference and forward-difference operator.

3. by the described method of claim 1, it is characterized in that described zoning R _iAnd R _jSimilarity undertaken by following formula:

$S{m}_{i,j}=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}|S{F}_{i,k}-{\mathrm{SF}}_{j,k}|$

Wherein, Sm _{I, j}The expression region R _iAnd region R _jBetween similarity, m represents region R _iAnd region R _jThe statistical nature number, SF _{I, m}The expression region R _iStatistical nature, SF _{J, m}The expression region R _jStatistical nature.

Images