CN103886603B - Left ventricle nuclear magnetic resonance image segmentation and three-dimensional reconstruction method - Google Patents

Abstract

The invention provides a left ventricle nuclear magnetic resonance image segmentation and three-dimensional reconstruction method. According to the method, a variational level set evolution model is established to be used for segmenting the boundary of the inner membrane and the outer membrane of the left ventricle, and meanwhile, a three-dimensional curved surface and entity reconstruction scheme of the left ventricle and a corresponding mesh generation scheme are also provided. Due to the fact that a novel quasi local binary simulation item based on an edge detection operator is adopted, the method can be used for effectively segmenting images which are uneven in gray level; only the area item is considered, therefore, theoretical models involved in the method can be simplified, and the algorithm operation time of the method can also be shortened; the application of the convex hull algorithm enables the method to be capable of making up the influences such as leakage caused by fuzziness of the boundary of the outer membrane. According to the method, the three-dimensional reconstruction scheme based on the contour surface generation algorithm and secondary development and application technologies of relevant commercial software is further provided, and can help a user realize three-dimensional reconstruction of the left ventricle.

Description

A kind of Left Ventricle MR Images segmentation and the method for three-dimensionalreconstruction

Technical field

The present invention relates to a kind of image processing field, the method for especially a kind of segmentation of nuclear magnetic resonance image and reconstruct.

Background technology

Heart nuclear magnetic resonance CMR (Cardiac Magnetic Resonance) imaging can provide high-resolution figure Picture, is one of the study hotspot in medical image analysis field, is also the important supplementary means of heart disease diagnosis.In order to fully sharp With the information in CMR image, provide quantization, intuitively reference for clinical diagnosis, the accurate segmentation of the inside and outside film of left room film is first Want task.Left ventricle image segmentation problem mainly faces following challenge：(1) heterogeneity of gradation of image；(2) epicardial border is not Clearly；(3) the inside and outside membrane boundary shape in different sectioning images is changeable etc..Although it is existing at present many inside and outside with regard to left ventricle The cutting techniques of membrane boundary, but how accurately to split its inside and outside membrane boundary and remain a major challenge.

Level set (Level Set) method can implicitly follow the trail of the border of target：First solve given partial differential equation In level set function, then by zero isopleth of level set function come implicit expression follow the tracks of target.At present, Level Set Method is extensive Be applied to the fields such as medical image segmentation, wherein more famous method surely belongs to the Active contour models method of Chan and the nothing of Li Reinitialize method.But, the Active contour models method of Chan can not process the image with non-homogeneous gray scale well, and The gradient terms that the nothing of Li reinitializes in the corresponding Partial Differential Equation Model of method are more, this undoubtedly increased the complexity of model with And amount of calculation during numerical discretization solution.

Content of the invention

In order to overcome the deficiencies in the prior art, the present invention provides a kind of new variation level diversity method, for automatic, fast The inside and outside membrane boundary of fast, accurate Ground Split left ventricle.Meanwhile, in order to meet the demand of clinical practice and mechanical analysis, also provide The three-dimension curved surface of left ventricle, entity reconstruct and corresponding mesh generation scheme.

The technical solution adopted for the present invention to solve the technical problems is to be broadly divided into segmentation module and three-dimensionalreconstruction module. Segmentation module is the core of the present invention, mainly employs the image Segmentation Technology of new variation level diversity method, and three-dimensionalreconstruction Module is then Visualized Post Processing supplementary module.

The present invention comprises the following steps that：

First, initialization step, i.e. the input of CMR image：First by the CMR Image semantic classification of dicom form become png, bmp, The picture format of jpg or jpeg, pretreated image is read into scientific and engineering meter by general built-in imread order Calculate in software Matlab；

2nd, segmentation step：

1. the foundation of model：After the completion of initialization step, first set up and intend local binary simulated items, area item and penalty term Three, described plan local binary simulated items, area item and penalty term three are added and just constitute gross energy functional, recycles General variation extreme value analysis method, can get variation level set evolutionary model, this model of iterative, until meeting following two Any one in individual condition, then stop iteration：

1) gross energy functional, that is, intend local binary simulated items, three summations being added of area item and penalty term are less than 10^-6；

2) iterative steps reach the upper limit 200；

2. initial value is chosen：In the image that Matlab reads in initialization step, arbitrarily choose initial level set function Zero isopleth；

3. parameter setting：Corresponding model parameter in the variation level set evolutionary model of the setting present invention；

3rd, the convex closure result step of output segmentation：Based on the result of segmentation step, available algorithm of convex hull obtains corresponding Convex closure result；

