Microscopy Image Analysis: Blob Segmentation using Geodesic Active Contours

An Insight Toolkit (ITK) processing framework for blob segmentation with applications to microscopy images is presented in this paper. Our algorithm is a reﬁnement of the work of Mosaliganti et al. [3] for splitting cell clusters. The basic idea is to incorporate as many cues as possible into developing a suitable level-set speed function so that an evolving contour exactly segments a cell/nuclei blob. We use image gradients, distance maps, multiple channel information and a shape model to drive the evolution. The framework consists of a linear pipeline of 6 new ITK ﬁlters which are applied in succession to generate the segmentation. The ﬁlters extract the cell foreground, construct the speed image, ﬁnd seed points for evolution and collect cell statistics from the segmentation. We include 2D/3D example code, parameter settings and show the results generated on confocal images of the zebraﬁsh embryo.


Introduction
Most biological experiments that involve microscopic imaging require cell profiling, counting, segmentation and tracking of cells.Cell/nuclei segmentation is often the first step in any analysis protocol and is becoming an increasingly popular topic in the image analysis community.The nucleus/cell forms the fundamental biological entity of interest.The problem arises when nuclei appear as overlapping or touching each other.Identifying each nucleus separately in a biologically consistent fashion is non-trivial.While some biochemical stains provide viable clues in the form of sharp color-space gradients at the boundaries, others exhibit a narrow neck at the site of overlap between two nuclei.We implement an approach that elegantly incorporates available cues and shapes into a viable segmentation pipeline [3].We use geodesic active contours incorporating shape priors in a level-set framework to provide nuclei segmentations.
The framework is implemented in 2D/3D for segmenting blobs (cells or nuclei) and for multiple stains that may be available.For example, often researchers use a nuclear stain and a separate cell membrane stain both of which are captured in two separate channels.The membrane channel can then be used to provides additional help in separating overlapping blobs etc.The implementation consists of 6 ITK filters whose description and use is documented in the rest of the document.Note that for convenience, we use the term nuclei but it can be replaced by any repeating small blob.We also assume that at this point, relevant preprocessing such as deconvolution, filtering and channel unmixing have been carried out to increase the SNR of the data.

Foreground Detection
In order to separate the nuclei, we employ a divide-and-conquer strategy to implement the filter itk::CellForegroundExtraction.The idea is to first separate nuclear foreground from the background before splitting them into individual entities.We identify nuclear material from the background by drawing information from the following sources: nuclei stain, membrane stain (identified membrane voxels which separate conjoined nuclei) and finally, the convex shape of the nuclei intensity function.A multiple threshold strategy is applied in each case.First, note that the nuclei stain primarily delineates nuclei material from the rest of the image.However, the staining can be inhomogenous and sometimes photon shot noise is captured in the acquisition.An upper and lower threshold is applied in this case.The upper threshold (m ThresholdCellmax) detects a nuclei voxel.The lower threshold rules out a nuclei voxel (m ThresholdCellmin).The intermediate range is further resolved as follows.
m ThresholdCellmin -lower threshold for nuclei stain given by the user in [0, 255].m ThresholdCellmax -upper threshold for nuclei stain given by the user in [0, 255].
Other stains, such as the membrane stain may be available, which can also separate the membrane from being classified as the nuclei foreground.A single threshold is applied to classify membrane pixels as being the background.
m ThresholdMembrane -threshold for the membrane stain set by the user in [0, 255] We also detect nuclear foreground by detecting the convex shape of the nuclear intensity function by comparing with a Gaussian kernel.Let C i represent a Gaussian kernel in the image I at location p, its domain (or extent of support) (δ given by m SigmaForm) and peak intensity value equal to unity.It is described by the following equation.
m SigmaForm -Standard deviation of the Gaussian kernel set by the user.Typically, half the size of the nucleus diameter.Usually 3-5µm.
At all pixels in the intermediate range, a Gaussian kernel with typical parameters, is placed and the Pearsons correlation (R(p)) is computed.A high value of the correlation corresponds to the center of the nucleus, a medium value to the background and a negative value to the region in-between the nuclei.Hence, a threshold is again applied on the correlation value to resolve the intermediate voxels.
The correlation (R(p)) is computed in a neighborhood patch (δ (p)) that is rectangular.The size of the patch is determined by m LargestCellRadius.
The detected foreground might still contain small holes, which have escaped detection, but which are filled using a hole-filling filter internally.
In order to speed up computation, the images can also be under-sampled during the intermediate computations and resampled back.The ratios for sub-sampling are set using m Sampling.The output consists of the foreground image (filter->GetOutput()) as well as the correlation map R(p) (filter->GetGaussCorrImage()).The correlation map has uses later in the pipeline.

