Abstract:
Segmenting a digitized image includes providing a diffusion tensor field image, partitioning the image into 2 regions, each point being assigned to one of the regions, associating a level-set function with each region, the level-set function having a negative value in one region and a positive value in the other region, calculating a mean value of the diffusion tensor field over each of 2 regions, initializing an energy defined as a functional of level-set functions and diffusion tensor field, changing the region membership of each point in the image if the energy functional value decreases as a result of region membership change and updating said mean value of said diffusion tensor field over each of 2 regions, and obtaining a segmentation of the image when the magnitude of the change of the energy function value resulting from changing the region membership of a point is less then a predetermined threshold.

Description:
CROSS REFERENCE TO RELATED UNITED STATES APPLICATIONS 
     This application claims priority from “Fast Tensor Field Segmentation Algorithm”, U.S. Provisional Application No. 60/644,144 of McGraw and Wang, filed Jan. 14, 2005, the contents of which are incorporated herein by reference. 
    
    
     TECHNICAL FIELD 
     This invention is directed to the segmentation of objects in digital medical images. 
     DISCUSSION OF THE RELATED ART 
     Variational and PDE methods have been widely used in image processing in the past few years. Image segmentation is an issue in early vision and has been extensively studied through these techniques. The basic idea is to represent a 2D image as a R 2  function (or R 3  for a 3D image) that satisfies a time dependent PDE that characterizes the problem. The solution of the differential equation gives the processed image at the appropriate scale. The PDE can be derived by minimizing an energy functional, which can be formulated as:
 
Arg{min u F(u)},
 
where F is a given energy computed over the image u. Let F′(u) denote the Euler-Lagrange derivative. The necessary condition for u to be a minimizer of F is that F′(u)=0, where the (local) minima may be computed via the steady state solution of the equation
 
 u   t   =−F ′( u ),
 
where t is an artificial time-marching parameter. The most conventional example is the Dirichlet integral:
 
                   min   u     ⁢     F   ⁡     (   u   )         =         ∫   Ω     ⁢              ∇   u          2     ⁢           ⁢     ⅆ   x     ⁢     ⅆ   y         +       λ   2     ⁢     ∫         (     u   -     u   0       )     2     ⁢     ⅆ   x     ⁢     ⅆ   y               ,         
which is associated with equation
   u   t   =Δu −λ( u−u   0 ). 
     The use of variational and PDE methods in image analysis leads to modeling images in a continuous domain. This simplifies the formalism, which becomes grid independent and isotropic. Conversely, when the image is represented as a continuous signal, PDEs can be seen as the iteration of local filters with an infinitesimal neighborhood. This interpretation of PDE&#39;s allows one to unify and classify a number of the known iterated filters, as well as to derive new ones. Many of the PDEs used in image processing and computer vision are based on moving curves and surfaces with curvature-based velocities. In this area, level set methods have proven to be useful. The basic idea is to represent the curves or surfaces as the zero level set of a higher dimensional hypersurface. This technique not only provides accurate numerical implementations but can also handle topological change easily. 
     Many fields in science and engineering use tensors to describe physical quantities e.g., in solid mechanics, stress and strain are commonly characterized by tensors, in fluid mechanics, anisotropic diffusion is characterized by a tensor, in image processing, tensors have been used for describing local structure in images. More recently, in medical imaging, water diffusion in tissues has been depicted by a diffusion tensor characterizing the anisotropy within the tissue. In this context, processing the data may involve restoration from noisy data, segmentation, visualization and possibly even tensor valued image registration. One factor in tensor field analysis is a proper choice of a tensor distance that measures the similarity or dissimilarity between tensors and is particularly important in the aforementioned tasks. 
     In general, any kind of matrix norm can be used to induce a tensor distance. One such example is the tensor Euclidean distance obtained by using the Frobenius norm. Due to its simplicity, tensor Euclidean distance has been used extensively in tensor field restoration. As compared to other similarity measures for matching diffusion tensor images, the Euclidean difference measure yields the best results. Though not many sophisticated tensor distances have been proposed in tensor field analysis, there are quite a few in the context of machine learning. One interesting tensor distance measure uses information geometry in the space of positive definite matrices to derive a Kullback-Leibler divergence for two matrices and then uses it in to approximate an incomplete kernel. 
     Diffusion tensor magnetic resonance imaging is a relatively new imaging modality that is able to quantify the anisotropic diffusion of water molecules in highly structured biological tissues. Diffusion refers to the movement of molecules as a result of random thermal agitation. In the context of diffusion tensor magnetic resonance imaging (DT-MRI), it refers specifically to the random translational motion of water molecules in the part of the anatomy being imaged with MR. In three dimensions, water diffusivity can be described by a 3×3 symmetric positive definite matrix D, known as a diffusion tensor, which is related to the geometry and organization of the microscopic environment. The probability density function of the molecular motion about r∈R 3  can be modeled by a Gaussian function whose covariance matrix is given by the diffusion tensor D. Diffusion Tensor Imaging (DTI) then produces a volumetric image containing, at each voxel, a 3×3 symmetric positive definite tensor. The estimation of these tensors requires the acquisition of diffusion weighted images in different sampling directions. 
     Diffusion tensor MRI is particularly relevant to a wide range of clinical pathologies investigations, such as acute brain ischemia detection, stroke, Alzheimer&#39;s disease, schizophrenia, etc. It is also useful in order to identify the neural connectivity of the human brain. As of today, diffusion MRI is the only non-invasive method that allows one to distinguish the various anatomical structures of the cerebral white matter such as the corpus callosum, the arcuate fasciculus or the corona radiata. These are examples of commisural, associative and projective neural pathways, the three major types of fiber bundles, respectively connecting the two hemispheres, regions of a given hemisphere or the cerebral cortex with subcortical areas. In the past, many techniques have been proposed to classify gray matter, white matter and cephalo-spinal fluid from T1-weighted MR images, but few techniques address segmentation of the internal structures of white matter. 
     SUMMARY OF THE INVENTION 
     Exemplary embodiments of the invention as described herein generally include methods and systems for segmenting a tensor field such as that obtained from diffusion tensor MRI. According to an embodiment of the invention, each voxel in the image is assigned to a class based on the statistics of each region. A segmentation model according to an embodiment of the invention uses a symmetrized Kullback-Liebler divergence to define a tensor distance along with a variance model for each region. A segmentation model according to an embodiment of the invention also includes an additional energy term for each region to penalize the distance from the mean fractional anisotropy in each region where the segmented fractional anisotropy fields are modeled by Gaussian distributions. According to an embodiment of the invention, a region based active contour model involves the definition of a new tensor discriminant based on information theory, a new technique for the computation of the mean tensor of a tensor field in closed form that facilitates the efficient segmentation of the tensor field, and an extension of the region-based active contour model to handle matrix-valued images. According to an embodiment of the invention, a technique for the segmentation of an arbitrary probability density function (pdf) fields examines the statistics of the distribution of the Kullback-Leibler distances between these pdfs to perform direct segmentation of internal structures of white matter. According to an embodiment of the invention, a fast levelset technique is adapted to solve the region based active contour model. 
     According to an aspect of the invention, there is provided a method for segmenting a digitized medical image, including providing a digitized image of an anatomical region, said image comprising a diffusion tensor field corresponding to a domain of points on an 3-dimensional grid, partitioning said image into 2 regions, wherein each point is initially assigned to one of said regions, associating a level-set function with each region, wherein said level-set function has a negative value in one region and a positive value in the other region, calculating a mean value of said diffusion tensor field over each of said 2 regions, initializing an energy defined as a functional of said level-set functions and said diffusion tensor field, changing the region membership of each point in said image if the energy functional value decreases as a result of said region membership change and updating said mean value of said diffusion tensor field over each of said 2 regions, and obtaining a segmentation of said image when the magnitude of the change of said energy function value resulting from changing the region membership of a point is less then a predetermined threshold. 
     According to a further aspect of the invention, the energy functional includes the terms
 