4th, obtain the convex closure result of inside and outside film segmentation：The CMR sectioning image of left ventricle is directed respectively into segmentation module, can Obtain the convex closure result of a series of inside and outside film image segmentation of this left ventricle；

5th, three-dimension curved surface, entity reconstruct：Using contour surface generating algorithm, the three-dimension curved surface reconstruct knot of inside and outside film can be obtained Really；The plane convex closure result of inside and outside film segmentation is imported the macro document in three-dimensional Cartography designing software SolidWorks, then generates 3D solid reconstruction result based on Spline Interpolating Surfaces；

6th, 3D solid mesh generation step：The 3D solid reconstruction result that SolidWorks is generated imports finite element Analysis software Abaqus, obtains 3D solid mesh generation by calling Mesh functional module, that is, obtain corresponding node and sit Mark, the information of node adjacency chained list, unit and its composition node.

The establishment step of the variation level set evolutionary model described in segmentation step of the present invention is as follows：

I sets up and intends local binary simulated items, for processing gray scale non-uniform image：

Gray value is taken to be I₀Image in any one initial open region D, its border isAnd D_inAnd D_out It is respectively the inside and outside region of open region D, introduce level set function φ (x), meet equation below：

Introduce edge detection operator g (x), construct new plan local binary simulated items F_New：

F_New=α ∫ g (x) | I₀(x)-c₁|²H(φ(x))dx+β∫g(x)|I₀(x)-c₂|²(1-H(φ(x)))dx (2)

c₁=∫ g (x) I₀(x)H(φ(x))dx/∫H(φ(x))dx (3)

c₂=∫ g (x) I₀(x)(1-H(φ(x)))dx/∫(1-H(φ(x)))dx (4)

Wherein α, β are model parameter, c₁、c₂For two value parameters, I₀X () is image intensity value at x, H () is Heaviside function, g (x) is then following edge detection operator

Wherein K_σ*I₀It is initial pictures I₀With Gaussian coreConvolution, parameter σ=1.5, n=10；

II introduces geometry item and penalty term

Introduce geometry item F_Geom：

F_Geom=γ ∫ g (x) H (φ (x)) dx (6)

Using following weighting penalty term F_P：

Wherein, γ, λ are weighted factor；

III constructs gross energy functional F_T, see below formula：

IV sets up new variation level set evolutionary model：

Preset parameter α, β, γ and λ, solve above-mentioned gross energy functional F_TVariation extreme value with regard to level set function φ (x) is asked Topic, derives the Euler-Lagrange equation with regard to φ (x)：

Wherein, δ () is Dirac Delta function, formula (9) the variation level set evolution mould that i.e. present invention is set up Type, can the inside and outside film of left ventricle effectively in Ground Split CMR image.

The invention has the beneficial effects as follows：Due to employing the new plan local binary simulated items based on edge detection operator, The present invention is enable effectively to split gray scale non-uniform image；Only consider area item, theoretical mould involved in the present invention can be simplified Type, can reduce the Riming time of algorithm of the present invention again；The application of algorithm of convex hull allows the present invention to repair due to epicardial border Smudgy leakage causing etc. affects.The present invention gives secondary based on contour surface generating algorithm and relative commercial software The three-dimensionalreconstruction scheme of development and application technology, can help user to realize the three-dimensionalreconstruction of left ventricle.

Brief description

Fig. 1 is the flow chart of initialization and segmentation module.

Fig. 2 is the flow chart of three-dimensionalreconstruction module.

Fig. 3 is the selection of initial value, and left figure is original image, and right figure white line is current initial zero isopleth chosen.

Fig. 4 is inner membrance segmentation result, and left and right in figure straight line is original segmentation result, and in right figure, the straight line with No. * is segmentation Convex closure result.

Fig. 5 is adventitia segmentation result：Left and right in figure straight line is original segmentation result, and in right figure, the straight line with No. * is segmentation Convex closure result.

Fig. 6 is a series of corresponding inside and outside film segmentation result of CMR images.

Fig. 7 is the reconstruct image being obtained by contour surface generating algorithm, and left figure is left ventricle inner membrance surface reconstruction, and right figure is adventitia Surface reconstruction.

Fig. 8 is three-dimensionalreconstruction result, and left figure reconstructs for SolidWorks 3D solid, and right figure is Abaqus 3D solid net Lattice subdivision.

Specific embodiment

The present invention is further described with reference to the accompanying drawings and examples.