Feature Extraction
A key step in using explicit level-set methods is the creation of a suitable feature or speed map f to drive the evolution.This function is designed to be minimized at nucleus boundaries and remains high elsewhere.The function is based on three ideas from literature that were developed independently namely, (i) inter-nuclear boundary gradients and (iii) morphological shape cues.We now describe the implementation of the filter itkCellFeatureGenerator.
Latest  The filter itk::CellForegroundExtraction is templated over the Input, Foreground and Correlation Image Types.It has a total of 6 parameters whose settings are described below.NOTE: The setting of parameters takes into account the spacing in the images.Hence, if the images are in µm, then the relevant parameters also need to be in the same units.
Let D(I f ) and D(I b ) represent the unsigned distance fields emanating from the nuclei contours that exist outside and within the contour respectively.Consider the image formed by P = D(I f ) + D(I b ).D(I b ) causes the appearance of a directional minimum in the neck region as represented by a white line segment.The appearance of a directional minimum is a direct consequence of the neck shape resulting in lower magnitudes of the distance field.Similarly, D(I f ) causes the appearance of a valley region around the nuclei foregroundbackground interface.Hence, P is a suitable speed function in such cases to control the evolution of a level-set to segment individual cells.Now, consider cases where image gradients are significant.
Some nuclei exhibit high gradients at the site of overlap as shown in Figure 2d.As in standard practice, we take the gradient information, g(I) into account.The gradient information g(I) is the sigmoid function of the derivative-of-Gaussian (doG) gradient magnitudes with an appropriate σ .The σ is automatically determined from the size of the cells (m LargestCellRadius).
m SigmaCell -doG filter σ for nuclear stain.Finally, we also use the correlation map (R(p)) computed as output in itkCellForegroundExtraction as yet another source.These three sources are assigned weights to yield the final speed function.
m DistanceMapWeight, m GradientMagnitudeWeight, m GaussCorrWeight -Weights for the three individual speed functions to create a cumulative speed function.Usually, we use a higher weight for the distance map (0.5) and the gradients (0.3) compared to the correlation map (0.2).
In order to speed up computation, the images can also be under-sampled during the intermediate computations and resampled back.The ratios for sub-sampling are set using m Sampling.

Seed Generation
The next logical step in using level-set methods is to generate a suitable initialization (seed points).In our approach, we generate a list of potential seed points that are used to generate level-set functions.The seed points are chosen as the maxima points in the distance map that are located in the nuclei foreground.
Maxima points are located by dilating the distance map and finding those points that do not change in value.
It is easy to see that dilation by a mask does not change the value at only the maxima points.The size of the mask determines the minimum separation distance between any two seed-points.It is set depending on the largest cell diameter supplied by the user (m LargestCellRadius).The filter itk::SeedExtraction is templated over the Feature and Foreground Image Types.It has a single parameter.NOTE: The setting of parameters takes into account the spacing in the images.Hence, if the images are in µm, then the relevant parameters also need to be in the same units.

