New Variational Level-Set Classes with Region Fitting Energy in ITK
The level-set methodology involves numerically evolving a contour or surface according to a rate change partial differential equation (PDE). The contour or surface is embedded as the zero level-set of a higher dimensional implicit function, also called the level-set function Φ(x,t). Other techniques in the literature with the same objective include snakes and parametric surfaces. The main advantage of using level-set techniques over other methods is that they can handle complex shapes and topological changes, such as merging and splitting, with ease. Different constraints on the contour smoothness, speed, size and shape are easily specified.

Level-sets are very popular in the image analysis community for use in image segmentation, smoothing and tracking problems. There are two main classes of the level-set methodology - explicit and implicit methods. The explicit methods have a generic PDE to determine the temporal rate of change in the level-set function. It is expressed as a sum of speed functions that contain propagation, smoothing, and advection terms. Each of these terms is constructed from the image intensity function and involves the usage of image gradients. Segmentation is accomplished by detecting edges in the object of interest. Hence, these techniques are also referred as edge-based level-sets.

Implicit level-sets are driven by PDE that is derived by using a variational formulation. Here, an energy functional is first constructed and then minimized. An example of such a variant is that proposed by Chan and Vese [1]. In this variant, objects are segmented in images where edge information may be poor or entirely absent. The core idea is that the region statistics of the foreground region is much different from the background region. The region mean intensity is a popularly used statistic for segmentation as shown in Figure 1. These techniques are also referred to as region-based or active contours without edges. The current level-set framework in ITK is completely explicit with no support for the implicit techniques; for more information, see our submissions to the Insight Journal [4], [5] and [6].

(a) Dense filter iterations at 0, 10, 20 and 30

(b) Sparse filter iterations at 0, 10, 20 and 30.

Figure 1: Implicit segmentation of cells using mean intensity as a region statistic to separate foreground and background (a) Top row showing dense filter iterations and (b) bottom row showing the sparse filter iterations. Note the difference in the structures local to initialization in the sparse case.

Normally, the solution of the level-set update equation can be done on the entire image domain. The dense filter iteratively updates every single pixel in the image regardless of its position with respect to the zero level-set. As a result, the convergence is much quicker and structures far from the zero level-set spawn contours for segmentation. The sparse implementation of the level-set solver maintains a band of pixels around the zero level-set which are iteratively updated. The PDE is solved exactly on those pixels that are on the zero level-set or are immediate neighbors of the zero level-set. A user-defined bandwidth of pixels is used in maintaining the level structure and for calculating the derivatives. Pixels in the image beyond this band around the zero level-set are not considered. As a result, the evolution gradually proceeds from the initialization and converges on local structures. Structures far off from the zero level-set may remain unaffected.

Multi Object Segmentations with Level-Sets and KD-Tree Optimizations
In biomedical image analysis, we are often interested in segmenting more than a single object from a given image. This especially happens when the objects to be segmented are adjacent to each other and the delineation of one object automatically affects the neighboring object. In such situations, it makes sense to concurrently process their segmentation in order to optimally segment the objects. As a simple example, there is a significant interest in segmenting nuclei or cells in microscopy images as shown in Figures 5 and 6. Each image is acquired at high resolution and could contain thousands of cells. These cells often cluster in regions and appear to overlap. The challenge is to split these cells into individual components. An iterative (or linear) cell extraction procedure using level-sets can cause inconsistent splits since each level-set function does not compete with the neighboring cells. There could be an overlap of the level-set functions. In such cases, it is necessary to use multiphase methods for segmentation.

There are several research papers devoted to multiphase methods that optimize the number of level-set functions used for a generic case of N phases or objects [2]. We are largely motivated by microscopy applications where cells are segmented and tracked with constraints placed on the area, volume and shape of each individual cell. Each cell has a unique fluorescence intensity level, therefore it is best to have a unique level-set function per cell. This is the approach adopted in the popular work by Dufour, et al. [3].

Figure 2: Inheritance model of the dense solver

Figure 3: Inheritance model of the sparse solver

Figure 4: Inheritance model of the region-based level-set function

Figure 5: Space optimizations. (a) ROI defined around individual cells. (b) kd-tree structure constructed from ROI centroid.

Memory Optimization Using KD-Trees
Computationally, it is memory intensive to have N level-set functions defined on the image domain. For a large image with many small objects (such as cells in microscopy images), it becomes a memory-intensive problem. Hence, we make the implementation robust by defining region-of-interest (ROI), using spatial data structures such as the kd-trees and cached lookup tables. Each level-set function is first defined in a region-of-interest (ROI) within the image domain. This is illustrated in Figure 5(a). The ROI should encompass the object to be segmented and its extent is specified by the attributes’ origin and size. The spacing is the same as the feature or raw intensity image. This saves us considerable computational memory space. The centroid of each ROI region is then placed in a kd-tree structure. In the update for each level-set function, the overlaps of ROI regions are calculated by querying the kd-tree for the k-nearest neighbors as illustrated in Figure 5(b). This saves considerable computational time. Note that there is a cost associated with building the kd-tree that can be avoided for a small number of phases. We only instantiate the kd-tree mechanism of search upon user-request, or when there are more than 20 phases involved.