Fig. 1 is the flow chart of initialization and segmentation module, and Fig. 2 is the flow chart of three-dimensionalreconstruction module, concrete steps explanation As follows：

1st, " initialize " step：

" image input "：Using RadiAnt DICOM Viewer software, the CMR image of dicom form is changed into png lattice The image of formula, reads in this image afterwards.Fig. 3 left figure is the png image after reading in.

2nd, " split " step：

1. first set up new variation level set model

I sets up and intends local binary simulated items, for processing gray scale non-uniform image

Consideration gray value is I₀Image in any one initial open region D, its border isAnd D_inWith D_outInside and outside region for D.In order to implicit expression follows the trail of object boundary, level set function φ (x) can be introduced, meet equation below

In the present invention, so that method can be applied to the segmentation of non-uniform image, it has been specifically incorporated rim detection Operator, constructs new plan local binary simulated items F_New：

F_New=α ∫ g (x) | I₀(x)-c₁|²H(φ(x))dx+β∫g(x)|I₀(x)-c₂|²(1-H(φ(x)))dx (2)

c₁=∫ g (x) I₀(x)H(φ(x))dx/∫H(φ(x))dx (3)

c₂=∫ g (x) I₀(x)(1-H(φ(x)))dx/∫(1-H(φ(x)))dx (4)

Wherein α, β are model parameter, c₁、c₂For two value parameters, I₀X () is image intensity value at x, H () is Heaviside function, g (x) is then following edge detection operator

Wherein K_σ*I₀It is initial pictures I₀With Gaussian coreConvolution, parameter σ=1.5, n=10.

II introduces geometry item and penalty term, for reducing Riming time of algorithm

For the motion of acceleration level set function zero isopleth, reduce Riming time of algorithm, need using following geometry item F_Geom：

F_Geom=γ ∫ g (x) H (φ (x)) dx (6)

When weighted factor γ be on the occasion of when, zero isopleth of φ (x) can be made to gradually taper up；Conversely, this activity zero is equivalent Line will be expanded.Geometry item (6) is only related to area, does not utilize other geometry items, such as length item, decreases part partial derivative Amount of calculation, also simplify final model simultaneously.

Finally, reinitialize, in order to eliminate, the time-consuming impact that level set function brings, also can be using following weighting punishment Item F_P：

Wherein, λ is weighted factor.

III constructs gross energy functional F_T, it is used for setting up final mask and whether evaluation algorithm terminates

Consider the impact of (2), (6) and (7), finally can be defined as follows gross energy functional F_T：

In formula (8), F_New+F_GeomLevel set function zero isopleth can be ordered about gradually move to object boundary, and F_P Then being to try to maintenance level set function is symbolic measurement, thus avoiding reinitializing of level set function.It is based on (8), following variation level set evolutionary model can be derived；In addition, when gross energy functional F in (8)_TValue be less than threshold values 10^-6When, Iterative algorithm can be terminated, complete cutting procedure.

IV finally sets up new variation level set evolutionary model, is used for realizing dynamic partition process

Preset parameter α, β, γ and λ, solve above-mentioned gross energy functional F_TVariation extreme value with regard to level set function φ (x) is asked Topic, can derive the Euler-Lagrange equation with regard to φ (x)：

Wherein, δ () is Dirac Delta function, and model (9) is exactly that the new variation level set used in the present invention is drilled Change model, can the inside and outside film of left ventricle effectively in Ground Split CMR image.

2. " initial value selection "：In order to quickly obtain left ventricle inside and outside film segmentation result it is proposed that choosing between inside and outside film Initial level set function zero isopleth, as shown in Fig. 3 right figure.For experienced user, initial value might as well be generally disposed at Near inside and outside film, so then can faster obtain and be satisfied with segmentation result.

Table 1 parameter declaration and its span

3. " parameter setting "：The parameter of setting variation level set model, parameter declaration and suggestion span are shown in Table 1.Under In mask body implementation process, when splitting left ventricle inner membrance, have chosen following parameter value：Δ t=5, Loops=200, α= 0.001, β=0.001, γ=1.5, λ Δ t=0.2；When splitting left ventricular epicardium it is only necessary to adjust γ=- 1.5, other specification Value is constant.

4. iterative variation level set evolutionary model (9) again

Partial derivative for the direction in space being related in model (9)WithEmploy Second-Order Central Difference near Seemingly；Partial derivative for time orientationSingle order forward difference is then adopted approximate.The finite difference of solution (9) thus can be arrived Form, the increase with time step number will form iterative process.

5. following criterion is utilized to control iterative process

