CN106780518B - A kind of MR image three-dimensional interactive segmentation method of the movable contour model cut based on random walk and figure - Google Patents
A kind of MR image three-dimensional interactive segmentation method of the movable contour model cut based on random walk and figure Download PDFInfo
- Publication number
- CN106780518B CN106780518B CN201710073688.XA CN201710073688A CN106780518B CN 106780518 B CN106780518 B CN 106780518B CN 201710073688 A CN201710073688 A CN 201710073688A CN 106780518 B CN106780518 B CN 106780518B
- Authority
- CN
- China
- Prior art keywords
- segmentation
- image
- region
- point
- energy
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10088—Magnetic resonance imaging [MRI]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30096—Tumor; Lesion
Abstract
The present invention discloses a kind of MR image three-dimensional interactive segmentation method of movable contour model cut based on random walk and figure comprising step: S1, and selected seed point interception includes the partial 3 d MR image data of hypophysoma, and then obtains the initial boundary curved surface of segmentation；S2 is based on initial boundary curved surface, establishes hybrid activity skeleton pattern, obtain model energy function；S3 carries out discretization to energy function；S4 carries out the building of figure using 6 neighborhoods of each tissue points using each tissue points of the partial 3 d MR image of interception as the node of figure；Initial value is assigned respectively to the tissue points inside and outside initial boundary curved surface；And corresponding weight is assigned respectively to the side connected between node and meeting point between node, between node and source point according to the energy function after discrete；S5, the figure based on building carry out figure and cut calculating, obtain segmentation result, and the boundary surface extracted in segmentation result is segmentation contour；S6 is repeated step S2 to S5, until segmentation result is restrained, is obtained final segmentation contour with the initial boundary curved surface in current segmentation contour substitution S1.Using the three-dimensional segmentation of the achievable hypophysoma MR image of the present invention, the accuracy of hypophysoma image segmentation is improved.
Description
Technical field
The present invention relates to technical field of image segmentation, especially carry out region segmentation to the human brain image of magnetic resonance imaging
, the MR image three-dimensional interactive segmentation method based on the movable contour model that random walk and figure are cut.
Background technique
Magnetic resonance imaging (Magnetic resonance imaging, MRI) has excellent soft tissue resolving power, can be right
Partes corporis humani position multi-angle, multiple plane imaging, the application of a variety of imaging methods and pulse train technology facilitate to hypophysoma
Positioning, etiologic diagnosis.Usually clinically medical image lesion segmentation can by doctor on bidimensional image by sequence hand drawing
Hypophysoma boundary is completed.According to solid tumor therapeutic evaluation (RECIST) standard, radiating doctor will recognise that brain on bidimensional image
The high brightness in portion/low-light level tissue measures hypophysoma lesion one, the diameters of both direction, but due to the shape of hypophysoma
Changeable, single diameter can not indicate the volume of hypophysoma.Huge, the people additionally, due to the data volume of brain scans sequence image
Work point cuts time and effort consuming.And manually divide and image information analysis result by doctor professional knowledge and skilled operation degree shadow
Ring, often vary with each individual, subjectivity it is larger, it is easy to generate uncertain result.Multi-modal hypophysoma magnetic resonance
(Magnetic resonance, MR) image processing and analysis, integrated use image procossing, image analysis, pattern-recognition, machine
The present computer technologies such as study, can greatly improve the status of clinical hypophysoma medical image processing and analysis, provide higher
It imitates, is accurate, objectively analyzing result.
Hypophysoma belongs to one of brain tumor type, and the medical image tool of multi-modal Pituitary adenoma is found in clinical practice
There is complicated characteristic.Its MR image and other brain tumors MR imaging characteristic have similarity, are mainly presented with: the structure of hypophysoma
It is complicated；Its space occupy-place, shapes and sizes are uncertain；Pathological tissues, the intensity profile of surrounding soft tissue are uneven, Density Distribution
There is aliasing, boundary is relatively fuzzyyer；There are larger differences for the institutional framework of Different Individual；Again since MRI technique principle and image are adopted
The influence of various factors during collection, there are noise, artifacts etc. for image.Therefore to realize that the correct segmentation to hypophysoma is very tired
It is difficult.The segmentation task of brain tumor include normal cerebral tissue's (white matter of brain, ectocinerea and cerebrospinal fluid), pathological tissues (oedema, tumour and
Its internal necrosis and capsule change etc.).Since the 1970s, various types of brain tumor dividing methods are emerged, mainly
The dividing method of the dividing method based on region, the dividing method based on edge and bond area and boundary can be divided into, also
Neural network, fuzzy clustering, Markov random field model, other dividing methods based on graph theory and based on map.
Due to the complexity of brain medical image, all kinds of dividing methods still can only specifically be solved with regard to particular problem, not had still
A kind of general dividing method is formed.Following feature is presented in current research tendency:
(1) a variety of partitioning algorithms for merging respective advantage combine.Such as combine the image segmentation of fuzzy clustering and deformation model
Method, the dividing method based on region growing and neural network；There are also what is developed in recent years, the partitioning algorithm based on graph theory,
Such as figure, which is cut, to be cut with active shape model joint partitioning algorithm, figure and combines partitioning algorithm with deformation model.
(2) automanual cutting techniques are still the emphasis direction of research.Full-automatic dividing algorithm can greatly improve brain
The segmentation efficiency of tumour, but also need that doctor is allowed to carry out special cutting operation according to specific case in clinic.Interactive half
Automatic division method allows doctor to assess segmentation result, manual intervention segmentation, makes up certainly on the basis of automatic segmentation result
Dynamic deficiency of the partitioning algorithm when carrying out special segmentation task reaches precisely segmentation, quantitative analysis.Such as apply the brain of random walk
Tumour and lung neoplasm partitioning algorithm；Also some Interactive Segmentation softwares, such as 3D Slicer are developed.
The brain tumor partitioning algorithm of automatic, the semi-automatic interactive mode occurred in recent years, can be roughly divided into based on engineering
Two major class practised: generating algorithm and distinguished number.Generating algorithm is from a series of prior images (known tissue classify image)
The middle gray space distribution for learning healthy brain tissue type, using study to distribution function go " to predict " health tissues, find out
Boundary segmentation tumour.Distinguished number directly learns the gray difference between tumour and normal cerebral tissue, with the differentiation letter learnt
The unlabelled pixel/voxel of number classification.Patient's training image that some manual markings had both can be used in learning procedure carries out offline
It executes, seed region can also be marked to complete online on current patient image by user.Generating algorithm is for normal brain activity group
It is advantageous but then difficult for the lesion segmentation of arbitrary shape to knit segmentation.The advantages of distinguished number is to have merged normal brain activity group
Knit with tumour of both intensity profile, be conducive to the segmentation of tumour；But local image characteristics are used only in it, to multi-modal image
The image grayscale of different tissues requires relatively uniform in data, this is extremely difficult to, for MR image.By two kinds
Method combines, and merges the brain Anatomical Structure knowledge of background and the gamma characteristic of lesion, may improve the correctness of segmentation.
The multi-modal image division method of brain tumor based on deformation model is had also appeared in recent years.Such method is from image
The priori knowledges such as constraint information, target position, size and shape are obtained, the bottom-up information of image itself had not only been utilized but also has combined
The upper layer information of human anatomic structure is suitable for the three-dimensional data segmentation of the multi-modal image of brain tumor.Water is based in deformation model
The theoretical geometric deformation model of flat collection, due to topologies change problem when it can solve curve evolvement and other be better than
The characteristics of parameter deformation model, has obtained extensive attention and development.Single level set function can only divide the image into target area
Domain and the part of background area two, that is, image two-phase segmentation produce for the needs for meeting multiple target area segmentations in image
Multi-phase horizontal set (Multiphase Level Set) method.By the development of recent two decades, multi-phase horizontal set thought is gradually answered
Use Mumford-Shah model, on Chan-Vese model.The LBF model that Li et al. people proposes introduces gaussian kernel function, enhances
Level set function is used only in certain multiphase images, solves multiphase for capture ability of the Level Set Method to edge
Level set algorithm need to be using computational complexity problem brought by multiple level set functions；And developed by Peng et al. as normalization
LBF model.But there are still method complexity, computationally intensive problem for the solution of Level Set Method at present.
Using J.Egger as the research team of representative propose using the figure based on graph theory cut with graph search method, design
Automanual hypophysoma dividing method based on template.But from the point of view of the research paper delivered from them, the disease that mainly uses
Example is regular shape, the hypophysis structure of approximate ellipse shape, and segmentation difficulty is relatively little.Algorithm and deformation model are cut in figure
The common dividing method that shape is also the brain tumor that segmented shape is relatively fixed, variation is little is added in energy function.But
Most hypophysoma shape does not have systematicness, thereby increases and it is possible to have to the appearance of the Infiltratings of other brain tissues, so, template or
Shape acts on limited played in the hypophysoma cutting procedure of irregular shape.
Summary of the invention
The technical problem to be solved in the present invention are as follows: propose that one kind is cut based on random walk (Random walk, RW) and figure
Movable contour model (the Random walks and graph cuts based active contour of (Graph cut)
Model, RGCACM) MR image three-dimensional interactive segmentation method, realize hypophysoma MR image three-dimensional segmentation, improve hypophysoma figure
As the accuracy of segmentation.
A kind of technical solution that the present invention takes specifically: MR figure of the movable contour model cut based on random walk and figure
As three-dimension interaction dividing method, comprising the following steps:
S1 obtains initial boundary curved surface, comprising steps of
S11, seed point needed for choosing Random Walk Algorithm, the point using the seed point as three-dimensional center, from brain MR tri-
Interception includes the partial 3 d MR image data of hypophysoma in dimension data；
S12 is handled using partial 3 d MR image data of the three-dimensional random migration algorithm to interception, obtains each voxel
Point reaches the probability three-dimensional figure of seed point；
S13, using maximum expected value algorithm picks probability threshold value, to be partitioned into suspected target region, by suspected target area
Foreground point of the voxel as movable contour model algorithm in domain, the boundary surface, that is, initial boundary curved surface in suspected target region；
S2 is based on initial boundary curved surface, establishes hybrid activity skeleton pattern: to the energy function of standard actions skeleton pattern
Middle boundary energy item uses standard Geodesic active contour, and region energy item is general using the gray scale maximum a posteriori of Local Gaussian Model description
Rate, and then obtain the energy function of hybrid activity skeleton pattern；
S3, model discretization, comprising:
S31 defines a binary variable to each tissue points, and assignment 1 and 0, which respectively corresponds, represents foreground point and background dot,
With by region energy item function discretization；
S32 is estimated using cutting by boundary energy item function i.e. Geodesic active contour discretization；
S33, the energy function after obtaining discretization；
S4, structure figures: for the partial 3 d MR image data of interception, using each tissue points therein as the section of figure
Point carries out the building of figure using 6 neighborhoods of each tissue points；To the tissue points in initial boundary curved surface and outside initial boundary curved surface
Initial value is assigned respectively；And the side being connect to the side connected between node, node with source point according to the energy function after discrete, and section
Point assigns corresponding weight with the side that meeting point connects respectively；
S5, the figure based on step S4 building carry out figure and cut calculating: segmentation result is obtained using max-flow and minimal cut algorithm,
The boundary surface extracted in segmentation result is segmentation contour；
S6, using current segmentation contour as initial boundary curved surface, iteration step S2 to step S5, until segmentation result
Convergence, exports final segmentation result.
Further, the invention also includes step S7, step S5 is obtained using three-dimensional median filtering evolution curved surface into
Row processing, segmentation result after being filtered；Step S6 is using the segmentation contour after current filter as initial boundary curved surface.
Preferably, in step S11, the selecting step of seed point is, in the two-dimensional display image of sagittal plain, Selection Center
Position slice, approximate center position is chosen a bit in hypophysoma in a slice, that is, is used as seed point.The approximate center position is
Near center location, this selecting step are to choose manually, therefore judge whether approximate center location-dependent query chooses people to selected seed point
It carries out, on the basis of as close as possible to center, general error range is within the 30% of hypophysoma diameter.
Preferably, in step S11, the size of the partial 3 d MR image data comprising hypophysoma is 71 × 71 × 41
A voxel.
It in step S13, is determined according to the grey value profile of probability graph, value range is 0~1.
Preferably, in step S2, given image domainI (x): Ω → R is two dimensional image；Closed curve C is by image
It is divided into two isolated areas, L (C) indicates the length of C；If curve C inner region is Ω1, curve C exterior domain is outside Ω2, for figure
Each pixel p of image field Ω considers that its radius is the circle shaped neighborhood region of ρ, is defined as Ox=y:| x-y | and≤ρ }, x is circle shaped neighborhood region
Central pixel point, y are any pixel in circle shaped neighborhood region in addition to x；
Then region Ox∩ΩiInterior pixel gray level probability density function indicates are as follows:
U in formulai(x) and σiIt (x) is respectively region Ox∩ΩiThe local gray level mean value and variance, I (y) of interior pixel x be
The gray value of pixel y,
Boundary energy item EEdgeUsing standard Geodesic active contour, region energy item EReigonIt is described using Local Gaussian Model
Gray scale maximum a posteriori probability form energy function E indicate are as follows:
β is any normal number in formula, and function w is a weighting function.
Further, in step S2 of the present invention, weighting function w is indicated using truncation Gauss kernel form are as follows:
α is so that the constant that ∫ w (x, y)=1 is set up in formula.That is, the present invention is to neighborhood OxIn close to central point x pixel y
Gray value I (y) assign big weighted value.
It is in order to which the energy function for cutting algorithm with figure blends that step S3 of the present invention, which carries out discretization to model,.Preferably,
In step S31, a binary variable x is defined to each tissue points pp, discretization is carried out with the region energy item to energy function:
Object indicates that foreground point, that is, segmentation object, Background indicate background dot；
It is indicated after region energy item function discretization are as follows:
P, q refers to voxel different in image, xp、xqThe bi-values (0 or 1) of different voxels are referred to, N (p) refers to voxel p's
Neighborhood.
us(p)、ut(p) and σs(p)、σtIt (p) is respectively ui(x) and σi(x) discrete form, is expressed as follows:
In step S32, indicated after boundary energy item function discretization are as follows:
ωpqFor non-negative side right value, it is expressed as three dimensional form i.e.:
ΔΦpq=Δ φpq·Δψpq=π2/ 4, the angle interval of unit circle in 6 neighborhoods is corresponded to, δ is the ruler of unit lattice
It is very little, can value δ=1, | epq| it is boundary epqUnit length, the present invention in value | epq|=1, D are that constant Riemann estimates.It is excellent
Choosing, ω in the present inventionpq=π/4.
Based on the area item and border item after discrete, can be obtained it is discrete after hybrid activity model energy function.
Preferably, to avoid figure from cutting the equilibrium problem in region and border item in the standard energy function of algorithm, and region
It will affect figure and cut algorithm and tend to be partitioned into zonule.Energy function after discrete uses multiplication form, that is, by region energy item
As the weight term of boundary energy item, it is discrete after energy function indicate are as follows:
E=EEdge(p,q)·E'Region(p) (10)
If neighboring voxels p, q are located in or beyond segmentation object simultaneously, have | Es(p)+Et(q) | 0 He of > | Es(q)+Et
(p) | > 0 is set up；If p, q are located at the both sides of object boundary, then have | Es(p)+Et(q) | ≈ 0 or | Es(q)+Et(p) | ≈ 0 at
It is vertical.So only working as neighboring voxels p, when q is located at the both sides of object boundary, i.e. p, q are belonging respectively to target and background region
When, formula (10) can just obtain minimum value.
Preferably, in step S4, defining n-link weight in side interconnected between tissue points is EEdge(p,q)·
E'Region(p)；If | Es(p) | > | Et(p) | and xp=1, then tissue points and meeting point T (refer to segmentation object)It is connected
The weight of side t-link is set as infinite；If | Es(p) | < | Et(p) | and xp=0, then tissue points and source point S (are referred to and are carried on the back
Scape) weight of side t-link that is connected is set as infinite.Finally i.e. available max-flow/minimal cut algorithm solves segmentation result,
The energy function of figure is solved, obtains the label (0 or 1) of each node in figure, this is the prior art.
Beneficial effect
The method of the present invention is movable contour model (the Graph cuts based active contour cut based on figure
Model, GCACM) algorithm improvement and expansion, using movable contour model solution figure cut algorithm segmenting edge ladderization lack
Point, while the max-flow also cut using figure/minimal cut algorithm solves the computational complexity of movable contour model Level Set Method
Problem, and it is easy the shortcomings that sinking into partial points.The boundary seed point select permeability in algorithm is cut for solution figure, using random trip
It walks algorithm and probably determines target area, and (nose rear, eyes center, view are handed over using the relatively fixed anatomical position of hypophysoma
Fork lower section etc.), for Random Walk Algorithm, interaction provides seed point manually.Binding site priori knowledge and figure cut, active contour mould
Type realizes the three-dimensional segmentation of hypophysoma, can be improved in the positioning of such as hypophysoma, is associated with the pathology of surrounding tissue and volume
Utilization validity in the quantitative analyses such as measurement.
I.e. the present invention provides a kind of semi-automatic hypophysoma MR image segmentation algorithm with feasibility and validity for the first time,
The two quasi-representative image segmentation algorithms based on graph theory and based on movable contour model are combined, regular and irregular may be implemented and hang down
The segmentation of body tumor and volumetric estimate especially also have good accuracy to the segmentation of infiltrative type hypophysoma.Interactive operation is simply easy
Row, remaining segmentation task can be automatically performed by only needing user's mouse to choose two points.It is retouched in algorithm using local Gaussian distribution
Image-region gray scale is stated, can effectively solve the problems, such as that MR image hypophysoma lesion region gray scale is non-uniform.
Detailed description of the invention
Fig. 1 show the method for the present invention flow diagram；
Fig. 2 show cutting procedure of the present invention and result schematic diagram, and wherein a- user chooses a foreground point (yellow dots)
It is superimposed upon on original sagittal plain two dimensional image with a background dot (green point), red boxes indicate the region of interest of interception；b-
The initial profile being superimposed upon on probability graph；The segmentation contour of c- superposition on the original image；The three-dimensional rendering mould of d- segmentation result
Type；
Fig. 3 show the circle shaped neighborhood region schematic diagram of pixel in step S2, and dashed circle indicates the neighborhood of pixel x.
Specific embodiment
It is further described below in conjunction with the drawings and specific embodiments.
The present invention is based on the MR image three-dimensional interactive segmentation methods for the movable contour model that random walk and figure are cut, including with
Lower step:
S1 obtains initial boundary curved surface, comprising steps of
S11, seed point needed for choosing Random Walk Algorithm, the point using the seed point as three-dimensional center, from brain MR tri-
Interception includes the partial 3 d MR image data of hypophysoma in dimension data；
S12 is handled using partial 3 d MR image data of the three-dimensional random migration algorithm to interception, obtains each voxel
Point reaches the probability three-dimensional figure of seed point；
S13, using maximum expected value algorithm picks probability threshold value, to be partitioned into suspected target region, by suspected target area
Foreground point of the voxel as movable contour model algorithm in domain, the boundary surface, that is, initial boundary curved surface in suspected target region；
S2 is based on initial boundary curved surface, establishes hybrid activity skeleton pattern: to the energy function of standard actions skeleton pattern
Middle boundary energy item uses standard Geodesic active contour, and region energy item is general using the gray scale maximum a posteriori of Local Gaussian Model description
Rate, and then obtain the energy function of hybrid activity skeleton pattern；
S3, model discretization, comprising:
S31 defines a binary variable to each tissue points, and assignment 1 and 0, which respectively corresponds, represents foreground point and background dot,
With by region energy item function discretization；
S32 is estimated using cutting by boundary energy item function i.e. Geodesic active contour discretization；
S33, the energy function after obtaining discretization；
S4, structure figures: for the partial 3 d MR image data of interception, using each tissue points therein as the section of figure
Point carries out the building of figure using 6 neighborhoods of each tissue points；To the tissue points in initial boundary curved surface and outside initial boundary curved surface
Initial value is assigned respectively；And the side being connect to the side connected between node, node with source point according to the energy function after discrete, and section
Point assigns corresponding weight with the side that meeting point connects respectively；
S5, the figure based on step S4 building carry out figure and cut calculating: segmentation result is obtained using max-flow and minimal cut algorithm,
The boundary surface extracted in segmentation result is segmentation contour；
S6, using current segmentation contour as initial boundary curved surface, iteration step S2 to step S5, until segmentation result
Convergence, exports final segmentation result.
Further, the present invention may also include step S7, implement between S5 and S6, i.e., using three-dimensional median filtering to step
The evolution curved surface that rapid S5 is obtained is handled, segmentation result after being filtered；Step S6 is made with the segmentation contour after current filter
For initial boundary curved surface.
Embodiment 1
Refering to what is shown in Fig. 1, the step of the present embodiment method, sketches are as follows: image segmentation seed point is chosen and initialization boundary is bent
Face acquisition, the three-dimensional median filtering post-processing of the iterative segmentation of improved GCACM algorithm, segmentation result.Namely it is following initial
Change step, segmentation step and post-processing step, is described in detail below.
1, initialization step.Predominantly seed point needed for the interactive Random walk algorithm of selection manually of user, and
Initial boundary curved surface needed for obtaining partitioning algorithm.
To reduce data computation complexity, and using the medical ground knowledge of user, cut in protocerebrum archicerebrum portion MR three-dimensional data
Take the data cube comprising hypophysoma.Specific practice is the Selection Center position slice in the two-dimensional display image of sagittal plain,
Approximate center position mouse is chosen a bit in hypophysoma again.Using the point as center data intercept cube, (this experiment selects 71
The size of × 71 × 41 sizes), while the point is also used as the foreground point of Random walk algorithm, then in the sectioning image
Arbitrarily choose the background dot for a little serving as Random walk algorithm in background area.
Three-dimensional Random walk algorithm process is carried out to the data cube of interception, each voxel is obtained and reaches foreground seeds
The probability three-dimensional figure of point.Using EM algorithm picks probability threshold value, primary segmentation goes out the region of suspected target, and the voxel in region will
It can be used as the foreground point of GCACM, the boundary surface in region is initial boundary curved surface.
Process such as Fig. 2 signal of initial boundary curved surface is obtained, with a width sagittal plain two dimension of approximate hypophysoma center
For image, a- user chooses a foreground point (yellow dots) first and a background dot (green point) is superimposed upon original sagittal plain
On two dimensional image, red boxes indicate the region of interest of interception；It is the initial profile being superimposed upon on probability graph in b figure, is in c figure
It is superimposed segmentation contour on the original image, is the three-dimensional rendering model of segmentation result in d.
2, segmentation step, the building and solution of design, model discretization, figure including hybrid activity skeleton pattern.
(1) foundation of hybrid activity skeleton pattern
Refering to what is shown in Fig. 3, given image domainI (x): Ω → R is two dimensional image.Closed curve C is divided the image into
Two isolated areas, L (C) indicate the length of C.If curve C inner region is Ω1, curve C exterior domain is outside Ω2.For image area
Each pixel p of Ω considers that its radius is the circle shaped neighborhood region of ρ, is defined as Ox={ y:x-y≤ρ }.
Region Ox∩ΩiInterior pixel gray level probability density function can be expressed as follows:
U in formulai(x) and σiIt (x) is respectively region Ox∩ΩiThe local gray level mean value and variance of interior pixel x.
Energy function combines boundary energy item and region energy item, wherein border item EEdgeUsing standard geodesic curve mould
Type, area item EReigonUsing gray scale Local Gaussian Model and it is expressed as maximum a posteriori probability (Maximum a posteriori
Probability, MAP) form.Energy function can be written as follow form:
β is any normal number in formula.Function w is a weighting function, to neighborhood OxIn close to central point x pixel
Y gray value I (y) assigns big weighted value.Using the w of truncation Gauss kernel form in this algorithm, it is expressed as follows:
Wherein α is constant, so that ∫ w (x, y)=1 is set up.
(2) model discretization
In order to which the energy function for cutting algorithm with figure merges, a binary variable x is defined to tissue points p each in figurep, assignment
1 and 0 respectively represents foreground point (segmentation object) and background dot (background), by area item energy function discretization, such as formula (4) institute
Show.
It recycles cutting for Boykov et al. proposition to estimate (Cut metric), Geodesic active contour discretization is indicated, then can
The border item of flow function indicates are as follows:
ω in formulapqIt is non-negative side weight, it is as follows can be expressed as three dimensional form:
ΔΦ in above formulapq=Δ φpq·Δψpq=π2The angle interval of unit circle, δ are unit lattice in/4 corresponding 6 neighborhoods
Size, δ=1, | epq| it is boundary epqUnit length, | epq|=1, D are that constant Riemann estimates, and take ω in the present inventionpq=π/
4。
The area item discrete form of energy function can be expressed as follows:
P, q refer to voxel different in image, x in formulap、xqThe bi-values (0 or 1) of different voxels are referred to, N (p) refers to voxel
The neighborhood of p；us(p)、ut(p) and σs(p)、σtIt (p) is respectively ui(x) and σi(x) discrete form, as follows:
To avoid figure from cutting the equilibrium problem in region and border item in the standard energy function of algorithm, and area item will affect
Figure cuts algorithm and tends to be partitioned into zonule.This algorithm uses the energy function of multiplication form, i.e., serves as border item using area item
Weight term, as shown in formula (10), (11):
E=EEdge(p,q)·E'Region(p) (10)
If neighborhood territory pixel p, q are located at simultaneously in or beyond segmentation object (i.e. hypophysoma region), have | Es(p)+Et(q)|
0 He of > | Es(q)+Et(p) | > 0 is set up；If p, q are located at the both sides of object boundary, then have | Es(p)+Et(q) | ≈ 0 or | Es
(q)+Et(p) | ≈ 0 is set up.So only working as neighborhood territory pixel p, when q is located at the both sides of object boundary, formula (10) can just be obtained
Minimum value.Formula (10) obtains minimum value and represents the energy of figure as minimum value, i.e., energy function obtains optimal result, and energy reaches
To it is optimal when segmentation complete, the tissue points for adhering to target and background region separately at this time have been assigned different mark values.
(3) building and solution of figure
The present invention serves as the node of figure with each voxel in MR data volume, constructs three-dimensional figure using 6 neighborhoods.Voxel in figure
Between n-link weight in side interconnected be set as EEdge(p,q)·E'Region(p)；If | Es(p) | > | Et(p) | and xp
=1, then voxel is set as infinite with meeting point T (corresponding target) the side t-link weight being connected；If | Es(p) | < | Et(p)|
And xp=0, then voxel is set as infinite with source point S (corresponding background) the side t-link weight being connected.Finally with max-flow/
Minimal cut algorithm solves segmentation result, is herein the prior art, does not repeat them here.
3, it post-processes
To the resulting evolution curved surface of segmentation step again with three-dimensional median filtering post-process, then according to segmentation result into
Row iteration segmentation obtains final segmentation result until result convergence
Iterative segmentation process description of the invention is as follows:
(1) initial boundary curved surface S is obtained using Random walk algorithmur；Initial assignment SurInterior voxel p corresponding two
Metavariable xpIt is 1, SurThe outer corresponding x of voxel ppIt is 0；
(2) mean value and variance u of local gray level are calculated with formula (8), (9)s(p)、ut(p) and σs(p)、σt(p)；
(3) the energy function structure figures indicated with formula (10)；
(4) segmentation result is solved with max-flow/minimal cut algorithm, updates the corresponding binary of voxel p further according to segmentation result
Variable xpValue, in fact, solving while obtaining segmentation result namely having updated the corresponding binary variable value of voxel p；
(5) smooth post-processing is carried out to segmentation result with three-dimensional median filter；
(6) initial profile is replaced with current segmentation contour, i.e. it is straight to repeat step (2)-(5) for initial boundary curved surface in S2
It is restrained to segmentation result, it is to restrain that the segmentation result acquired, which no longer changes,.
4, experimental result
Split-run test, result schematic diagram such as Fig. 2 are carried out in the MR T1W data of 23 hypophysoma patients using the present invention
Shown in middle c, d.Table 1 provides 23 case segmentation result Dice coefficient values (Dice similarity coefficient, DSC)
Minimum value, maximum value and mean value and standard deviation square value.DSC value is defined as experimental result VEWith goldstandard VGRelative superposition
Volume ratio, as shown in formula (12):
It furthermore is the superiority for showing the method for the present invention segmentation effect, by this algorithm and other two kinds of segmentations delivered
Algorithm (loading on the GrowCut algorithm and general GCACM algorithm based on region growing in software 3DSlice) is split knot
The comparison of fruit DSC value.Comparison result is refering to table 1, it can be seen that the segmentation result DSC value mean value of this algorithm is 88.36%, is better than
The 70.88% of 81.61% and general GCACM algorithm of GrowCut algorithm.Illustrate that this algorithm has in hypophysoma segmentation application
Higher validity and accuracy.
1 algorithm of table and two kinds have delivered the segmentation result DSC value of partitioning algorithm
| Results | Proposed method | GrowCut | GCACM |
| min | 74.86% | 67.62% | 49.45% |
| max | 93.35% | 89.07% | 78.55% |
| mean±SD | 88.36% ± 5.61% | 81.61% ± 7.37% | 70.88% ± 8.08% |
The above is only a preferred embodiment of the present invention, it is noted that for the ordinary skill people of the art
For member, without departing from the technical principles of the invention, several improvement and deformations can also be made, these improvement and deformations
Also it should be regarded as protection scope of the present invention.
Claims (10)
1. a kind of MR image three-dimensional interactive segmentation method of movable contour model cut based on random walk and figure, characterized in that
The following steps are included:
S1 obtains initial boundary curved surface, comprising steps of
S11, seed point needed for choosing Random Walk Algorithm, the point using the seed point as three-dimensional center, from tri- dimension of brain MR
It include the partial 3 d MR image data of hypophysoma according to middle interception；
S12 is handled using partial 3 d MR image data of the three-dimensional random migration algorithm to interception, is obtained each tissue points and arrive
Up to the probability three-dimensional figure of seed point；
S13, will be in suspected target region to be partitioned into suspected target region using maximum expected value algorithm picks probability threshold value
Foreground point of the voxel as movable contour model algorithm, the boundary surface, that is, initial boundary curved surface in suspected target region；
S2 is based on initial boundary curved surface, establishes hybrid activity skeleton pattern: to side in the energy function of standard actions skeleton pattern
Boundary's energy term uses standard Geodesic active contour, and region energy item uses the gray scale maximum a posteriori probability of Local Gaussian Model description,
And then obtain the energy function of hybrid activity skeleton pattern；
S3, model discretization, comprising:
S31 defines a binary variable to each tissue points, and assignment 1 and 0, which respectively corresponds, represents foreground point and background dot, will
Region energy item function discretization；
S32 is estimated using cutting by boundary energy item function i.e. Geodesic active contour discretization；
S33, the energy function after obtaining discretization；
Structure figures: S4 adopts the partial 3 d MR image data of interception using each tissue points therein as the node of figure
The building of figure is carried out with 6 neighborhoods of each tissue points；To the tissue points difference in initial boundary curved surface and outside initial boundary curved surface
Assign initial value；And the side that is connect with source point according to the energy function after discrete to the side connected between node, node and node with
The side of meeting point connection assigns corresponding weight respectively；
S5, the figure based on step S4 building carry out figure and cut calculating: obtaining segmentation result using max-flow and minimal cut algorithm, extract
Boundary surface in segmentation result is segmentation contour；
S6, using current segmentation contour as initial boundary curved surface, iteration step S2 to step S5, until segmentation result is received
It holds back, exports final segmentation result.
2. according to the method described in claim 1, it is characterized in that, step S5 further includes, using three-dimensional median filtering to obtaining
Segmentation contour is handled, segmentation result after being filtered；Step S6 is using the segmentation contour after current filter as initial boundary
Curved surface.
3. method according to claim 1 or 2, characterized in that in step S11, the selecting step of seed point is, in sagittal
In the two-dimensional display image of position, Selection Center position slice, approximate center position is chosen a bit in hypophysoma in a slice, that is, makees
For seed point.
4. method according to claim 1 or 2, characterized in that in step S11, the data cube comprising hypophysoma
Body is having a size of 71 × 71 × 41 voxels.
5. method according to claim 1 or 2, characterized in that in step S2, given image domainI(x):Ω→
R is two dimensional image；Closed curve C divides the image into two isolated areas, and L (C) indicates the length of C；If curve C inner region is
Ω1, curve C exterior domain is Ω2, for each pixel p of image area Ω, consider that its radius is the circle shaped neighborhood region of ρ, be defined as Ox=
Y:| x-y | and≤ρ }, x is the central pixel point of circle shaped neighborhood region, and y is any pixel in circle shaped neighborhood region in addition to x；
Then region Ox∩ΩiInterior pixel gray level probability density function indicates are as follows:
U in formulai(x) and σiIt (x) is respectively region Ox∩ΩiThe local gray level mean value and variance of interior pixel x, I (y) are pixel
The gray value of y,
Boundary energy item EEdgeUsing standard Geodesic active contour, region energy item EReigonThe gray scale described using Local Gaussian Model
The energy function E of maximum a posteriori probability form is indicated are as follows:
β is any normal number in formula, and function w is a weighting function.
6. according to the method described in claim 5, it is characterized in that, in step S2, weighting function w using truncation Gauss kernel form,
It indicates are as follows:
α is so that the constant that ∫ w (x, y)=1 is set up in formula.
7. according to the method described in claim 5, it is characterized in that, in step S31, define binary to each tissue points p and become
Measure xp:
Object indicates that foreground point, that is, segmentation object, Background indicate background dot；
Discretization is carried out to the region energy item of energy function, is indicated after region energy item function discretization are as follows:
P, q refers to voxel different in image, xp、xqThe bi-values of different voxels are referred to, the bi-values are 0 or 1, and N (p) refers to
The neighborhood of voxel p；
us(p)、ut(p) and σs(p)、σtIt (p) is respectively ui(x) and σi(x) discrete form is expressed as follows:
In step S32, indicated after boundary energy item function discretization are as follows:
ωpqFor non-negative side right value, it is expressed as three dimensional form i.e.:
ΔΦpq=Δ φpq·Δψpq=π2/ 4, the angle interval of unit circle in 6 neighborhoods is corresponded to, δ is the size of unit lattice, | epq
| it is boundary epqUnit length, D is that constant Riemann estimates.
8. according to the method described in claim 7, it is characterized in that, ωpq=π/4.
9. according to the method described in claim 7, it is characterized in that, it is discrete after energy function use multiplication form, that is, by region
Weight term of the energy term as boundary energy item, it is discrete after energy function indicate are as follows:
E=EEdge(p,q)·E'Region(p) (10)
If neighborhood territory pixel p, q are located in or beyond segmentation object simultaneously, have | Es(p)+Et(q) | 0 He of > | Es(q)+Et(p) | >
0 sets up；If p, q are located at the both sides of object boundary, then have | Es(p)+Et(q) | ≈ 0 or | Es(q)+Et(p)|≈0。
10. according to the method described in claim 8, it is characterized in that, define tissue points between n-link weight in side interconnected
For EEdge(p,q)·E'Region(p)；If | Es(p) | > | Et(p) | and xp=1, then the side that tissue points are connect with meeting point T-phase
The weight of t-link is set as infinite, and the meeting point T refers to segmentation object；If | Es(p) | < | Et(p) | and xp=0, then
The weight for the side t-link that tissue points are connected with source point S is set as infinite, and the source point S refers to background.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710073688.XA CN106780518B (en) | 2017-02-10 | 2017-02-10 | A kind of MR image three-dimensional interactive segmentation method of the movable contour model cut based on random walk and figure |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710073688.XA CN106780518B (en) | 2017-02-10 | 2017-02-10 | A kind of MR image three-dimensional interactive segmentation method of the movable contour model cut based on random walk and figure |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106780518A CN106780518A (en) | 2017-05-31 |
| CN106780518B true CN106780518B (en) | 2019-06-25 |
Family
ID=58956878
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710073688.XA Active CN106780518B (en) | 2017-02-10 | 2017-02-10 | A kind of MR image three-dimensional interactive segmentation method of the movable contour model cut based on random walk and figure |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106780518B (en) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107292896A (en) * | 2017-08-15 | 2017-10-24 | 电子科技大学 | Contour extraction method based on Snake models |
| EP3668408A1 (en) * | 2017-08-17 | 2020-06-24 | Koninklijke Philips N.V. | Ultrasound system with extraction of image planes from volume data using touch interaction with an image |
| CN107818566A (en) * | 2017-10-27 | 2018-03-20 | 中山大学 | A kind of image partition method based on partial structurtes information around random walk combination pixel |
| CN108257118B (en) * | 2018-01-08 | 2020-07-24 | 浙江大学 | Fracture adhesion segmentation method based on normal corrosion and random walk |
| CN110415252B (en) * | 2018-04-26 | 2022-08-05 | 北京连心医疗科技有限公司 | CNN-based periocular organ segmentation method, CNN-based periocular organ segmentation equipment and CNN-based periocular organ segmentation storage medium |
| CN108460783B (en) * | 2018-05-09 | 2019-03-12 | 电子科技大学 | A kind of cerebral magnetic resonance image organizational dividing method |
| CN109191397A (en) * | 2018-08-28 | 2019-01-11 | 广州智美科技有限公司 | Voxel-based hole repairing method and device |
| CN109741439B (en) * | 2018-12-07 | 2023-12-15 | 广州医科大学 | Three-dimensional reconstruction method of two-dimensional MRI fetal image |
| US11423544B1 (en) | 2019-11-14 | 2022-08-23 | Seg AI LLC | Segmenting medical images |
| US10762629B1 (en) | 2019-11-14 | 2020-09-01 | SegAI LLC | Segmenting medical images |
| CN112801851B (en) * | 2021-01-28 | 2021-12-17 | 太原科技大学 | Hardware system for cutting field plant leaves and cutting method thereof |
| CN116759052B (en) * | 2023-06-20 | 2024-04-05 | 华平祥晟(上海)医疗科技有限公司 | Image storage management system and method based on big data |
| CN116630358B (en) * | 2023-07-25 | 2023-09-26 | 潍坊医学院附属医院 | Threshold segmentation method for brain tumor CT image |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103208114A (en) * | 2013-01-25 | 2013-07-17 | 西安电子科技大学 | Stomach adipose tissue extraction method based on interactive segmentation |
| CN104835154A (en) * | 2015-05-03 | 2015-08-12 | 华东理工大学 | Color image object acquisition method based on random walk |
| CN105701832A (en) * | 2016-01-19 | 2016-06-22 | 苏州大学 | PET-CT lung tumor segmentation method combining three dimensional graph cut algorithm with random walk algorithm |
| CN105957063A (en) * | 2016-04-22 | 2016-09-21 | 北京理工大学 | CT image liver segmentation method and system based on multi-scale weighting similarity measure |
-
2017
- 2017-02-10 CN CN201710073688.XA patent/CN106780518B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103208114A (en) * | 2013-01-25 | 2013-07-17 | 西安电子科技大学 | Stomach adipose tissue extraction method based on interactive segmentation |
| CN104835154A (en) * | 2015-05-03 | 2015-08-12 | 华东理工大学 | Color image object acquisition method based on random walk |
| CN105701832A (en) * | 2016-01-19 | 2016-06-22 | 苏州大学 | PET-CT lung tumor segmentation method combining three dimensional graph cut algorithm with random walk algorithm |
| CN105957063A (en) * | 2016-04-22 | 2016-09-21 | 北京理工大学 | CT image liver segmentation method and system based on multi-scale weighting similarity measure |
Non-Patent Citations (2)
| Title |
|---|
| Graph cuts based active contour model with selective local or global segmentation;Zheng Q et al.;《Electronics Letters》;20120930;第48卷(第9期);第961-966页 |
| 基于随机游走的交互式图像分割算法研究;程伟;《中国优秀硕士学位论文全文数据库信息科技辑（月刊）》;20150515(第05期);论文全文 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106780518A (en) | 2017-05-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106780518B (en) | A kind of MR image three-dimensional interactive segmentation method of the movable contour model cut based on random walk and figure | |
| Mharib et al. | Survey on liver CT image segmentation methods | |
| CN107808377B (en) | The positioning device of lesion in a kind of lobe of the lung | |
| CN106204587B (en) | Multiple organ dividing method based on depth convolutional neural networks and region-competitive model | |
| CN104809723B (en) | The three-dimensional CT image for liver automatic division method of algorithm is cut based on super voxel and figure | |
| Gonçalves et al. | Hessian based approaches for 3D lung nodule segmentation | |
| US8571278B2 (en) | System and methods for multi-object multi-surface segmentation | |
| US8358819B2 (en) | System and methods for image segmentation in N-dimensional space | |
| CN107230206A (en) | A kind of 3D Lung neoplasm dividing methods of the super voxel sequence lung images based on multi-modal data | |
| CN102385751B (en) | Liver tumor region segmentation method based on watershed transform and classification through support vector machine | |
| CN102048550B (en) | Method for automatically generating liver 3D (three-dimensional) image and accurately positioning liver vascular domination region | |
| CN107590809A (en) | Lung dividing method and medical image system | |
| Ma et al. | Two graph theory based methods for identifying the pectoral muscle in mammograms | |
| CN103294883B (en) | For the method and system for through the implantation of conduit aorta petal intervene planning | |
| CN102509286B (en) | Target region sketching method for medical image | |
| CN104915950B (en) | A kind of region based on energy constraint increases ultrasonoscopy automatic division method | |
| CN103310458A (en) | Method for elastically registering medical images by aid of combined convex hull matching and multi-scale classification strategy | |
| CN105389811A (en) | Multi-modality medical image processing method based on multilevel threshold segmentation | |
| CN102938027A (en) | Realization method of computer-assisted liver transplantation operation planning system | |
| CN106780497B (en) | A kind of organ vascular tree extraction method based on statistical information | |
| CN109300113A (en) | A kind of Lung neoplasm assisted detection system and method based on improvement Convex Hull Method | |
| CN105139377A (en) | Rapid robustness auto-partitioning method for abdomen computed tomography (CT) sequence image of liver | |
| CN103473805A (en) | Method for measuring size of three-dimensional reconstruction liver model on basis of improved region growing algorithm | |
| CN104616289A (en) | Removal method and system for bone tissue in 3D CT (Three Dimensional Computed Tomography) image | |
| CN107067393A (en) | A kind of three-dimensional medical image segmentation method based on user mutual and shape prior knowledge |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20190605 Address after: 215011 Bamboo Garden Road, Suzhou high tech Zone, Jiangsu Province, No. 209 Applicant after: Suzhou were Medical Technology Co. Ltd. Address before: 215123 No. 199 Renai Road, Suzhou Industrial Park, Suzhou City, Jiangsu Province Applicant before: Soochow University |
|
| TA01 | Transfer of patent application right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20220315 Address after: 341000 Building 2, Ganzhou national high-level talent science and Innovation Park, No. 1, Wudangshan Road, high tech Zone, Zhanggong District, Ganzhou City, Jiangxi Province Patentee after: Jiangxi Bigway Medical Technology Co.,Ltd. Address before: 215011 No. 209 Chuk Yuen Road, hi tech Zone, Jiangsu, Suzhou Patentee before: SUZHOU BIGVISION MEDICAL TECHNOLOGY Co.,Ltd. |
|
| TR01 | Transfer of patent right |