NLM Home Page VHP Home Page

Next: References Previous: Title Page Contents: Conference page

Hybrid Segmentation of the Visible Human Data

 

     

Introduction

      In this paper we develop and test new hybrid methods for segmenting the Visible Human color cryosection and radiological (patient) data (e.g. CT, MRI, PET). The novelty stems from the integration of region-based and deformable model-based segmentation methods with a variety of region-based and statistical methods which aims toward the development of segmentation methods that yield high precision, accuracy and efficiency.This is a collaborative project between the University of Pennsylvania and Columbia University, a part of a larger effort to provide a fully implemented and tested Visible Human Project Segmentation and Registration Toolkit (Insight).

     We have already developed powerful region and contour-based segmentation tools that we plan to further extend to meet the above goals ( [5], [7], [8], [9], [17], [18]). The Visible Human data require a revision of the existing segmentation methods to accommodate both the color and unusual detail that is present in the dataset [13]. Also, new segmentation methods tested on the Visible Human data might improve segmentation and understanding of the clinical radiological images.

     Radiological image data of internal organs and in particular of the lungs, the heart and the brain are routinely used at the University of Pennsylvania, where the Interactive Virtual Environment for Modeling Anatomy and Physiology has been developed (Project funded by NLM, used the Visible Human Datasets in addition to Physiology to go beyond the anatomy of the cardiopulmonary system for educational purposes.), in addition to 3DVIEWNIX [16], a Unix-based software system for the visualization, manipulation, and analysis of multidimensional multi-parametric, multi-modality images. This latter system is freely available via Internet and there are hundreds of sites worldwide where the system is used for various medical and non- medical applications.

     Highly detailed segmentations and 3D visualizations of the Visible Human dataset are used in interactive applications, such as the Vesalius™ Project (Trademark held by Columbia University. The Vesalius Project is named after Andreas Vesalius, a sixteenth century anatomist whose work laid the foundation for all subsequent anatomical research. The Columbia University Health Sciences Library owns several first editions of his work.) at Columbia University, an interdisciplinary effort to create a network-based, platform-independent electronic learning environment for anatomy and other computer based medical training applications [6].

     Integration of boundary-based and region-based methods amplifies the strengths but reduces the weaknesses of both approaches. Our preliminary results are promising. In the next section we will give an overview of our current efforts towards the development of robust and efficient segmentation methods suitable for segmenting the Visible Human Data as well as for use in clinical applications

Approach and Methodology

     Automatic internal organ segmentation from various imaging modalities, including color Visible Human cryosection data, is a very important yet open research problem. Over the past several years, a variety of segmentation methods have been developed. Boundary-based techniques such as snakes [11] start with a deformable boundary and attempt to align this boundary with the edges in the image. The solution to these systems generally involves minimizing an energy functional which quantifies the shape of the model and image information near the boundary of the model. To avoid becoming stuck in local minima, most model-based techniques require that the model be initialized near the solution or supervised by an interactive interface or some higher-level reasoning entity. User steered methods such as live wire that are in day- to-day clinical research use that find a global optimum through dynamic programming have also been developed. Region- based or statistical techniques such as region growing ( [3], [7]) or MAP-based methods [2] assign membership to objects based on homogeneity statistics. The advantage here is that image information inside the object is considered as well as on the boundaries. However, there is no provision in the region based framework for including the shape of the region in the decision making process, which can lead to noisy boundaries and holes in the interior of the object.

     Like several other recent approaches ( [4], [20]), our design integrates the boundary and region-based techniques into a hybrid framework. By combining these techniques, hybrid approaches offer greater robustness than either technique alone. However, most previous work still requires significant initialization to avoid local minima. Furthermore, most of the earlier approaches use prior models for their region-based statistics, which we would rather avoid to increase usefulness in situations where a comprehensive set of priors may not be available.

     We have recently developed a new approach to internal organ segmentation that is based on the integration of region-based and physics-based boundary estimation methods ( [8], [9]). Starting from a single pixel within the interior of an object, we make an initial estimate of the object's boundary using the techniques of fuzzy connectedness [17] and clustering. A deformable surface model is then fitted to the extracted boundary data to fill in the missing boundary data and to override the spurious boundary data due to image noise. This is achieved by generalizing the formulation of our deformable models [12] to incorporate simple domain-specific knowledge.

     In this paper we will further develop this approach by integrating deformable models with fuzzy connectedness, and the region-based color segmentation method developed at Columbia [7] which we tested on the color cryosection Visible Human data. In addition, we will present a method based on an integration of Markov Random Field (Gibbs prior) and deformable model. In particular we develop two types of automatic segmentation algorithms: 1) Integration of fuzzy connectedness , Voronoi Diagram, and deformable models; and 2) Integration of Gibbs prior and deformable models.

     Hybrid Method I: Integration of Fuzzy Connectedness, Voronoi Diagram and Deformable Models

     We present a hybrid segmentation method which requires minimal manual initialization, where we integrate the fuzzy connectedness, Voronoi Diagram, and deformable models. We will start with fuzzy connectedness algorithm to generate a region with a fragment of tissue which we plan to segment. From the sample region, we will generate automatically homogeneity statistics for the Voronoi Diagram-based algorithm which will produce an estimation of the boundary in a number of iterations. The final boundary of the segmented region can be smoothed over with the deformable models method. We will give an overview of the component algorithms used in the hybrid segmentation algorithm below.

          Fuzzy Connectedness Algorithm: This method developed by Udupa and Samarasekera ,and described in detail in [17], uses the fact that medical images are inherently fuzzy. We define affinity between two elements in an image (e.g. pixels, voxels, spels [17]) via a degree of adjacency and the similarity of their intensity values. The closer the elements are and more similar their intensities are, the greater is the affinity between them. There are two important characteristics of a medical image. First, it has graded composition coming from material heterogeneity (e.g. bone consists of hard cortical tissue and softer cancellous tissues), blurring, noise and background variation. Second, the image elements that constitute an anatomical object hang together in a certain way (Figure 1). Both these properties, hanging togetherness and graded composition are fuzzy properties. The aim of fuzzy connectedness is to capture the global hanging togetherness using image-based local fuzzy affinity.

      Let us define a scene over a digital space (Zn , a) as a pair W=(C,f), where C is an n-dimensional array  of spels (elements)  and  f: C ® [0,1]. Fuzzy affinity k is any reflexive and symmetric fuzzy relation in C , that is:

                                           k = {((c,d), mk(c,d)) | c,d Î C}
                                           mk::  C ´ C ® [0,1]
                                           mk(c,c) = 1, for all c Î C
                                           mk(c,d) = mk(d,c) , for all c,d Î C.  