Iterative variation level set evolutionary model, until gross energy functional F_T<10^-6Or iterative steps<Greatest iteration Step number Loops=200.Based on this criterion, dynamic zero isopleth by automatic, rapid desufflation or is expanded to the inside and outside film of left ventricle Near, complete cutting procedure.The segmentation result of inside and outside film asks for an interview the straight line in Fig. 4, Fig. 5.

3rd, " the convex closure result of output segmentation " step：

Result Ji Yu " segmentation " step, the algorithm of convex hull of available Preparata and Hong (1977) is acted upon, energy Obtain more rational convex closure result, as shown in the straight line with No. * in Fig. 4, Fig. 5.Especially, because epicardial border obscures not Clearly, this can cause algorithm can not converge to epicardial border exactly, and the application of algorithm of convex hull then can eliminate due to adventitia significantly The unclear impact causing of obscurity boundary, can be referring to the right figure in Fig. 5：Convex closure result (straight line with No. *) is than original segmentation result (straight line) is more accurate.

4th, " the convex closure result of inside and outside film segmentation " step：

After the CMR sectioning image of left ventricle is directed respectively into above-mentioned segmentation module, just can obtain the one of this left ventricle is Arrange the convex closure result (closed curve) of inside and outside film image segmentation, Fig. 6 can be refer to.

5th, " three-dimension curved surface, entity reconstruct " step：

The convex closure result of the inside and outside film image segmentation 1. being obtained based on step 4, using Hansen and Johnson (2004) Contour surface generating algorithm, the three-dimension curved surface of the inside and outside film of left ventricle can be generated, refer to Fig. 7.Can see, at this moment generate Inside and outside face is rough.Next visual effect more preferably entity reconstruction result will be generated based on SolidWorks.

2. the 3D solid based on Spline Interpolating Surfaces reconstructs：The convex closure result of inside and outside film segmentation is imported SolidWorks macro document, then run this macro document in SolidWorks, you can insertion closing convex closure curve.Repeat this behaviour Make, a series of left ventricle inside and outside film outline line of closings can be obtained；Next, being inserted using the spline surface in SolidWorks Value function and hypostazation function, can obtain the ideal three-dimensional entity reconstruction result as shown in Fig. 8 left figure.

6th, " 3D solid mesh generation " step：

First the 3D solid reconstruction result that SolidWorks generates is derived and be saved as the file of x_t form, afterwards by it Import Abaqus, 3D solid mesh generation result be can get by the Mesh functional module calling Abaqus, such as Fig. 8 right figure Shown tetrahedron mesh generation.

Above-mentioned embodiment is merely to illustrate the present invention, and wherein implementation steps of each method etc. are all to be varied from , every equivalents carrying out on the basis of technical solution of the present invention and improvement, all should not exclude the protection model in the present invention Outside enclosing.

Claims (1)

1. a kind of method of Left Ventricle MR Images segmentation and three-dimensionalreconstruction is it is characterised in that segmentation and three-dimensionalreconstruction include Following steps：

First, initialization step, i.e. the input of CMR image：First by the CMR Image semantic classification of dicom form become png, bmp, jpg or The picture format of jpeg, pretreated image is read into scientific and engineering computing software by general built-in imread order In Matlab；

2nd, segmentation step：

1. the foundation of model：After the completion of initialization step, first set up and intend local binary simulated items, area item and penalty term three , described plan local binary simulated items, area item and penalty term three are added and just constitute gross energy functional, recycles and lead to Variation extreme value analysis method, can get variation level set evolutionary model, this model of iterative, until meeting following two Any one in condition, then stop iteration：

1) gross energy functional, that is, intend local binary simulated items, three summations being added of area item and penalty term are less than 10^-6；

2) iterative steps reach the upper limit 200；

2. initial value is chosen：In the image that Matlab reads in initialization step, arbitrarily choose initial level set function zero etc. Value line；

3. parameter setting：Corresponding model parameter in the variation level set evolutionary model of the setting present invention；

3rd, the convex closure result step of output segmentation：Based on the result of segmentation step, available algorithm of convex hull obtains corresponding convex closure Result；

4th, obtain the convex closure result of inside and outside film segmentation：The CMR sectioning image of left ventricle is directed respectively into segmentation module, can obtain The convex closure result of a series of inside and outside film image segmentation of this left ventricle；

5th, three-dimension curved surface, entity reconstruct：Using contour surface generating algorithm, the three-dimension curved surface reconstruction result of inside and outside film can be obtained； The plane convex closure result of inside and outside film segmentation is imported the macro document in three-dimensional Cartography designing software SolidWorks, then generates base 3D solid reconstruction result in Spline Interpolating Surfaces；