(1−β)∫ Ω (log p d,1 (d 2 (T(x),T 1 ))H ε (φ)+log p d,2 (d 2 (T(x),T 2 ))(1−H ε (φ)))dx
 
wherein φ:Ω→R 3  denotes the level set function defined on the image domain Ω whose zero isosurface coincides with a curve enclosing one of said 2 regions, T 1  and T 2  are the mean values of the tensor field in each of the two regions, H ε (φ) is the regularized version of the Heaviside function of the level set function, p d,1  and p d,2  are probability distribution functions of the tensor distance in each of the two regions, β is a weighting parameter, and d 2 (.,.) is a tensor distance.
 
     According to a further aspect of the invention, the tensor distance between tensors T 1  and T 2  is defined by 
               d   ⁡     (       T   1     ,     T   2       )       =       1   2     ⁢         tr   ⁡     (         T   1     -   1       ⁢     T   2       +       T   2     -   1       ⁢     T   1         )       -     2   ⁢   n                 
wherein tr(•) is the matrix trace operator, and n is the size of the square matrices T 1  and T 2 .
 
     According to a further aspect of the invention, the mean value  M  of a diffusion tensor field T(x) in a region R is defined by
 
  M =√{square root over (B −1 )}[√{square root over (√{square root over (B)}A√{square root over (B)})}]√{square root over (B −1 )},
 
wherein A=∫ R T(x)dx, and B=∫ R T −1 (x)dx, and R is the region of integration.
 
     According to a further aspect of the invention, the level-set function has a value of −1 in one region and a value of +1 in the other region. 
     According to a further aspect of the invention, the probability distribution of the distance from said diffusion tensor field to a mean value of said field is defined as 
               p     d   ,   i       =       1       2   ⁢   π   ⁢           ⁢     σ     i   ,   d     2           ⁢     exp   ⁡     (       -       d   2     ⁡     (       T   ⁡     (   x   )       ,     T   i       )           2   ⁢     σ     i   ,   d     2         )               
where T i  is the mean value of said tensor field in one of said 2 regions, T(x) is said tensor field, σ i,d  is the variance of said distribution in one of said 2 regions.
 
     According to a further aspect of the invention, the method comprises calculating a fractional anisotropy defined by 
                 A   ⁡     (     T   ⁡     (   x   )       )       =             (       λ   1     -     λ   2       )     2     +       (       λ   2     -     λ   3       )     2     +       (       λ   1     -     λ   3       )     2             2     ⁢         λ   1   2     +     λ   2   2     +     λ   3   2               ,         
wherein the λ 1 , λ 2 , λ 3  are eigenvalues of the diffusion tensor field, and calculating a mean and variance of said fractional anisotropy over each of said 2 regions.
 
     According to a further aspect of the invention, the method comprises updating the mean value and variance of said fractional anisotropy of the diffusion tensor in each of said 2 regions if the energy functional value decreases as a result of said region membership change. 
     According to a further aspect of the invention, the energy function further includes terms
 
−β∫ Ω (log p α,1 (A(T(x)))H ε (φ)+log p α,2 (A(T(x)))(1−H ε (φ)))dx,
 
