WO2000079481A1 - Mra segmentation using active contour models - Google Patents

Mra segmentation using active contour models Download PDF

Info

Publication number
WO2000079481A1
WO2000079481A1 PCT/US2000/017282 US0017282W WO0079481A1 WO 2000079481 A1 WO2000079481 A1 WO 2000079481A1 US 0017282 W US0017282 W US 0017282W WO 0079481 A1 WO0079481 A1 WO 0079481A1
Authority
WO
WIPO (PCT)
Prior art keywords
volume
distance function
segmentation
image data
dimensional image
Prior art date
Application number
PCT/US2000/017282
Other languages
French (fr)
Other versions
WO2000079481A9 (en
Inventor
Liana Lorigo
W. Eric L. Grimson
Olivier Faugeras
Renaud Keriven
Carl-Fredrik Westin
Ronald Kikinis
Original Assignee
Massachusetts Institute Of Technology
Institut National De Recherche En Informatique Et En Automatique
Ecole Nationale Des Ponts Et Chaussees
The Brigham And Women's Hospital, Inc.
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Massachusetts Institute Of Technology, Institut National De Recherche En Informatique Et En Automatique, Ecole Nationale Des Ponts Et Chaussees, The Brigham And Women's Hospital, Inc. filed Critical Massachusetts Institute Of Technology
Priority to EP00943090A priority Critical patent/EP1208535A1/en
Publication of WO2000079481A1 publication Critical patent/WO2000079481A1/en
Publication of WO2000079481A9 publication Critical patent/WO2000079481A9/en

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/12Edge-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/149Segmentation; Edge detection involving deformable models, e.g. active contour models
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/755Deformable models or variational models, e.g. snakes or active contours
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10088Magnetic resonance imaging [MRI]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20112Image segmentation details
    • G06T2207/20161Level set
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30101Blood vessel; Artery; Vein; Vascular

Definitions

  • the invention relates to the field of volumetric three-dimensional image data segmentation, and in particular to MRA segmentation.
  • MRA magnetic resonance angiography
  • MRA magnetic resonance angiography
  • blood vessels appear in MRA images as bright curve-like patterns which may be noisy and have gaps. What is shown is a "maximum intensity projection".
  • the data is a stack of slices where most areas are dark, but vessels tend to be bright. This stack is collapsed into a single image for viewing by performing a projection through the stack that assigns to each pixel in the projection the brightest voxel over all slices. This image shows projections along three orthogonal axes.
  • Thresholding is one possible approach to this segmentation problem and works adequately on the larger vessels.
  • the problem arises in detecting the small vessels. Thresholding cannot be used for the small vessels for several reasons.
  • the voxels may have an intensity that is a combination of the intensities of vessels and background if the vessel is only partially inside the voxel. This sampling artifact is called partial voluming.
  • Other imaging conditions can cause some background areas to be as bright as other vessel areas, complicating threshold selection.
  • the images are often noisy, and methods using local contextual information can be more robust.
  • Mean curvature evolution schemes for segmentation, implemented with level set methods, have become an important approach in computer vision. This approach uses partial differential equations to control the evolution.
  • the invention presents the first implementation of geodesic active contours in 3D.
  • the system and method of the invention use these techniques for automatic segmentation of blood vessels in MRA images.
  • the dimension of the manifold is 1, and its co-dimension is 2.
  • the invention utilizes the fact that the underlying structures in the image are indeed 3D curves and evolves an initial curve into the curves in the data (the vessels).
  • the segmentation techniques of the invention are based on the concept of mean curvature flow, or curve-shortening flow, from the field of differential geometry.
  • the proposed MRA segmentation method uses a mathematical modeling technique that is well-suited to the complicated curve-like structure of blood vessels.
  • the segmentation task is defined as an energy minimization over all 3D curves and uses a level set method to search for a solution.
  • the approach is an extension of previous level set segmentation techniques to higher co- dimension.
  • the process proceeds by iteratively updating v according to
  • v, ⁇ (Vv(x,t), ⁇ 7 v(x,t)) + — Vv(x,t) -Hi — r , the updating terminates at convergence or as g
  • S' is then defined to be the zero level set of the current distance function v' and reinitialize v' to be a distance function to S'.
  • the volume is continually iteratively updated such that a final distance function v is obtained.
  • the first output obtained from this volume is a segmentation of vessels in the MRA data, obtained by computing the zero level set of v.
  • FIG. 1 is a maximum intensity projection of a phase-contrast MRA image of blood vessels in a brain
  • FIG. 2 illustrates level sets of an embedding function u, for a closed curve in R 2 ;
  • FIG. 3 illustrates a single sequence showing eight successive stages of a tubular object undergoing mean curvature flow;
  • FIG. 4A shows a curve having a tangent to C at p, the normal plane, the image-based vector, and its projection onto the normal plane
  • FIG. 4B shows a curve using the ⁇ -level set method
  • FIGs. 5 A-C illustrate an evolving helix under mean curvature flow with additional vector field: target curve, initial level set, level set after evolution with endpoints constrained, respectively;
  • FIG. 6 is an operational block diagram of a MR segmentation system utilizing the invention
  • FIG. 7 is a flowchart of the segmentation algorithm in accordance with the invention
  • FIG. 8 is a flowchart showing the details of the surface and volume initialization portion of the algorithm in accordance with the invention.
  • FIGs. 9A-D illustrate a vertical bar evolving into a segmentation of a first dataset
  • FIG. 10 is a single 3D dataset, the first image in each row is the maximum intensity projection of the raw data, and the second image is the segmentation result from three orthogonal viewpoints;
  • FIG. 11 illustrates an image of a partial segmentation of the first dataset in FIG. 10, the colorscale is continuous from darkest to lightest intensities, with darkest indicating a radius of curvature ⁇ 1mm and lightest indicating a radius of curvature >2mm.
  • mean curvature flow refers to some curve evolving in time so that at each point, the velocity vector normal to the curve is equal to the mean curvature vector.
  • This concept is normally defined for arbitrary generic surfaces, but only curves are necessary for the invention, so the definition has been restricted for purposes of illustration. More formally, let C Constant t > 0 be a family of curves in R 2 or R 3 , N the normal for a given orientation. That is, C is a curve, and t represents the time parameter or the index into the family of curves, not position.
  • the mean curvature flow equation is then given by the vector equation
  • FIG. 2 illustrates level sets of an embedding function u, for a closed curve in R 2 .
  • C is a curve in 3D.
  • FIG. 3 demonstrates this evolution by illustrating evolving curves under mean curvature flow.
  • FIG 3 illustrates a single sequence showing eight successive stages of a tubular object undergoing curve-shortening flow (mean curvature flow), where the curve is the centerline of the tubular shape. The bumps are first smoothed out until the shape approximates a torus, then the torus shrinks to a point.
  • FIG. 4A shows a curve having a tangent to C at p, the normal plane, the image-based vector, and its projection onto the normal plane.
  • FIGs. 5A-C show how underlying image information can attract the evolving tube.
  • FIGs. 5A-C illustrate an evolving helix under mean curvature flow with additional vector field: target curve, initial level set, level set after evolution with endpoints constrained, respectively.
  • the underlying volumetric image data is shown, as a maximum intensity projection in FIG. 5 A. This volume was generated by drawing a cosine curve in the volume, then smoothing with a Gaussian filter.
  • FIG. 5B shows the initial curve, a helix.
  • FIG. 6 is an operational block diagram of a MR segmentation system 600 utilizing the invention.
  • the system includes a conventional MR scanner 602 running in conjunction with a MR computer that generates and stores raw MRA image slices.
  • the MR data is then segmented by a computer 606.
  • a flowchart of the segmentation algorithm in accordance with the invention is shown in FIG. 7.
  • the invention produces 3D surface models of blood vessels based on magnetic resonance angiography (MRA) data.
  • MRA magnetic resonance angiography
  • the patient's head is imaged in a magnetic resonance scanner 602.
  • the image produced is a three-dimensional image. This means that it is actually a stack of many (often
  • the head is this region.
  • the invention utilizes a computer 606 to generate a 3D surface model of the blood vessels in the head, based on this 3D image.
  • This surface model could be displayed and manipulated on a standard computer.
  • the surface model can be viewed by clinicians, radiologists, and other persons.
  • the surface model is often preferable to the raw 3D image in the areas of ease of interpretation, ease of further measurements, incorporation with other anatomical information, and other areas.
  • Blood vessels appear in MRA images as bright curve-like patterns that may be noisy and have gaps.
  • the data is a stack of slices where most areas are dark, but vessels tend to be bright. This stack is collapsed into a single image for viewing by performing a projection through the stack that assigns to each pixel in the projection the brightest voxel over all slices.
  • MRA segmentation is complicated by the presence of imaging artifacts which appear visually similar to true vessel structures and also to partial voluming, the case of a small imaging area having a brightness value that is a combination of the brightness values of vessels and of background because the vessel is only partially inside the area imaged.
  • This specific segmentation problem is part of the high-level problem of developing computerized techniques for the analysis of medical images.
  • Automatic and semi-automatic techniques can potentially assist clinicians and radiologists, saving them much of the time required to manually segment large datasets, or more generally facilitating measurements and interpretation of the images.
  • banding will now be described. Instead of evolving the entire volume, only the portion of the volume within a narrow band of the zero level set is evolved (the current surface). Normally, the band is set to include voxels that are up to 6 voxels away from the surface.
  • the advantage of this technique is efficiency, and the disadvantage is that structures that are outside the band may be missed if the potential function g does not have a large enough capture range to attract the segmentation to these structures.
  • the interpretation of banding is different from that in previous level set methods; they propagate image information from the zero level set to the rest of the bands, while the invention uses image information at each point.
  • the algorithm used in the method of the invention is described.
  • the 3D MRA volume is loaded into the computer.
  • the vessels appear bright, background appears dark.
  • An initial surface S is generated either by thresholding the inputted data set or by using a previously generated surface.
  • a signed distance function to S is then generated.
  • This distance function is a 3D volume, v.
  • FIG. 8 The detailed flowchart of the surface and volume initialization is shown in FIG. 8.
  • the method continues by iteratively updating v according to Equation 5.
  • the algorithm may also incorporate the image scaling term previously described and/or an orientation term.
  • the process terminates at convergence or as determined by the user.
  • S is redefined periodically or when needed, to be the zero level set of the current distance function v' and reinitialize v' to be the distance function to S'.
  • the process continues by updating the volume according to the previous step.
  • the loop above yields a final distance function v.
  • the first output obtained from this volume is a segmentation of the blood vessels, obtained by computing the zero level set of v.
  • the centerlines are obtained as the local minima of the distance function.
  • estimates of vessel diameter are obtained as a by-product of the computation of ⁇ in
  • Equation 5 To control the trade-off between fitting the surface to the image data and enforcing the smoothness constraint on the surface, an image scaling term imscale is added to Equation
  • v, ⁇ (Vv(x ),V 2 v(x,t)) + imscale *(VvV I) * ⁇ Vv(x,t)- H ⁇ -. (7) g
  • For example, if the two vectors point in the same direction, then the brighter region is inside the surface and the darker region is outside. The angle between the vectors is 0, whose cosine is 1 , so the image term is fully counted. However, if they point in opposite directions, the negative weighting prevents the evolving vessel walls from being attracted to image gradients that point in the opposite direction. As customary in level set segmentation methods, the volume v is periodically reinitialized to be a distance function.
  • each point in the volume is set to be its distance to S.
  • this reinitialization is itself a level set method.
  • the surface is propagated outward at constant speed of 1, and the distance at each point is determined to be the time at which the surface crossed that point.
  • a second step propagates the surface inward to obtain the negative distances analogously.
  • FIG. 8 shows a flowchart of a more detailed portion of the algorithm used to generate the initial surface.
  • This initial surface (and thus the initial volume) is normally generated by thresholding the MRA dataset. However, the method does not require that the initial surface be near the target surface but may use any initial surface.
  • FIGs. 9A-D illustrate a vertical bar evolving into the segmentation of the first dataset in FIG. 11.
  • the datasets may be pre-processed to reduce noise and smooth.
  • the segmentations are post-processed to remove any surface patches whose surface area is less than some threshold (a parameter of the method) to eliminate patches corresponding to noise in the original data.
  • the larger principal curvature can be useful in measuring the radii of the vessels for a particular application, since radius is the inverse of curvature. This curvature can be easily computed when the smaller principal curvature is computed for the segmentation.
  • the option to color-code the segmentations can be added based on vessel radii, as estimated from the local larger principal curvature of the tubular surface.
  • FIG. 10 is a single 3D dataset, the first image in each row is the maximum intensity projection of the raw data, and the second image is the segmentation result from three orthogonal viewpoints.
  • a cleaning threshold c indicated the minimum surface area of connected components of the surface to be retained in the post-processing cleaning step.
  • FIG. 10 shows a single 3D data set (in maximum intensity projection) and segmentation result from three orthogonal viewpoints. The result shown is obtained by thresholding the raw data set. Finally, the capability to color-code the vasculature surface is demonstrated based on local curvature. With reference to FIG. 11, it will be appreciated that for a ribbon-like vessel, the flatter sides shows a large radius, and the sharply curved edges show a small radius. In this image of a partial segmentation of the first dataset in FIG. 10, the colorscale is continuous from darkest to lightest intensities, with darkest indicating a radius of curvature ⁇ 1mm and lightest indicating a radius of curvature >2mm. The curvatures output by the evolution have been smoothed by a 3x3x3 filter prior to coloring the surface.

Abstract

A method and system of providing segmentation of volumetric three-dimensional image data set such as MRA images. Initially, a three-dimensional MRA volume is input. An initial surface S is then generated by thresholding the input. A signed distance function v to S is then generated, where v=v(x,t) and S is the zero level set of v(.,0). The process proceeds by iteratively updating v according to Formula (I) , the updating terminates at convergence or as determined by an operator. S' is then defined to be the zero level set of the current distance function v' and reinitialize v' to be a distance function to S'. The volume is continually iteratively updated such that a final distance function v is obtained. The first output obtained from this volume is a segmentation of vessels in the MRA data, obtained by computing the zero level set of v.

Description

MRA SEGMENTATION USING ACTIVE CONTOUR MODELS
PRIORITY INFORMATION
This application claims priority from provisional application Ser. No. 60/140,609 filed June 23, 1999 and incorporated herein by reference in its entirety.
BACKGROUND OF THE INVENTION
The invention relates to the field of volumetric three-dimensional image data segmentation, and in particular to MRA segmentation. Automatic and semi-automatic magnetic resonance angiography (MRA) segmentation techniques can potentially save radiologists large amounts of time required for manual segmentation and can facilitate further data analysis. It is desirable to develop computer vision techniques for the segmentation of medical images. Specifically, the segmentation of volumetric vasculature images is considered, such as the magnetic resonance angiography (MRA) image pictured in FIG. 1.
As shown in FIG. 1, blood vessels appear in MRA images as bright curve-like patterns which may be noisy and have gaps. What is shown is a "maximum intensity projection". The data is a stack of slices where most areas are dark, but vessels tend to be bright. This stack is collapsed into a single image for viewing by performing a projection through the stack that assigns to each pixel in the projection the brightest voxel over all slices. This image shows projections along three orthogonal axes.
Thresholding is one possible approach to this segmentation problem and works adequately on the larger vessels. The problem arises in detecting the small vessels. Thresholding cannot be used for the small vessels for several reasons. The voxels may have an intensity that is a combination of the intensities of vessels and background if the vessel is only partially inside the voxel. This sampling artifact is called partial voluming. Other imaging conditions can cause some background areas to be as bright as other vessel areas, complicating threshold selection. Finally, the images are often noisy, and methods using local contextual information can be more robust. Mean curvature evolution schemes for segmentation, implemented with level set methods, have become an important approach in computer vision. This approach uses partial differential equations to control the evolution. The fundamental concepts from mathematics from which mean curvature schemes derive were explored several years earlier when smooth closed curves in 2D were proven to shrink to a point under mean curvature motion. Mean curvature flow of any hypersurface was framed as a level set problem. For application to image segmentation, a vector field was induced on the embedding space, so that the evolution could be controlled by an image gradient field or other image data. The same results of existence, uniqueness, and stability of viscosity solutions were obtained for the modified evolution equations for the case of planar curves, and experiments on real- world images demonstrated the effectiveness of the approach.
Curves evolving in the plane became surfaces evolving in space, called "minimal surfaces". Although the theorem on planar curves shrinking to a point could not be extended to the case of surfaces evolving in 3D, the existence, uniqueness, and stability results of the level set formalism held analogously to the 2D case. Thus, the method was feasible for evolving both curves in 2D and surfaces in 3D. Beyond elegant mathematics, spectacular results on real- world data sets established the method as an important segmentation tool in both domains. One fundamental limitation to these schemes has been that they describe only the flow of hypersurfaces, i.e., surfaces of co-dimension 1.
The problem of curve-shortening flow for 3D curves has been studied, and the level set technique has been generalized to arbitrary manifolds in arbitrary dimension. Analogous results were provided and extend the level set evolution equation to account for an additional vector field induced on the space.
SUMMARY OF THE INVENTION
Accordingly, the invention presents the first implementation of geodesic active contours in 3D. Specifically, the system and method of the invention use these techniques for automatic segmentation of blood vessels in MRA images. The dimension of the manifold is 1, and its co-dimension is 2.
The invention utilizes the fact that the underlying structures in the image are indeed 3D curves and evolves an initial curve into the curves in the data (the vessels). In particular, the segmentation techniques of the invention are based on the concept of mean curvature flow, or curve-shortening flow, from the field of differential geometry. The proposed MRA segmentation method uses a mathematical modeling technique that is well-suited to the complicated curve-like structure of blood vessels. The segmentation task is defined as an energy minimization over all 3D curves and uses a level set method to search for a solution.
The approach is an extension of previous level set segmentation techniques to higher co- dimension.
The invention thus sets forth a method of providing segmentation of volumetric three-dimensional image data such as MRA images. Initially, a three-dimensional MRA volume is provided. An initial surface S is then generated by thresholding input. A signed distance function v to S is then generated, where v=v(x,t) and S is the zero level set of v(-,0).
The process proceeds by iteratively updating v according to
v, = λ(Vv(x,t),\7 v(x,t)) + — Vv(x,t) -Hi — r , the updating terminates at convergence or as g |W| determined by an operator. S' is then defined to be the zero level set of the current distance function v' and reinitialize v' to be a distance function to S'. The volume is continually iteratively updated such that a final distance function v is obtained. The first output obtained from this volume is a segmentation of vessels in the MRA data, obtained by computing the zero level set of v.
BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a maximum intensity projection of a phase-contrast MRA image of blood vessels in a brain;
FIG. 2 illustrates level sets of an embedding function u, for a closed curve in R2; FIG. 3 illustrates a single sequence showing eight successive stages of a tubular object undergoing mean curvature flow;
FIG. 4A shows a curve having a tangent to C at p, the normal plane, the image-based vector, and its projection onto the normal plane; FIG. 4B shows a curve using the ε-level set method; FIGs. 5 A-C illustrate an evolving helix under mean curvature flow with additional vector field: target curve, initial level set, level set after evolution with endpoints constrained, respectively;
FIG. 6 is an operational block diagram of a MR segmentation system utilizing the invention; FIG. 7 is a flowchart of the segmentation algorithm in accordance with the invention;
FIG. 8 is a flowchart showing the details of the surface and volume initialization portion of the algorithm in accordance with the invention;
FIGs. 9A-D illustrate a vertical bar evolving into a segmentation of a first dataset;
FIG. 10 is a single 3D dataset, the first image in each row is the maximum intensity projection of the raw data, and the second image is the segmentation result from three orthogonal viewpoints;
FIG. 11 illustrates an image of a partial segmentation of the first dataset in FIG. 10, the colorscale is continuous from darkest to lightest intensities, with darkest indicating a radius of curvature < 1mm and lightest indicating a radius of curvature >2mm.
DETAILED DESCRIPTION OF THE INVENTION
Intuitively, mean curvature flow refers to some curve evolving in time so that at each point, the velocity vector normal to the curve is equal to the mean curvature vector. This concept is normally defined for arbitrary generic surfaces, but only curves are necessary for the invention, so the definition has been restricted for purposes of illustration. More formally, let C„ t > 0 be a family of curves in R2 or R3, N the normal for a given orientation. That is, C is a curve, and t represents the time parameter or the index into the family of curves, not position. The mean curvature flow equation is then given by the vector equation
C, = κN (1) with given initial curve C(0)=C0, K the curvature of the curve, and Ct the time derivative of the curve. It will be appreciated that only ID curves are considered here, as opposed to evolving surfaces, the mean curvature is just the usual curvature of the curve. This motion is also called "curve-shortening flow" since it is the solution, obtained by Euler-Lagrange equations, to the problem of minimizing curve length:
Figure imgf000006_0001
where p is the spatial parameter of the curve.
The basic idea of the level set method is given to evolve a planar curve C. Define a function u: R2 → R so that C is a level-set of u. The convention that C is, in particular, the zero level set of u, although this choice is not necessary for the method. The function u is now an implicit representation of the curve C. The advantages of this representation are that it is intrinsic (independent of parameterization) and that it is topologically flexible since different topologies of C are represented by the constant topology of u. Let C0 be the initial curve.
Evolving C according to
C, = βN (2) with initial condition C(-,0) = C0(-) for any function β, is equivalent to evolving u according to
Figure imgf000007_0001
with initial condition u(-,0) = Uo(-) and u0(C0) = 0.
This result is independent of the choice of function u. As is customary in the literature, u to be the signed distance function to the curve C as shown in FIG. 2. FIG. 2 illustrates level sets of an embedding function u, for a closed curve in R2.
The level set evolution equations that follow enable one to evolve space curves, with evolution driven by both mean curvature and image information. In the following discussion, C is a curve in 3D.
Let v:R3—»[0,∞) be an auxiliary function whose zero level set is identically C, that is smooth near C, and such that Vv is non-zero outside C. For a nonzero vector qeR", define
Figure imgf000007_0002
as the projector onto the plane normal to q. Further define λ(Vv(x,t)N2v(x,t)) as the smaller nonzero eigenvalue of PVvV2v PVv. The level set evolution equation is then vt=λ(Vv(x,t)N2v(x,t)). That is, this evolution is equivalent to evolving C according to Ct=κN in the sense that C is the zero level set of v throughout the evolution.
FIG. 3 demonstrates this evolution by illustrating evolving curves under mean curvature flow. FIG 3 illustrates a single sequence showing eight successive stages of a tubular object undergoing curve-shortening flow (mean curvature flow), where the curve is the centerline of the tubular shape. The bumps are first smoothed out until the shape approximates a torus, then the torus shrinks to a point.
The situation where there is an underlying vector field driving the evolution, in combination with the curvature term will now be described. Assume the desired evolution equation is of the form Ct=κN-πd where π is the projection operator onto the normal space of C (which is a vector space of dimension 2) and d is a given vector field in R3. The evolution equation for the embedding space then becomes vt=λ(Vv,V2v)+Vvd.
For the case of ID structures in 3D images, it is desirable to minimize
Figure imgf000008_0001
where C(p):[0,l]→R3 is the ID curve, l:[0,a]χ [0,b]χ[0,c]→[0,∞) is the image, and g:[0,∞)-»R+ is a strictly decreasing function such that g(r)-»0 as r→∞. For the current implementation, g(r)=exp(-r) is used because it works well in practice. Another common choice is g(|VI|)=l/[l+|V I|2. By computing the Euler-Lagrange equations, it is found that the curve evolution equation is
c, = κN-^u(H^- (4)
where Η is the Hessian of the intensity function. The second term in the above equation is illustrated in FIG. 4A, which shows a curve having a tangent to C at p, the normal plane, the image-based vector, and its projection onto the normal plane. FIG. 4B is a curve using the ε-level set method. That is, s' VI d = £-H- g lV/l so the equation for the embedding space is
Figure imgf000008_0002
Thus, the mean curvature flow and level set methods have now been used to segment ID structures in 3D. FIGs. 5A-C show how underlying image information can attract the evolving tube. FIGs. 5A-C illustrate an evolving helix under mean curvature flow with additional vector field: target curve, initial level set, level set after evolution with endpoints constrained, respectively. The underlying volumetric image data is shown, as a maximum intensity projection in FIG. 5 A. This volume was generated by drawing a cosine curve in the volume, then smoothing with a Gaussian filter. FIG. 5B shows the initial curve, a helix.
The result of the evolution is shown in FIG. 5C.
Segmentation refers to the task of separating foreground regions of an image from background regions. The system and method for segmentation of vessels from MRA using the described level set method in accordance with the invention will now be described. FIG. 6 is an operational block diagram of a MR segmentation system 600 utilizing the invention.
The system includes a conventional MR scanner 602 running in conjunction with a MR computer that generates and stores raw MRA image slices. The MR data is then segmented by a computer 606. A flowchart of the segmentation algorithm in accordance with the invention is shown in FIG. 7.
The invention produces 3D surface models of blood vessels based on magnetic resonance angiography (MRA) data. A scenario in which one could use the invention is illustrated in FIG. 6. Consider a patient who is suspected to have abnormalities in blood vessels in the brain, such as an aneuryism or other. Or perhaps one would like to view the blood vessels for study purposes.
The patient's head is imaged in a magnetic resonance scanner 602. The image produced is a three-dimensional image. This means that it is actually a stack of many (often
60) two-dimensional images, each of which pictures a "slice" or cross-section of the anatomical region imaged and stored by MR computer 604. In this example, the head is this region.
The invention utilizes a computer 606 to generate a 3D surface model of the blood vessels in the head, based on this 3D image. This surface model could be displayed and manipulated on a standard computer. The surface model can be viewed by clinicians, radiologists, and other persons. The surface model is often preferable to the raw 3D image in the areas of ease of interpretation, ease of further measurements, incorporation with other anatomical information, and other areas.
Blood vessels appear in MRA images as bright curve-like patterns that may be noisy and have gaps. The data is a stack of slices where most areas are dark, but vessels tend to be bright. This stack is collapsed into a single image for viewing by performing a projection through the stack that assigns to each pixel in the projection the brightest voxel over all slices.
The task of MRA segmentation is complicated by the presence of imaging artifacts which appear visually similar to true vessel structures and also to partial voluming, the case of a small imaging area having a brightness value that is a combination of the brightness values of vessels and of background because the vessel is only partially inside the area imaged.
This specific segmentation problem is part of the high-level problem of developing computerized techniques for the analysis of medical images. Automatic and semi-automatic techniques can potentially assist clinicians and radiologists, saving them much of the time required to manually segment large datasets, or more generally facilitating measurements and interpretation of the images.
The ε-level set method is now described. Since the projection operator Pq is defined only for non-zero vectors q, the method is undefined at Vv=0, which is the curve itself, and is numerically unstable near the curve. For this reason, v is regarded as a distance function to a "tube" of small radius ε around the curve, instead of extracting the true ID curve. That is, ε-level set is evolved instead of evolving the true curve as shown in FIG. 4B. It will be appreciated that ε does not denote a fixed value here. The evolving shape is a "tubular" surface of some (unspecified and variable) nonzero width. In addition to being more robust, this method better captures the geometry of blood vessels, which have nonzero diameter.
Banding will now be described. Instead of evolving the entire volume, only the portion of the volume within a narrow band of the zero level set is evolved (the current surface). Normally, the band is set to include voxels that are up to 6 voxels away from the surface. The advantage of this technique is efficiency, and the disadvantage is that structures that are outside the band may be missed if the potential function g does not have a large enough capture range to attract the segmentation to these structures. The interpretation of banding is different from that in previous level set methods; they propagate image information from the zero level set to the rest of the bands, while the invention uses image information at each point.
A description of the use of curvature instead of eigenvalues is now provided. For computational efficiency and because of numerical instability of the gradient computations and thus the evolution equation near Vv=0, the level sets of the function v flow in the direction of the normal with velocity equal to the sum of their smaller principal curvature and the dot product of Vv with the image-based vector field d. Therefore, the smaller curvature is computed directly from v instead of as an eigenvalue of PVvV2v PVv, and the invention uses λ=| Vv| p2, where p2 is the smaller curvature.
With reference now to FIGs. 7 and 8, the algorithm used in the method of the invention is described. Initially, the 3D MRA volume is loaded into the computer. The vessels appear bright, background appears dark. An initial surface S is generated either by thresholding the inputted data set or by using a previously generated surface. A signed distance function to S is then generated. This distance function is a 3D volume, v. The volume v is a function of spatial location x. Since it will be repeatedly updated over time, it can also be considered a function of time. That is, v=v(x,t) and S is the zero level set of v(-,0). The detailed flowchart of the surface and volume initialization is shown in FIG. 8.
The method continues by iteratively updating v according to Equation 5. The algorithm may also incorporate the image scaling term previously described and/or an orientation term. The process terminates at convergence or as determined by the user.
S is redefined periodically or when needed, to be the zero level set of the current distance function v' and reinitialize v' to be the distance function to S'. The process continues by updating the volume according to the previous step. The loop above yields a final distance function v. The first output obtained from this volume is a segmentation of the blood vessels, obtained by computing the zero level set of v. Second, the centerlines are obtained as the local minima of the distance function. Third, estimates of vessel diameter are obtained as a by-product of the computation of λ in
Equation 5. To control the trade-off between fitting the surface to the image data and enforcing the smoothness constraint on the surface, an image scaling term imscale is added to Equation
5 to obtain σ' \7I v, = λ(Vv(x,t),V2v(x,t)) + imscale * — Vv(x,t) - H (6) g VI imscale being set by the user or pre-set to a default value. With respect to gradient directionality, because vessels appear brighter than the background, the image term is weighted by the cosine of the angle between the normal to the surface and the gradient in the image. This cosine is given by the dot product of the respective gradients of v and I, so the update equation becomes
v, = λ(Vv(x ),V2v(x,t)) + imscale *(VvV I) *^Vv(x,t)- H^-. (7) g |V/| For example, if the two vectors point in the same direction, then the brighter region is inside the surface and the darker region is outside. The angle between the vectors is 0, whose cosine is 1 , so the image term is fully counted. However, if they point in opposite directions, the negative weighting prevents the evolving vessel walls from being attracted to image gradients that point in the opposite direction. As customary in level set segmentation methods, the volume v is periodically reinitialized to be a distance function. The zero level set S is extracted, then each point in the volume is set to be its distance to S. For the implementation of the invention, this reinitialization is itself a level set method. To obtain the positive distances, the surface is propagated outward at constant speed of 1, and the distance at each point is determined to be the time at which the surface crossed that point. A second step propagates the surface inward to obtain the negative distances analogously. With respect to the initial surface, FIG. 8 shows a flowchart of a more detailed portion of the algorithm used to generate the initial surface. This initial surface (and thus the initial volume) is normally generated by thresholding the MRA dataset. However, the method does not require that the initial surface be near the target surface but may use any initial surface. FIGs. 9A-D illustrate a vertical bar evolving into the segmentation of the first dataset in FIG. 11.
As shown in FIG. 8, the datasets may be pre-processed to reduce noise and smooth. For the results presented here, the raw datasets were convolved with an isotropic Gaussian ofσ=0.5.
The segmentations are post-processed to remove any surface patches whose surface area is less than some threshold (a parameter of the method) to eliminate patches corresponding to noise in the original data.
With respect to vessel radii estimation, the larger principal curvature can be useful in measuring the radii of the vessels for a particular application, since radius is the inverse of curvature. This curvature can be easily computed when the smaller principal curvature is computed for the segmentation. The option to color-code the segmentations can be added based on vessel radii, as estimated from the local larger principal curvature of the tubular surface.
The results of the segmentation on four datasets in accordance with the invention are demonstrated in FIG. 10. FIG. 10 is a single 3D dataset, the first image in each row is the maximum intensity projection of the raw data, and the second image is the segmentation result from three orthogonal viewpoints.
The following parameters were used in these experiments. All settings were chosen empirically. For the method, imscale varied across the datasets depending on the noise present. A cleaning threshold c indicated the minimum surface area of connected components of the surface to be retained in the post-processing cleaning step.
FIG. 10 shows a single 3D data set (in maximum intensity projection) and segmentation result from three orthogonal viewpoints. The result shown is obtained by thresholding the raw data set. Finally, the capability to color-code the vasculature surface is demonstrated based on local curvature. With reference to FIG. 11, it will be appreciated that for a ribbon-like vessel, the flatter sides shows a large radius, and the sharply curved edges show a small radius. In this image of a partial segmentation of the first dataset in FIG. 10, the colorscale is continuous from darkest to lightest intensities, with darkest indicating a radius of curvature < 1mm and lightest indicating a radius of curvature >2mm. The curvatures output by the evolution have been smoothed by a 3x3x3 filter prior to coloring the surface.
Although the present invention has been shown and described with respect to several preferred embodiments thereof, various changes, omissions and additions to the form and detail thereof, may be made therein, without departing from the spirit and scope of the invention.
What is claimed is:

Claims

A method of providing segmentation of volumetric three-dimensional image data set comprising: (a) inputting a volume of said volumetric three-dimensional image data set; (b) generating an initial surface S by thresholding said input volume; (c) generating a signed distance function v to S, where v=v(x,t) and S is the zero level set of v(-,0);
(d) iteratively updating v according to v, = λ( Vv(x, t)N v(x,t)) -\ — Vv(x, t) • H g |W| said updating terminates at convergence or as determined by an operator; (e) defining S' to be the zero level set of the current distance function v' and reinitialize v' to be a distance function to S'; and (f) continuing updating the volume according to step (d) such that a final distance function v is obtained, wherein the first output obtained from said volume is a segmentation of vessels in the three- dimensional image data set, obtained by computing the zero level set of v.
2. The method of claim 1 , wherein the function g in the equation of step (d) can be any function such that g:[0,∞)— »R+ is strictly decreasing and g(r)→0 as r— »∞.
3. The method of claim 1, wherein an image scaling term imscale is added to the
equation of step (d) so that v, = λ(Vv(x,t),V2v(x,t)) + imscale *
Figure imgf000014_0001
4. The method of claim 3, wherein an orientation term costerm = -Vv-VI, so that
v, = λ(Vv(x,t),V2v(x,t))- imscale *(Vv I) *^-Vv(x,t) - H^- g V/
5. The method of claim 4 further comprising varying the frequency of reinitializing v to be the signed distance function to its zero level set.
6. The method of claim 5, wherein any distance function algorithm can be used to reinitialize v during the reinitialization step.
7. The method of claim 6, wherein a predefined shape with a predefined distance function is chosen instead of thresholding the volumetric three-dimensional image data set.
8. The method of claim 7 further comprising extracting the centerlines of the segmented volumetric three-dimensional image data set, said extraction is performed by: extracting the set of points that are inside the segmented surface and that are also local minima of the distance function v; and defining the hierarchical shortest path.
9. The method of claim 8, wherein said step of defining the hierarchical shortest path comprises: (a) defining the set of primary nodes of a graph to be the extracted local minima and also secondary nodes to be the other voxels inside the vessels of the volumetric three- dimensional image data set; and (b) connecting nodes that are adjacent in the three-dimensional volume with the Euclidean distance between the nodes; (c) finding, for each primary node, the shortest path to another primary node; (d) combining connected paths into a single primary node and again finding the shortest path to another primary node; and (e) repeating steps (a)-(d) until entire structure is connected or no more connections are possible within the given segmentation.
10. The method of claim 7 further comprising inferring local diameter information about the vessels of the volumetric three-dimensional image data set, said information being visualized by coloring the segmented surface according to local diameter, said diameter information being obtained directly in the algorithm since it is the inverse of the larger principal curvature of the surface and this curvature is found during the computation of λ in the update rule for v.
11. The method of claim 7 further comprising removing from the segmented surface any surface patches whose area or volume is beneath a predetermined threshold.
12. The method of claim 8 further comprising removing from the segmented surface any surface patches whose area or volume is beneath a predetermined threshold prior to extraction of the centerlines.
13. The method of claim 9 further comprising removing from the segmented surface any surface patches whose area or volume is beneath a predetermined threshold.
14. The method of claim 1, wherein said volumetric three-dimensional image data set comprises MRA data.
15. A system of providing segmentation of volumetric three-dimensional image data set comprising: (a) means for inputting a volume of said volumetric three-dimensional image data set; (b) means for generating an initial surface S by thresholding said input volume; (c) means for generating a signed distance function v to S, where v=v(x,t) and S is the zero level set of v(-,0); Cd) means for iteratively updating v according to
v, = λ(Vv(x,t),V 7 v(x,t)) + , said updating terminates at convergence or
Figure imgf000016_0001
as determined by an operator; (e) means for defining S' to be the zero level set of the current distance function v' and reinitialize V to be a distance function to S'; and (f) means for continuing updating the volume according to (d) such that a final distance function v is obtained, wherein the first output obtained from said volume is a segmentation of vessels in the three- dimensional image data set, obtained by computing the zero level set of v.
PCT/US2000/017282 1999-06-23 2000-06-23 Mra segmentation using active contour models WO2000079481A1 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
EP00943090A EP1208535A1 (en) 1999-06-23 2000-06-23 Mra segmentation using active contour models

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US14060999P 1999-06-23 1999-06-23
US60/140,609 1999-06-23

Publications (2)

Publication Number Publication Date
WO2000079481A1 true WO2000079481A1 (en) 2000-12-28
WO2000079481A9 WO2000079481A9 (en) 2002-06-06

Family

ID=22492009

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US2000/017282 WO2000079481A1 (en) 1999-06-23 2000-06-23 Mra segmentation using active contour models

Country Status (2)

Country Link
EP (1) EP1208535A1 (en)
WO (1) WO2000079481A1 (en)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002052509A1 (en) * 2000-12-22 2002-07-04 Koninklijke Philips Electronics N.V. Method of analyzing a data set comprising a volumetric representation of an object to be examined
WO2003034337A2 (en) * 2001-10-16 2003-04-24 Koninklijke Philips Electronics N.V. Method for branch selection for probe alignment
WO2004079654A2 (en) * 2003-03-07 2004-09-16 Martin, Philip, John Image processing apparatus and methods
US6845260B2 (en) 2001-07-18 2005-01-18 Koninklijke Philips Electronics N.V. Automatic vessel indentification for angiographic screening
WO2005027053A1 (en) * 2003-09-18 2005-03-24 Politecnico Di Milano Method for determining the three-dimensional surface of an object
WO2006121410A1 (en) * 2005-05-11 2006-11-16 Agency For Science, Technology And Research Method, apparatus and computer software for segmenting the brain from mr data
CN102651130A (en) * 2012-03-30 2012-08-29 宋怡 Level set image processing method
US8922552B2 (en) 2003-01-15 2014-12-30 Koninklijke Philips N.V. Image processing method for automatic adaptation of 3-D deformable model onto a substantially tubular surface of a 3-D object
CN106373097A (en) * 2016-08-29 2017-02-01 合肥康胜达智能科技有限公司 Image processing method
CN112085723A (en) * 2020-09-09 2020-12-15 哈尔滨市科佳通用机电股份有限公司 Automatic detection method for spring jumping fault of truck bolster

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5920319A (en) * 1994-10-27 1999-07-06 Wake Forest University Automatic analysis in virtual endoscopy
US6058218A (en) * 1997-11-10 2000-05-02 General Electric Company Enhanced visualization of weak image sources in the vicinity of dominant sources

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5920319A (en) * 1994-10-27 1999-07-06 Wake Forest University Automatic analysis in virtual endoscopy
US6058218A (en) * 1997-11-10 2000-05-02 General Electric Company Enhanced visualization of weak image sources in the vicinity of dominant sources

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
KIM D K ET AL: "BOUNDARY SEGMENTATION AND MODEL FITTING USING AFFINE ACTIVE SURFACEMODEL", IEEETENCON - DIGITAL SIGNAL PROCESSING APPLICATIONS,US,NEW YORK, NY: IEEE, 26 November 1996 (1996-11-26), pages 147 - 150, XP000781671, ISBN: 0-7803-3680-1 *
SNELL J W ET AL: "MODEL-BASED SEGMENTATION OF THE BRAIN FROM 3-D MRI USING ACTIVE SURFACES", PROCEEDINGS OF THE NORTHEAST BIOENGINEERING CONFERENCE,US,NEW YORK, IEEE, vol. CONF. 19, 18 March 1993 (1993-03-18), pages 164 - 165, XP000399775 *
VERDOCK B ET AL: "ACCURATE SEGMENTATION OF BLOOD VESSELS FROM 3D MEDICAL IMAGES", PROCEEDINGS OF THE INTERNATIONAL CONFERENCE ON IMAGE PROCESSING (ICIP),US,NEW YORK, IEEE, 16 September 1996 (1996-09-16), pages 311 - 314, XP000704040, ISBN: 0-7803-3259-8 *

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002052509A1 (en) * 2000-12-22 2002-07-04 Koninklijke Philips Electronics N.V. Method of analyzing a data set comprising a volumetric representation of an object to be examined
US6845260B2 (en) 2001-07-18 2005-01-18 Koninklijke Philips Electronics N.V. Automatic vessel indentification for angiographic screening
WO2003034337A2 (en) * 2001-10-16 2003-04-24 Koninklijke Philips Electronics N.V. Method for branch selection for probe alignment
WO2003034337A3 (en) * 2001-10-16 2004-11-11 Koninkl Philips Electronics Nv Method for branch selection for probe alignment
US8922552B2 (en) 2003-01-15 2014-12-30 Koninklijke Philips N.V. Image processing method for automatic adaptation of 3-D deformable model onto a substantially tubular surface of a 3-D object
WO2004079654A2 (en) * 2003-03-07 2004-09-16 Martin, Philip, John Image processing apparatus and methods
WO2004079654A3 (en) * 2003-03-07 2004-10-28 Martin Weber Image processing apparatus and methods
WO2005027053A1 (en) * 2003-09-18 2005-03-24 Politecnico Di Milano Method for determining the three-dimensional surface of an object
US7626584B2 (en) 2003-09-18 2009-12-01 Politecnico Di Milano Method for determining the three-dimensional surface of an object
WO2006121410A1 (en) * 2005-05-11 2006-11-16 Agency For Science, Technology And Research Method, apparatus and computer software for segmenting the brain from mr data
CN102651130A (en) * 2012-03-30 2012-08-29 宋怡 Level set image processing method
CN106373097A (en) * 2016-08-29 2017-02-01 合肥康胜达智能科技有限公司 Image processing method
CN112085723A (en) * 2020-09-09 2020-12-15 哈尔滨市科佳通用机电股份有限公司 Automatic detection method for spring jumping fault of truck bolster

Also Published As

Publication number Publication date
WO2000079481A9 (en) 2002-06-06
EP1208535A1 (en) 2002-05-29

Similar Documents

Publication Publication Date Title
Lorigo et al. Co-dimension 2 geodesic active contours for MRA segmentation
Tsai et al. Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification
Cremers et al. Motion competition: A variational approach to piecewise parametric motion segmentation
JP3902765B2 (en) Various approaches to left ventricular segmentation in MR heart images
Lorigo et al. Curves: Curve evolution for vessel segmentation
US7095890B2 (en) Integration of visual information, anatomic constraints and prior shape knowledge for medical segmentations
Montagnat et al. Anisotropic filtering for model-based segmentation of 4D cylindrical echocardiographic images
Suri Two-dimensional fast magnetic resonance brain segmentation
US20030095121A1 (en) Vessel detection by mean shift based ray propagation
US6718054B1 (en) MRA segmentation using active contour models
US7672492B2 (en) Method of incorporating prior knowledge in level set segmentation of 3D complex structures
Tai et al. Wavelet frame based multiphase image segmentation
JP2006302291A (en) Method and system for segmentation of object in image data set
Deschamps et al. Minimal paths in 3D images and application to virtual endoscopy
Zhang et al. A local information based variational model for selective image segmentation
EP1208535A1 (en) Mra segmentation using active contour models
Yang et al. Neighbor-constrained segmentation with 3d deformable models
Holtzman-Gazit et al. Hierarchical segmentation of thin structures in volumetric medical images
Park et al. Medical image segmentation using level set method with a new hybrid speed function based on boundary and region segmentation
Yang et al. 3D image segmentation of deformable objects with shape-appearance joint prior models
Luo Automated medical image segmentation using a new deformable surface model
Rudra et al. 3D Graph cut with new edge weights for cerebral white matter segmentation
Slabaugh et al. A variational approach to the evolution of radial basis functions for image segmentation
Wang et al. Level set segmentation based on local Gaussian distribution fitting
Suri et al. Modeling segmentation via geometric deformable regularizers, pde and level sets in still and motion imagery: a revisit

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A1

Designated state(s): JP

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE

121 Ep: the epo has been informed by wipo that ep was designated in this application
DFPE Request for preliminary examination filed prior to expiration of 19th month from priority date (pct application filed before 20040101)
WWE Wipo information: entry into national phase

Ref document number: 2000943090

Country of ref document: EP

WWP Wipo information: published in national office

Ref document number: 2000943090

Country of ref document: EP

AK Designated states

Kind code of ref document: C2

Designated state(s): JP

AL Designated countries for regional patents

Kind code of ref document: C2

Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE

COP Corrected version of pamphlet

Free format text: PAGES 1/11-11/11, DRAWINGS, REPLACED BY NEW PAGES 1/9-9/9; DUE TO LATE TRANSMITTAL BY THE RECEIVINGOFFICE

NENP Non-entry into the national phase

Ref country code: JP

WWW Wipo information: withdrawn in national office

Ref document number: 2000943090

Country of ref document: EP