The general form of  mk can be written as follows. For all c,d Î C,                         

                               mk(c,d)  =  g(ma(c,d) , my(c,d),mf(c,d) , c, d)

where:

ma(c,d) represents the degree of adjacency of c and d;
my(c,d) represents the degree of intensity homogeneity of c and d;
mf(c,d) represents the degree of similarity of  the intensity features of c and d to expected object features.  

      Fuzzy k-connectedness K is a fuzzy relation in C, where mK(c,d ) is the strength of the strongest path between c and d, and the strength of a path is the smallest affinity along the path, Figure 2.

      To define the notion of a fuzzy connected component, we need the following hard binary relation Kq based on the fuzzy relation K. Let W=(C,f)  be a membership scene over a fuzzy digital space (Zn , a), and let k be a fuzzy spel affinity in W. We define a (hard) binary relation Kq in C as

 

                                              

Let Oq be an equivalence class [14, Chap.10] of the relation Kq in C. A  fuzzy k-component  of C of strength q is a fuzzy subset of C defined by the membership function

 

                                             

The equivalence class Oq  Ì  C, is shown in Figure 3,  such that for any c,d Î C, mk(c,d) ³ q, q Î [0,1], and for any e Î C - Oq, mk(c,d) < q. We use the notation to denote the equivalence class of  Kq    that contains  for any  Î C. The fuzzy k-component of C that contains  , denoted , is a fuzzy subset of C whose membership function is

 

                                              

A fuzzy kq-object of W is a fuzzy  k-component of  W of strength q. For any spel  Î C, a fuzzy kq-object of  W that contains  is a fuzzy k-component of  W of strength q that contains .Given k, , q, and W,  a fuzzy kq-object  of W of  strength q Î [0,1]  containing , for any  Î C,  can be computed via dynamic programming [17]. This means that for finding the fuzzy kq-object containing  Î C, it is not necessary to compute mK(c,d) for each possible pair (c,d) of spels in C. It is sufficient to compute mK (,c) for each spel c Î C, and achieve a vast computational reduction.

           Voronoi Diagram-based Algorithm:  This algorithm, which is described in detail in [7], is based on repeatedly dividing an image into regions using Voronoi diagram and clasifying the Voronoi regions based on a selected homogeneity classifier for the segmented anatomical tissue. We will use the algorithm as a component in the hybrid method where the classifiers for different tissue type will be generated automatically from the region segmented by the fuzzy connectedness method.

      Voronoi diagram and Delaunay triangulation play a central role in the algorithm. In 2D, the Voronoi diagram for a set V of points  is a partition of the Euclidean plane into Voronoi regions of points closer to one point of V than to any other seed point [14], see Figure 4(a).

For any ,

Similarly, we define the Voronoi diagram in 3D:

Two Voronoi regions are adjacent if they share a Voronoi edge. The Delaunay triangulation, DT, of V is a dual graph of  the Voronoi diagram of V, obtained by joining two points whose Voronoi regions are adjacent, see Figure 4(b).

The following is a pseudocode for the Voronoi diagram-based algorithm. For a given color Visible Human 2D cryosection image:

  1. Input a set of n points.
  2. Compute Voronoi diagram.
  3. For each Voronoi region, classify it as interior or exterior using a homogeneity classifier.
  4. Label a subset of the exterior regions, adjacent to interior regions, as boundary regions.
  5. Compute Delaunay triangulation and display segments which connect Voronoi seed points of the boundary regions.
  6. Add seeds to the Voronoi edges of boundary regions.
  7. GoTo Step2 until the algorithm converges to a stable state or until the user chooses to quit.
           Hybrid Method I. The hybrid algorithm integrates three methods, the fuzzy connectedness, Voronoi diagram-based algorithm and the deformable models. We will outline the algorithm first, and explain the component steps later. The fuzzy connectedness algorithm is used to segment a fragment of the target tissue. From the sample, a set of statistics is generated automatically, in RGB and HVC color spaces, to define the homogeneity operator. The homogeneity operator will be used as a multi-channel classifier for the Voronoi diagram-based algorithm. We will use the deformable model, described in the next section, to determine the final (3D) boundary. We must note that the deformable model method can be used, in the future , not only for smoothing of the final boundary but also for "tuning up" the homogeneity classifier. Below, we outline the hybrid method:
Step 1. We run the fuzzy connectedness algorithm to segment a sample of the target tissue, and generate statistics, average and variance, in three "strong" color channels, in two color spaces, RGB and HVC.
Step 2. Run the Voronoi diagram-based algorithm using multiple color channels, until it converges:
(a) For each Voronoi region, classify it as interior/exterior/boundary region using multi-channel homogeneity operator.
(b) Compute Delaunay triangulation and display segments which connect boundary regions.
(c) Add seeds to Voronoi edges of Voronoi boundary regions.
(d) GoTo Step 2(a) until the algorithm converges to a stable state or until the user chooses to quit.
Step 3. Use the deformable model to determine the final (3D) boundary.

      Implementation of Step 1. To initialize the fuzzy connectedness algorithm and establish the mean and standard deviation of spel values and their gradient magnitudes, the user collects the pixels within the region of interest, by clicking on the image and selecting at each time a square region with 7x7 pixels, see Figure.5. Then an initial seed spel    is selected to compute the fuzzy connectedness, using the dynamic programming approach [17]. We determine q, q Î [0,1], by letting the user to select interactively its threshold value, such that the initially segmented sample of the target tissue resembles the shape of the underlying image, see Figure 6.
      For a binary image with a roughly segmented sample of a tissue, we generate the "strongest" three channels in two color spaces, RGB (Red, Green, Blue) and HVC (Hue, Value, Chroma), for average and variance, respectively. The HVC system represents images with hue, value and chroma (saturation) levels, which in effect separate the luminant and chromatic constituents of a color [15]. First, we define for the binary image, the smallest enclosing rectangle, a region of interest (ROI), in which we identify the segmented image and its background. Within the ROI, we calculate the mean and variance in each of the six color channels (R,G,B,H,V,C) for the object and its background, respectively. Then three channels with the largest difference (relative difference) in mean value between the object and its background are selected. Another three channels with largest (relative) difference in variance value are also selected. The homogeneity operator for the Voronoi diagram-based algorithm uses the expected mean/variance values of the object together with tolerance values, computed for each selected channel, for classifying the internal and external region. This set of statistics plays a role of a homogeneity operator, see Figure 7.

      Implementation of Step 2. We build an initial Voronoi diagram by generating some number of random seed points (Voronoi points) and then run the QuickHull, [1], to calculate the Voronoi diagram. Once the initial Voronoi diagram has been generated, the program visits each region to accumulate classification statistics and makes a determination as to the identity of the region. For each Voronoi region, the mean/variance value for the preselected channels (strongest 3 for mean and another strongest for variance) are computed, if they are similar with the object ( within the tolerance range) it is marked as internal, otherwise external. Those external regions that have at least one internal neighbor are marked as boundary. Each boundary region is divided for next iterations until the total number of pixels within it is less than 20, see Figure 7, Figure8 and Figure 9.