wherein φ:Ω→R 3  denotes the level set function defined on the image domain Ω whose zero isosurface coincides with a curve enclosing one of said 2 regions, H ε (φ) is the regularized version of the Heaviside function of the level set function, p α,1  and p α,2  are probability distribution functions of the fractional anisotropy in each of the two regions, and β is a weighting parameter.
 
     According to a further aspect of the invention, the probability distribution function of said fractional anisotropy is defined as 
                 p     a   ,   i       =       1       2   ⁢   π   ⁢           ⁢     σ     i   ,   A     2           ⁢     exp   ⁡     (                A   ⁡     (   T   )       -     A   _            2       2   ⁢     σ     i   ,   A     2         )           ,         
wherein Ā is the mean value of A(T) over the region i, and σ i,A  is the variance of p α,i  over the region i.
 
     According to a further aspect of the invention, the method comprises calculating variances of the distance of said diffusion tensor to the mean value over each of said two regions, and updating said variances when the energy functional value decreases as a result of said region membership change. 
     According to a further aspect of the invention, a change in energy ΔE ij  from changing the membership of a point from region i to region j includes 
               Δ   ⁢           ⁢     E     ij   .   d         =       (     1   -   β     )     ⁢       (             d   2     ⁡     (       T   ⁡     (   x   )       ,     T   j       )         2   ⁢     σ     j   ,   d     2         ⁢       n   j         n   j     +   1         -           d   2     ⁡     (       T   ⁡     (   x   )       ,     T   i       )         2   ⁢     σ     i   ,   d     2         ⁢       n   i         n   i     -   1           )     .             
wherein the value of the point being changed is T(x), T i  and T j  are the mean values of the tensor field in each region, respectively, σ i,d  and σ j,d  are the variances of p d,i , p d,j  over the regions i and j, respectively, and n i  and n j  are the respective number of points for each region.
 
     According to a further aspect of the invention, a change in energy ΔE ij  from changing the membership of a point from region i to region j includes 
               Δ   ⁢           ⁢     E     ij   ,   A         =     β   ⁡     (                    A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     j            2       2   ⁢     σ     j   ,   A     2         ⁢       n   j         n   j     +   1         -                  A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     i            2       2   ⁢     σ     i   ,   A     2         ⁢       n   i         n   i     -   1           )             
wherein Ā i  and Ā j  are the mean values of A(T) over the regions i and j, respectively, σ i,A  and σ j,A  are the variances of p α,i , p α,j , over the regions i and j, respectively, and n i  and n j  are the respective number of points for each region.
 
     According to another aspect of the invention, there is provided a program storage device readable by a computer, tangibly embodying a program of instructions executable by the computer to perform the method steps for segmenting a digitized medical image, 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIGS. 1 and 2  depict an exemplary segmentation of the Corpus Callosum, with and without an FA slice, respectively, according to an embodiment of the invention. 
         FIGS. 3 and 4  depict an exemplary segmentation of the Corona Radiata, with and without an FA slice, respectively, according to an embodiment of the invention. 
         FIGS. 5 and 6  depict an exemplary segmentation of the ventricle, with and without an FA slice, respectively, according to an embodiment of the invention. 
         FIG. 7  is a flow chart of an exemplary method for tensor field segmentation, according to an embodiment of the invention. 
         FIG. 8  is a block diagram of an exemplary computer system for implementing a tensor field segmentation algorithm according to an embodiment of the invention. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     Exemplary embodiments of the invention as described herein generally include systems and methods for segmenting a tensor field such as that obtained from diffusion tensor MRI. A tensor field segmentation method according to an embodiment of the invention incorporates a new discriminant for tensors into an active contour model without edges. The discriminant employed according to an embodiment of the invention is based on information theory, which follows from the physical phenomena of diffusion, is affine invariant and allows for the computation of the mean of the tensor field in closed form. A model according to an embodiment of the invention is implemented in a level set framework to take advantage of the easy ability of this framework to change topologies when desired. A technique for the segmentation of probability density fields according to an embodiment of the invention can extract anatomical structures in anisotropic biological tissues such as the brain white matter. 
     In the context of DT-MRI, diffusion of water molecules may be characterized by a 2-tensor T which is positive definite. This T is related to the displacement r of water molecules at each lattice point in the image at time t via 
               p   ⁡     (       r   ❘   t     ,   T     )       =       1           (     2   ⁢   π     )     n     ⁢          2   ⁢           ⁢   tT                ⁢     exp   ⁡     (         -     r   T       ⁢     T     -   1       ⁢   r       4   ⁢   t       )               
where n is the dimension of the tensor. Thus it is natural to use the distance measure between Gaussian distributions to induce a distance between these tensors. The most frequently used information theoretic distance measure is the Kullback-Leibler (KL) divergence defined as:
 
             KL   (       p   ⁢        q   )       =     ∫       p   ⁡     (   x   )       ⁢   log   ⁢       p   ⁡     (   x   )         q   ⁡     (   x   )         ⁢     ⅆ   x                 
for two given densities p(x) and q(x). The KL divergence is not symmetric and a popular way to symmetrize it is given by:
 
               J   ⁡     (     p   ,   q     )       =       1   2     ⁢     (     KL   (         p   ⁢        q   )       +     KL   ⁡     (     q   ⁢        p   )       )         ,                 
which is called the J-divergence. According to an embodiment of the invention, a definition of tensor distance for symmetric positive definite (SPD) tensors is the square root of the J-divergence:
   d ( T   1   ,T   2 )=√{square root over ( J ( p ( r|t,T   1 ), p ( r|t,T   2 )))}{square root over ( J ( p ( r|t,T   1 ), p ( r|t,T   2 )))}. 
     It is known that twice the KL divergence and thus twice the J-divergence is the square distance of two infinitesimally nearby points on a Riemannian manifold of parameterized distributions. Thus, taking the square root of the above definition is justified. Furthermore, this definition has a very simple form given by: 
               d   ⁡     (       T   1     ,     T   2       )       =       1   2     ⁢         tr   ⁡     (         T   1     -   1       ⁢     T   2       +       T   2     -   1       ⁢     T   1         )       -     2   ⁢   n                 