Level-set Segmentation
We use the active contour formulation with shape priors developed by Leventon et al. [2] for nuclei segmentation.In this formulation, the evolution is given by the following pde: where ψ is the level-set function and the zero level-set corresponds to the segmentation curve/surface.The first terms include propagation, curvature and advection terms.The final term (ψ * − ψ) is a recent addition by Leventon et al. [2] wherein they incorporated principal components of the segmentation shape model to drive the update equation.The surface ψ * is the maximum a-posteriori shape given the current segmentation and image information.The parameters α, β , γ are user-defined settings for the relative scaling of the three speeds.Note that we typically want the level-set to converge and not escape out of the object.Hence, we specify higher values of γ and β relative to α.
Latest version available at the Insight Journal [ http://hdl.handle.net/1926/1]Distributed under Creative Commons Attribution License Since, we have a large number of cells and each cell has a small finite size, it is encapsulated in a ROI.The size of the ROI, once again, depends on the largest cell diameter(m LargestCellRadius).

Shape Model:
Training data consisting of manually segmented nuclei from 3D images are used in estimating a PCA-based shape model with user-chosen number of PCA modes of variation.The advantage of using this model is that a good number of cells within a phenotype share similar sizes, shapes and even orientations.A further discussion of the shape parameters is available in the section on level-sets in the Insight Toolkit [1].
m ShapePriorScaling -Shape scaling parameters γ m SeedValue -Seed value for initializing level-set.Usually 1/3rd the typical nucleus diameter.m Iterations -User chosen number of maximum iterations.Typically, 100 iterations are sufficient.m MaxRMSChange -Maximum permissible RMS change to halt evolution.m NumberOfPCAModes -Number of PCA modes chosen, say 2 or 3 depending on the shape model.

Usage
const unsigned int Dimension = 2; typedef itk::ShapeLevelSetBasedCellSegmentation< FeatureImage,CompImage > FilterType; ... 6 Voronoi Refinement The level-set evolution partitions the foreground into individual nuclei.Their evolution is guided by the speed function.The affiliation of minute fragments of the nuclear foreground remains unresolved.In order to resolve the affiliation, we assign pixels in a given connected component to the nearest cell available in that component.Some of the unresolved pixels may also come from components that were too small to represent cells.These components are eliminated by using a size threshold (m MinComponentSize).

Results
The results in this example can be obtained by using VoronoiBasedCellSegmentation2D.cxx on the level-set image 30-levelset.mhausing the parameter settings given below.The output is written out in the images 30-voronoi.mhawhich is the final segmentation.Potential users of this code can execute the example with the following command line: ./voronoiSegment2D 30-fg.mha30-levelset.mha30-voronoi.mha//Filter settings: filter->SetMinComponentSize( 400 );

Cell Statistics Computation
The filter itk::CellStatistics computes relevant statistical/shape parameters of each individual segmented cell such as centroid, intensity-weighted centroid, standard deviation of intensities, average intensity, principal directions of the shape etc.It also relabels the cells to have consecutive labels.The filter is templated over the Segmentation Image Type and has no parameters.

Results
The results in this example can be obtained by using CellStatistics2D.cxx on the segmentation image 30-voronoi.mhausing the parameter settings given below.The tabbed output is written out on standard output which can be redirected to a text file.The output image (30-segment.mha)consists of relabeled cells after deleting cells below a threshold size.Potential users of this code can execute the example with the following command line: //Filter settings: filter->SetMinComponentSize( 400 ); This contribution also use some of the class in the Review directory of ITK.You need to compile ITK with ITK USE REVIEW ON to be able to compile our code.

Figure 1 :Figure 2 :
Figure 1: (a) Confocal image showing nuclei of the zebrafish ear.(b) Extracted foreground (c) Gaussian correlation map showing convex regions and suppressing non-convex ones.

7. 3 Acknowledgments
These classes use the Insight Journal contribution from Lehmann G. named "Binary Attribute Morphology".A version of this contribution is included in the code distributed along with this paper.Only the CMakeFiles have been modified to fit into our build framework.More recent version of this contribution might be found here: http://hdl.handle.net/1926/584.
It has a total of 6 parameters whose settings are described below.NOTE: The setting of parameters takes into account the spacing in the images.Hence, if the images are in µm, then the relevant parameters also need to be in the same units.
The filter itk::CellForegroundExtraction is templated over the Input, Foreground and Correlation Im-Latest version available at the Insight Journal [ http://hdl.handle.net/1926/1]Distributed under Creative Commons Attribution License age Types.
The filter itk::SeedExtraction is templated over the Feature and Foreground Image Types.It has a single parameter.