
1.IntroductionImage segmentation is a fundamental problem in computer vision and image analysis. In the image segmentation community, level setbased approaches are important tools, because they are able to handle nuclei shapes and contours with complex variations. Chan and Vese initiated a regionbased active contour model with a level set formulation based on MumfordShah’s functional.^{1}^{,}^{2} This model does not depend on gradient information and thus can detect nuclei contours with weak edges. However, Chan–Vese model does not include prior shape knowledge to restrain shapes in appearance and may, therefore, detect meaningless shapes. This drawback becomes more serious when nuclei are partially occluded, corrupted, or represented in low contrast imaging data. In the literature, some research work emerged to mitigate this problem by incorporating shape prior information into the level set formulation so that the detected shape can be regulated by a selected reference shape.^{3}^{–}^{5} For example, Chan and Zhu^{6} proposed a variational model introducing a shape difference term with a shape prior as the reference shape.^{1} With prior information in place, shapes similar to the shape prior can be successfully segmented. Meanwhile, nuclei presenting meaningless shapes are restrained. Based on the shape prior segmentation model, Yan et al.^{7} introduced a geodesic length term to drive contours to nuclei edges. Ali and Madabhushi^{8} proposed an active contour model to accommodate complex shapes by learning shape priors through principal component analysis. Nevertheless, shape priorbased level set segmentation methods are still challenged by multiple problems. (1) It is difficult to segment nuclei from raw images without knowing the number and position of nuclei. Thus, one important step in this class of segmentation approaches is to appropriately identify number and positions of nuclei of interest for initialization purpose. (2) Thus far, shape priorbased segmentation approaches exploit libraries consisting of a small number of shape priors as reference to nuclei in similar shapes. However, in most realworld scenarios, such as nuclei segmentation, it could be very complex to model shape variations explicitly. Therefore, it is practically difficult to find a limited number of shape priors that could represent all shapes reasonably well. (3) Simultaneous segmentation of mutually occluded nuclei remains challenging. Recent development in level set methods introduced a repulsion term to prevent two nuclei from becoming identical through evolution. However, it is common to have more than two nuclei involved in mutual occlusion. As a result, we conceive that larger penalties should be given to regions shared by more nuclei. Nuclei initialization is an important step that guides the following modules in segmentation. The geometric centers of nuclei, socalled seeds, are natural nuclei indicators in practice, as nuclei contours can be appropriately initialized with them. Parvin et al.^{9} proposed an iterative voting method for detecting seeds of overlapping nuclei. AlKofahi et al.^{10} proposed a multiscale Laplacian of Gaussian (LoG) filtering method for automatic detection and segmentation of nuclei. By distancemapbased adaptive scale selection, nuclear seed points were detected to perform an initial segmentation. However, the multiscale LoG method is sensitive to minor peaks on the distance map, resulting in oversegmentation and detection. Qi et al.^{11} improved the method by applying a shifted Gaussian kernel and mean shift to a singlepass voting algorithm. The basic idea is to search for the voting direction over a coneshaped voting area. The twodimensional (2D) Gaussian kernel was designed to assign larger weights for area centers in the voting process. Finally, mean shift was applied to determine the seed points. To accommodate the variation in nuclei shapes, a large set of shape priors can be used as a reference library. However, it is challenging to model the relationship between the nuclei of interest to be segmented and the shape priors. Representation learning methods, such as graph learning and sparse coding, were successfully applied in many areas of computer vision research, such as face recognition,^{12} image restoration,^{13} graph learning,^{14} and medical image analysis.^{15} Wright et al.^{12} proposed a robust face recognition algorithm based on sparse representation. The method can effectively handle occlusion and corruption errors by enforcing sparse distribution with respect to the pixel/sample basis. Zhang et al.^{15} proposed a sparse shape composition model to incorporate shape prior information for detection and segmentation tasks in radiologic images. Cheng et al.^{14} introduced a process to build ${\mathcal{L}}_{1}$ graph and multiple algorithms for data clustering, subspace learning, and semisupervised learning based on ${\mathcal{L}}_{1}$ graph. Inspired by theories in graph learning, spectral clustering, and sparse representation, an ${\mathcal{L}}_{1}$ graphbased shape proximity was learned to cluster the shape priors with which a compact dictionary was created for sparse coding. Given a shape of interest, it can be approximated by a sparse representation based on the compact dictionary within the level set framework. The resulting reconstruction error can be included as a shape penalty term in the variational level set model. To date, numerous investigators have devised methods to segment overlapping nuclei. These methods extended the shape prior segmentation approaches to segment multiple similar nuclei under mutual occlusion. In Refs. 1617.–18, various repulsive force terms were proposed and included in the cost functional to penalize two overlapping nuclei contours and prevent them from becoming identical. In our approach, we extend this idea by introducing an adaptive occlusion penalty term that penalizes regions occupied by multiple nuclei based on the scope of occlusion. Following this principle, we assign a larger penalty to a region overlapped by more nuclei. Recently, deep learningbased approaches greatly improved the stateoftheart in research fields, such as computer vision, speech recognition, and natural language processing. Deep learning methods require largescale annotated data to preserve generative power and prevent overfitting. However, collecting both data and annotation for medical imaging tasks is timeconsuming and requires much domain expertise. Specifically, a large number of bounding boxes or shape masks of nuclei are required to train the network for deep learningbased nuclei detection and segmentation methods. Compared to the formidable data scale required by deep learning methods, the only prior information used in our method is a set of representative shapes in vector representations with only nuclei contour coordinates. In this paper, we propose a level setbased variational model that allows simultaneous segmentation of multiple nuclei with mutual occlusion. The main contributions of our work are summarized as follows:
As shown in Fig. 1, our method framework consists of a seed detection algorithm for nuclei contour initialization, and an integrated contour deformable model that incorporates region, shape, and boundary information. The final nuclei contours are converged through an iterative contour evolution process. This journal paper extends our earlier work^{19}^{,}^{20} through substantial method improvement on shape prior library generation, more comprehensive experiments, more indepth parameter sensitivity analysis. In detail, these important extensions include:
2.Shape RepresentationLet us consider an image that contains $N$ nuclei of interest ${\mathcal{O}}_{1},{\mathcal{O}}_{2},\dots ,{\mathcal{O}}_{N}$, where each nuclei ${\mathcal{O}}_{i}$ has a closed and bounded shape contour $\partial {\mathcal{O}}_{i}$ in image domain $\mathrm{\Omega}\subset {\mathbb{R}}^{2}$. The basic idea of the level set framework is to implicitly represent $\partial {\mathcal{O}}_{i}$ as a zero level set of a Lipschitz function ${\varphi}_{i}:\mathrm{\Omega}\to \mathbb{R}$. ${\varphi}_{i}(\mathbf{x})$ has positive and negative value when $\mathbf{x}$ is inside and outside ${\mathcal{O}}_{i}$, respectively. Note that due to memory limit, we use image patches instead of whole slide images in our experiments. 2.1.Distance Map FunctionA distance map function $\mathrm{\Gamma}[{\varphi}_{i}(\mathbf{x})]$ represents the shortest distance $d(\mathit{x},\partial {\mathcal{O}}_{i})$ from a current pixel $\mathbf{x}$ in the image domain $\mathrm{\Omega}$ to a given nuclei contour $\partial {\mathcal{O}}_{i}=\{\mathbf{x}{\varphi}_{i}(\mathbf{x})=0,\forall \text{\hspace{0.17em}\hspace{0.17em}}\mathbf{x}\in \mathrm{\Omega}\}$: In this way, any nuclei in the image domain can be represented as a distance map or vice versa. Given a distance map, the corresponding nuclei contour can be recovered by searching for the zerovalued pixel locations. As 2D images only capture 2D projections of 3D nuclei, it is not uncommon to observe overlapped nuclei. Instead of partitioning an image into disjoint regions, we allow a pixel to be associated with multiple nuclei in an image with intersecting contours. Therefore, the corresponding level set functions of all those intersecting nuclei have positive values over the overlapped regions.2.2.Shape Prior ModellingInstead of a single shape prior by previous shape prior segmentation methods, we use a large set of shape priors to deal with the complex shape variation observed in most realworld applications. Training shape priors are manually extracted from raw images and aligned with generalized Procrustes analysis,^{21} as illustrated in Fig. 2. The resulting shapes are represented by a set of vectors of uniformly sampled local landmarks: $\{{\gamma}_{1},{\gamma}_{2},\dots ,{\gamma}_{K}\}$, where $K$ is the size of the shape prior set. For better computational efficiency, the shape priors are partitioned into $M$ clusters, where $M\ll K$. We cluster the shape priors with a learned ${\mathcal{L}}_{1}$ graph learned by solving the minimization problem for each ${\mathit{\gamma}}_{i}$:^{14} Eq. (2)$${\widehat{\mathit{z}}}_{i}={\Vert {\mathit{\gamma}}_{i}B{\mathit{z}}_{i}\Vert}_{2}^{2}+\lambda {\Vert {\mathit{z}}_{i}\Vert}_{1},\phantom{\rule[0.0ex]{1em}{0.0ex}}i=1,\dots ,K\phantom{\rule[0.0ex]{1em}{0.0ex}}\mathrm{s.t.}\text{\hspace{0.17em}\hspace{0.17em}}{z}_{ii}=0,$$Eq. (3)$${W}_{ij}=({\widehat{z}}_{ij}+{\widehat{z}}_{ji})/2,\phantom{\rule[0.0ex]{1em}{0.0ex}}i,j=1,\dots ,K.$$2.3.Shape AlignmentTo represent a nuclei ${\mathcal{O}}_{i}$ with shape priors in the training set $\mathrm{\Psi}$, we first align the associated distance map $\mathrm{\Gamma}[{\varphi}_{i}(\mathit{x})]$ to shape priors by rotation and translation. The transformation from a pixel $\mathit{x}={[{x}_{1}\text{\hspace{0.17em}}{x}_{2}]}^{T}$ to the corresponding point ${\tilde{\mathit{x}}}_{i}$ in shape prior after alignment is formulated as follows: Eq. (4)$${\tilde{\mathit{x}}}_{i}={s}_{i}\left[\begin{array}{cc}\mathrm{cos}({\theta}_{i})& \mathrm{sin}({\theta}_{i})\\ \mathrm{sin}({\theta}_{i})& \mathrm{cos}({\theta}_{i})\end{array}\right]\mathit{x}+{\mathit{t}}_{i},$$2.4.Sparse Shape RepresentationGiven the set of distance maps derived by mapping the $i$’th nuclei to shape prior set $\mathrm{\Psi}({\tilde{\mathit{x}}}_{i})$, we assume that the distance map $\mathrm{\Gamma}({\varphi}_{i}(\mathit{x}))$ of ${\mathcal{O}}_{i}$ can be approximately represented as a linear composition of shape priors in $\mathrm{\Psi}({\tilde{\mathit{x}}}_{i})$. The distance maps and shape priors are vectorized and denoted as $\overline{\mathrm{\Gamma}}[{\varphi}_{i}(\mathit{x})]$ and $\overline{\mathrm{\Psi}}({\tilde{x}}_{i})=({\overline{\psi}}_{1},{\overline{\psi}}_{2},\dots ,{\overline{\psi}}_{M})$, respectively. By the linearity assumption, $\overline{\mathrm{\Gamma}}[{\varphi}_{i}(\mathit{x})]$ can be represented as follows: Eq. (5)$$\overline{\mathrm{\Gamma}}[{\varphi}_{i}(x)]=\overline{\mathrm{\Psi}}({\tilde{x}}_{i}){\alpha}_{i}+{e}_{i},$$To avoid the curse of dimensionality problem, we reduce $\overline{\mathrm{\Gamma}}$ dimensionality, leading to a less computational cost. The process of dimension reduction can be modeled as left multiplication by a nonzero projection matrix $R\in {\mathbb{R}}^{d\times m}$, $d\ll m$, where $m$ is the total number of pixels in an image. As proved in 12, the choice of matrix $R$ does not critically affect the ability to recover the sparse solution. For computation simplicity, we compose a matrix $R$ with entries randomly generated from standard Gaussian distribution and subject to unit length for all rows. The corresponding low dimension representations are denoted as $\stackrel{\u02c7}{\mathrm{\Gamma}}({\varphi}_{i})=R\overline{\mathrm{\Gamma}}({\varphi}_{i})$, $\stackrel{\u02c7}{\mathrm{\Psi}}=R\overline{\mathrm{\Psi}}$, and ${\stackrel{\u02c7}{\mathit{e}}}_{i}=R{\mathit{e}}_{i}$. Thus, we have Eq. (6)$$\stackrel{\u02c7}{\mathrm{\Gamma}}[{\varphi}_{i}(x)]=\stackrel{\u02c7}{\mathrm{\Psi}}({\tilde{x}}_{i}){\alpha}_{i}+{\stackrel{\u02c7}{e}}_{i}=\stackrel{\u02c7}{\mathrm{\Psi}}({\tilde{x}}_{i}){\alpha}_{i}+I{\stackrel{\u02c7}{e}}_{i},$$Eq. (7)$${\widehat{\alpha}}_{i},{\widehat{e}}_{i}=\underset{{\alpha}_{i},{\stackrel{\u02c7}{e}}_{i}}{\mathrm{arg}\mathrm{min}}{\Vert \stackrel{\u02c7}{\mathrm{\Gamma}}[{\varphi}_{i}(x)]\stackrel{\u02c7}{\mathrm{\Psi}}({\tilde{x}}_{i}){\alpha}_{i}I{\stackrel{\u02c7}{e}}_{i}\Vert}_{2}^{2}\phantom{\rule{0ex}{0ex}}+{\lambda}_{1}{\Vert {\alpha}_{i}\Vert}_{0}+{\lambda}_{2}{\Vert {\stackrel{\u02c7}{e}}_{i}\Vert}_{0}.$$However, the problem of finding the sparsest solution in an underdetermined linear system such as Eq. (7) is proved to be nondeterministic polynomialtime hard.^{24} Due to recent developments in sparse theory, it has been revealed that such a problem can be solved via ${\mathcal{L}}_{1}$ relaxation:^{25}^{–}^{27} Eq. (8)$${\widehat{\alpha}}_{i},{\widehat{e}}_{i}=\underset{{\alpha}_{i},{\stackrel{\u02c7}{e}}_{i}}{\mathrm{arg}\mathrm{min}}{\Vert \stackrel{\u02c7}{\mathrm{\Gamma}}[{\varphi}_{i}(x)]\stackrel{\u02c7}{\mathrm{\Psi}}({\tilde{x}}_{i}){\alpha}_{i}I{\stackrel{\u02c7}{e}}_{i}\Vert}_{2}^{2}\phantom{\rule{0ex}{0ex}}+{\lambda}_{1}{\Vert {\alpha}_{i}\Vert}_{1}+{\lambda}_{2}{\Vert {\stackrel{\u02c7}{e}}_{i}\Vert}_{1}.$$3.Seed DetectionIn this section, we present an algorithm to recognize nuclei spatial locations by nuclei seed detection. Nuclei seed detection is essential for the follow up nuclei segmentation, as it decides the number and initial contour locations. The proposed seed detection algorithm utilizes joint information of spatial connectivity, distance constraint, image edge map, and a shapebased voting map derived from eigenvalue analysis of Hessian matrix across multiple scales. 3.1.Initialization and PreprocessingThe pathology images of tumor specimens for routine diagnostics usually contain two primary chemical stains: hematoxylin and eosin (H&E). Hematoxylin presents a dark blue or violet color, which positively binds to nuclei. We decompose signal intensity components of H&E stains in original images and use only signals of hematoxylin stain for analysis. Based on Lambert–Beer’s law, the relationship between the intensity of light entering a specimen ${L}_{i}$ and that through a specimen ${L}_{o}$ can be characterized as ${L}_{o}={L}_{i}{e}^{(Ab)}$, where $A$ and $b$ are the amount of stain and the absorption factor, respectively. By this law, we can compute the stainspecific strength values with a predefined stain decomposition matrix $M$ and the observed optical density (OD): $Y=\mathrm{log}({L}_{o}/{L}_{i})$.^{28} The stain composition matrix $M$ consists of three normalized column vectors representing OD values of red, green, and blue channels from hematoxylin, eosin, and null stain. Therefore, stain specific optical signals can be computed by $X={M}^{1}Y$. We retain the first entry of the resulting stain specific signal vector at each pixel and denote the decomposed hematoxylin channel signal as ${I}_{H}$. Note that ${I}_{H}$ is often coupled with noise produced during the slide preparation process, including uneven tissue slide cut, heterogeneous histology process, staining artifacts, and digital scanning noise. We can normalize the background noise in ${I}_{H}$ by morphological reconstruction with two morphological components, namely mask image ${I}_{\cap}$ and marker image ${I}_{\cup}$. Initially, we set ${I}_{\cap}$ as complementary image of ${I}_{H}$ and ${I}_{\cup}$ as ${I}_{H}\oplus {\rho}_{r}$, where $\oplus $ denotes the morphological dilation operator and ${\rho}_{r}$ is a circular structural element with radius $r$. With marker and mask image, an image $R({I}_{\cup},{I}_{\cap})$ can be reconstructed by iterative dilation and intersection until convergence.^{29} The difference image $h({I}_{\cup},{I}_{\cap},{\rho}_{r})={I}_{\cap}R({I}_{\cup},{I}_{\cap})$ consists of a near zerolevel background and a group of enhanced foreground peaks, each representing an object of interest. In Fig. 4(a), we present a typical background normalization result with morphological reconstruction. 3.2.Voting Map ConstructionAssuming a typical nuclei shape is close to either a circle or an ellipse, we proceed with a shapebased voting process by analyzing nuclei structures within local image regions. Before the voting process, we first enhance nuclei regions by eigenvalue analysis with the Hessian matrix convolved with Gaussian filters at multiple scales.^{30} With this approach, we search for circular structures based on geometric structures characterized by the neighboring pixel intensity profiles. Specifically, for a pixel at $({i}_{0},{j}_{0})$, its local image intensity change can be represented by Taylor expansion: Eq. (9)$${h}_{\sigma}({i}_{0}+\delta i,{j}_{0}+\delta j)={h}_{\sigma}({i}_{0},{j}_{0})+(\delta i,\delta j)\mathcal{D}({h}_{\sigma}({i}_{0},{j}_{0}))\phantom{\rule{0ex}{0ex}}+\frac{1}{2}(\delta i,\delta j){\mathcal{D}}^{2}({h}_{\sigma}({i}_{0},{j}_{0})){(\delta i,\delta j)}^{T}\phantom{\rule{0ex}{0ex}}+\mathcal{O}[({\delta}^{3})],$$3.3.Dynamic Seeds Detection and MergingGiven the derived voting map, we dynamically adjust a seed list based on candidate spatial connectivity, distance constraint, and image edge map. The proposed method can produce robust and accurate seed detection result, especially for overlapped nuclei. As any pixel on the voting map is no larger than the number of scales $S$, we consider the voting map as a surface in a threedimensional space, as shown in Fig. 4(b). Those strong voting peaks, representing consistent negative eigenvalue pairs from Hessian matrix at different smoothing scales due to radially outward decreasing pixel intensity profiles in local regions, suggest the presence of nuclei. The strong peaks on the voting surface can be detected as we gradually slide down an imaginary horizontal plane (e.g., from the blue to green plane) intersecting with the voting surface. We can generate a binary image from the original voting map with threshold $v$ for each intersection plane. We begin with an empty seed list $L$ and append to it the resulting centroids of voting peaks satisfying all the following conditions: (a) such voting peak centers do not exist in the peak list $L$ yet. (b) They come from strong peaks that suggest consistent local intensity change property, a key to a reliable nuclei detection. Such strong peaks can be identified with their sizes no less than area threshold $A(v)=\eta +\mathrm{exp}[\rho f(v)]$, where $\rho $ and $\eta $ are predefined scalars that determine the lower area threshold for detected peaks in the voting map at each step and $f(v)=\frac{{v}_{\mathrm{max}}v}{{v}_{\mathrm{max}}{v}_{\mathrm{min}}}$ is a function normalizing the voting value $v$; specifically, $\eta +1$ is the lowest peak crosssection size expectation when the imaginary horizontal plane is the highest. $\rho $ determines how peak areas are expected to grow as the imaginary horizontal plane slides down. Their choices depend on the nuclei texture homogeneity and nuclei size. We choose $\eta $ to be small enough not to miss any strong peaks, and $\rho $ big enough to represent the rapid growth of strong voting peak crosssections as the imaginary horizontal plane slides down. (c) Such points are within the foreground region in the binary mask detected by the adaptive thresholding method. With a list of seed candidates, we perform seed pruning to eliminate false positive seeds. We compute the pairwise distances for all peaks in $L$ and iteratively merge adjacent centroids within distance threshold ${D}_{1}$. This is followed by a second round of distancebased merging with a relaxed threshold ${D}_{2}$ (${D}_{1}<{D}_{2}$). In the second round, two peaks are merged only if the following conditions are true: (a) distance between these two points is less than ${D}_{2}$ and (b) the path connecting a pair of points does not intersect with any canny edge derived from the original image. The edge map is used to prevent centroids of closely clumped nuclei from being merged. With this seed merging process, we can retain seeds from true clumped nuclei in histological images without tedious parameter tuning process. 4.Occluded Nuclei SegmentationWith nuclei detected in Sec. 3, we further develop an improved level set based segmentation method to drive initialized nuclei contours to true nuclei edges. Our development is driven by the use case of automatic analysis of occluded nuclei in wholeslide histopathologic images. In our method, each nuclei is described by a level set function and we aim to simultaneously obtain all $N$ level set functions by optimizing a variational model. 4.1.Prior Information IntegrationAs described in Sec. 2.4, given a large enough set of shape priors $\mathrm{\Psi}=\{{\psi}_{j}\}$, a shape of interest ${\varphi}_{i}$ can be encoded as sparse linear superposition of shape priors. Thus, we define a shape term as follows: Eq. (10)$${E}_{1}=\sum _{i=1}^{N}{\int}_{\mathrm{\Omega}}{[\mathrm{\Gamma}[{\varphi}_{i}(\mathit{x})]\sum _{j=1}^{M}{\psi}_{j}({\tilde{x}}_{i}){\widehat{\alpha}}_{ij}]}^{2}\mathrm{d}x,$$4.2.Adaptive Penalty on Mutual OcclusionIt is common to have mutually occluded nuclei in 2D histopathologic images due to the fact that these images represent 2D projected signals of tissue nuclei in 3D space. It is a challenging task to segment mutually occluded nuclei and identify hidden boundaries across each other. This problem becomes exponentially complicated when there are more than two nuclei involved in occlusion. In the level set framework, intersecting nuclei may all have positive function values after contour evolutions, making it difficult to differentiate them from each other. To address this problem, we introduce an adaptive occlusion penalty term to dynamically suppress nuclei intersection events. The occlusion penalty term is determined by the number of nuclei that are overlapped. Meanwhile, this term prevents deformable contours from becoming identical after iterations of evolution. Specifically, we define the adaptive occlusion penalty term as follows: Eq. (11)$${E}_{2}={\int}_{\mathrm{\Omega}}\sum _{i=1}^{N}H[{\varphi}_{i}(x)][1\prod _{k=1}^{N}(1H({\varphi}_{k}(\mathit{x})))]\mathrm{d}x,$$4.3.Edge Guided Contour and Evolution StabilityTo further encourage contour convergence and retain contour smoothness, we define an edgeguided contour regularity term as follows: Eq. (12)$${E}_{3}=\sum _{i=1}^{N}{\int}_{\mathrm{\Omega}}Q(x)\nabla H[{\varphi}_{i}(x)]\mathrm{d}x,$$Additionally, we include an evolutionary stability term^{31} to regulate the property of level set function as follows: where $R(x)$ is defined as a doublewell potential function: This term in Eq. (13) is used to retain signed distance property for stable level set function computation.4.4.Improved Variational Level Set FormulationBy combining terms in Eqs. (10)–(13), we formulate our method with a variational level set framework.^{1} To evolve nuclei contours to desired nuclei boundaries, we minimize the following functional: Eq. (15)$$E(C,\mathrm{\varphi},\mathrm{\Theta},T)={\int}_{\mathrm{\Omega}}\mathcal{L}(x,\mathrm{\varphi},C,\mathrm{\Theta},T)\mathrm{d}x={E}_{cv}+\sum _{i=1}^{4}{E}_{i}\phantom{\rule{0ex}{0ex}}={\lambda}_{o}\sum _{i=1}^{N}{\int}_{\mathrm{\Omega}}{(I{c}_{i})}^{2}H[{\varphi}_{i}(x)]\mathrm{d}\mathit{x}\phantom{\rule{0ex}{0ex}}+{\lambda}_{b}{\int}_{\mathrm{\Omega}}{(I{c}_{b})}^{2}\prod _{i=1}^{N}[1H({\varphi}_{i}(x))]\mathrm{d}\mathit{x}\phantom{\rule{0ex}{0ex}}+\nu \sum _{i=1}^{N}{\int}_{\mathrm{\Omega}}{[\mathrm{\Gamma}({\varphi}_{i}(x))\sum _{j=1}^{M}{\psi}_{j}({x}_{i}){\widehat{\alpha}}_{ij}]}^{2}\mathrm{d}x\phantom{\rule{0ex}{0ex}}+\omega {\int}_{\mathrm{\Omega}}\sum _{i=1}^{N}H[{\varphi}_{i}(x)][1\prod _{k=1}^{N}(1H({\varphi}_{k}(x)))]\mathrm{d}x\phantom{\rule{0ex}{0ex}}+\mu \sum _{i=1}^{N}{\int}_{\mathrm{\Omega}}Q(x)\nabla H[{\varphi}_{i}(x)]\mathrm{d}x+\xi \sum _{i=1}^{N}{\int}_{\mathrm{\Omega}}R(\nabla {\varphi}_{i}(x))\mathrm{d}x,$$5.Numerical ComputationsIn this section, we present numerical computations to minimize the functional in Eq. (15). We optimize the functional iteratively by updating functions $\mathrm{\varphi}$ and variables $\{C,\mathrm{\Theta},T\}$ alternatively. We begin with updating functions $\mathrm{\varphi}$ first. Parameterizing iteration as an artificial time variable $t>0$, we minimize the functional by solving Euler–Lagrange equation based on theory of calculus of variations: Eq. (16)$$\frac{\partial {\varphi}_{i}}{\partial t}=\frac{\partial E}{\partial {\varphi}_{i}}=(\frac{\partial \mathcal{L}(\mathbf{x},{\varphi}_{i})}{\partial {\varphi}_{i}}\sum _{j=1}^{2}\frac{\partial}{\partial {x}_{j}}(\frac{\partial \mathcal{L}(\mathbf{x},{\varphi}_{i})}{\partial \frac{\partial {\varphi}_{i}}{\partial {x}_{j}}}\left)\right)\phantom{\rule{0ex}{0ex}}={\lambda}_{o}{(I{c}_{i})}^{2}\delta ({\varphi}_{i})\phantom{\rule{0ex}{0ex}}+{\lambda}_{b}{(I{c}_{b})}^{2}\prod _{k=1,k=i}^{N}(1H({\varphi}_{k}(\mathit{x})))\delta ({\varphi}_{i})\phantom{\rule{0ex}{0ex}}+\mu \delta ({\varphi}_{i})(\nabla Q(\mathit{x})\xb7\frac{\nabla {\varphi}_{i}}{{\varphi}_{i}}+Q(\mathit{x})div(\frac{\nabla {\varphi}_{i}}{\nabla {\varphi}_{i}}\left)\right)\phantom{\rule{0ex}{0ex}}+\xi div(\frac{{R}^{\prime}(\nabla {\varphi}_{i})}{\nabla {\varphi}_{i}}\nabla {\varphi}_{i})\phantom{\rule{0ex}{0ex}}\omega \delta ({\varphi}_{i})(1\prod _{k=1,k=i}^{N}(1H({\varphi}_{k}(\mathit{x}))))\phantom{\rule{0ex}{0ex}}2\nu (\mathrm{\Gamma}({\varphi}_{i}(\mathit{x}))\sum _{j=1}^{M}{\psi}_{j}({\tilde{\mathit{x}}}_{i}){\widehat{\alpha}}_{ij}),$$Next, we fix $\mathrm{\varphi}$ and update transformation parameters. We derive updating equations for transformation parameters by computing gradient descent of functional $E(\mathrm{\Theta},T)$: Eq. (17)$$\frac{\partial {\theta}_{i}}{\partial t}=2\nu {\int}_{\mathrm{\Omega}}(\mathrm{\Gamma}({\varphi}_{i}(\mathit{x}))\sum _{j=1}^{M}{\psi}_{j}({\tilde{\mathit{x}}}_{i}){\widehat{\alpha}}_{ij})\phantom{\rule{0ex}{0ex}}(\sum _{j=1}^{M}({\nabla}_{{\tilde{\mathit{x}}}_{i}}{\psi}_{j}\xb7{\nabla}_{{\theta}_{i}}{\tilde{\mathit{x}}}_{i}){\widehat{\alpha}}_{ij})\mathrm{d}\mathbf{x},$$Eq. (18)$$\frac{\partial {t}_{ik}}{\partial t}=2\nu {\int}_{\mathrm{\Omega}}(\mathrm{\Gamma}({\varphi}_{i}(x))\sum _{j=1}^{M}{\psi}_{j}({\tilde{x}}_{i}){\widehat{\alpha}}_{ij})\phantom{\rule{0ex}{0ex}}(\sum _{j=1}^{M}({\nabla}_{{\tilde{\mathit{x}}}_{i}}{\psi}_{j}\xb7{\nabla}_{{t}_{ik}}{\tilde{x}}_{i}){\widehat{\alpha}}_{ij})\mathrm{d}x,$$Eq. (19)$${\nabla}_{{t}_{ik}}{\tilde{x}}_{i}=\{\begin{array}{cc}{[\begin{array}{cc}1& 0\end{array}]}^{T},& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}k=1\\ {[\begin{array}{cc}0& 1\end{array}]}^{T},& \text{if}\text{\hspace{0.17em}\hspace{0.17em}}k=2\end{array}.$$Finally, we fix $\mathrm{\varphi}$, and $\{\mathrm{\Theta},T\}$, and update ${c}_{i}$ and ${c}_{b}$ by setting the partial derivative of $E({c}_{i},{c}_{b})$ with respective to ${c}_{i}$ and ${c}_{b}$ and solving these equations, respectively. The optimal values turn out to be the average image intensities in corresponding area:^{1} Eq. (20)$$\frac{\partial E({c}_{i},{c}_{b})}{\partial {c}_{i}}=0\Rightarrow {c}_{i}=\frac{{\int}_{\mathrm{\Omega}}I(x)H({\varphi}_{i})\mathrm{d}x}{{\int}_{\mathrm{\Omega}}H({\varphi}_{i})\mathrm{d}x},$$Eq. (21)$$\frac{\partial E({c}_{i},{c}_{b})}{\partial {c}_{b}}=0\Rightarrow {c}_{b}=\frac{{\int}_{\mathrm{\Omega}}I(x)\prod _{i=1}^{N}[1H({\varphi}_{i}(x))]\mathrm{d}x}{{\int}_{\mathrm{\Omega}}\prod _{i=1}^{N}[1H({\varphi}_{i}(x))]\mathrm{d}x}.$$In this way, this minimization problem is solved by iterative computation of Euler–Lagrange equation and gradient descent approach until the functional is converged. The zero level sets of the converged functions $\mathrm{\varphi}$ indicate the final nuclei contours. 6.Experimental Results6.1.Dataset and Parameter SettingWe present experimental results of our algorithm for analysis of nuclei within histopathologic images of glioblastoma brain tumor specimens. The effectiveness of our algorithm is verified with two datasets of H&E stained GBM specimens captured at $40\times $ magnification, namely GBM40 dataset and TCGA FFPE dataset: http://cancergenome.nih.gov/. These images are obtained after glioblastoma brain tissues are processed with a tissue preparation protocol. Image patches with size $512\times 512$ are used for experiments due to memory limit. For GBM40 dataset, there are 5396 and 1849 manually annotated nuclei for seed detection and boundary delineation, respectively. For TCGA dataset, the total number of annotated nuclei for boundary segmentation is 4961. Shape profiles of 27,000 glioblastoma brain nuclei are manually extracted from the dataset to form the set of shape priors. All shape priors are aligned with generalized Procrustes analysis.^{21} As discussed in Secs. 2.1 and 2.2, we represent shapes as distance maps, cluster shape priors with an ${L}_{1}$ graph into $M$ groups and use one representative shape from each group to form a training shape dictionary $\mathrm{\Psi}$ for sparse representation. We apply the proposed method to datasets with the following parameter setup: $r=10$, ${\sigma}_{i}=\{3.3,3.6,\dots ,10\}$, ${D}_{1}=15$, ${D}_{2}=25$, ${\lambda}_{o}={\lambda}_{b}=1$, $\mu =5000$, $\xi =2$, $\omega =2000$, $\nu =3000$, $M=100$, $\eta =10$, $\rho =10$. Note that our approach is an image data driven process. Therefore, the scanner setting such as magnification factor does not have significant impact on the parameter setting. As seeds are detected in the voting map that is produced by eigenvalue analysis of image local Hessian matrix associated with a set of Gaussian filter scales, transforming from the original image to voting map can partially help take care of nuclei scale variations. In our study, this set of Gaussian filter scales is chosen to cover varying nuclei sizes, with $3\text{\hspace{0.17em}}{\mathrm{max}}_{i}({\sigma}_{i})$ approximately equal to the radius of a typical large nucleus. For the weighting parameters for nuclei segmentation, we assign parameter values so that all terms in Eq. (15) are appropriately balanced in numerical values. In our experiment, following three parameters have similar value and scale: the coefficient for edgegradient weighted contour length $\mu $, the coefficient for the occlusion penalty term $\omega $, and the coefficient for squared fitting error of shapederived distance map $\nu $. The weight $\xi $ is numerically set by referring to typical value of the doublewell potential function. We use the l1ls MATLAB toolbox^{32} to solve the ${\mathcal{L}}_{1}$minimization problem. 6.2.Seed DetectionFor quantitative analysis, we assess performance of seed detection method with reference to annotations by a pathologist. Note that we evaluate our approach with nontouching and occluded nuclei in each image separately. Five metrics are computed from each image to show seed detection performance: nuclei number error ($C$), miss detection ($M$), false recognition ($F$), over segmentation ($O$), under segmentation ($U$), and count error (CE)%. Nuclei number error is used to demonstrate the absolute difference between the number of nuclei detected by machine and that reported by human expert. Miss and false detection represent the number of missing and false recognition events when machinebased seed detection method is used to detect individual nuclei with no occluded neighbors. Meanwhile, we use over and undersegmentation to record events where the number of machinedidentified nuclei from a nuclei clump is more and less than the true number marked by human expert, respectively. Finally, we use the count error% to represent the nuclei number error in reference to the true nuclei number. The resulting outputs are shown in Tables 1 and 2. Additionally, we compare our method with a singlepass voting method based on mean shift^{11} and an oriented kernelbased iterative voting method.^{9} By contrast to our method, both the singlepass voting method and the iterative voting method tend to miss detecting nuclei. Note that GBM40 dataset is produced by our local laboratory. Compared to the public TCGA dataset, GBM40 dataset is more carefully prepared with better staining and contrast, leading to higher SNR. This results in better performances of all methods for comparison consistently. We also illustrate seed detection results on several typical image patches in Fig. 5. Table 1Imagewise seed detection performance on GBM40 dataset.
Note: Bold numbers indicate the best performance. Table 2Imagewise seed detection performance on TCGA dataset.
Note: Bold numbers indicate the best performance. 6.3.Nuclei SegmentationWe model the nuclei segmentation problem by the variational level set framework mathematically described in Eq. (15) and solve it by the numerical computations in Sec. 4. The level set functions do not stop evolving until they either reach convergence or exceed iteration number limit. We present the evolution results of zero level sets at iteration 10, 20, and 30 in Fig. 6. The detected nuclei shapes are well defined, as shown in Fig. 6. Notably, those overlapping nuclei can also be correctly segmented, with occluded boundaries reasonably recovered. It is also observed that most zero level sets can rapidly converge to true nuclei boundaries within 10 iterations. After that, zero level sets fine tune themselves to better fit to nuclei contour details, especially over those overlapped nuclei. With our experiment data, we only observe minor improvements in nuclei segmentation from 10 to 30 iterations. In general, the optimal number of iterations depends on data properties. End users can determine the accuracyspeed tradeoff and select the best number of iterations based on individual experiment scenarios. In addition, we present experimental results of three images presenting the best, median, and worst overall segmentation performance in Fig. 7. The bar charts in Fig. 7 present our method efficacy measured by the Jaccard coefficient, precision, recall, and ${F}_{1}$ score. The small variation in these metrics over the best, median, and worst images suggests a good consistency and generalization ability of our method. To quantitatively assess our method’s performance, we use human annotations from pathologist as ground truth. By comparing the experimental result $B$ with ground truth annotation $A$, we evaluate the performance of our method with multiple metrics:^{36}^{,}^{37} Jaccard coefficient $J(A,B)=\frac{A\cap B}{A\cup B}$, precision rate $P(A,B)=\frac{A\cap B}{B}$, recall rate $R(A,B)=\frac{A\cap B}{A}$, F1 score ${F}_{1}(A,B)=2\frac{P(A,B)R(A,B)}{P(A,B)+R(A,B)}$, and Hausdorff distance $H(A,B)=\mathrm{max}(h(A,B),h(B,A))$, where $h(A,B)=\underset{a\in A}{\mathrm{max}}\text{\hspace{0.17em}}\underset{b\in B}{\mathrm{min}}\Vert ab\Vert $. The testing results evaluated with these metrics are presented in Table 3. Our proposed method achieves better performance in most metrics than the other methods. A qualitative performance comparison across several compared methods and our proposed method is shown in Fig. 8. We also demonstrate the results of the proposed nuclei segmentation method with different seed detection methods in Fig. 10. It is notable that our proposed method is able to better capture nuclei boundaries than other methods for comparison. In particular, only our method can recover boundaries of overlapped nuclei by comparisons. This property makes our method superior to the other methods when analytics of occluded nuclei is crucial in investigations. Table 3The performance comparison of different nuclei segmentation methods measured by region and distancebased metrics on GBM40 dataset and TCGA dataset.
Note: Bold numbers indicate the best performance. 6.4.Parameter Sensitivity Analysis for SegmentationWe further investigate the segmentation contribution from individual terms in our variational model by testing with different parameters and comparing their associated segmentation results. The segmentation result with our default parameter setting is shown in Fig. 9(a). When $\nu =0$, the shape prior fitting term ${E}_{1}$ does not take effect in the contour deformation process. The resulting segmentation outcome is presented in Fig. 9(b), where finalized nuclei contours are less regulated. Similarly, we can remove the occlusion penalty term ${E}_{2}$ from the variational model by setting $\omega =0$. The associated result is illustrated in Fig. 9(c). Under this setting, the detected nuclei present a strong inclination to overlap with each other. When the shape prior fitting term ${E}_{1}$, the dynamic occlusion penalty term ${E}_{2}$, and the evolutionary stability term ${E}_{4}$ are all removed ($\nu =\omega =\xi =0$), the resulting nuclei contours become significantly degraded, as shown in Fig. 9(d). Note that shapes appear to be less regulated in Figs. 9(b)–9(d) without one or more terms in the variational model. In addition to the final results with only a subset of terms in Figs. 9(b)–9(d), we also investigate the sensitivity of parameters to final results. Our investigations with parameter deviation from our proposed value set suggest that results remain similar even when we change $\mu $ by 22% ($\mu =3900$), 10% ($\mu =4500$), and $\nu $ by 67% ($\nu =5000$), as presented in Figs. 9(e)–9(g). Figure 9(h) presents results when $\mu $ and $\nu $ are simultaneously changed by 22% and 67%. Overall, a larger $\mu $ leads to a better contour convergence to true nuclei boundary by energy term ${E}_{3}$; a larger $\nu $ forces contours to look more from the reconstructed sparse shape priors by ${E}_{1}$; a larger $\omega $ tends to prevent nuclei more from overlapping with each other by ${E}_{2}$. To illustrate the effect of seed detection, we demonstrate the segmented nuclei contours in Fig. 10 by replacing the seed detection algorithm in our pipeline with other seed detection methods. Since the number of detected nuclei depends on the seedbased initialization, miss detection of seeds may result in undersegmentation of nuclei, as shown in Figs. 10(b) and 10(c). We also test the impact of the shape prior library size on segmentation by changing the number of shape prior clusters $M$. The segmented nuclei contours in Fig. 11 show that the proposed method does not degrade significantly with less shape priors. In practice, we carefully choose the number of clusters $M=100$ so that $M$ is large enough to cover shape variations, but small enough to avoid high computation burden. 6.5.Limitation and Future WorkAlthough the proposed method can achieve good quantitative and visual results, it is limited by the following factors: (a) Nuclei contour evolution depends on accurate detection of seeds, therefore, the segmentation performance may degrade when seeds are not correctly detected. In cases where seeds are missing, our proposed method, like other level set based methods in literature, does not produce correct segmentation result as no level set function is correctly initialized for deformation. However, we have shown in our experimental tests that our proposed seed detection has a good performance with a very low missing rate, as suggested in Tables 1 and 2. When nuclei seeds are slightly shifted within nuclei regions, our approach is able to produce similar and correct segmented nuclei contours for overlapping nuclei. (b) Application of our method to whole slide images can be timeconsuming compared to other stateoftheart methods. To address the problem, we will further develop a MapReduce based high performance image analysis framework to make this process scalable and cost effective.^{43}^{,}^{44} 7.ConclusionIn this paper, we propose a nuclei segmentation method based on the level set framework, aiming to simultaneously identify contours of multiple nuclei with mutual occlusion. We present our method with its application to nuclei segmentation using histopathologic images of glioblastoma brain tumors. First, a seed detection method is developed to automatically initialize nuclei contours. For better nuclei contour deformation, we incorporate into the model a set of typical nuclei shapes as prior information through spectral clustering on ${L}_{1}$ graph. Meanwhile, an adaptive nuclei occlusion penalty term is designed to dynamically penalize nuclei contour occlusion based on the number of overlapped nuclei. It also prevents recognized contours of overlapped nuclei from being identical. A numerical optimization algorithm is used to iteratively search for the desired level set functions. Experiments on glioblastoma brain histopathologic images produce better results as compared with other stateoftheart experiments, suggesting the effectiveness of our detection and segmentation approach for nuclei analysis with glioblastoma histopathologic images. AcknowledgmentsThis research is supported in part by grants from National Institute of Health under Grant No. K25CA181503, National Science Foundation under Grant Nos. ACI 1443054 and IIS 1350885, and CNPq. ReferencesT. F. Chan and L. Vese,
“Active contours without edges,”
IEEE Trans. Image Process., 10
(2), 266
–277
(2001). https://doi.org/10.1109/83.902291 IIPRE4 10577149 Google Scholar
D. Mumford and J. Shah,
“Optimal approximations by piecewise smooth functions and associated variational problems,”
Commun. Pure Appl. Math., 42
(5), 577
–685
(1989). https://doi.org/10.1002/cpa.3160420503 CPMAMV 00103640 Google Scholar
A. Tareef et al.,
“Automatic segmentation of overlapping cervical smear cells based on local distinctive features and guided shape deformation,”
Neurocomputing, 221 94
–107
(2017). https://doi.org/10.1016/j.neucom.2016.09.070 NRCGEO 09252312 Google Scholar
F. Xing, Y. Xie and L. Yang,
“An automatic learningbased framework for robust nucleus segmentation,”
IEEE Trans. Med. Imaging, 35
(2), 550
–566
(2016). https://doi.org/10.1109/TMI.2015.2481436 ITMID4 02780062 Google Scholar
Z. Lu, G. Carneiro and A. P. Bradley,
“An improved joint optimization of multiple level set functions for the segmentation of overlapping cervical cells,”
IEEE Trans. Image Process., 24
(4), 1261
–1272
(2015). https://doi.org/10.1109/TIP.2015.2389619 IIPRE4 10577149 Google Scholar
T. Chan and W. Zhu,
“Level set based shape prior segmentation,”
in IEEE Comput. Soc. Conf. Comput. Vision and Pattern Recognit.,
1164
–1170
(2005). https://doi.org/10.1109/CVPR.2005.212 Google Scholar
P. Yan et al.,
“Automatic segmentation of highthroughput RNAi fluorescent cellular images,”
IEEE Trans. Inf. Theory. Biomed., 12
(1), 109
–117
(2008). https://doi.org/10.1109/TITB.2007.898006 Google Scholar
S. Ali and A. Madabhushi,
“An integrated region, boundary, shapebased active contour for multiple object overlap resolution in histological imagery,”
IEEE Trans. Med. Imaging, 31
(7), 1448
–1460
(2012). https://doi.org/10.1109/TMI.2012.2190089 ITMID4 02780062 Google Scholar
B. Parvin et al.,
“Iterative voting for inference of structural saliency and characterization of subcellular events,”
IEEE Trans. Image Process., 16
(3), 615
–623
(2007). https://doi.org/10.1109/TIP.2007.891154 IIPRE4 10577149 Google Scholar
Y. AlKofahi et al.,
“Improved automatic detection and segmentation of cell nuclei in histopathology images,”
IEEE Trans. Biomed. Eng., 57
(4), 841
–852
(2010). https://doi.org/10.1109/TBME.2009.2035102 IEBEAX 00189294 Google Scholar
X. Qi et al.,
“Robust segmentation of overlapping cells in histopathology specimens using parallel seed detection and repulsive level set,”
IEEE Trans. Biomed. Eng., 59
(3), 754
–765
(2012). https://doi.org/10.1109/TBME.2011.2179298 IEBEAX 00189294 Google Scholar
J. Wright et al.,
“Robust face recognition via sparse representation,”
IEEE Trans. Pattern Anal. Mach. Intell., 31
(2), 210
–227
(2009). https://doi.org/10.1109/TPAMI.2008.79 ITPIDJ 01628828 Google Scholar
J. Mairal, M. Elad and G. Sapiro,
“Sparse representation for color image restoration,”
IEEE Trans. Image Process., 17
(1), 53
–69
(2008). https://doi.org/10.1109/TIP.2007.911828 IIPRE4 10577149 Google Scholar
B. Cheng et al.,
“Learning with ${\mathcal{L}}_{1}$graph for image analysis,”
IEEE Trans. Image Process., 19
(4), 858
–866
(2010). https://doi.org/10.1109/TIP.2009.2038764 IIPRE4 10577149 Google Scholar
S. Zhang et al.,
“Sparse shape composition: a new framework for shape prior modeling,”
in IEEE Comput. Soc. Conf. Comput. Vision and Pattern Recognit.,
1025
–1032
(2011). https://doi.org/10.1109/CVPR.2011.5995322 Google Scholar
Q. Zhang and R. Pless,
“Segmenting multiple familiar objects under mutual occlusion,”
in Proc. Int. Conf. Image Process.,
197
–200
(2006). https://doi.org/10.1109/ICIP.2006.312454 Google Scholar
C. Le Guyader and L. A. Vese,
“Selfrepelling snakes for topologypreserving segmentation models,”
IEEE Trans. Image Process., 17
(5), 767
–779
(2008). https://doi.org/10.1109/TIP.2008.919951 IIPRE4 10577149 Google Scholar
A. Dufour et al.,
“Segmenting and tracking fluorescent cells in dynamic 3d microscopy with coupled active surfaces,”
IEEE Trans. Image Process., 14
(9), 1396
–1410
(2005). https://doi.org/10.1109/TIP.2005.852790 IIPRE4 10577149 Google Scholar
J. Kong et al.,
“Robust cell segmentation for histological images of glioblastoma,”
in Proc. IEEE 13th Int. Symp. Biomed. Imaging,
1041
–1045
(2016). https://doi.org/10.1109/ISBI.2016.7493444 Google Scholar
P. Zhang et al.,
“Automated level set segmentation of histopathologic cells with sparse shape prior support and dynamic occlusion constraint,”
in Proc. IEEE 14th Int. Symp. Biomed. Imaging,
718
–722
(2017). https://doi.org/10.1109/ISBI.2017.7950620 Google Scholar
C. Goodall,
“Procrustes methods in the statistical analysis of shape,”
J. R. Stat. Soc., 53
(1), 285
–339
(1991). https://doi.org/10.2307/2345744 Google Scholar
B. Mohar et al.,
“The Laplacian spectrum of graphs,”
Graph Theory, Combinatorics and Applications, 2 12 Springer, Berlin
(1991). Google Scholar
B. Mohar,
“Some applications of Laplace eigenvalues of graphs,”
Graph Symmetry, Springer, Dordrecht
(1997). Google Scholar
D. L. Donoho and M. Elad,
“Optimally sparse representation in general (nonorthogonal) dictionaries via l1 minimization,”
Proc. Natl. Acad. Sci. U. S. A., 100
(5), 2197
–2202
(2003). https://doi.org/10.1073/pnas.0437847100 Google Scholar
E. J. Candes, J. K. Romberg and T. Tao,
“Stable signal recovery from incomplete and inaccurate measurements,”
Commun. Pure Appl. Math., 59
(8), 1207
–1223
(2006). https://doi.org/10.1002/cpa.20124 CPMAMV 00103640 Google Scholar
E. J. Candes and T. Tao,
“Nearoptimal signal recovery from random projections: universal encoding strategies?,”
IEEE Trans. Inf. Theory, 52
(12), 5406
–5425
(2006). https://doi.org/10.1109/TIT.2006.885507 IETTAW 00189448 Google Scholar
D. L. Donoho,
“For most large underdetermined systems of linear equations the minimal l1norm solution is also the sparsest solution,”
Commun. Pure Appl. Math., 59
(6), 797
–829
(2006). https://doi.org/10.1002/cpa.20132 CPMAMV 00103640 Google Scholar
A. Ruifrok and D. A. Johnston,
“Quantification of histochemical staining by color deconvolution,”
Anal. Quant. Cytol. Histol., 23
(4), 291
–299
(2001). AQCHED 08846812 Google Scholar
L. Vincent,
“Morphological grayscale reconstruction in image analysis: applications and efficient algorithms,”
IEEE Trans. Image Process., 2
(2), 176
–201
(1993). https://doi.org/10.1109/83.217222 IIPRE4 10577149 Google Scholar
A. F. Frangi et al.,
“Multiscale vessel enhancement filtering,”
Lect. Notes Comput. Sci., 1496 130
–137
(1998). https://doi.org/10.1007/BFb0056181 LNCSD9 03029743 Google Scholar
C. Li et al.,
“Distance regularized level set evolution and its application to image segmentation,”
IEEE Trans. Image Process., 19
(12), 3243
–3254
(2010). https://doi.org/10.1109/TIP.2010.2069690 IIPRE4 10577149 Google Scholar
S. Kim et al.,
“A method for largescale l1regularized least squares problems with applications in signal processing and statistics,”
IEEE J. Sel. Top. Signal. Process., 1
(4), 606
–617
(2007). https://doi.org/10.1109/JSTSP.2007.910971 Google Scholar
H. Xu, C. Lu and M. Mandal,
“An efficient technique for nuclei segmentation based on ellipse descriptor analysis and improved seed detection algorithm,”
IEEE J. Bio. Health Info., 18
(5), 1729
–1741
(2014). https://doi.org/10.1109/JBHI.2013.2297030 Google Scholar
C. Jung and C. Kim,
“Segmenting clustered nuclei using hminima transformbased marker extraction and contour parameterization,”
IEEE Trans. Biomed. Eng., 57
(10), 2600
–2604
(2010). https://doi.org/10.1109/TBME.2010.2060336 IEBEAX 00189294 Google Scholar
F. Xing et al.,
“Automatic ki67 counting using robust cell detection and online dictionary learning,”
IEEE Trans. Biomed. Eng., 61
(3), 859
–870
(2014). https://doi.org/10.1109/TBME.2013.2291703 IEBEAX 00189294 Google Scholar
T. Fawcett,
“An introduction to roc analysis,”
Pattern. Recognit. Lett., 27
(8), 861
–874
(2006). https://doi.org/10.1016/j.patrec.2005.10.010 Google Scholar
D. P. Huttenlocher et al.,
“Comparing images using the Hausdorff distance,”
IEEE Trans. Pattern Anal. Mach. Intell., 15
(9), 850
–863
(1993). https://doi.org/10.1109/34.232073 ITPIDJ 01628828 Google Scholar
D. Comaniciu and M. Peter,
“Mean shift: a robust approach toward feature space analysis,”
IEEE Trans. Pattern Anal. Mach. Intell., 24
(5), 603
–619
(2002). https://doi.org/10.1109/34.1000236 ITPIDJ 01628828 Google Scholar
L. Grady and E. Schwartz,
“Isoperimetric graph partitioning for image segmentation,”
IEEE Trans. Pattern Anal. Mach. Intell., 28
(3), 469
–475
(2006). https://doi.org/10.1109/TPAMI.2006.57 ITPIDJ 01628828 Google Scholar
K. Zhang et al.,
“Active contours with selective local or global segmentation: a new formulation and level set method,”
Image. Vision Comput., 28
(4), 668
–676
(2010). https://doi.org/10.1016/j.imavis.2009.10.009 IVCODK 02628856 Google Scholar
H. Xu et al.,
“Automatic nuclear segmentation using multiscale radial line scanning with dynamic programming,”
IEEE Trans. Biomed. Eng.,
(99), 1–1
(2017). https://doi.org/10.1109/TBME.2017.2649485 IEBEAX 00189294 Google Scholar
O. Ronneberger, P. Fischer and T. Brox,
“Unet: convolutional networks for biomedical image segmentation,”
Lect. Notes Comput. Sci., 9351 234
–241
(2015). https://doi.org/10.1007/9783319245744 LNCSD9 03029743 Google Scholar
H. Vo et al.,
“Cloudbased whole slide image analysis using MapReduce,”
in Proc. VLDB Workshop Data Manange. Anal. Med. Healthcare,
62
–77
(2016). https://doi.org/10.1007/9783319577418_5 Google Scholar
G. Teodoro et al.,
“Algorithm sensitivity analysis and parameter tuning for tissue image segmentation pipelines,”
Bioinformatics, 33
(7), 1064
–1072
(2017). https://doi.org/10.1093/bioinformatics/btw749 BOINFP 13674803 Google Scholar
BiographyPengyue Zhang is a PhD candidate at the Department of Computer Science at Stony Brook University. His research interests include medical image analysis, computer vision, and machine learning. Fusheng Wang is an associate professor at the Department of Biomedical Informatics and Department of Computer Science at Stony Brook University. He received his PhD in computer science from the University of California, Los Angeles, in 2004. Prior to joining Stony Brook University, he was an assistant professor at Emory University. His research interest crosscuts data management and biomedical informatics. George Teodoro received his PhD in computer science from the Universidade Federal de Minas Gerais (UFMG), Brazil, in 2010. Currently, he is an assistant professor in the Computer Science Department at the University of Brasilia (UnB), Brazil. His primary areas of expertise include high performance runtime systems for efficient execution of biomedical and datamining applications on distributed heterogeneous environments. Yanhui Liang is works at Google Brain as a software engineer. She received her PhD in biomedical informatics from Stony Brook University. Her research interests include 2D/3D medical imaging, computer vision, machine learning, and largescale spatial data analytics. Mousumi Roy is a PhD student at Stony Brook University with research focus on computer vision, data analysis, and machine learning modeling for biomedical research. Daniel Brat received his MD and PhD degrees from Mayo Medical and Graduate Schools. He completed residency and a fellowship at Johns Hopkins Hospital. He is Magerstadt professor and chair of pathology at the Northwestern University Feinberg School of Medicine and pathologistinchief of Northwestern Memorial Hospital. He directs a basic and translational research lab that investigates mechanisms of glioma progression, including the contributions of hypoxia, genetics, tumor microenvironment, and stem cells. Jun Kong received a PhD in electrical and computer engineering from the Ohio State University. Currently, he is an associate professor in the Department of Mathematics and Statistics at Georgia State University. He is an affiliated faculty and a member of the Winship Cancer Institute in Emory university. His primary areas of expertise include biomedical image analysis algorithms, computeraided diagnosis systems, and largescale integrative analysis for quantitative biomedical research. 