3D Cell Nuclei Segmentation with Balanced Graph Partitioning
33D Cell Nuclei Segmentation with BalancedGraph Partitioning
Julian Arz , Peter Sanders , Johannes Stegmaier , and Ralf Mikut Institute for Applied Computer Science, KIT, Karlsruhe, Germany.
Abstract
Cell nuclei segmentation is one of the most important tasks in the analysis of biomedicalimages. With ever-growing sizes and amounts of three-dimensional images to be processed,there is a need for better and faster segmentation methods. Graph-based image segmenta-tion has seen a rise in popularity in recent years, but is seen as very costly with regard tocomputational demand. We propose a new segmentation algorithm which overcomes theselimitations. Our method uses recursive balanced graph partitioning to segment foregroundcomponents of a fast and efficient binarization. We construct a model for the cell nucleito guide the partitioning process. Our algorithm is compared to other state-of-the-art seg-mentation algorithms in an experimental evaluation on two sets of realistically simulatedinputs. Our method is faster, has similar or better quality and an acceptable memoryoverhead.
The ultimate objective of developmental biology is to understand the development of a singlesymmetric cell to a complex three-dimensional organism. A big step towards this goal arerecent advances in light microscopy technology, which have paved the way for the large-scale invivo investigation of biological structure. For an example, the experiments of Tomer et al. [49]on the embryonic development of
Drosophilia result in more than one million images with atotal size of about 10 TB per embryo.Already the amounts and sizes of these images make it difficult for them to be manuallyanalyzed. Another challenge is posed by the imaging quality. The images suffer from thepresence of noise, low contrast and a blurring effect induced by the physical properties of theoptical system. When combining the image stacks to three-dimensional images, a degradationof image quality with increasing axial depth can be observed. This phenomenon is due to thescattering of light in living tissue and shadowing of tightly clustered objects. Another obstaclewhich has to be overcome is the anisotropic resolution. The resolution in axial direction istypically lower than in lateral direction, by a factor of three to five. This renders it difficult toreconstruct a three-dimensional boundary of an object.To overcome these challenges, automated methods have to be developed which process andanalyze these images. The automation also has the advantage that it frees the image analysis i.e. the number of voxels per unit of length a r X i v : . [ c s . C V ] F e b rocess from the researcher’s subjective bias and guarantees the reproducibility of the results.Typically, image processing methods consist of a pipeline of multiple tools: from multiviewimage restoration to object segmentation, annotation and tracking to the modeling of entireprocesses. Each stage relies on the results of the previous ones. It is important that each toolhave results of high quality, despite the aforementioned image deficiencies. Furthermore, it isnatural to demand that the total time to analyze the images is in the order it takes for themto be taken, such that computation time does not evolve into a bottleneck of developmentalbiology research. Therefore, new algorithms have to be efficient, both in quality and in time.One of the most important steps in nuclei image analysis is segmentation . The goal ofsegmentation is to accurately delineate the boundaries of all observed cell nuclei. It translatesthe pixel-based representation of the data to an object-based one [28] and is therefore crucialfor downstream tasks such as cell state annotation, cell tracking or lineaging.Many multi-purpose image segmentation algorithms are based on a graph representation ofthe image [42, 3]. By mapping each pixel to a node in a graph, the image segmentation problemcan be reduced to graph theoretic algorithms. The drawback of the employment of the graph-based techniques is their high memory overhead, which, in the context of cell segmentation,makes them suffer most from the steady increase in the size and the amount of images to beprocessed.In our work, we develop a new graph-based 3D cell nuclei segmentation algorithm whichovercomes the current limitations of these methods. We extract foreground regions from thebackground with a fast and effective binarization. We perform graph-based segmentation suc-cessively on each connected foreground component. This approach solves the problem of thegraph representation induced memory overhead, because the size of these components does notdirectly depend on the image size. The splitting phase uses a probabilistic model based ondomain knowledge to make decisions. We perform an experimental evaluation on two sets ofrealistically simulated inputs. The experiments confirm that our algorithm is both efficient andallows a high quality segmentation even of tightly clustered cells. Confocal Fluorescence Microscopy.
Biological structures at cell level are mostly transparent.Therefore, fluorescent markers, such as the green fluorescent protein (GFP), are introduced intothe structures to be observed. When illuminated with an excitation wavelength, the fluorophoreemits light of a longer wavelength. In conventional fluorescent microscopes, the entire specimenis illuminated. This poses a problem because emitted light from out-of-focus areas reduces thesignal contrast from the in-focus plane [24].Confocal laser-scanning microscopy (CLSM) [34] limits the field of illumination to the regionin focus. A point of light scans the specimen by moving over all planes, lines and points. Thisallows for optical sectioning of an organism, i.e. creating a 3D image. Confocal fluorescencemicroscopy was among the most frequently used techniques in the last decades. However, thistechnique still suffers from the limited accessible depths, because the laser has to penetratethe entire specimen. Another problem is the effect of photobleaching, a process in which themarkers are left unable to fluoresce. These limitations hinder its service for long-time detailedimaging of developing organisms.In contrast, in selective plane illumination microscopy (SPIM) [21], the entire in-focus planeof the sample is illuminated orthogonally to the detection axis. This technique increases thefeasible depth of the sample, as the illumination light does not have to penetrate the sample upto the in-focus plane. It also increases the attainable axial resolution. Furthermore, only the2urrent plane is affected by photobleaching. A further reduction in phototoxicity is achievedwith digital scanned laser light sheet fluorescence microscopy (DSLM) [27]. These advancesmake it possible to image organisms over long periods of time, with high spatiotemporal res-olutions, enabling new studies of cell behavior in complex developing specimens [28]. For anexample, Keller et al. [27] record entire zebrafish and analyze cell nuclei positions and move-ment over the first 24 hours of embryogenesis. The study produces stacks of about 400 imageswith a size of 8 MP and a temporal resolution of up to 60 s. Tomer et al. [49] combine light-sheet microscopy with multi-view imaging. Their experiments on the embryonic developmentof
Drosophilia result in more than one million images with a total size of about 10 TB perembryo.
Due to the difficulties of the cell nuclei segmentation task, most algorithms apply a sequence ofprocessing methods. We explain some of these methods, their goals and their relation to eachother.
A straightforward approach for graph-based image segmentation is to define edges with weightsbased on the similarity between neighboring pixels and then minimize the cut (maximize dis-similarity) between two regions. However, this approach is biased towards unnaturally smallregions, because these have less outgoing edges and hence often a smaller cut. Therefore, asecond important graph-based image segmentation technique is the minimization of the normal-ized cut [42]. This criterion maximizes both the similarity inside a region and the dissimilaritybetween regions.Another popular method is the graph cut framework (see Section 1.2.4). Al-Kofahi et al. [1]apply it both for image binarization and splitting of nuclei clusters. In recent experimentalevaluations, graph cut-based cell nuclei segmentation methods have ranked best regarding thesegmentation quality [45].The graph cut multi-label techniques’ running times grow linearly with the number of labels,i.e. the number of image objects. This drawback makes additional work necessary to renderthem practicable for more than 20 labels [1].
Often a first step in image analysis is a binarization of the image, with the aim to distinguishthe foreground from the background.Image binarization algorithms such as Otsu’s method [36] or adaptive iterative thresholding[27] can work well if the cells are isolated and foreground and background intensity ranges aresufficiently far apart. Combined with morphological filtering, they enable a fast binarization[29].Graph cut based binarization in the context of cell nuclei has been applied several times.Daněk [8] heuristically determines hard constraints for the data term, the boundary termemploys a sophisticated computation of the distance between pixels in a gradient inducedRiemannian manifold. This technique has high computational demands. Al-Kofahi et al. [1]perform minimum error thresholding to gain probability distributions for the data term, anduse a simple gradient-based smoothness term. 3f the cells are isolated, a single image binarization can be used to successfully segment thenuclei. However, for touching or overlapping objects, the task is more complex. The binarizationstep computes an under-segmentation of the image, resulting in clusters of cell nuclei. Thesehave to be split in a cluster separation step.
In terms of image analysis, detection refers to the task of placing a marker, or a seed point, oneach desired object in the image. A basic insight of image analysis is that detection is easierthan segmenting, i.e. delineating the border of objects. Therefore, an actual segmentationis often preceded by a seed detection phase. This especially applies to the field of cell nucleisegmentation, where the desired objects are roundish blobs of similar sizes.The popular multiscale Laplacian of Gaussian (LoG) provides an effective blob detector.Stegmaier et al. [45] find peaks in the maximum intensity projection of a range of predefinedscales, while Al-Kofahi et al. [1] implement a computationally more expensive automatic scaledetection, exploiting cues in the distance transform of the binary foreground mask. Daněk andMatula [9] search for peaks in the distance transform. While this method is fast, it is errorprone for tightly clustered nuclei, especially in 3D. Lou et al. [33] use a blob detector based oneigenvalues of the Hessian matrix.The success rate of later stages depend on the accuracy of the seeds, as each seed will laterform a final object.
When faced with an under-segmentation, object clusters have to be split to distinguish each cellnuclei. The classical approach is the watershed algorithm. However, it often results in over-segmentation, making a post-processing step necessary [29]. The graph cut framework withits multi-label variant has been used as well. Daněk and Matula [9] use the above describedseed points as initial labels and combine the boundary term from the binarization phase witha distance transform induced shape term. In [1], the data term is modeled with a Gaussianmixture for each cell, and the ad hoc smoothness term is reused. The α -expansion techniqueis sped up by reducing the number of labels. Instead of introducing one label per seed, thelabels are found by coloring a graph which is built on the initial labels. In [33] the probabilitydistribution is gained with supervised learning, by training a random forest on a set of localfeatures. This approach requires a set of images with ground truth – a fact which makes thismethod difficult to employ for huge 3D images with thousands of cell nuclei. The authors alsointroduce shape priors: they augment the graph cut formulation with a term which encouragescuts to be aligned to a vector field centered at the initial seeds.He et al. [19] employ the normalized cut algorithm [42] to segment 2D cell nuclei images.The image is represented with a fully connected graph. The edge weights encode gradientinformation, and weights of edges which cross binarization contours are set to zero. The optimalnormalized cut is computed by solving an eigenvector system. Clusters consisting of more thantwo nuclei are split recursively, with a heuristically fixed number of recursions. Depending onthis number, the image might be over- or under-segmented. The complexity of solving theeigenvector system is O ( n ) , where n is the number of image pixels. It is unclear if this methodcan be extended to big 3D images. 4 .1.5 Cluster Merging Instead of performing a top-down approach, where an under-segmentation is split up in a post-processing step, it is also possible to segment the image in a bottom-up fashion. The idea isto use an algorithm which results in an over-segmentation and then merge the components tointact nuclei. Lin et al. [29] apply a 3D watershed algorithm on a gradient-weighted distancetransform, which allows them to incorporate both geometric and intensity cues. In a postpro-cessing step, they recursively merge touching objects. Merging blocks of pixels is conceptuallysimpler than splitting them, because the cut is already defined. The difficulty is to determinewhich segments are to be merged. In [29], a score is computed which signifies the coherenceof the potential merged objects to a nuclei model. The highest score determines which objectsare merged next. However, future potential merging steps are not considered. An improvementLin et al. [30] is to make merge decisions by finding the best solution in merge trees built onthe nodes of a region adjacency graph. This graph encodes all possible sequences of mergingsteps. Liu and Sclaroff [32] combine splitting and merging of objects in one algorithm. Thealgorithm recursively splits regions by straight lines between candidate endpoints.
Best segmentation results cannot be obtained with an all-purpose algorithm. The assumptionthat the information required for a perfect segmentation is provided by the image only is wrong.Hence, segmentation algorithms have to incorporate domain knowledge to be successful. Thisknowledge can be of qualitative or quantitative nature. In the context of nuclei segmentation,for example, the qualitative knowledge that cell nuclei are blob-shaped is exploited by using aLoG seed detector. Quantitative knowledge is introduced by making use of known values forcertain object-specific features. In Stegmaier et al. [45], for example, the multiscale LoG-filterworks on a range of user-defined scales, which reflects the expected diameter range of the nuclei.The knowledge on a feature or a set of features can be expressed in a probabilistic objectmodel . This model can guide a segmentation process by reflecting the probability that an imageobject is a desired object or not. Instead of predefining the model parameters, it is also possibleto learn them from the data. The merge decision process in [29, 30], for example, makes use of amodel built on eight features, including texture, volume and several quantifiable shape features.The parameters of a multi-dimensional Gaussian distribution over these features are estimatedfrom heuristically chosen objects in the image. This distribution enables the computation of aconfidence score for any object. Liu and Sclaroff [32]The split and merge algorithm in [32] is model-driven as well, the parameters for the threefeatures color, shape and area overlap are gained by supervised learning.
This section introduces a set of terms and their relation to each other. For a deeper explanationof these terms, we refer to [43, 16, 48]. A scene is a segment of the real world. It is composedof three-dimensional objects . A real image is a continuous function f : R → R . It mapspoints from the two-dimensional image plane to the intensity . The 2D real image is the resultof a projection of the scene onto the plane [43]. To be processed by a computer, the realimage is digitized . A digital image I is a function from a bounded discrete rectangle P = , , . . . , S x } × { , , . . . , S y } to the discrete brightness levels V = { , , . . . , K } . An element p = ( i, j ) ∈ P is called a pixel . The rectangle P is called the pixel space . For the value of apixel p , we also write I p instead of I ( p ) .The digitization is the conversion of the real image to the digital image. The discretization ofthe plane is done by sampling with a discrete set of sampling points in the plane. Quantization splits the intensity range into K intervals. Usually, K is of the form b , where b is the bit depth .The histogram h : V → N of a digital image is the absolute frequency of discrete brightnesslevels in an image.Digital images formed by the described process are two-dimensional, because the imagesensor used in the digitization is two-dimensional. Images of higher dimensions are formedby grouping a set of images to an image stack : a three dimensional image is a stack of two-dimensional images. Multiple 3D images which represent different points in time can be stackedto a 4D (or 3D+t) image. This work, only deals with 3D images. In this case, the pixel spaceis a discrete cuboid: P = { , , . . . , S x − } × { , , . . . , S y − } × { , , . . . , S z − } . A 3D pixel might also be called voxel . The size S of an three-dimensional image is the vec-tor ( S x , S y , S z ) . In the following, when we speak about images, we will refer to 3D digitalimages.There are two notions of distances in digital images: physical distance or pixel spacing , and pixel distance . To be able to measure the physical distance between two real image featuresusing the digital image, the pixel spacing has the be known. It is defined as the shortestdistance between two real-world points which can be distinguished with the combination ofthe optical system and the employed camera, along each of the three axes. Given the imagespacing δ ∈ R , δ = ( δ x , δ y , δ z ) , the image resolution denotes the pixel per unit distance r =( r x , r y , r z ) = (1 /δ x , /δ y , /δ z ) . Then the physical distance d ( p, q ) between two pixels p, q ∈ P is computed by interpreting P as the euclidean three-dimensional space and determining theeuclidean distance between p and q .The pixel distance is a dimensionless measure for the amount of basic steps required to movebetween discrete pixels. We define two variants. The distance D is the 3D extension of thecity-block distance and is defined as D (( x , y , z ) , ( x , y , z )) = | x − x | + | y − y | + | z − z | .D is known as the checker-board distance: D (( x , y , z ) , ( x , y , z )) = max( | x − x | , | y − y | , | z − z | ) . Two pixels p , p are k -adjacent , k ∈ { , } , if D k ( p, q ) = 1 . Then, p is a k -neighbor of q . The k -neighborhood of a pixel is the set of pixels which are k -adjacent to it. When clear from thecontext, or not necessary, we omit the specifier in the following. A neighborhood system N isthe set of adjacent pairs of pixels.A path from a pixel p to a pixel q is a sequence of pixels A , A , . . . , A n where p = A , q = A n and A i +1 is a neighbor of A i for ≤ i < n . Two pixels are connected if there exists a pathbetween them. A connected component or simply component is a set of pixel which are pairwiseconnected. The border of a component C is the set of pixels of the component with one ormore neighbors outside of C . 6 segmentation family of an image is a finite set of subregions R , . . . , R l , where • R i is a connected component, • P = (cid:83) li =1 R i , and • R i ∩ R j = ∅ for i (cid:54) = j .It is possible to relax the second demand by introducing a special not necessarily connectedsubregion R b for the background defined as R b := R \ (cid:83) li =1 R i .For a given segmentation family, a segmentation is a function which maps pixels to thesubregions they are contained in. The subregions R i are then also called labels . A special caseof a segmentation is a binarization , where the number of subregions is l = 2 , and the regionsare foreground and background, denoted as R f and R b , respectively.Informally, each component of a segmentation corresponds to a real world object in the scenerepresented by the digital image. Hence, each pixel is mapped to either a real world object,or background. Another notion of segmentation is soft segmentation , which takes into accountthat a particular pixel can incorporate signals of multiple objects. This is accomplished bymaking use of fuzzy set theory [50]. Each pixel is mapped to a fuzzy set over R , . . . , R s .This technique should not be confused with the uncertainty propagation discussed in the nextsection, where fuzzy sets are used to measure the uncertainty of the entire (hard) segmentation.We only consider hard segmentation in our work. Most operators in an image processing pipeline are unable to give precise results for all inputs.On the one hand, this is due to the inherent deficiency of the imaging quality. The low signal-to-noise ratio might it make hard to distinguish small nuclei from background, or aberrationsmight be wrongly detected as nuclei. On the other hand, many tasks in image processing areill-posed problems [28, 8]: the solution to such a problem are subject to ambiguities. Theremight not be exactly one correct answer regarding a certain image feature, e.g. a cell which isin the mitosis during the imaging process can be detected as a single cell or two cells, dependingon how far the mitosis is already visible. In order to prevent errors made by one operator fromaffecting downstream tasks, but also to avert the loss of improbable but viable solutions, astage should be able to inform its successors of the uncertainty involved in its output. This ismade possible by introducing a measure for the validity for every extracted piece of informationand propagating it together with the results to the next pipeline stage [44].We follow Stegmaier et al. [44] in making use of fuzzy set theory to measure uncertainty,based on a priori knowledge. They introduce a fuzzy set membership function µ ijm : R → [0 , for every operator i , feature m and linguistic term j . One way to define such a function isby using the following set of four terms: “Feature is. . . { . . . as expected”;. . . too small, but use-ful”;. . . too large, but useful”;. . . not useful” } . This set lets us model the uncertainty with atrapezoidal membership function µ m ( x, θ m ) defined by a quadruple θ m = ( a, b, c, d ) as Most definitions[16] further demand that Q ( R i ) = TRUE and Q ( R i ∪ R j ) = FALSE for a logical, domain-specific predicate Q . While this is certainly true for all segmentations, it is rather vague not necessary for adefinition. m ( x, θ m ) = x − ab − a if a ≤ x < b ; if b ≤ x < c ; d − xd − c if c ≤ x < d ; otherwise.Following the rules of fuzzy set theory, multiple uncorrelated features can be combined bymultiplying their respective membership functions. Given an undirected graph G = ( V, E ) with non-negative edge weights c : E −→ R + , and anumber k ∈ N , k > , the graph partitioning problem is to find k sets of nodes V , . . . , V k ,with (cid:83) ki =1 V i = V and V i ∩ V j = ∅ ∀ i (cid:54) = j . Given an additional imbalance factor ε ∈ R + ,the balanced graph partitioning problem further demands that the sets are balanced, i.e. that | V i | ≤ (1 + ε ) (cid:100)| V | /k (cid:101) . The goal is to minimize an objective function. In this work, we will aimto minimize the sum of the weights of the set of cut edges {{ u, v } ∈ E | u ∈ V i , v ∈ V j , i (cid:54) = j } .The graph partitioning problem has been shown to be NP-complete [22, 14]. Practical graphpartitioning tools rely on heuristics to find good partitions. Most tools use the multilevelapproach [20], where the graph is first iteratively coarsened, i.e. a hierarchy of graphs isconstructed which keeps the original structure but reduces the input size. A common methodfor coarsening is contracting edges in a matching in G [20]. Then a partitioning is computed onthe coarsened graph. One way to do this is to search for k seed nodes which are far away fromeach other and then alternately run k BFS, one from each seed [12]. Finally the graph withits partitioning is unpacked, while employing local search methods, such as max-flow min-cutbased search [40] and FM [13]. As we cannot cover all methods here, we refer the reader to[2] for a deeper insight. From the many available tools and libraries for graph partitioning, wetested METIS [26] and KaHIP [41], and achieved best results with a fine-tuned configurationfor KaHIP.compute exact solutions to the two-label segmentation problem, i.e. image binarization.
The graph cut framework is a powerful graph-based segmentation method. It has been shown[17, 6, 3] that it is possible to efficiently minimize an energy function of the form E ( A ) for abinarization A defined as E ( A ) = λE data ( A ) + E smooth ( A ) . A minimum of this function corresponds to an exact solution to the two-label segmentationproblem, i.e. image binarization.The data term is typically E data ( A ) = (cid:88) p ∈P D p ( A p ) , where D p measures the pixel-wise penalty of mapping pixel p to label A p , i.e. fore- or back-ground. The smoothness term E smooth ( A ) penalizes neighboring pixels with different labels.Using a neighborhood system N , it is defined as E smooth ( A ) = (cid:88) ( p,q ) ∈N B p,q · A p (cid:54) = A q . A p (cid:54) = A q is an indicator function which is iff A p (cid:54) = A q . The function B p,q is the boundaryenergy. The energy function E ( A ) is minimized by computing a minimum cut on a specialgraph.Boykov et al. [6] extended the method to multi-label energy minimization, which is NP-hardin general. They find approximate solutions by introducing two approximation techniques, α -expansion, and α/β -swap. The former performs iterative binary segmentations between onelabel and every other label. The latter computes the minimum cut between all pairs of labels.The graph cut segmentation method is only a framework. To employ it in practice, the terms D p and B p,q , and the scalar λ have to be parameterized. Boykov and Funka-Lea [3] suggest todetermine models for the pixel intensity distributions of the objects to be segmented. The dataterm can then be parameterized using the fit of the pixel value to the models. Hard constraintscan be incorporated as well.Due to the extensive use of the graph cut framework for segmentation problems, maximumflow algorithms which exploit the special structure of the graph have been developed [5, 15]and parallelized [10, 31, 25]. Due to the discrete nature of digital images, computing the surface area of a real-world objectusing a digital 3D image is a nontrivial problem. Interpreting each voxel as a cuboid andsumming up the area of their surfaces is only a bad approximation, since with that method,every object has the same surface area as its bounding box. In the context of adopting graphcut to compute geodesics in Riemannian spaces, Boykov and Kolmogorov [4] have introducedthe notion of cut metrics . The authors have shown that it is possible to compute edge weightsfor a regular grid graph such that the length of a contour is approximated by the sum of theedge weights cut by that contour. Their method is based on the Cauchy-Crofton formula. Thisformula relates the Euclidean length |C| of a curve C to the function n C ( l ) which measures thenumber of intersections between the curve C and a line l ∈ L : |C| = 12 (cid:90) L n C ( l ) d l = 12 (cid:90) π (cid:90) + ∞−∞ n C ( φ, ρ ) d φ d ρ, (1)using that the space of lines can be represented as [0 , π ] × [ −∞ , + ∞ ] . Boykov and Kol-mogorov [4] approximate the integral with partial sums by partitioning the space of lines amongthe edges of the grid graph’s neighborhood system. This method yields edge weights ω k = δφ k δρ k , (2)where δφ k is the angular partition of edge k and δρ k is the distance of the line induced by edge k of the neighborhood system to the nearest parallel edge induced line.The method can be formulated for 3D grids as well. However, for that case Boykov andKolmogorov [4] do not specify how to partition the space of angular orientations among theedges of the neighborhood system. Daněk and Matula [9] proposed a solution both in 2D and3D, as well as for grids with anisotropic node spacing. They compute the Voronoi diagram ofthe intersections of the neighborhood system with a unit hypersphere. The angular partitioningof an edge is equal to the fraction of its Voronoi cell. or the length of a contour in 2D Our Method
This section deals with our segmentation method. Our method performs a top-down approach.We first perform a fast and robust foreground detection, detailed in Section 2.1. The connectedcomponents in the binarized image are then recursively split. We make use of balanced graphpartitioning to cope with the bias towards small partitions when minimizing a cut (cf. Sec-tion 2.2). Our segmentation framework computes a score for each component which is based onprior knowledge. The nucleus model is presented in Section 2.3 and the segmentation frame-work Section 2.4. The advantage of our cluster segmentation approach over the approachesdescribed in Section 1.1 is that our method does not depend on a heuristic seed detection,which improves the detection and segmentation quality for highly clustered cell nuclei.
The first step of our segmentation algorithm is to binarize the image. We aim to classifythe image pixels into two groups, the foreground and the background. This step immediatelyreduces the input, as all background pixels do not have to be further examined. For our purpose,it is important that foreground objects not be wrongly classified as background. An error of thiskind leads to undetected nuclei because our subsequent steps are not capable of detecting theseobjects. Another binarization mistake are non-foreground areas which are wrongly detected asforeground. This case might arise due to defective camera sensors, e.g. isolated white pixels.To some extend, our splitting framework is able to detect and discard these erroneous regionsespecially in the aforementioned case.A third binarization error is an inexact delineation of the foreground object’s border. Theobject is correctly detected as foreground, but its shape or size differs. This error commonlyarises in binarized images of confocal fluorescence microscopy. On the one hand, unstained partsof a nucleus may appear as dark as the background [47]. On the other hand, an enlargement ofthe nuclei can be observed for binarized cell nuclei images. The diameter of binarized nuclei isseveral pixels larger than the real extent of the nuclei. The reason for the inflated objects is theoptical diffraction, which can be expressed by the point spread function (PSF) of the imagingsystem. The PSF represents the response of a pixel sensor to an ideal point light source [48]. Apoint source appears blurred in the image and might be surrounded by concentric rings. Thearea directly adjacent to the foreground objects is then still brighter than the dark background.This effect can add up if several objects are close to each other. In these cases, clustering-basedthresholding approaches which do not use domain knowledge result in inflated objects.This phenomenon effects our nuclei segmentation method both in quality and in runningtime. Firstly, the quality suffers. Our method relies on splitting the foreground regions of theimage. It only adds cuts inside the foreground, but does not change the delineation betweenfore- and background. Furthermore, the inflation of the nuclei can also cause well-separatedobjects to be connected in the binarization. These connecting bridges will then be part of onethe segmented objects. These circumstances would have a negative impact on the quality of thesegmentation results. Secondly, the segmentation speed would be slowed down. Our methoduses graph partitioning to split the foreground objects. The larger number of merged objectsrequires more partitioning steps, and their inflated size results in overall bigger graphs to bepartitioned.The PSF depends on the imaging system. It is a combination of the blur induced by theoptical system and the finite integration area of a chip sensor [48]. Hence, the extend of the10oreground blur can differ across distinct datasets. In one dataset, the blur might be strong,and measures have to be taken to detect and cope with it. In another set, simple thresholdingalready leads to good results. In this work, we experimented with three different binarizationmethods. The first is Otsu’s method [36], a simple thresholding algorithm. The other two arebased on work by Restif [37], who models the gray level histogram of a blurred image. Themodel is described in the next section. We either use this model to compute a threshold, orparameterize the graph cut framework with its probability distributions.For all algorithms, binarization quality could be improved by convolving the image with a low-pass Gaussian filter with small standard deviation σ s . This method reduces background noiseand eliminates one-pixel errors. Another image degradation specific to confocal microscopy isuneven illumination, a brightness gradient in the image, often along the axial direction [35].We overcome this deficiency by dividing the image into sets of consecutive slices of size m andbinarizing each set separately. We compute the normalized gray level histogram h ( i ) with the goal to infer two distributions,one for the foreground and one for the background. Al-Kofahi et al. [1] assumed both to bePoisson distributed, citing the image formation process and their findings of the histogramsto be bimodal. However, for some of our data we observed unimodal histograms, due to thesmall fraction of foreground compared to the entire image, the high intensity variance in thefluorescing nuclei, and the system-specific point spread function, which blurs the image andresults in actual background regions to be brighter around the foreground objects. We thereforeadopt a model inspired by [38]. The background is subdivided into non-illuminated (NB) andilluminated background (IB). The former is a normal distribution with parameters µ b and σ b ,and a priori probability p b : N B ( i ) = p b √ πσ b exp (cid:18) − ( i − µ b ) σ b (cid:19) . (3)The IB accounts for the fluorescence microscopy-inherent blurring effect: an increased graylevel intensity of the background around the fluorescing nuclei. This blur is modeled with anexponential decline, from which the model for the IB is derived (see [38]) as IB ( i ) = 2 αAi − µ b log ( I f − µ b )( i − µ b ) , (4)where I f is the intensity at the border of the object and A is the object’s area. The factor α defines the blur’s rate of decline. Restif [38] introduces a foreground normal distribution andan IB term for each nuclei, which is practical as the images used in his work on average onlycontain about seven foreground objects. For our purpose, we model the nuclei with one normaldistribution F ( i ) with parameters µ f , σ f , and p f , and the entire illuminated background withone term IB ( i ) . We set I f to µ f , and A to p f , the foreground’s a priori probability. Becausethe IB encroaches into the NIB, we only compute the IB ( i ) in the range [ µ b + 2 σ b , µ f ) . Then,the histogram model is h model ( i ) = N B ( i ) + IB ( i ) + F ( i ) (5) = p b √ πσ b exp − ( i − µ b ) σ b + 2 αp f i − µ b log ( µ f − µ b )( i − µ b ) + p f √ πσ f exp − ( i − µ f ) σ f (6)11ith the set of missing parameters { p b , µ b , σ b , p f , µ f , σ f , α } .Like Restif [38], we determine the parameters with the expectation–maximization (EM) algo-rithm [11]. The EM algorithm iterates between an expectation step (E-step) and a maximizationstep (M-step). The former computes weights which define the proportion of each summand inEquation 5, given an estimate of the parameters. The latter estimates the parameters, given aset of weights.We introduce weights ω nbi , ω ibi , ω fi for the non-illuminated background, illuminated back-ground, and foreground, respectively. In the E-step, we use Equation 5 to compute the weightsfor each histogram bin i , e.g. ω bi = N B ( i ) /h model ( i ) (7)In the M-step, the parameters for the normal distributions are estimated based on theweighted histogram: p b = (cid:88) i ω bi h ( i ) , (8) µ b = (cid:88) i iω bi h ( i ) , and (9) σ b = (cid:115)(cid:88) i ( i − µ b ) ω bi h ( i ) , (10)for the background distribution and equivalently for the foreground distribution. The parameter α is subsequently determined, analogous to [38], as α = 1log ε (cid:118)(cid:117)(cid:117)(cid:116) p f µ f − (cid:88) i = µ b +2 σ b ω ibi h ( i ) , (11)where ε = 2 σ b / ( µ f − µ b ) is a correction factor for the truncation of the illuminated backgroundmodel.Alternating E-step and M-steps, the algorithm converges after a small number of iterations(about 5 in our settings). We found the outcome of the EM algorithm to be robust and largelyindependent on the initialization. Therefore, we initialize the parameters with a simple iterativethresholding algorithm [39], where a threshold is defined as the histogram mean, and a newthreshold is iteratively computed as the average between the means of the histogram belowand above the old threshold. We use the final threshold to compute the means and a prioriprobabilities of the foreground and background distribution, and set the variances to 1 and α to 0.01.For subsequent computations which require probability distributions for foreground and back-ground pixels, we will combine the illuminated and non-illuminated background into one term B ( i ) : B ( i ) = N B ( i ) + B ( i ) . (12)Hence, B ( i ) ( F ( i ) ) can be interpreted as the probability that a pixel in a background (fore-ground) area has gray level i .We implemented two binarization methods using the above distributions. The first is tocompute an intensity threshold t (cid:63) as t (cid:63) = min i ∈ V F ( i ) < B ( i ) . (13)12ll pixels with values above this threshold are determined as foreground, all others are back-ground.The second method uses these distributions to parameterize the graph cut framework. Withthe notations in Section 1.2.4, we set the pixel-wise data term to D p ( A p ) = (cid:40) − log F ( I p ) if A p = R f ; − log B ( I p ) if A p = R b . (14)For the smoothness term, we use B p,q := exp (cid:18) − ( I p − I q ) σ (cid:19) d ( p, q ) . (15)This function was suggested in [3] and also adopted in [1]. However, it introduces a parameter, σ bin . It results in a high penalty for discontinuities when the intensity difference is below σ bin .A the function resembles a Gaussian probability distribution, the parameter can be interpretedas the standard deviation of the image noise. This parameter, and the scalar parameter λ , haveto be set by hand. The binarized foreground in cell nuclei images decomposes into a large number of small con-nected components. These can either form single nuclei, clusters of nuclei or falsely detectedregions due to imaging defects. We perform graph-based segmentation successively on eachcomponent.We use balanced graph partitioning to split the binarized foreground. The graph partitioningaims to minimize the cut between partitions, i.e. the sum of weights of the cut edges. Thisallows the method to take advantage of intensity or gradient information of the original image,but also of shape cues of the binary foreground mask.In this perspective, the approach is related to the multi-label cut methods from the graphcut framework, namely α -expansion, and α/β -swap. However, these methods rely on a properparameterization to show good performance. The data term of the labeling energy graph cutframework reflects the attractions of a voxel to all labels. These have to be determined in ad-vance. With graph partitioning, the number and, more importantly, the positions of the labelsdo not have to be known. Another difference is that the graph cut methods perform multipleiterations of optimization steps. The α -expansion, for example, iterates over all partitions andcomputes a minimum cut to the rest of the graph. It repeats this process multiple times untilconvergence. In our method, each graph partitioner call is only done once, and its outcome iseither kept or discarded.The balance constraint serves multiple purposes. Most importantly, it counteracts the biastowards unnaturally small components. Without this measure, the process would almost alwayscut only a small set of pixels from the main component, as these are connected by a smallnumber of edges. Figure 1a gives an example for this problem.The balance constraint forces the graph partitioner to select edges in the green area of Fig-ure 1b and thus finds the cut which correctly segments the image. This image also illustratesanother reason to employ balanced graph partitioning for cell nuclei segmentation. The objectsare often roughly equally sized and separated only by tight bottlenecks. Therefore, the combi-nation of minimizing a cut and balancing the two resulting partitions makes adequate use of13 a) Minimum cut bipartition of the component. (b) Minimum cut balanced bipartition of the compo-nent. To obey the maximum imbalance factor,cuts have to lie in the green area. Figure 1: Purpose of the balance factor for graph partitioning.the shape cues of the binary foreground mask: cuts are smaller when they contain less edges,which favors bottleneck regions to be cut.The balancing also improves the speed of our algorithm for big components which have tobe split multiple times. While the number of required partitioner calls remains asymptoticallyunchanged, the subcomponents sizes decrease faster, and the recursion depth is decreased.For each connected component in the binarized image we construct a graph with nodesrepresenting voxels and edges connecting neighboring voxels. A 6-neighborhood was sufficientin our experiments. As the number k the foreground components have to be split into is notknown a priori, we set k = 2 and compute recursive bipartitions. The model-based stopping-criterion will be explained in the next section. We describe the edge weights for the isotropic case. We implemented two methods of settingthe edge weights. One ( grad ) sets weights based on gradient, the other ( prob ) on the probabilityof a pixel intensity belonging to background.For the gradient-based edge weights, we reuse the method for the graph cut binarization (seeSection 2.1). The intuition behind this choice is that neighboring pixels with high intensitydifference are less likely to belong to same label. We compute edge weights c grad ( { p, q } ) betweentwo pixels p, q with intensities I p and I q as c grad ( { p, q } ) = exp (cid:32) − ( I u − I v ) σ (cid:33) . The parameter σ grad is chosen smaller than the σ bin in the binarization.However, in our case all pixels are already detected as foreground. We experienced problemswith this method for images with high variations of the fluorescent marker. The resultingirregular inner-cell texture of the nuclei can lead to wrong cuts if the inner-nuclei gradient istoo high.In contrast to 2D images, tightly clustered cell nuclei in 3D images do not overlap. In ourinstances, we observed that even between very close or touching nuclei there was still a tightband of darker pixels. Therefore, a second method is to induce the edge weights directly fromthe gray levels. The cut should tend to be near the darker area between the bright nuclei. In14rder to be able to quantify the notion of “dark areas”, we reuse the histogram image modeldescribed in Section 2.3.We maximize the probability that the nodes a found cut S is adjacent to represent backgroundpixels. To this end, we disregard the non-independence of neighboring pixel values, and aimto maximize the product of the background probabilities of the adjacent pixels. Formally, wemaximize (cid:89) { u,v }∈ S min { P( B | I ( u )) , P( B | I ( v )) } , where we use the conditional probability P( B | I ( u )) to denote the probability that pixel u with intensity I ( u ) is in the background. Maximizing the above expression is equivalent tominimizing (cid:88) { u,v }∈ S − log min { P( B | I ( u )) , P( B | I ( v )) } . From the last formula, we derive the edge weights c prob ( { u, v } ) as c prob ( { u, v } ) = − log min { P( B | I ( u )) , P( B | I ( v )) } . (16)The required probabilities that a given intensity represents background can be computed with P( B | I ( u )) = p b P( I ( u ) | B )) p b P( I ( u ) | B ) + p f P( I ( u ) | F ) . (17)for our model described in the previous section, this is equal to P( B | I ( u )) = B ( i ) h model ( i ) . (18)For experimental purposes, we also tried setting all edge weights to the same constant. Thenthe cut is independent of the intensity levels and only depends on the shape of the foregroundobjects.Once the graph is constructed, we partition it into two blocks with a balanced graph parti-tioning solver. In general, these blocks do not have to be connected. This circumstance occursin particular if a component consists of more than two nuclei of differing sizes, as depicted inFigure 2a. Neither the red nor the green cut are eligible. One of the smaller cell nuclei wouldform one block, and the other block would consist of the big nuclei paired with the other smallnuclei, which is not feasible if the imbalance factor is too tight. Instead, the balanced graphpartitioning results in the green cut shown in Figure 2b. It partitions the graph into two blocks,one (lilac) of which is not connected. To cope with this phenomenon, we perform a breadth-firstsearch after the partitioner call to identify the connected components and their sizes. In thecase shown in the example, the component is partitioned into three subcomponents.In order for the cut in Figure 2b still being valid, even for equally sized nuclei, the balanceconstraint has to be relaxed. Hence, we set ε = 0 . to cope with this case. In most confocal fluorescence microscopy images, the resolution along the lateral axes ( r x and r y ) is smaller that the axial resolution ( r z ). When setting the edge weights, we take theanisotropic resolution into account by dividing the weight of an edge between two pixels u and v by their distance d ( u, v ) . 15 a) Neither the red nor the green cut obey the imbal-ance factor for a bipartitioning. (b) The minimum cut balanced bipartitioning. Oneof the two blocks (lilac) is not connected. A BFSresults in three subcomponents. Figure 2: Unconnected blocks of the graph partitioning.This method is based on the observation that the cut of an object should be invariant of itsorientation in the image. The number of edges cut by a 2D plane in the 3D image dependson its orientation. For an example, a cut of the form of an xy -plane of area will cut (i.e. isperpendicular to) r x r y edges, while a xz -plane of the same size will cut r x r z edges. By dividingthe edges by their length the sum of these cut edges will be the equal. We construct a mathematical model of a cell nucleus. The goal is twofold: on the one hand,the model should decide whether an image component represents a single object such that therecursion can be stopped, a cluster of multiple objects to be split, or background noise to bediscarded. On the other hand, we want to compute a score s ( C i ) for each component C i whichsignifies the probability that it represents a single nucleus.A splitting step results in a number of subcomponents, each of which is smaller in size than theoriginal component. If a connected component is too big to be a nuclei, we can safely performa splitting operation. If it is too small, we can be certain that it is not a nuclei. Therefore, ina first approach, we base the scoring function only on the volume V of the components. Thevolume V i of an object C i is approximated by multiplying its voxel number n with the spacing: V i = nδ x δ y δ z . (19)The voxel number n is counted during the breadth-first search.From user-defined minimal and maximal cell volumes ( V min and V max , respectively) a trape-zoidal fuzzy set membership function µ V indicates whether a component has a plausible volume.We define µ V with the parameter vector θ V = ( V min , (1 + λ ) V min , (1 − λ ) V max , V max ) . The pa-rameter λ gives us a fine control over the shape of the lateral sides of the trapezoid. In ourexperiments, it is set to . .A component with volume below V min is discarded. A component is repartitioned if it is bigenough to contain a minimum sized nucleus. Due to the imbalance factor ε of balanced graphpartitioning, the reduction in size of a subcomponent compared to its original component is atleast (1 + ε ) / . Therefore, if V i is greater than ε V min , the component C i could be constitutedof one minimum sized nucleus C c . The remaining part, C i \ C c , does not have to represent anucleus. 16he volume can be used as a single criterion if the objects’ sizes have a small variance.However, we observed that the volume range of cell nuclei allowed both a parent componentand its subcomponents to be inside the user-given size limits.To improve the decision making process and give better estimates for the uncertainty values,we use the knowledge that our desired objects are cell nuclei, which can be seen as roundish,blob-like objects. To distinguish between an ellipsoid-like nucleus and one cut into two nearlyequally sized partitions (i.e. similar to an ellipsoidal dome), we compute the sphericity Ψ c defined [51] as the ratio between the surface of a sphere with volume V i and the surface A i ofthe object C i . Ψ i = π (6 V i ) A i . (20)where A i is the surface area of a C i .In theory, values for the sphericity Ψ are in a range between and , where is the sphericityof a perfect sphere. Depending on organism and cell type, different values of sphericity can beexpected for cell nuclei in practice.We assume an ideal nucleus to be an ellipsoid with near-equal semi-axes and deduced asphericity of Ψ ideal = 0 . and Ψ min = 0 . . The sphericity membership function is thendefined as µ Ψ ( x ) = if x ≤ Ψ min ; (cid:16) x − Ψ min Ψ ideal − Ψ min (cid:17) if Ψ min < x < Ψ ideal ; if Ψ ideal ≤ x ; (21)We compute the surface area A with the cut metric results described in Section ?? . The set ofcut metric edge weights in a regular grid are the same for each node. We use a 26-neighborhoodfor the cut metric evaluation because the approximation gained with a 6-neighborhood is toinaccurate.We multiplicatively combine the two membership function values to gain the score s ( C i ) ofa component: s ( C i ) := µ V ( V i ) µ Ψ (Ψ i ) . We recursively split the connected components with balanced graph bipartitioning, computingthe above score at each step. A way to see this splitting process is a depth-first search througha tree, a set of potential desired objects. Each inner node in the tree represents a component.The initial component forms the root of the tree, the children of a node are the outcomes ofa partitioning step. The leaves are the desired objects. Listing 1 gives a top-level overview ofthis depth-first search framework.The framework’s core is a scoring function (
ScoreFunction ), which determinates the nextprocessing step for a given component. A component can either be kept, discarded, or reparti-tioned. Keeping a component means that the search is stopped at that point and the componentis not partitioned. In the DFS analogy the component is marked as a leaf. A component whichis repartitioned might, but does not have to, consist of multiple desired objects. For a discardedcomponent the search is stopped as well, but the component is deleted from the search tree.The scoring function also computes a score for each component.17 : function RecursiveSplit ( C ) returns ( k , R ) // k : number of components R ←
PartitionerCall ( C , k ) // R : partitioning of C for each C i ∈ R do (decision , s ( C i )) ← ScoreFunction ( C i ) if decision = repartition then ( k i , R i ) ← RecursiveSplit ( C i ) // recursive call if k i = 0 then // backtrack if s ( C i ) > then decision ← keep else decision ← discard else merge R i into C i of R if decision = keep then k i ← else if decision = discard then remove C i from R k i ← k ← (cid:80) i k i Algorithm 1: The cell partitioning framework.The score s ( C i ) measured for a component C i is interpreted as the probability for an event X i , that the component C i is a desired object, under the condition that its parent component C p is not a desired object: s ( C i ) = P( X i | ¯ X p ) . (22)A component which is kept must have a score value greater than . If both the score s ( C p ) of a component C p and the score s ( C c ) of one of its children is strictly positive, we keep thecomponent which has a higher probability to be a nuclei. Using the fact that for the event X c to occur, the event X p cannot occur, we deduce: P( X c ) > P( X p ) ⇔ P( X c ∩ ¯ X p ) > P( X p ) ⇔ P( X c | ¯ X p ) · P( ¯ X p ) > P( X p ) ⇔ s ( C c ) · (1 − s ( C p )) > s ( C p ) (23)Thus, for each child component of a component C p , we check whether it is discarded or notwith Equation 23. As an implication, components with a score greater than . are not split.Figure 3: An exemplary execution of our pipeline. The input image (A) is first binarized (B),then recursively split (C-E). The split nuclei in the last image (E) is detected by ourcell nucleus model. Image D is the final segmentation.18e examine each child separately such that discarded components have no influence on theirsiblings. That way, our algorithm is allowed to discard areas which do not contain any nuclei,reflecting the possible existence of wrongly binarized areas. In the backtracking step of ourdepth-first search, we check if any of the child components were kept. Otherwise, the current(parent) component is discarded or kept if its score is greater or equal to , respectively. Figure 3gives an example. The purple nuclei in image D has a low score of . due to its upper-rangedsize. It is split into two components (E), with scores . and . . After comparing both tothe original score, they are discarded during backtracking. Image D is the final segmentation. We implemented our algorithm in C ++ using the Insight Toolkit [23] for image operations.The test machine is an Intel Core i7-920, equipped with four cores clocked at 2.67 GHz. Themachine has 12 GiB of RAM. All experiments were performed on a single core. The graphs arepartitioned with the open source graph partitioner KaHIP [41] with a customized configuration.We compare our algorithm to TWANG [45] and the graph cut (GC) segmentation described in[1]. We also include our binarization method (Otsu) in the evaluation.We measured the quality of a segmentation by counting the number of added, split, mergedand missing objects. Added and split are false positives, areas detected as nuclei where thereis none. Merged and missing are false negatives. We divide by the number of labeled nucleiin the respective image. The percentages are then averaged over all images in the set. Thetimes measured are averaged as well and do not include image I/O. None of the implementationswere explicitly optimized for memory consumption. Nevertheless, we display values for memoryconsumption by extracting the peak resident set size of the respective process. Ground truth of real world 3D microscopy images which contain thousands of cell nuclei doesnot exist. Manual segmentation of these images is not an option due to the large amount ofobjects on the one hand and the researcher’s bias on the other hand. We therefore benchmarkthe quality of our algorithm on two sets of synthetic cell nuclei images. We downloaded a setof 30 HL60 cell line images generated with the CytoPacq toolbox [47], with high noise andclustering probability of 75% (High75). They contain 20 nuclei each and have an isotropicresolution. We scaled the images down by a factor of two in each dimension to let nucleisize resemble those of real world images. Each downscaled image consists of 64 slices with aresolution of × . Based on this toolbox, [46] have conceived a more realistic set of 3D+tbenchmark nuclei images. We use a set of 10 images of different time points (SBDE) whichcontain between 316 and 1016 nuclei. They have a size of × × and the physicaldistance between two slices is five times greater than the resolution of each slice. On the synthetic datasets, the best binarization results were achieved with applying Otsu’sthresholding method after convolving the input image with a low-pass Gaussian filter withparameter σ s . This procedure showed robust results on our image sets, even on those without19trictly bimodal intensity histograms. Computing a sophisticated image histogram model (Sec-tion 2.1.1) and using the parameters to threshold the image had worse results. Parameterizingthe graph cut framework exhibited a bad trade-off between binarization quality and speed.For High75, the partitioner edge weights using on the intensity probability ( c prob ) had slightlybetter segmentation result than the gradient-based edge weights ( c grad ). For SBDE however,the situation was inversed, with the edge weights computed by c prob failing to capture thetrue object boundaries. The reasons are that the inner-nuclei variance is very low for SBDE,while the nuclei are very tightly clustered, such that there is no “background area” between theobjects. We therefore use c grad in our comparison.Table 1: Parameters for High75 and SBDE, respectivelyParameter High75 SBDE σ s m σ grad
15 100 ( V min , V max ) (20000 , , Table 1 contains values for the parameters we used in our comparisons. The standard devi-ation σ s is used for the Gaussian filter preceding the binarization and depends on the size ofthe objects to detect. The parameter m is the number of regions the image is split into for thebinarization and reflects the change of the intensity distribution along the visual axis. The edgeweights are parameterized with the parameter σ grad , which separates the inner-nucleus intensityvariance from the inter-nucleus intensity variance. The minimum and maximum volumes V min and V max are the limits used by our model. The difference in volume between the two imagesets stems from the dimensionless nature of the synthetic images. Table 2: Segmentation results for High75.Algo [s] [MiB] Added Missed Merged SplitOtsu 0.4 0 0 71.3 0GC 244 174 0 0 0.5 1.2Ours 8 432 0 0 0 0TWANG 14 251 0 0 0 0Table 2 shows results for the High75 dataset. Both TWANG and our algorithm have nosegmentation error on this dataset. GC has a low ratio of merged nuclei and split nuclei. Ouralgorithm is the fastest. With the unscaled images, quality is unchanged. GC takes threeand a half hours per image, TWANG about 4 minutes and our algorithm less than 90 s, whichshows that our algorithm scales best with increasing resolution. This result is due to the fastbinarization, which can quickly (415 ms) reduce the amount of voxels to look at by about 91 %.Table 3 lists results on the SBDE dataset. None of the tested algorithms exhibits perfectsegmentation quality. The seed-based methods have very low rate of false positives, however20able 3: Segmentation results for SBDEAlgo [s] [GiB] Added Missed Merged SplitOtsu 4 0 0 56.2 0GC 147 1.02 0 6.0 3.3 0.2Ours 56 0.83 0.4 0.8 1.3 1.8TWANG 140 1.35 0.1 4.9 1.3 0fail to detect at least 5% of the nuclei. The missing nuclei of GC are mainly at the far end ofthe visual axis, because GC has no mechanism to cope with the variations of illumination. Thequality of our algorithm is similar or better than the other methods, only the number of splitnuclei is elevated.Our algorithm is again best regarding time. However, the running time varied between 10 sand 150 s. Running times of a segmentation algorithm can depend on several factors. Therunning time of our method is primarily affected by the number of clustered nuclei and theimage resolution. The former is reflected in the number of partitioner calls, the latter increasesthe time needed per single partitioner call due to a larger graph. The running time variance ofour method is caused by the variance in clustered objects for this data set, which ranges from200 to 700 nuclei. Errors in an image analysis pipeline are inevitable. An approach to copewith this fact is to quantify the uncertainty in the outcome of an operator [44]. Our algorithmcould be augmented by interpreting the score of our nuclei model as an uncertainty value. Theerror rate could then be reduced for narrow decisions of the nuclei model. A manual analysisof the errors in one image of the set showed that this was the case for 25% of the errors made.15% of the errors could be attributed to the graph partitioner. The remaining errors wouldrequire a more exact nuclei model to be eliminated.
In this work we have presented a novel algorithm to split clustered nuclei by combining graphpartitioning with a simple object model. We have shown the performance of our method ontwo sets of synthetic cell nuclei images. Our method is very fast and performs qualitativelycomparable or better to other cell nuclei segmentation techniques.Future work is necessary for a better automation. The required user input should be reduced,e.g. by unsupervised learning of the model parameters. This would allow the model to includemore features, which could further improve the quality. Another approach to reducing the userinput, but also to allow more sophisticated binarization methods, would be to find a unifyingmodel which can explain all used parameters. For an example, the parameters σ grad and σ bin ,used in the computation of the edge weights, should be deduced from the parameters of theautomatically computed image model. We tried to find such a relationship, but it did not workin our experiments.The anisotropic resolution was one of the biggest challenges for this work. The ratio betweenspacing in axial direction to the lateral directions was as high as in one of our datasets.This negatively affected several aspects of our method. For example, the approximation ofsphericity is hard to obtain, as a nucleus extends only about slices in the axial direction. Thiscan make it difficult to distinguish between a sphere, and a half-sphere. One solution for this21ssue would be to include more features in the nucleus model, for example the extent of theobject in z -direction.The anisotropic resolution is especially problematic for the gradient-based edge weightingmethod. Here, the edge weights are a measure for the probability of a change in intensity.However, such a change has a higher probability when moving along the z -axis, because pixelsalong this axis have a higher spacing. As a result, the weights for the respective edges are toosmall, favouring cuts of these edges. One idea to cope with this problem could be to adjustthe parameter σ grad accordingly. For example, one could use separate parameters dependingon the edge orientation.This work has dealt with image segmentation of 3D images with anisotropic spacing. Anotherway to consider the images is as a stack of 2D images. This perception relates closely to theimage formation process, as the image sensor which digitizes the real image is a 2D arrayof capacitors. One approach could be to segment all slices separately and then combine thedetected 2D objects. This also has the advantage of an easy parallelization. Another ideais to separate the slices for the sphericity measure in the nucleus model: instead of using athree-dimensional sphericity, we could compute a two-dimensional roundness measure for eachinvolved slice.From a computer scientist’s point of view, the problem in the development of image segmen-tation algorithms is that the problem is ill-posed [8, 28]. For a well-posed problem, a solutionexists, is unique, and is stable, i.e. it depends continuously on the data [18]. The image seg-mentation problem does not have a unique solution, because the digitization process is notinjective: multiple different images can be mapped to the same digital image. This is especiallythe case for the segmentation of live cell nuclei in 3D fluorescence microscopy. The low imagingquality and the anisotropic resolution lead to ambiguities. Even human experts can differ ontheir choice of the detection and/or the boundary of an object [7, 8].It is difficult to design a segmentation algorithm when the correct solution is not defined.This problem results from the lack of ground truth for large-scale image datasets. How shouldthe quality of an algorithm be evaluated? How can the quality be optimized? How can itbe compared to other algorithms? Some authors resort to subjective quality comparison [45].Most either compare against hand-segmented datasets or synthetic benchmark datasets [46].The latter choice is also the path taken in this work. We have optimized our algorithm andthe parameters with two sets of synthetic images. While our method shows superior results onthese instances, its performance on real-world images has yet to be evaluated.22 eferences [1] Yousef Al-Kofahi, Wiem Lassoued, William Lee, and Badrinath Roysam. Improved auto-matic detection and segmentation of cell nuclei in histopathology images. IEEE Transac-tions on Biomedical Engineering , 57(4):841–852, 2010.[2] Charles-Edmond Bichot and Patrick Siarry.
Graph partitioning . John Wiley & Sons, 2013.[3] Yuri Boykov and Gareth Funka-Lea. Graph cuts and efficient N-D image segmentation.
International Journal of Computer Vision , 70(2):109–131, 2006.[4] Yuri Boykov and Vladimir Kolmogorov. Computing geodesics and minimal surfaces viagraph cuts. In ,pages 26–33. IEEE Computer Society, 2003.[5] Yuri Boykov and Vladimir Kolmogorov. An experimental comparison of min-cut/max-flowalgorithms for energy minimization in vision.
IEEE Transactions on Pattern Analysis andMachine Intelligence , 26(9):1124–1137, 2004.[6] Yuri Boykov, Olga Veksler, and Ramin Zabih. Fast approximate energy minimization viagraph cuts.
IEEE Transactions on Pattern Analysis and Machine Intelligence , 23(11):1222–1239, 2001.[7] Luís Pedro Coelho, Aabid Shariff, and Robert F. Murphy. Nuclear segmentation in micro-scope cell images: A hand-segmented dataset and comparison of algorithms. In
Proceedingsof the 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro ,pages 518–521. IEEE, 2009.[8] Ondřej Daněk.
Graph Cut Based Image Segmentation in Fluorescence Microscopy . PhDthesis, Masaryk University, 2012.[9] Ondřej Daněk and Pavel Matula. On euclidean metric approximation via graph cuts. InPaul Richard and José Braz, editors,
Computer Vision, Imaging and Computer Graphics.Theory and Applications , volume 229 of
Communications in Computer and InformationScience , pages 125–134. Springer, 2010.[10] Andrew Delong and Yuri Boykov. A scalable graph-cut algorithm for N-D grids. In . IEEEComputer Society, 2008.[11] Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin. Maximum likelihood fromincomplete data via the EM algorithm.
Journal of the Royal Statistical Society. Series B(Methodological) , pages 1–38, 1977.[12] Ralf Diekmann, Robert Preis, Frank Schlimbach, and Chris Walshaw. Shape-optimizedmesh partitioning and load balancing for parallel adaptive FEM.
Parallel Computing , 26(12):1555–1581, 2000.[13] Charles M. Fiduccia and Robert M. Mattheyses. A linear-time heuristic for improvingnetwork partitions. In , pages 175–181. ACM/IEEE,1982. 2314] Michael R. Garey, David S. Johnson, and Larry J. Stockmeyer. Some simplified np-complete problems. In Robert L. Constable, Robert W. Ritchie, Jack W. Carlyle, andMichael A. Harrison, editors,
Proceedings of the 6th Annual ACM Symposium on Theoryof Computing , pages 47–63. ACM, 1974.[15] Andrew V. Goldberg, Sagi Hed, Haim Kaplan, Robert Endre Tarjan, and Renato Fon-seca F. Werneck. Maximum flows by incremental breadth-first search. In Camil Deme-trescu and Magnús M. Halldórsson, editors,
Proceedings of the 19th Annual European Sym-posium on Algorithms , volume 6942 of
Lecture Notes in Computer Science , pages 457–468.Springer, 2011.[16] Rafael C. Gonzalez and Richard E. Woods. Image processing.
Digital image processing , 2,2007.[17] Dorothy M. Greig, Bruce T. Porteous, and Allan H. Seheult. Exact maximum a pos-teriori estimation for binary images.
Journal of the Royal Statistical Society. Series B(Methodological) , 51(2):271–279, 1989.[18] Jacques Hadamard. Sur les problèmes aux dérivées partielles et leur signification physique.
Princeton university bulletin , 13(49-52):28, 1902.[19] Yong He, Hui Gong, Benyi Xiong, Xiaofeng Xu, Anan Li, Tao Jiang, Qingtao Sun, SiminWang, Qingming Luo, and Shangbin Chen. icut: an integrative cut algorithm enablesaccurate segmentation of touching cells.
Scientific Reports , 5, 2015.[20] Bruce Hendrickson and Robert W. Leland. A multi-level algorithm for partitioning graphs.In Sidney Karin, editor,
Supercomputing ’95 , page 28. IEEE Computer Society / ACM,1995.[21] Jan Huisken, Jim Swoger, Filippo Del Bene, Joachim Wittbrodt, and Ernst HK Stelzer.Optical sectioning deep inside live embryos by selective plane illumination microscopy.
Science , 305(5686):1007–1009, 2004.[22] Laurent Hyafil and Ronald L. Rivest. Graph partitioning and constructing optimal de-cision trees are polynomial complete problems. Technical report, IRIA – Laboratoire deRecherche en Informatique et Automatique, 1973.[23] L. Ibanez, W. Schroeder, L. Ng, and J. Cates.
The ITK Software Guide . Kitware, Inc.,second edition, 2005. ISBN 1-930934-15-7.[24] Shinya Inoué. Foundations of confocal scanned imaging in light microscopy. In
Handbookof biological confocal microscopy , pages 1–19. Springer, 2006.[25] Ondrej Jamriska, Daniel Sýkora, and Alexander Hornung. Cache-efficient graph cuts onstructured grids. In ,pages 3673–3680. IEEE Computer Society, IEEE Computer Society, 2012.[26] George Karypis and Vipin Kumar. A fast and high quality multilevel scheme for parti-tioning irregular graphs.
SIAM Journal on scientific Computing , 20(1):359–392, 1998.2427] Philipp J. Keller, Annette D. Schmidt, Joachim Wittbrodt, and Ernst H.K. Stelzer. Re-construction of zebrafish early embryonic development by scanned light sheet microscopy.
Science , 322(5904):1065–1069, 2008. doi: 10.1126/science.1162493.[28] Khaled Khairy and Philipp J. Keller. Reconstructing embryonic development.
Genesis ,49(7):488–513, 2011.[29] Gang Lin, Umesh Adiga, Kathy Olson, John F. Guzowski, Carol A. Barnes, and BadrinathRoysam. A hybrid 3d watershed algorithm incorporating gradient cues and object modelsfor automatic segmentation of nuclei in confocal image stacks.
Cytometry Part A , 56A(1):23–36, 2003.[30] Gang Lin, Monica K. Chawla, Kathy Olson, John F. Guzowski, Carol A. Barnes, andBadrinath Roysam. Hierarchical, model-based merging of multiple fragments for improvedthree-dimensional segmentation of nuclei.
Cytometry Part A , 63(1):20–33, 2005.[31] Jiangyu Liu and Jian Sun. Parallel graph-cuts by adaptive bottom-up merging. In
TheTwenty-Third IEEE Conference on Computer Vision and Pattern Recognition , pages 2181–2188. IEEE Computer Society, 2010.[32] Lifeng Liu and Stan Sclaroff. Region segmentation via deformable model-guided split andmerge. In
ICCV , pages 98–104. IEEE Computer Society, 2001.[33] Xinghua Lou, Ullrich Köthe, Jochen Wittbrodt, and Fred A. Hamprecht. Learning tosegment dense cell nuclei with shape prior. In , pages 1012–1018. IEEE Computer Society, IEEE ComputerSociety, 2012.[34] Marvin Minsky. Microscopy apparatus, December 19 1961. US Patent 3,013,467.[35] Michiel Muller.
Introduction to confocal fluorescence microscopy , volume 69. SPIE press,2006.[36] Nobuyuki Otsu. A threshold selection method from gray-level histograms.
Automatica , 11(285-296):23–27, 1975.[37] Christophe Restif.
Segmentation and evaluation of fluorescence microscopy images . PhDthesis, Oxford Brookes University, 2006.[38] Christophe Restif. Towards safer, faster prenatal genetic tests: Novel unsupervised, auto-matic and robust methods of segmentation of nuclei and probes. In Ales Leonardis, HorstBischof, and Axel Pinz, editors, , vol-ume 3954 of
Lecture Notes in Computer Science , pages 437–450. Springer, 2006.[39] T.W. Ridler and S. Calvard. Picture thresholding using an iterative selection method.
IEEE transactions on Systems, Man and Cybernetics , 8(8):630–632, 1978.[40] Peter Sanders and Christian Schulz. Engineering multilevel graph partitioning algorithms.In Camil Demetrescu and Magnús M. Halldórsson, editors,
Proceedings of the 19th AnnualEuropean Symposium on Algorithms , volume 6942 of
Lecture Notes in Computer Science ,pages 469–480. Springer, 2011. 2541] Peter Sanders and Christian Schulz. Think locally, act globally: Highly balanced graphpartitioning. In
Proceedings of the 12th International Symposium on Experimental Al-gorithms , volume 7933 of
Lecture Notes in Computer Science , pages 164–175. Springer,2013.[42] Jianbo Shi and Jitendra Malik. Normalized cuts and image segmentation.
IEEE Trans.Pattern Anal. Mach. Intell. , 22(8):888–905, 2000.[43] Milan Sonka, Vaclav Hlavac, and Roger Boyle.
Image processing, analysis, and machinevision . Cengage Learning, 2014.[44] Johannes Stegmaier, Arif ul Maula Khan, Markus Reischl, and Ralf Mikut. Challengesof uncertainty propagation in image analysis. In
Proceedings of the 22nd Workshop onComputational Intelligence , pages 55–69. KIT Scientific Publishing, 2012.[45] Johannes Stegmaier, Jens C. Otte, Andrei Kobitski, Andreas Bartschat, Ariel Garcia,G. Ulrich Nienhaus, Uwe Strähle, and Ralf Mikut. Fast segmentation of stained nuclei interabyte-scale, time resolved 3d microscopy image stacks.
PloS one , 9(2):e90036, 2014.[46] Johannes Stegmaier, Julian Arz, Benjamin Schott, Jens C. Otte, Andrei Kobitski, G. Ul-rich Nienhaus, Uwe Strähle, Peter Sanders, and Ralf Mikut. Generating semi-syntheticvalidation benchmarks for embryomics. In
Proceedings of the 13th IEEE InternationalSymposium on Biomedical Imaging , pages 684–688. IEEE, 2016.[47] David Svoboda, Michal Kozubek, and Stanislav Stejskal. Generation of digital phantomsof cell nuclei and simulation of image formation in 3d image cytometry.
Cytometry PartA , 75(6):494–509, 2009.[48] Richard Szeliski.
Computer vision: algorithms and applications . Springer Science & Busi-ness Media, 2010.[49] Raju Tomer, Khaled Khairy, Fernando Amat, and Philipp J. Keller. Quantitative high-speed imaging of entire developing embryos with simultaneous multiview light-sheet mi-croscopy.
Nature Methods , 9(7):755–763, 2012.[50] Jayaram K. Udupa and Supun Samarasekera. Fuzzy connectedness and object definition:Theory, algorithms, and applications in image segmentation.
CVGIP: Graphical Modeland Image Processing , 58(3):246–261, 1996.[51] Hakon Wadell. Volume, shape, and roundness of quartz particles.