where tr(•) is the matrix trace operator, and n is the size of the square matrices T 1  and T 2 .
 
     When a coordinate system undergoes an affine transformation, the tensor field will also be transformed. If the coordinate system undergoes an affine transform y=Ax+b, then the displacement of the water molecules will be transformed as {circumflex over (r)}=Ar. Since r has a Gaussian distribution with covariance matrix 2tT, the transformed displacement {circumflex over (r)} has a covariance matrix of 2tATA T . Thus, the transformed tensor field will be:
 
 {circumflex over (T)} ( y )= AT ( x ) A   T   , y=Ax+b.  
 
     The above definition of tensor distance is invariant to such transformations:
 
 d ( T   1   ,T   2 )= d ( AT   1   A   T   ,AT   2   A   T ).
 
     The mean value M(T,R) of a tensor field T over a region R is defined as: 
                   M   _     ⁡     (     T   ,   R     )       =       min     M   ∈     SPD   ⁡     (   n   )           ⁢       ∫   R     ⁢         d   2     ⁡     [     M   ,     T   ⁡     (   x   )         ]       ⁢           ⁢     ⅆ   x             ,         
where SPD(n) denotes the set of symmetric positive definite matrices of size n. It can be shown that this mean value can be computed according to the formula
   M =√{square root over (B −1 )}[√{square root over (√{square root over (B)}A√{square root over (B)})}]√{square root over (B −1 )}, 
where A=∫ R T(x)dx, and B=∫ R T −1 (x)dx. Since A and B are both SPD matrices, M is also an SPD matrix. The mean value of a tensor field over a region is used in the region-based contour model used in a segmentation according to an embodiment of the invention.
 
     According to an embodiment of the invention, a model for piecewise constant tensor field segmentation in R 2  is obtained by minimizing the following energy integral:
 
 E ( C,T   1   ,T   2 )=∫ Ω   d   2 ( T ( x ), T   1 ) dx+∫   Ω     c     d   2 ( T ( x ), T   2 ) dx+α|C|.   (1)
 
     Here, the curve C is the boundary of the desired unknown segmentation, Ω is the region enclosed by C and Ω C  is the region outside C, T 1  and T 2  are the mean values of the tensor fields in the regions Ω and Ω C  respectively, |C| is the arclength of the curve C, α is a regularization parameter, and d(.,.) is the tensor distance as defined above. 
     An active contour model according to an embodiment of the invention above can segment tensor fields with two piecewise constant regions, where each region type can have disconnected parts, and incorporates the above-defined tensor distance. The Euler-Lagrange equation for the above variational principle is given by:
 
(α k+d   2 ( T,T   1 )− d   2 ( T,T   2 )) N= 0,
 
where T 1 =M(T,Ω), T 2 =M(T,Ω C ), k is the curvature of the curve C at location x, and N is the outward normal to the curve C. In a two phase implementation, the curve evolution of the Euler-Lagrange is governed by:
 
                   ∂   C       ∂   t       =       -     (       α   ⁢           ⁢   k     +       d   2     ⁡     (     T   ,       T   1     ⁡     (   t   )         )       -       d   2     ⁡     (     T   ,       T   2     ⁡     (   t   )         )         )       ⁢   N       ,         
which can be easily implemented in a level set framework. The corresponding level set formulation is given by:
 
     
       
         
           
             
               
                 ∂ 
                 ϕ 
               
               
                 ∂ 
                 t 
               
             
             = 
             
               
                 [ 
                 
                   
                     α 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ∇ 
                       
                         · 
                         
                           
                             ∇ 
                             
                                 
                             
                             ⁢ 
                             ϕ 
                           
                           
                              
                             
                               ∇ 
                               ϕ 
                             
                              
                           
                         
                       
                     
                   
                   + 
                   
                     
                       d 
                       2 
                     
                     ⁡ 
                     
                       ( 
                       
                         T 
                         , 
                         
                           T 
                           1 
                         
                       
                       ) 
                     
                   
                   - 
                   
                     
                       d 
                       2 
                     
                     ⁡ 
                     
                       ( 
                       
                         T 
                         , 
                         
                           T 
                           2 
                         
                       
                       ) 
                     
                   
                 
                 ] 
               
               ⁢ 
               
                 
                    
                   
                     ∇ 
                     ϕ 
                   
                    
                 
                 . 
               
             
           
         
       
     
     Assuming a partition of the data Ω, one seeks an optimal separating surface Γ between a segment Ω 1  and the rest of the volume Ω 2 . The most representative mean of tensor fields that minimizes the tensor distance can be denoted by T 1  and T 2 . Furthermore, it is possible to model the distribution of the tensor distances to T 1  and T 2  in their respective domains by suitable densities p d,1 , p d,2 . It can be assumed that p d,1  and p d,2  are Gaussians of zero mean and variances σ 1,d   2 ,σ 2,d   2 . The mean distance to T 1  and T 2  should be as small as possible, while retaining a degree of freedom by considering the variances of those distributions. The following energy can be defined in order to maximize the likelihood of these densities on their associated domain: 
                 E   ⁡     (       Ω   i     ,     σ     i   ,   d     2     ,     T   i       )       =       ∑     i   =   1     2     ⁢       ∫     Ω   i       ⁢       -   log     ⁢           ⁢       p     d   ,   i       ⁡     (       d   2     ⁡     (       T   ⁡     (   x   )       ,     T   i       )       )       ⁢           ⁢     ⅆ   x             ,     
     ⁢   where                 p     d   ,   i       =       1       2   ⁢   π   ⁢           ⁢     σ     i   ,   d     2           ⁢       exp   ⁡     (       -       d   2     ⁡     (     T   ,     T   i       )           2   ⁢     σ     i   ,   d     2         )       .             