Usage Example
We consider an example for using multiphase level-set filters. After initialization, it is imperative that the user specify the number of level-set functions and set the feature image and initialization for each level-set function. We illustrate our example for N = 3.

Appropriate global settings of the level-set include the number of iterations, maximum permissible change in RMS values and whether to use image spacing.

Using a for-loop over all the level-set functions, we call the i-th difference function and set the corresponding attributes of that level-set function.

Tracking Multiple Objects Using Level-Sets
Many biological experiments that involve microscopic imaging require segmentation and temporal tracking of cells as part of the analysis protocol. For example, development biologists are interested in reconstructing cell lineages during embryonic development. Migratory behavior and rearrangement of cells is a fascinating topic of research. Cancer researchers track cells in colonies to determine growth kinetics and the effects of different chemical agents. The cell forms the fundamental biological entity of interest and its tracking is essential in these applications.

Dufour, et al. [3] proposed a solution to the tracking problem by using coupled active surfaces. In this method, each cell is represented by a unique level-set function. Energy functions involving the level-set contours are defined to partition the image into constant intensity background and constant intensity foreground components. The foreground components are regularized, in terms of their area and length, for smoothness. Several other properties such as continuity in volume and shape across time-points are maintained. Their solution as proposed is robust and elegant for small datasets since each cell only requires a unique level-set function of the same size as the image domain. In our implementation of the method in [6], we make use of performance optimization using ROI bounding boxes that now make the tracking filter scale up to larger datasets containing many cells.

Figure 6: 3D confocal images of a developing zebrafish embryo. (a-c) Raw images at 1, 5 and 10 time-points. (d-f) Tracking results at 1, 5 and 10 time-points.

Conclusion and Future Work
In this work, we developed region-based methods for variational contour evolution using the level-set strategy. We extended the methods to simultaneously segment and track multiple objects in the images and thereby use a kd-tree spatial portioning to be efficient. Our goal was to use these methods in cell segmentation and tracking analysis. We obtained very good results in our work. The main limitation of multi-object segmentations using the current level-set formulation is their interactions are not well defined. We often observe absorption of one level-set function by a neighboring one. We plan on further investigating these methods to make these methods robust.

References
[1]  T. Chan and L. Vese. An active contour model without edges. In Scale-Space Theories in Computer Vision, pages 141–151, 1999.
[2]  L. Vese and T. Chan. A multiphase level-set framework for image segmentation using the Mumford and Shah model. International
      Journal of Computer Vision, 50:271–293, 2002.
[3]  A. Dufour, V. Shinin, S. Tajbakhsh, N. Guillen-Aghion, J. C. Olivo-Marin, and C. Zimmer. Segmenting and tracking ?uorescent cells in
      dynamic 3-d microscopy with coupled active surfaces. IEEE Transactions on Image Processing, 14(9):1396–1410, 2005.
[4]  K. Mosaliganti, B. Smith, A. Gelas, A. Gouaillard, and S. Megason. Level-set segmentation: Active contours without edges. The Insight
      Journal, 2008.
[5]  K. Mosaliganti, B. Smith, A. Gelas, A. Gouaillard, and S. Megason. Segmentation using coupled active surfaces. The Insight Journal,
      2008.
[6]  K. Mosaliganti, B. Smith, A. Gelas, A. Gouaillard, and S. Megason. Cell Tracking using Coupled Active Surfaces for Nuclei and
      Membranes. The Insight Journal, 2008.

Acknowledgements
This work was funded by a grant from the NHGRI (P50HG004071-02) to found the Center for in-toto genomic analysis of vertebrate development. Benjamin Smith at Simon Fraser University, Vancouver, Canada first developed a working prototype of the code using ITK. It was significantly improved with bug corrections and enhanced with newer C++ constructs and optimized for its speed by the team comprising of Kishore Mosaliganti, Arnaud Gelas, Alexandre Gouaillard and Sean Megason at Harvard Medical School. We are currently working on the development of GoFigure2 – an open source application for biomedical image analysis, visualization and archival using ITK, VTK and Qt.

  Kishore Mosaliganti is a Research Fellow in the Megason Lab at Harvard Medical School where he He is currently developing algorithms for the extraction of zebrafish ear lineages using confocal microscopy images to be included in GoFigure2.

  Arnaud Gelas is a Research Associate in the Megason Lab at Harvard Medical School where he is in charge of the development of GoFigure2. Currently, he manages the software development process of GoFigure2.

  Alexandre Gouaillard is the President of CoSMo software which proposes software development services in medical image processing and modeling. Formerly he was a Research Associate at the Megason Lab at Harvard Medical School where he was in charge of the design and development of GoFigure2. He is now a PI at the Singaporean Immunology Network (SIgN) in a Systems Immunology research group using Complex Systems Modeling.

  Sean Megason is an Assistant Professor in the Department of Systems Biology at Harvard Medical School where he overlooks the development of GoFigure2. He is working in systems biology research on the developmental circuits of zebrafish.