Results

      In this section we present the results from experiments using Hybrid Method I. We tested different tissue types (muscle [Figure10], brain tissue [Figure11], and fat [Figure12] (including patient MRI fat [Figure13] data provided to us by Dr. S. Heymsfield)), to show that the algorithm can handle a variety of human tissue. In [Figure11], we segment brain gray matter using the Hybrid Method, and we would like to note that we have used fuzzy connectedness extensively, in 3D, to segment gray and white matter in patient image data [19]. Gray matter and white matter are fuzzy connected entities in 3D. With a few seed points they can be, and they are, routinely segmented. With multiple points, you can segment it well, even without using scale-based fuzzy connectedness.

      In this experiment we showed that the algorithm can generate automated statistical homogeneity classifiers for a variety of human tissue, and use the classifiers successfully in generating the outlines of the boundary of the target structures. We tested the method with the Visible Human data and as well with a sample of MRI patient data. These are just preliminary results and this experiment can lead to a rigorous hybrid segmentation methodology. The algorithm has to be validated, in the future, using ground truth anatomical structures obtained from a variety of medical image modalities.

     Hybrid Method II: Integration of Gibbs Prior and Deformable Models.

       Most medical images are Markov Random Field images, that is, the statistics of a pixel is related to the statistics of pixels in its neighborhood. According to the Equivalent Theorem proved by Hammersley and Clifford, the Markov Random Field is equivalent to the Gibbs field under certain conditions. Thus for medical images which are MRF images, their joint distribution can be written in the Gibbsian form, which is shown as follows:

where H(x) represents the energy function of image x, X is the set of all possible configuration of the image, Z is the normalization factor or partition function in the terminology of statistical mechanics. The local and global properties of images will be incorporated into the model by designing an appropriate energy function. The lower the value of the energy function, the higher the value of the Gibbs distribution, the better the image is fit to the prior distribution.

      We began the establishment of the Gibbs model by designing the energy function as:

              

 models the piecewise homogeneity statistics and  models the continuity of the boundary,  is the noise model term.

  1. Find the boundary according to the variance between objects and background. ()
  2. If the variance is not large enough or the object gradually turns into the background,  is capable of finding a smooth and continuous boundary estimation.
  3. Erase the noise using term .
      We use the iterated conditional mode (ICM) method to minimize the energy function. This will give an estimation of the object boundary which can be used by the deformable model. Our deformable model is a superellipsoid with local deformations. The global parameter of the model is its natural shape, and the local deformations determine its displacement from the natural shape for discrete nodes on the surface of the model. Given the reference shape s and displacement d, points on a model p are defined by:

     p = s + d

      To keep the continuity of the model surface, we impose a  continuous deformation strain energy on it.

,

where d is the node local deformation,  controls the local magnitude and ,  control the local variation of the deformation. We can calculate the stiffness matrix K of the model based on the strain energy. The model nodes will move under the influence of internal and external forces. The internal forces constrain the deformation of a model according to the stiffness matrix to maintain its natural shape, while the external forces make the model expand until they meet a balance with the internal forces. The dynamics can be described by the first order Lagrangian method:

,                                 

where f is the summation of external forces and boundary forces.

      The nodes of the model are forced outward in the direction of their normal vector. So the model will expand like a balloon being filled with gas. When the model is reaching the estimated boundary, it will also move under the influence of boundary forces, which is in the opposite direction of the balloon force. The Finite Element Analysis (FEA) method is used in our work to calculate the deformation on each node. The deformation is updated using Euler step:

,                          [10]

where t is the time step.

      When most nodes (>90%) on the surface stop moving, the model stops.

      The whole method can be summarized as follows:

  1. First create the Gibbs prior models using arbitrary parameters.
  2. Use the Gibbs prior models estimate the boundaries of objects.
  3. Fit the deformable models to this data according to the combination of balloon and the boundary forces.
  4. Recalculate the parameters for the Gibbs prior model according to the results of Step3. Then go back to Step 2.

      In this method we presented an new approach to segmentation of medical images by integrating a probabilistic model and a deformable model. By using the Gibbs prior, we generated a better estimation of the boundary and it was used as an initialization for the deformable model. Then, the deformable model provided updated parameters to the Gibbs prior model, and the iterative algorithm was applied recursively until we got a refined segmentation. The algorithm is efficient and is not sensitive to the selection of the initial parameters. See Figure 14, Figure 15, and Figure 16.

Acknowledgements

     This work was supported, in part, by a NLM contract on the "VHP Segmentation and Registration Toolkit" - NLM99-103/DJH. We also would like to thank Dr. S. Heymsfield for providing us with the MRI patient data, and Michael Downes for discussing with us implementation details of the Hybrid Method I.
Next: References Previous: Title Page Contents: Conference page
Office of High Performance Computing and Communications
Lister Hill National Center for Biomedical Communications
U.S. National Library of Medicine, 8600 Rockville Pike, Bethesda, MD 20894
National Institutes of Health
Department of Health & Human Services
Copyright and Privacy Policy
Last updated: 2 July 2001