The level set distance function whose zero isosurface coincides with Γ is denoted by φ:Ω→R 3 . One can define an energy functional incorporating the variance and with a regularity constraint on Γ, using H ε (z), the regularized version of the Heaviside function, as follows:
 ∫ Ω −log p d,1 (d 2 (T(x),T 1 ))H ε (φ)−log p d,2 (d 2 (T(x),T 2 ))(1−H ε (φ))+α|∇H ε (φ)|dx,  (2) 
which is equivalent to
 
                 -       ∫   Ω     ⁢     log   ⁢           ⁢       p     d   ,   1       ⁡     (       d   2     ⁡     (       T   ⁡     (   x   )       ,     T   1       )       )       ⁢           ⁢     ⅆ   x           -       ∫     Ω   c       ⁢     log   ⁢           ⁢       p     d   ,   2       ⁡     (       d   2     ⁡     (       T   ⁡     (   x   )       ,     T   2       )       )       ⁢           ⁢     ⅆ   x         +     α   ⁢        C            =     
     ⁢         ∫   Ω     ⁢       (           d   2     ⁡     (       T   ⁡     (   x   )       ,     T   1       )         2   ⁢     σ     1   ,   d     2         +     log   ⁢       2   ⁢   π   ⁢           ⁢     σ     i   ,   d     2             )     ⁢           ⁢     ⅆ   x         +       ∫     Ω   c       ⁢       (           d   2     ⁡     (       T   ⁡     (   x   )       ,     T   2       )         2   ⁢     σ     2   ,   d     2         +     log   ⁢       2   ⁢   π   ⁢           ⁢     σ     2   ,   d     2             )     ⁢           ⁢     ⅆ   x         +     α   ⁢          C        .               
The derivation of the Euler-Lagrange equations for this class of energy yields the following evolution for φ:
 
                   ∂   ϕ       ∂   t       =           δ   ɛ     ⁡     (     ϕ   ⁡     (   x   )       )       ⁡     [       α   ⁢           ⁢     ∇     ·       ∇           ⁢   ϕ            ∇   ϕ                  +       1   2     ⁢   log   ⁢       p     d   ,   2         p     d   ,   1             ]       ⁢     ∀     x   ∈   Ω           ,         
where δ ε (φ(x)) is the regularized version of the Dirac function.
 
     When considering the region terms, the initialization is important and in many cases, several seeding points have to be set manually to avoid having the surface evolving to and remaining in a local minima. This can be overcome by using a fractional anisotropy measure: 
                 A   ⁡     (     T   ⁡     (   x   )       )       =             (       λ   1     -     λ   2       )     2     +       (       λ   2     -     λ   3       )     2     +       (       λ   1     -     λ   3       )     2             2     ⁢         λ   1   2     +     λ   2   2     +     λ   3   2               ,         
where the λ&#39;s are eigenvalues of the tensors. Note that this anisotropy measure is computed for each pixel within a region. An additional term is then defined to impose a given distribution of the anisotropy inside each region. Let p α,1  and p α,2  be the probability distribution functions of the anisotropy inside and outside the region bounded by the curve C, approximated by Gaussian densities defined as
 
                 p     a   ,   i       =       1       2   ⁢   π   ⁢           ⁢     σ     i   ,   A     2           ⁢     exp   ⁡     (                A   ⁡     (   T   )       -     A   _            2       2   ⁢     σ     i   ,   A     2         )           ,         
where Ā is the mean value of A(T) over the region i, and σ i,A  is the variance of p i,A  over the region i. Then, the partitioning is obtained by minimizing:
 −∫ Ω log p α,1 (A(T(x)))H ε (φ)+log p α,2 (A(T(x)))(1−H ε (φ))dx, 
where Ω∈R 3  is the image domain, and H ε  is the Heaviside step function. This term is added to the objective function defined in equation (2), above. A new energy functional can be obtained for the level set function φ composed of three terms where the influence from the distribution of tensor distance and the fractional anisotropy can be controlled by adjusting a weight β between zero and one:
 
                     ∫   Ω     ⁢       [         -     (     1   -   β     )       ⁢   log   ⁢           ⁢       p     d   ,   1       ⁡     (       d   2     ⁡     (       T   ⁡     (   x   )       ,     T   1       )       )       ⁢       H   ɛ     ⁡     (   ϕ   )         -       (     1   -   β     )     ⁢   log   ⁢           ⁢       p     d   ,   2       ⁡     (       d   2     ⁡     (       T   ⁡     (   x   )       ,     T   2       )       )       ⁢     (     1   -       H   ɛ     ⁡     (   ϕ   )         )       -     β   ⁢           ⁢   log   ⁢           ⁢       p     a   ,   1       ⁡     (     A   ⁡     (     T   ⁡     (   x   )       )       )       ⁢       H   ɛ     ⁡     (   ϕ   )         -     β   ⁢           ⁢   log   ⁢           ⁢       p     a   ,   2       ⁡     (     A   ⁡     (     T   ⁡     (   x   )       )       )       ⁢     (     1   -       H   ɛ     ⁡     (   ϕ   )         )       +     α   ⁢          ∇       H   ɛ     ⁡     (   ϕ   )                  ]     ⁢     ⅆ   x               (   3   )               
