The present invention relates to a method and system for segmenting a space spanned by a set of medical data, such as a 2- or 3-dimensional MRI image of at least part of a brain or other structure consisting of a set of intensity values at corresponding locations in the space. The segmentation produces a model of an object having medical significance, that is an estimate of the portion of the space occupied by the object.
Since their introduction by Kass et al. (Kass 1988), a variety of deformable models have been proposed and used in a numerous areas such as image segmentation, shape modeling, and object tracking. Existing deformable models can be classified as either parametric deformable models (PDMs) (Kass, 1988; Cohen, 1991; Xu, 1998; McInery, 1996; Cottes, 1991, 2001; Shen, 2001; Liu, 2005) or geometric deformable models (GDMs) (Malladi, 1995; Caselles, 1997; Yezzi, 1997; Sethian, 1999; Han, 2003) according to their representation and implementation. In particular, PDMs are represented explicitly as parameterized curves or surfaces in a Lagrangian formulation. GDMS, on the other hand, are represented implicitly as level sets of higher-dimensional level set functions. That is, to model a contour within a 2-D space, or a surface within a 3-D space, we define a continuous function (a level set function) which takes a certain value (e.g. zero) on the contour or surface (referred to as the “level set”), and a non-zero value elsewhere. The level set function evolves according to an Eulerian formulation (Sethian 1999).
GDMs are superior to PDMs in the following aspects:
Model representation and construction. GDMs are completely intrinsic and independent of the parameterization of the evolving contour. Thus, there is no need to add or remove nodes from the initial parameterization or adjust the spacing of the nodes. For PDMs, particularly, in the 3D cases, it is not a trivial task to generate parameterized representations for surface models.
Model property calculation. The intrinsic geometric properties of a geometric deformable model, such as its normal and curvature, can be easily determined from the level set function.
Shape adaptation. In GDMs, the propagating contour or surface can automatically change its topology (e.g., merge or split) without requiring an elaborate mechanism to handle such changes.
Self-intersection detection. When the topology preserving transformations are required for a deformable model, self-intersections can become a problem. In GDMs, digital topology methods (Bertrand 1994) are available to detect self-intersections, while in PDMs, the computational demands related to self-intersection detection are usually very high (MacDonald 2000).
On the contrary, PDMs have the following merits over GDMs:
Model constraints. Since PDMs are represented globally, it is convenient to add high level constraints on shape and topology. These constraints are particularly important when local features of the model and image are liable to deform the model in unexpected ways.
Leakage prevention. Leakage from the boundary is less common in PDMs, whereas it is sensitive to boundary gaps and weak edges in GDMs.
Convergence detection. The deformation procedures in PDMs have an explicit stop condition which is the minimization of some kinds of energy function, whereas in GDMs, the procedure usually stops when the certain number of iterations reaches a predetermined maximum, but no simple and fast methods are available to automatically estimate a suitable maximum number of iterations.
The present invention aims to provide new and useful ways to estimate, within a space spanned by a set of medical data (e.g. MRI data consisting of intensity values at many corresponding locations in a two- or three-dimensional space), the portion of the space occupied by an object.
In general terms, the present invention proposes a boundary-based geometrical deformable model (BGDM). That is, a model is represented by a level set function, which is modified by force calculated using the medical data and the boundary of the model.
Thus, in a way similar to traditional geometrical deformable models, the invention employs a model which is embedded in a level set of a higher-dimensional function, while in a way similar to the traditional parametric deformable models, the deformation of the model is calculated by force on its boundary.
The convergence condition used for termination of the iteration may be the minimization of an energy function related to the force. Thus, the embodiments of the invention may have the convergence advantages of PDMs.
The force used to deform the level set function may be calculated at each point of a narrowband region (neighbourhood) of the level set, and the level set function modified only at these locations. This reduces the computational burden, compared to modifying the entire level set function in all iterations. Upon detecting that the level set approaching the boundary of the narrowband region, portions of level set function outside the narrowband region may be refreshed.
The narrowband may include two sets of locations (pixels in a 2D space or voxels in a 3D space): first locations at the boundary of the model (level set) and second locations in a neighbourhood of the boundary (“narrowband region”). The force at the first locations is calculated using the medical data, and the force at the second locations is calculated using the calculated force at the first locations. Specifically, around each of the second locations a neighbourhood may be defined, and the force at the second locations can be found by reference to the force at those first locations which are within the neighbourhood. It has been found that this use of neighbourhoods reduces (or even avoids) the leakage problem (in which the object may wrongly grow to other regions of the space).
High level knowledge about the object can be incorporated in any one or more of several ways. First, it may be used in the definition of the force which modifies the level set function (e.g. such that when it appears that a predetermined fact about the object is to be contravened by the modification, the force switches from being one based on the medical data to one which incorporates the prior knowledge). Secondly, it may be used to choose the exact form of the force which is calculated using the medical data. Thirdly, the modification of the level set function may only be performed when this does not violate some prior knowledge of the object.
The force which is calculated using the medical data may be a “gradient flux”-based force, which is a function of the gradient of intensity values of the medical data. The gradient flux-based force at each point on the level set may be obtained efficiently by searching along a line segment passing through that point and normal to the level set at that point, for a location which maximises the value of a gradient flux density which is a function of the gradient of the intensity values (or more generally, a location which is a stationary point of the gradient flux density). In a manner inspired by electromagnetic systems, the gradient flux density may vary as a linear function of the gradient of the intensity values divided by a linear function of the square of the distance of the location from the level set (i.e. analogous to an inverse square law). The energy function may be average of the gradient flux density over the level set.
One specific expression of the invention is a method of image segmentation and modelling, within a space spanned by a set of medical data, the portion of the space occupied by an object, the method comprising:
(a) defining a level set function, a level set of the level set function being the boundary of a deformable model of the object;
(b) iteratively performing at least once the set of steps:
(c) modelling the boundary of the said object as the level set of the modified level set function.
However, the invention may alternatively be expressed as a computer system arranged to carry out this method, or as a computer program product (i.e. software) which when run by a computer causes the computer to carry out the method.
An embodiment of the invention will now be described for the sake of example only with reference to the following figures, in which:
An embodiment of the invention will now be described with reference to
The embodiment employs an initial model which is defined in the space of the image. If the space of the image is 2D, the initial model is a contour represented polygon, or, if the space of the image is 3D, the model is a surface represented by triangle meshes. The contour or surface forms one or more closed inside regions which indicate object itself. The model is chosen manually based on the prior-knowledge of the object's shape and size, and initially located according the location of the object in the image. This model is then allowed to evolve, according to equations defined below, within the space of the image to form a final model which is the estimate of the position of the object.
Positions in the space are denoted by a vector x. The initial model is expressed as a zero level set of a higher dimensional level set function φ(x,t), defined as follows. The initial value of φ(x,t) at time t=0, that is φ(x,0), is defined as the signed Euclidian distance map (Maurer 2003) of the model: with negative value inside the model, and positive outide the model. The evolution of the function φ(x,t) is governed by the following equation:
φt(x,t)=−n(x,t)·f(x,t) (1)
where, φt(x,t) denotes the time derivative of φ(x, t), f(x,t) denotes the total force which is a function of the model at time t and of the image dataset (explicit expressions for f(x, t) are given below), and n(x, t) is the unit normal of the level set, i.e. the unit vector which is defined by:
n(x,t)=∇φ(x,t)/|∇φ(x,t)|. (2)
The numerical solution of Equation (1) is obtained by approximating the real time variable t by a series of discrete time ti=i·Δt, where i=0, 1, 2, . . . and the time step is denoted as Δt. The functions φt(x,t), n(x,t) and f(x,t) can be approximated as φ(x,ti), n(x,ti) and f(x,ti) respectively. The time derivative φt(x,t) is approximated as the difference between the respective values of φ(x,ti) at two successive time steps, divided by Δt. The spatial derivative on the right hand side of (2) is treated by upwind numerical schemes (Sethian, 1999), which gives:
φ(x,tt+1)=φ(x,ti)−Δt·n(x,ti)·f(x,ti)(i=0,1,2, . . . ). (3)
To describe the evaluation state of the deformable model at a given time t, we introduce several notations:
That is, the set of all locations x such that φ(x, t) is at least zero.
We now turn to the problem of defining f(x,ti). To do so, we first define a novel function called the “gradient flux”.
To begin with, according to the object we are trying to model, we define a “gradient characteristic number” g, which describes the gradient direction at a boundary of the object. In general terms, we select g according to whether the desired object is a valley of the intensity (g=1), or a peak of the intensity (g=−1). If the desired object is neither of these things, such as a saddle-point of the intensity, we select g=0. More specifically, if the object is such that all gradients on the object boundary point to the outside, g is set to 1. If all gradients on the objection boundary point to the inside, g is set to −1. Otherwise, g is set to 0. For example, in a T1 MR orbital image, g is 1 for orbital muscles; in T1 MR human brain images, g is 1 for the hippocampus, −1 for ventricles, and 0 for caudate nuclei.
Now we define the “gradient flux” as follows. This is intended to produce the result that the nearer the boundary of the model is to that of the object in the image, the larger the average gradient flux on the boundary of the model. Referring to
Consider a point 5 on the boundary 1. The point 5 has coordinate x, and the normal direction to the boundary at this point is shown on
To give a numerical value to the gradient flux, suppose that a point charge were located at point x+r·n(x,t). Thus, the electric flux density at point x would decrease with the distance r according to the function 3/(4πr2). By analogy, we will define the gradient flux density at any point along the line segment such that the gradient flux decreases with the distance to the boundary in the same way as the electric flux. In order to avoid the gradient flux becoming infinite for r=0, we define the gradient flux using the function 3/(4πr2+1) instead of 3/(4πr2), so that gradient flux is defined at r=0.
According to these heuristic rules, the embodiment calculates the gradient flux density e(x, n(x, t)) at any model boundary point 5 having location x as:
where n(x, t) denotes the outside normal of the model boundary 1 at point x; ∇I(x+r·n(x,t)) is the gradient of the image at x+r·n(x,t); d is the predefined maximum search length; the function ξ indicates the consistency requirement concerning the gradient direction, and is defined as:
Having now defined the gradient flux, we can now define the force f(x,t) used in Eqn. (1). Unlike traditional level set techniques, in which the force f(x,t) is directly calculated at the point x, we first calculate the force on the model boundary, and then expand the force to a narrowband of the boundary.
Specifically, for a point x on the boundary ∂(t), we define its neighborhood Nr(x,t) as the set of points that are on the boundary ∂(t) and which have distance to x less than a given threshold r. As in deformable surface models (Liu 2005, Shen 2001), a boundary point x of the model is considered as the vertex of its boundary patch Nr(x,t), on which the force at point x is calculated. The calculation combines the model constraint(s) with image features. This is done by defining the total force f(x,t) using (i) a gradient flux-based force fa(x,t) obtained using the gradient flux as defined above, which in turn is based on the intensity values, and (ii) an intensity based force fb(x,t) which ensures that model obeys the model constraints.
First we obtain the gradient flux-based force. As a first step, the total gradient flux q0(x,t) on the boundary patch Nr(x,t) is calculated as:
Then, the gradient fluxes q+(x,t) and q−(x,t) are calculated. These correspond respectively to the patch Nr(x,t) being shifted to the outside of the model region Ω(r) with unit translation n(x,t), and inside the model region Ω(t) with unit translation −n(x,t)).
where, x′+n(x′,t) and x′−n(x′,t) are points obtained by shifting the point x′ with translation n(x′,t) and −n(x′,t) respectively.
When the q+(x,t) is higher than q−(x,t) and q0(x,t), the model boundary ∂(t) is attracted to expand at point x; when q−(x,t) is higher than q+(x,t) and q0(x,t), the boundary ∂(t) is urged inward at point x; in other cases, the boundary should remain stationary. Based on these heuristic rules, the gradient flux-based force fa(x,t) is calculated as:
where, n− is the unit vector n(x,t) with an additional flag to indicate that the force is used to keep the boundary ∂(r) fixed at point x.
Now we calculate the intensity-based force. Let μ(t) and σ(t) denote the mean and variation of the image intensity within the region, and let and l(x,t) denote the mean of the image intensity on the boundary patch Nr(x,t), i.e.:
In the case that the gradient characteristic number g of an object is 1, i.e. the intensity inside the object is less than that outside the object, if the intensity l(x,t) at the boundary point x is greater than μ(t)+σ(t), it is expected that it will be necessary to remove point x from region Ω(t). Thus, we generate an intensity-based force at x which pushes the boundary ∂(t) toward the inside of the region Ω(t). Conversely, otherwise, if l(x,t)<μ(t)+σ(t), then we generate an intensity-based force to push x outwards. The cases for g=−1 and 0 are similar to g=1. Specifically, motivated by these heuristic rules, the embodiment generates an intensity-based force fb(x,t) at a boundary point x which is calculated as:
where the cut function φ is defined by:
φ(u)=min(1,max(−1,u)). (16)
and accordingly limits the force fb(x,t) within the interval [−1, 1].
Let fa(t) and fmax(t) denote the average and maximum magnitude of forces fa(x,t) on the boundary ∂(t), i.e.:
The total force f(x,t) is finally estimated from the gradient flux-based force fa(x,t) and the intensity-based force fb(x,t) with the force fa(x,t) as the priority:
We have now fully defined all elements of Eqn. (1) for a point on the boundary. The algorithm implements the level set evolution by a fast narrowband approach (Sethian 1999). We expand the boundary force f(x, t) into a narrowband br(t) of the boundary ∂(t).
In the 3D case, for any voxel x we use N26(x) to denote its 26-neighborhood (i.e. its 26 nearest neighbours, which are directly adjacent). At any given time t, we define a set of neighbourhoods {Ω(k)}, for k=0, 1, 2, . . . , with corresponding boundaries {∂(k)}. Ω(0) is defined to be equal to Ω(t), and its boundary ∂(0) is initially set to ∂(t). For k=0, 1, 2, etc. we define the region Ω(k+1) and its boundary ∂(k+1) iteratively from Ω(k) and ∂(k):
∂(k+1)=[∪xε∂
The force f(x,t) on the set of points ∂(0) is first calculated according to Eqn. (19). The embodiment then uses the force on each set of points ∂(k) to obtain the force on the set of points ∂(k+1) as follows. For any voxel which is ∂(k+1) and has position x, the force f(x,t) is a weighted average of the force on those of its 26 nearest-neigbors which are in the set ∂(k), i.e.:
where, d(x,x′) denotes the Euclidean distance of two points x and x′.
The force expansion in the 2D case is similar to that in the 3D case. It differs just by changing the 26-neighborhood in the 3D case into a 8-neighborhood in the 2D case.
The embodiment makes use of an “energy function” to determine when its iterations should terminate. The energy ē(t) of the model is estimated as the negative of the average gradient flux density on the model boundary ∂(t):
We are now in a position to define the complete algorithm. This will be explained with reference to
The embodiment was used to segment orbital muscles from MR images.
This application is a U.S. National Phase Application of PCT International Application No. PCT/SG2007/000275, filed on Aug. 24, 2007, which claims the benefit of U.S. Application Ser. No. 60/839,684, filed on Aug. 24, 2006, both of which are expressly incorporated by reference herein.
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/SG2007/000275 | 8/24/2007 | WO | 00 | 2/19/2009 |
Publishing Document | Publishing Date | Country | Kind |
---|---|---|---|
WO2008/024083 | 2/28/2008 | WO | A |
Number | Name | Date | Kind |
---|---|---|---|
6766043 | Zeng et al. | Jul 2004 | B2 |
6813373 | Suri et al. | Nov 2004 | B1 |
6845260 | Liu et al. | Jan 2005 | B2 |
6892089 | Prince et al. | May 2005 | B1 |
7359538 | Zeng et al. | Apr 2008 | B2 |
7391893 | Liang et al. | Jun 2008 | B2 |
7901357 | Boctor et al. | Mar 2011 | B2 |
8064664 | Suri et al. | Nov 2011 | B2 |
Number | Date | Country | |
---|---|---|---|
20100166275 A1 | Jul 2010 | US |
Number | Date | Country | |
---|---|---|---|
60839684 | Aug 2006 | US |