6th, 3D solid mesh generation step：The 3D solid reconstruction result that SolidWorks is generated imports finite element analysis Software Abaqus, obtains 3D solid mesh generation by calling Mesh functional module, that is, obtains corresponding node coordinate, section The information of point corresponding adjacent chained list, unit and its composition node；

Described in step 2 segmentation step, the method for building up of variation level set evolutionary model, comprises the steps：

I sets up and intends local binary simulated items, for processing gray scale non-uniform image：

Gray value is taken to be I₀Image in any one initial open region D, its border isAnd D_inAnd D_outRespectively Inside and outside region for open region D, introduces level set function φ (x), meets equation below：

$\left\{\begin{array}{cc}x\&Element;D_{in},& \mathrm{\&phi;}\left(x\right)>0,\\ x\&Element;C,& \mathrm{\&phi;}\left(x\right)=0,\\ x\&Element;D_{out},& \mathrm{\&phi;}\left(x\right)<0.\end{array}\right.---\left(1\right)$

Introduce edge detection operator g (x), construct new plan local binary simulated items F_New：

F_New=α ∫ g (x) | I₀(x)-c₁|²H(φ(x))dx+β∫g(x)|I₀(x)-c₂|²(1-H(φ(x)))dx (2)

c₁=∫ g (x) I₀(x)H(φ(x))dx/∫H(φ(x))dx (3)

c₂=∫ g (x) I₀(x)(1-H(φ(x)))dx/∫(1-H(φ(x)))dx (4)

Wherein α, β are model parameter, c₁、c₂For two value parameters, I₀X () is image intensity value at x, H () is Heaviside letter Number, g (x) is then following edge detection operator

G=1/ (1+ K_σ*I₀) (5)

Wherein K_σ*I₀It is initial pictures I₀With Gaussian coreConvolution, parameter σ=1.5, n=10；

II introduces area item and penalty term

Introduce area item F_Geom：

F_Geom=γ ∫ g (x) H (φ (x)) dx (6)

Using following weighting penalty term F_P：

$F_p=\mathrm{\&lambda;}\&Integral;\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}\left(x\right)|-1)}^{2}dx---\left(7\right)$

Wherein, γ, λ are weighted factor；

III constructs gross energy functional, defines gross energy functional F_TAs follows：

$\begin{array}{c}{F}_T=F_{New}+F_{Geom}+F_P\\ =\mathrm{\&alpha;}\&Integral;g\left(x\right)|I_0\left(x\right)-c_1{|}^{2}H\left(\mathrm{\&phi;}\left(x\right)\right)dx\\ +\mathrm{\&beta;}\&Integral;g\left(x\right)|I_0\left(x\right)-c_2{|}^2\left(1-H\left(\mathrm{\&phi;}\left(x\right)\right)\right)dx\\ +\mathrm{\&gamma;}\&Integral;g\left(x\right)H\left(\mathrm{\&phi;}\left(x\right)\right)dx+\mathrm{\&lambda;}\&Integral;\frac{1}{2}{\left(|\&dtri;\mathrm{\&phi;}\left(x\right)|-1\right)}^{2}dx\end{array}---\left(8\right)$

IV sets up new variation level set evolutionary model：

Preset parameter α, β, γ and λ, solve above-mentioned gross energy functional F_TWith regard to the variation extreme-value problem of level set function φ (x), lead Go out the Euler-Lagrange equation with regard to φ (x)：

$\frac{\&part;\mathrm{\&phi;}\left(x\right)}{\&part;t}=-\mathrm{\&alpha;}g{\left(I_0\left(x\right)-c_1\right)}^2+\mathrm{\&beta;}g{\left(I_0\left(x\right)-c_2\right)}^2+\mathrm{\&gamma;}g\mathrm{\&delta;}\left(\mathrm{\&phi;}\left(x\right)\right)+\mathrm{\&lambda;}\&lsqb;\mathrm{\&Delta;}\mathrm{\&phi;}\left(x\right)-\&dtri;\&CenterDot;\frac{\&dtri;\mathrm{\&phi;}\left(x\right)}{|\&dtri;\mathrm{\&phi;}\left(x\right)|}\&rsqb;---\left(9\right)$

Wherein, δ () is Dirac Delta function, formula (9) the variation level set evolutionary model that i.e. present invention is set up, energy The inside and outside film of left ventricle effectively in Ground Split CMR image.