In practice, a small weight on the anisotropy term is sufficient for the surface avoid a local minima.
 
     According to an embodiment of the invention, equation (3) can be efficiently solved by a fast level-set method. A level set formulation represents a front as the zero level set of a function defined in a higher dimensional space. Consider a closed moving interface Γ(t) in R 3 . Let Ω(t) be the region (possibly multi-connected) that Γ(t) encloses. One can associate with Ω(t) an auxiliary function φ(x,t), called the level set function, which satisfies: 
             {               ϕ   ⁡     (     x   ,   t     )       &gt;   0         for         x   ∈   Ω                 ϕ   ⁡     (     x   ,   t     )       =   0         for         x   ∈     ∂   Ω                   ϕ   ⁡     (     x   ,   t     )       &lt;   0         for         x   ∈     Ω   c                         
where x∈R 3 ,t∈R + .
 
     Conversely, knowing φ allows one to locate the interface by finding the zero level set of φ, which is Γ(t)={x:φ(x,t)=0}. So moving the interface is equivalent to updating φ, which can be done by solving a Hamilton-Jacobi equation, such as the time evolution equation for φ derived above. 
     But for some cases, one does not need the value of φ, only its sign. From the optimization point of view, this opens the possibility of other methods for solving a minimization problem directly and much more quickly. For instance, for a 2-phase image segmentation, one seeks a particular partition of a given image into two regions, one representing the objects to be detected and one representing the background. Assuming that the image u 0  is a 2-phase image with piecewise constant values u 0   i  and u 0   o  and that the object to be detected is represented by the value u 0   i , and that C 0  denotes the boundary of the object, then the fitting energy is defined as:
 
 E   1 ( C )+ E   2 ( C )=∫ inside(C)   |u   0   −c   1 | 2 +∫ outside(C)   |u   0   −c   2 | 2 ,
 
where C is any other variable curve, and the constants c 1 , c 2  are the averages of u 0  inside and outside of C respectively. The fitting energy will be minimized if C=C 0 . In this model according to an embodiment of the invention, there can also be a regularizing term, such as the length of C, or the area inside C, to control the smoothness of the boundary. Therefore, the energy E(C, c 1 , c 2 ) is defined by:
 
 E ( C,c   1   ,c   2 )=α·(length( C ))+λ 1 ∫ inside(C)   |u   0   −c   1 | 2 +λ 2 ∫ outside(C)   |u   0   −c   2 | 2 .
 
     If a level set represents C, that is, C is the zero level set of a Lipschitz function φ:R 2 →R, then φ can replace the unknown variable C, and the energy functional E(C, c 1 , c 2 ) can be written as:
 
 E ( H (φ), c   1   ,c   2 )=α(∫ Ω   |∇H (φ)|)+λ 1 ∫ Ω   |u   0   −c   1 | 2   H (φ) dx+λ   2 ∫ Ω   |u   0   −c   2 | 2 (1 −H (φ)) dx,  
 
where c 1 , c 2  are also functions of H(φ). In this equation, the two fitting terms are easy to compute directly, while ∫|∇H(φ)|dx can be approximated by:
 
                 ∑     i   ,   j       ⁢           (       H   ⁡     (     ϕ       i   +   1     ,   j       )       -     H   ⁡     (     ϕ     i   ,   j       )         )     2     +       (       H   ⁡     (     ϕ     i   ,     j   +   1         )       -     H   ⁡     (     ϕ     i   ,   j       )         )     2           ,         
where φ i,j  is the value of φ at the i,j th  pixel. The summand can only take the values 0, 1 or √{square root over (2)}, depending on whether the 3 distinct pair of points from the set {φ i,j , φ i+1,j , φ i,j+1 } belong to the same or different regions. Thus, the length term can be easily computed knowing only H(φ), and there is no need to know φ. This computed value can be interpreted as the discretized length of the zero level set.
 
     Consider a two phase image, where one can represent an object by A (possibly multiconnected), the background by B, and the corresponding values for A and B by a and b. An initial partition is given by φ 1 &gt;0 and φ 2 &lt;0, where there are m points in φ 1  and n points in φ 2 , and c i  is the average for φ i , i=1, 2. Then, if |c 1 −c 2 |&gt;c for some constant c, and if the condition 
                   n   -   1     n     ⁢     m     m   +   1         ≤         (     a   -     c   2       )     2         (     a   -     c   1       )     2       ≤       m     m   -   1       ⁢       n   +   1     n             
is not satisfied, then an algorithm (with α=0) according to an embodiment of the invention can be shown to converge in one sweep for a 2-phase image, using either Jacobi or Gauss-Seidel iterations. In practice, when applying the algorithm, one can directly apply the algorithm to the model with the length term included, or one can first consider α=0, then followed by α&gt;0 to have the full effect of regularization. Another choice is to consider α=0, followed by a PDE-based algorithm.
 
     An outline of a fast algorithm for solving the above model for scalar valued image is as follows: 
     Step 1. Initialize: given an initial partition of the image, set φ=1 for one part and φ=−1 for another part, and compute the value of the energy E according to φ. 
     Step 2. Advance: assume that the value of current pixel x is u, c 1  and c 2  are the pixel value averages for φ=1 and φ=−1, respectively, m and n are the number of pixels for φ=1 and φ=−1. If φ(x)=1, then compute the difference between the new and the old energy: 
               Δ   ⁢           ⁢     E   12       =                u   -     c   2            2     ⁢     n     n   +   1         -              u   -     c   1            2     ⁢       m     m   -   1       .               
If ΔE 12 &lt;0, then change φ(x) from 1 to −1. Similarly, if φ(x)=−1 case, compute ΔE 21 :
 
                 Δ   ⁢           ⁢     E   21       =                u   -     c   1            2     ⁢     m     m   +   1         -              u   -     c   2            2     ⁢     n     n   -   1             ,         
and if ΔE 21 &lt;0, change φ(x) from −1 to 1. If there is no change in ΔE, φ(x) remains unchanged. In 2D, if the length term is also considered, then the change of the length is easy to compute since only four neighbor points will be affected by a change of point value. The pixels can be looped through in any prescribed order, and can be updated using any iteration technique as is known in the art, such as Gauss-Seidel or Jacobi iteration.
 
     Step 3. Repeat step 2 until the change in energy E is sufficiently small or zero. 
     Now, when applying the above fast level-set method to equation (3), some care is needed to handle the tensors and variance. Note that there are three terms in equation (3). 
     The first term in equation (3) represents the contribution of the tensor distance distribution where the tensor distance is computed between the tensor field and the segmentation mean values T 1  and T 2 . If φ(x)=1, then the difference between the new and the old energy caused by this term is computed as: 
               Δ   ⁢           ⁢     E     12   ,   d         =       (     1   -   β     )     ⁢       (             d   2     ⁡     (       T   ⁡     (   x   )       ,     T   2       )         2   ⁢     σ     2   ,   d     2         ⁢     n     n   +   1         -           d   2     ⁡     (       T   ⁡     (   x   )       ,     T   1       )         2   ⁢     σ     1   ,   d     2         ⁢     m     m   -   1           )     .             
Similarly, if φ(x)=−1 case, ΔE 21,d  can be computed as:
 
     
       
         
           
             
               Δ 
               ⁢ 
               
                   
               
               ⁢ 
               
                 E 
                 
                   21 
                   , 
                   d 
                 
               
             
             = 
             
               
                 ( 
                 
                   1 
                   - 
                   β 
                 
                 ) 
               
               ⁢ 
               
                 
                   ( 
                   
                     
                       
                         
                           
                             d 
                             2 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 T 
                                 ⁡ 
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                               , 
                               
                                 T 
                                 1 
                               
                             
                             ) 
                           
                         
                         
                           2 
                           ⁢ 
                           
                             σ 
                             
                               1 
                               , 
                               d 
                             
                             2 
                           
                         
                       
                       ⁢ 
                       
                         m 
                         
                           m 
                           + 
                           1 
                         
                       
                     
                     - 
                     
                       
                         
                           
                             d 
                             2 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 T 
                                 ⁡ 
                                 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                               , 
                               
                                 T 
                                 2 
                               
                             
                             ) 
                           
                         
                         
                           2 
                           ⁢ 
                           
                             σ 
                             
                               2 
                               , 
                               d 
                             
                             2 
                           
                         
                       
                       ⁢ 
                       
                         n 
                         
                           n 
                           - 
                           1 
                         
                       
                     
                   
                   ) 
                 
                 . 
               
             
           
         
       
     
     The second term in equation (3) represents the contribution of the fractional anisotropy distribution. If φ(x)=1, then the difference between the new and the old energy caused by this term is computed as: 
               Δ   ⁢           ⁢     E     12   ,   A         =       β   ⁡     (                    A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     2            2       2   ⁢     σ     2   ,   A     2         ⁢     n     n   +   1         -                  A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     1            2       2   ⁢     σ     1   ,   A     2         ⁢     m     m   -   1           )       .           
Similarly, if φ(x)=−1 case, the change ΔE 21,A  is computed as:
 
               Δ   ⁢           ⁢     E     21   ,   A         =       β   ⁡     (                    A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     1            2       2   ⁢     σ     1   ,   A     2         ⁢     m     m   +   1         -                  A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     2            2       2   ⁢     σ     2   ,   A     2         ⁢     n     n   -   1           )       .           
The third term represents the smoothness of the boundary. In 2D, suppose the location of the point being changed is x=(x 1 ,x 2 ), then the change in the energy caused by this term using a calculation using forward difference in φ is:
   ΔE   ij,r =α(√{square root over ((φ( x   1 +1 ,x   2 )−φ j ) 2 +(φ( x   1   ,x   2 +1)−φ j ) 2 )}{square root over ((φ( x   1 +1 ,x   2 )−φ j ) 2 +(φ( x   1   ,x   2 +1)−φ j ) 2 )}−√{square root over ((φ( x   1 +1 ,x   2 )−φ i ) 2 +(φ( x   1   ,x   2 +1)−φ i ) 2 )}{square root over ((φ( x   1 +1 ,x   2 )−φ i ) 2 +(φ( x   1   ,x   2 +1)−φ i ) 2 )}) 
where φ i  is the level set value for region i, φ j  is the level set value for region j, either φ i =1 and φ j =−1 or φ i =−1 and φ j =1. The case for 3D is similar: suppose the location point being changed is x=(x 1 , x 2 , x 3 ), then
 
               Δ   ⁢           ⁢     E     ij   ,   r         =     α   ⁡     (                   (       ϕ   ⁡     (         x   1     +   1     ,     x   2     ,     x   3       )       -     ϕ   j       )     2     +                   (       ϕ   ⁢     (       x   1     ,       x   2     +   1     ,     x   3       )       -     ϕ   j       )     2     +       (       ϕ   ⁡     (       x   1     ,     x   2     ,       x   3     +   1       )       -     ϕ   j       )     2               -                 (       ϕ   ⁡     (         x   1     +   1     ,     x   2     ,     x   3       )       -     ϕ   i       )     2     +       (       ϕ   ⁡     (       x   1     ,       x   2     +   1     ,     x   3       )       -     ϕ   i       )     2     +                 (       ϕ   ⁡     (       x   1     ,     x   2     ,       x   3     +   1       )       -     ϕ   i       )     2               )             
and the above calculation can also based on central difference in φ.
 
     The total change of the energy function obtained by combining all the three terms decides how to change the sign of φ(x): that is, if ΔE Tot &lt;0, the total energy decreases, then the region membership of the point in question is changed. 
     A flow chart of an exemplary method for tensor field segmentation, according to an embodiment of the invention, is presented in  FIG. 7 . At step  70 , a diffusion tensor image acquires through DT-MRI is provided. In DT-MRI, what is measured is the diffusion weighted echo intensity image (DWI) S l  for different directions l. These directions are related to the diffusion tensor T through the following equation:
 
 S   l   =S   0  exp(− b   l   :T )= S   0  exp(−Σ i=1   3 Σ j=1   3   b   l,ij   T   ij ),
 
where b l  is the diffusion weighting of the l-th magnetic gradient, and “:” denotes the generalized inner product for matrices. Given several non-collinear diffusion weighted intensity measurements, T can be estimated via multivariate regression techniques and a diffusion tensor image is constructed. At step  71 , an initial partition of the image into separate regions is provided. According to one embodiment of the invention, the image is partitioned into two regions, with the auxiliary function φ being set to 1 inside an initial partition curve, and being set to −1 outside the curve. At step  72 , the mean value of the tensor fields and the mean value of the fractional anisotropy field are then calculated over the two regions, and an initial value of the energy is calculated, using the energy functional defined above. At step  73 , the iteration is advanced by looping through all pixels. For each pixel at location x near the region boundary, if φ(x)=1, then change it so φ(x)=−1, and vice versa if φ(x)=−1, and then compute the change in energy. In general, the change in the tensor distribution energy ΔE ij,d  from changing the membership of a pixel from region i to region j is
 
               Δ   ⁢           ⁢     E     ij   ,   d         =       (     1   -   β     )     ⁢       (             d   2     ⁡     (       T   ⁡     (   x   )       ,     T   j       )         2   ⁢     σ     j   ,   d     2         ⁢       n   j         n   j     +   1         -           d   2     ⁡     (       T   ⁡     (   x   )       ,     T   i       )         2   ⁢     σ     i   ,   d     2         ⁢       n   i         n   i     -   1           )     .             
and the corresponding change in the factional anisotropy distribution energy ΔE ij,A  is
 
               Δ   ⁢           ⁢     E     ij   ,   A         =     β   ⁡     (                    A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     j            2       2   ⁢     σ     j   ,   A     2         ⁢       n   j         n   j     +   1         -                  A   ⁡     (     T   ⁡     (   x   )       )       -       A   _     i            2       2   ⁢     σ     i   ,   A     2         ⁢       n   i         n   i     -   1           )             
where n i  and n j  are the respective number of pixels for each region. In this exemplary method which achieves a piecewise constant segmentation, all the variances can be set to 1 and the energy change can be computed accordingly. If the total energy change, which includes the above two energy changes and the changes in boundary smoothness, is decreased by changing the region membership of a pixel, then the region membership of the pixel should be changed at step  74 . At step  75 , the mean value of the tensor fields and the mean value of the fractional anisotropy field are updated when the region changes. At step  76 , if the magnitude of the change in energy is greater than a predetermined value, the process returns to step  73  to perform another iteration. An appropriate threshold can be determined experiment by one of ordinary skill in the art. Once the iterations fail to produce a change in energy, the resulting segmentation is obtained at step  77 , and the process is finished.
 
     An implementation of an embodiment of the invention not incorporating variance has been tested on a DELL precision 670 workstation with a Xeon™ 2.8 Ghz CPU with 2.0 GB of RAM, and running under the Microsoft Windows XP Professional Version 2002, Service Pack 2 operating system. Segmentation results of three brain structures are presented in  FIGS. 1-6 , according to an embodiment of the invention. Each segmentation is presented in two figures, one being paired with an FA slice as a context, the other being simply the segmentation itself.  FIGS. 1 and 2  depict segmentation of the Corpus Callosum  10 ,  FIGS. 3 and 4  depict segmentation of the Corona Radiata  30 , and  FIGS. 5  and  6  depict segmentation of the ventricle  50 . All the above three segmentation were achieved within 5 seconds. 
     It is to be understood that the present invention can be implemented in various forms of hardware, software, firmware, special purpose processes, or a combination thereof. In one embodiment, the present invention can be implemented in software as an application program tangible embodied on a computer readable program storage device. The application program can be uploaded to, and executed by, a machine comprising any suitable architecture. 
       FIG. 8  is a block diagram of an exemplary computer system for implementing a level-set based tensor field segmentation according to an embodiment of the invention. Referring now to  FIG. 8 , a computer system  81  for implementing the present invention can comprise, inter alia, a central processing unit (CPU)  82 , a memory  83  and an input/output (I/O) interface  84 . The computer system  81  is generally coupled through the I/O interface  84  to a display  85  and various input devices  86  such as a mouse and a keyboard. The support circuits can include circuits such as cache, power supplies, clock circuits, and a communication bus. The memory  83  can include random access memory (RAM), read only memory (ROM), disk drive, tape drive, etc., or a combinations thereof. The present invention can be implemented as a routine  87  that is stored in memory  83  and executed by the CPU  82  to process the signal from the signal source  88 . As such, the computer system  81  is a general purpose computer system that becomes a specific purpose computer system when executing the routine  87  of the present invention. 
     The computer system  81  also includes an operating system and micro instruction code. The various processes and functions described herein can either be part of the micro instruction code or part of the application program (or combination thereof) which is executed via the operating system. In addition, various other peripheral devices can be connected to the computer platform such as an additional data storage device and a printing device. 
     It is to be further understood that, because some of the constituent system components and method steps depicted in the accompanying figures can be implemented in software, the actual connections between the systems components (or the process steps) may differ depending upon the manner in which the present invention is programmed. Given the teachings of the present invention provided herein, one of ordinary skill in the related art will be able to contemplate these and similar implementations or configurations of the present invention. 
     While the present invention has been described in detail with reference to a preferred embodiment, those skilled in the art will appreciate that various modifications and substitutions can be made thereto without departing from the spirit and scope of the invention as set forth in the appended claims.