A Novel Evolution Strategy of Level Set Method for the Segmentation of Overlapping Cervical Cells

: Development of an accurate and automated algorithm to completely segment cervical cells in Pap images is still one of the most challenging tasks. The main reasons are the presence of overlapping cells and the lack of guiding mechanism for the convergence of ill-deﬁned contours to the actual cytoplasm boundaries. In this paper, we propose a novel method to address these problems based on level set method (LSM). Firstly, we proposed a morphological scaling-based topology ﬁlter (MSTF) and derived a new mathematical toolbox about vector calculus for evolution of level set function (LSF). Secondly, we combine MSTF and the mathematical toolbox into a multifunctional ﬁltering algorithm 2D codimension two-object level set method (DCTLSM) to split touching cells. The DCTLSM can morphologically scale up and down the contour while keeping part of the contour points ﬁxed. Thirdly, we design a contour scanning strategy as the evolution method of LSF to segment overlapping cells. In this strategy, a cutting line can be detected by morphologically scaling the union LSF of the pairs of cells. Then, we used this cutting line to construct a velocity ﬁeld with an effective guiding mechanism for attracting and repelling LSF. The performance of the proposed algorithm was evaluated quantitatively and qualitatively on the ISBI-2014 dataset. The experimental results demonstrated that the proposed method is capable of fully segmenting cervical cells with superior segmentation accuracy compared with recent peer works.


Introduction
Cervical cancer, which is primarily caused by the infection of some types of Human Papilloma Virus (HPV), is one of the most common gynecological cancer in the world and has a very high probability of fatality if it is left untreated [1][2][3]. However, if the precursor lesions caused by HPV can be detected, it is quite possible to cure cervical cancer. Fortunately, cervical tumor has a precancerous condition that can be recognized and treated early, so cervical cancer has a great potential for prevention and cure. Over the past three decades, cervical cancer incidence and mortality rates have decreased by more than 50%, and most of the reduction is attributed to the Papanicolaou (Pap) smear screening technique introduced by [4] for detecting the precursor lesions of cervical cancer [5]. However, it still leaves much to be desired.
The Pap smear screening work in current clinical practice is mainly performed manually by cytotechnologists who must be critically trained and well-experienced. This manual process is very difficult and time consuming, and cytotechnologists are susceptible to make erroneous decisions [6] mainly caused by intra-and interobserver variability. Therefore, as a preparatory work of automated diagnosis, automated slide analysis has been a research focus. As the most essential task in automated slide analysis, the automated segmentation of cervical cells in the extended depth-of-field (EDF) images captured in different fields of view (FOVs) [7] is an ongoing concern [8][9][10]. Although some work has been done for cervical cell segmentation, such as segmentation of isolated cells (i.e., free-lying cells) and nuclei, it is still a challenge to develop an accurate and complete segmentation algorithm for overlapping cells in Pap images [10]. Due to the influence of poor contrast, overlapping, folding, polymorphic nonconvex shapes and the impurity introduced from the cell deposition process, accurate and complete segmentation of cytoplasm from cellular clumps in Pap images is the biggest obstacle for this challenge. Therefore, it is essential to develop a new system that can commendably carry out automated segmentation of overlapping cytoplasm for cervical cells in the Pap images.
The literature reviews about related studies on segmentation of cervical cells are in the following. Prior to the first attempt [11] to segment the overlapping cervical cells completely, the systems could be divided into three groups, including (1) segmentation of the nuclei [12][13][14][15][16][17][18] of cells, (2) segmentation of the isolated cells [19][20][21] and (3) segmentation of clusters of cells crowded together (i.e., cellular clumps) [22]. Then, a few complete segmentation systems for overlapping cervical cells emerged in recent years. In general, the main algorithms used in these systems can be roughly classified into the following algorithm groups: (1) graph-based methods [22][23][24][25][26], (2) morphological processing-based methods [27][28][29] and (3) level set method (LSM)-based approaches [30][31][32][33]. Graph-based methods, such as Voronoi diagram method [23], provide a new idea for cervical cells segmentation; however, the segmentations look unnatural. Morphological processingbased methods are simply to be used but rely on the parameter configuration. The LSMbased approach is more suitable for completely automated cervical cell segmentations, because it is less dependent on seeds initialization and parameter configuration. As an effective method for object delineation, LSM was first proposed by [34]. In the development of LSM, scientists contributed some classic works. For instance, after the original snake model was proposed in [35], a geometrically flexible and topologically changeable snake model was introduced independently [36,37] by partial differential equation (PDE) and has been successfully applied in two-and three-dimensional spaces. After Mumford [38] first proposed a functional framework for image segmentation based on the variation method and Zhao [39] first introduced the Heaviside and Delta functions to represent the length and area of the region delineated by level set function (LSF), Chan [40] proposed a successful region-based LSM known as the C-V model. Moreover, based on the previous works [38,41], Chen [42] first proposed an edge-based LSM with a prior shape. Li [43] proposed a new LSM (DRLSE) without periodic reinitialization LSF as signed distance function (SDF). DRLSE is used in the baseline method [33] to speed up the contour evolution. LSM can be used to detect the textural edges like Canny edge detector [44], because its adaptive topology changes can fit multiple target contours at the same time. The LSM can even control non-closed curves in three-dimensional space by means of the codimension twolevel set method (CTLSM), which was first carried out by [45] and developed by [46]. In addition, the LSM has been successfully applied in image segmentation tasks [25,[47][48][49]. Lu [33] proposed an algorithm for cervical cells segmentation based on multiple level set functions, which has been regarded as the baseline and comparative algorithm by a lot of literature [26,28,[50][51][52][53][54]. However, it has a problem that the poorly initialized level set function never converges to the cell boundaries.
With the rapid improvement of the graphic processing unit (GPU), CNNs (convolutional neural networks) have been widely used in the field of medical image segmentation. However, these methods mainly focus on clinical parameters, such as the volume and shape obtained from organ and structure segmentation for quantitative analysis and organ substructure obtained from lesion segmentation for histopathological diagnosis [55]. Af-ter the U-net was introduced by Ronneberger [56], U-net-like architecture convolutional neural networks have achieved good results in semantic segmentation tasks, such as the segmentations of organs, tissues or lesions areas in CT, X-ray and MR medical images. In contrast, cervical cell segmentation is more like an instance segmentation task that attempts to delineate all cells from Pap images to separated instances. Although CNN-based methodologies, such as Mask R-CNN [57], have good performances for instance segmentation, they are not suitable for the automated and complete cervical cells segmentation, because the overlapping regions in pairs of overlapping cells cannot be separated. In the latest proposed panoptic segmentation task, which combines semantic with instance segmentation, nonoverlapping is clearly declared in the seminal work [58]. Therefore, CNN-based segmentation methods are difficult to use to achieve end-to-end cervical cell segmentation. This means classical image processing strategies are essential for CNNassisted methods [22,59] that are successfully used in the segmentation of overlapping cervical cells. Moreover, the other obstacle for cervical cell segmentation by CNNs is a lack of datasets, and the CNN-based segmentation methods are particularly dependent on datasets.
In this paper, we propose a novel method to automatically segment the cytoplasm of touching and overlapping cells in Pap images. Our proposed method is capable of fully segmenting cervical cells compared with the method that can only segment isolated cells or segment regions of cellular clumps [60,61]. Our method is unsupervised, which provides the possibility for cervical cell segmentation with limited data, and is potentially more valuable in guiding automated slide analysis for cervical pathology. The main goal of the proposed method is to alleviate nonconvergence of the level set function (LSF) in delineating the overlapping regions of cells caused by poor initialization of the cellular contours in [33].
Our main contributions lie in the following aspects: (1) We designed a morphological scaling-based topology filter (MSTF) to filter out the false positive fragments caused by improper allocation of the initial contour points of touching cells (i.e., misallocation). In MSTF, we constructed the signed distance function (SDF) as the LSF of the initialized cytoplasm contour based on the linear time Euclidean distance transform (LTEDT) algorithm [62], which is denoted by LTEDT-SDF. (2) We theoretically derived a new mathematical toolbox about vector calculus for evolution of the LSF as the supplementary for the codimension two-level set method (CTLSM) [63], aiming to keep the initialized contours of the nonoverlapping region fixed. Our proposed evolution method of partially fixed contour is called the 2D codimension two-object level set method (DCTLSM), which can alleviate the accuracy loss of a MSTF. (3) We proposed a novel evolution strategy of LSF inspired by the watershed method [64].
In this strategy, we provided an effective guidance mechanism for attracting and repelling the LSF to converge towards its actual cell boundary. (4) We used the dataset published by the First Overlapping Cervical Cytology Image Segmentation Challenge held in conjunction with the IEEE International Symposium on Biomedical Imaging (ISBI-2014 challenge) to evaluate our proposed method. The experimental results showed that cellular clumps consisting of two to 10 cells under an overlap ratio less than 0.2 can be accurately segmented. Furthermore, the segmentation of cellular clumps consisting of two to four cells can be effectively segmented with an overlap ratio less than 0.5. By qualitive and quantitative comparisons, our method outperformed the other segmentation methods.

Methodology
As shown in Figure 1, the architecture of the proposed method consists of contour initialization, touching cell splitting (pink arrow) and overlapping cell segmentation (purple arrow). In Section 2.1, the contour initialization includes nuclei detection and segmentation, as well as the segmentation of the contours formed by the boundary of the cellular clumps. After contour initialization, we obtained the initial contour of cytoplasm for each cervical cell. Based on this initialization, we proposed a new method for the touching cell splitting described in Section 2.2 and overlapping cell segmentation described in Section 2.3. As a major task, the proposed touching cell splitting (pink arrows) can alleviate the effect of misallocation of the pixels on the clump boundary, which may cause the cells slight touching or the isolated ones misjudged as overlapping. After the isolation judgement of the isolated cells, the algorithm proceeds to the overlapping cell segmentation, which is the other major task, aiming to segment the overlapping cell from its neighbor cells iteratively. Figure 1. The architecture of the proposed method takes the contour initialization as the input to perform touching cell splitting (pink arrow) and overlapping cell segmentation (purple arrow). MSTF: morphological scaling-based topology filter and DCTLSM: 2D codimension two-object level set method. The " Fig." and "Eq." denote " Figure" and "Equation" respectively.

Cellular Component Segmentation
The goal of cellular component segmentation is to achieve contour initialization by nuclei detection and cellular clump segmentation. In this paper, we used the contour initialization method [33], because it provides an efficient architecture for cellular clumps segmentation, which can be easily achieved by using the high grayscale contrast between the cellular components and the background. As illustrated in Figure 2, the contour initialization method can be divided into the following steps. Firstly, a super-pixel map (see Figure 2b, which is supposed as the preparation for the generation of a clean edge map (see Figure 2c) can be acquired after a quick shift algorithm [65]. Then, a convex hull (see Figure 2d) is generated after a connected components analysis [66] based on the clean edge map. In addition, a Gaussian mixture model (GMM) is fitted iteratively to cluster pixels in the region of the convex hull into the segmentation of cellular clumps (see Figure 2e). Afterwards, the nuclei detection can be completed (see Figure 2f) by using the MSER algorithm [67] based on the low grayscale values and homogeneous texture of the nuclei within the clump region. At last, cytoplasm contour is initialized as a shape prior by the extrapolation of the allocation that clusters each pixel of the clump boundary to its owner nucleus based on the Euclidean distance, illustrated as Figure 2g. Lu [33] used this initial contour to construct the LSF and a distance map as the velocity terms for contour evolution. In this way, the coarse-to-fine segmentation was achieved after the minimization of the multiple energy functionals. However, this nucleus-based allocation caused the misallocation of the initial contour points, as illustrated in Figure 3a.

Touching Cell Spliting
Touching cell splitting serves two purposes. One is to alleviate the impact of the misallocation that exists among the touching cells and produces the false positive fragments. If the misallocation exists, it will cause the free-lying cells to be misidentified as overlapping ones (see Figure 3a). The other purpose is to make use of the accurate segmentation of cellular clumps, in which the high contrast between the background and the cellular clumps can be used to accurately delineate the contours of cellular clumps. If the misallocation can be alleviated and the accurate cellular component segmentation can be preserved, the performance of the final segmentation will be improved. In this section, to achieve these two purposes, we designed a morphological scaling-based topology filter (MSTF) and 2D codimension two-object level set method (DCTLSM) to filter out the false positive fragments.

Morphological Scaling-Based Topology Filter
Instead of the conventional crossing time [68] and fast marching method (FMM) [69], we used a linear time Euclidean distance transform (LTEDT) algorithm [62] based on dimensionality reduction and partial Voronoi diagram construction to calculate the L 2 norm distance, which was used to construct the signed distance function (SDF) as the LSF of the initial contour (see Figure 4a). The evolution of LSF φ can be represented by Equation (1): where v represents the velocity. The numerical discretization of (1) can be expressed by Equation (2): Given the velocity v, N isocontours with approximate shapes (see Figure 4d) can be obtained step-by-step for N times from the temporal iteration denoted by the superscript n in Equation (2). This means the SDF constructed by LTEDT (LTEDT-SDF) can be used to zoom in and out of the contour φ = 0 while preserving its morphological shape.
By setting v = −1, the contour (φ = 0) can keep shrinking until the topological connectivity of the region (contour) is broken down, illustrated as the green isocontour in Figure 4d. During the shrinking process, false positive fragments caused by misallocation can naturally disappear or be abandoned if the topological connectivity changes (see Figure 3b). We constructed a LTEDT-SDF for φ based on the shrunken contour (see Figure 4b), and the instance contour (φ = 0) can be expanded by setting v= 1 in Equation (2) (see Figure 4e). When the contour marches to the boundary of the target cell (see Figure 3c), a MSTF completes the topological filtering process. The false negative pixels appear in Figure 3d, because some points on the cytoplasm boundary are missing during the shrinking process of the MSTF, and they cannot be found during the boundary reconstruction by expansion. Therefore, we proposed the DCTLSM to solve this problem.

2D Codimension Two-Object Level Set Method
To reduce the false negative pixels caused by expansion in the MSTF, we used the accurate segmentation of cellular clumps to determine the fixed contour, shown as the green contour in Figure 3e. We defined the LSF φ 0 and φ 1 generated by LTEDT for the fixed and active contours. φ 0 can be associated with φ 1 to produce a union LSF φ, represented by Equation (3), which is often used to combine two LSF in traditional level set methods [36][37][38][39][40][41][42][43].
To make the fixed contour inactive and perform morphological scaling on the active contour, we proposed a novel evolution strategy DCTLSM incited by the idea of CTLSM [45]. Since it is not easy to derive the functional variation about φ based on Equation (3), we express Equation (3) as a new form, Equation (4): where H 0 is defined in Equation (5).
where H is the Heaviside function defined by Equation (6).
where H(φ) is the Heaviside function from Equation (6) is the Dirac delta function [63] acting on the value of φ at each two-dimensional point → x = (x, y . H 1 is defined by Equation (5), and δ 0 is calculated by applying the Delta function from Equation (7) on φ 0 − φ 1 . The magnitude of the gradient vector of φ 0 and φ 1 is denoted by the operator |·| in Equation (11). Based on this mathematical toolbox, the energy functionals of φ are designed to implement the DCTLSM, represented as Equations (12) and (13).
where M(φ) denotes the morphological term that is used as the attraction force by setting v = −1 in the shrinking process of the DCTLSM (see Figure 3f) and the repulsion force by setting v= 1 in the expanding process of the DCTLSM (see Figure 3g). E(φ) denotes the edge term, which is used in the refinement process of the DCTLSM (see Figure 3h). The Canny mountain map g is obtained by applying the LTEDT-SDF on the edge map generated by the Canny detector. By performing the Gâteaux derivative of φ 1 , we can obtain the first-order variations (14) and (15) of the energy functionals (12) and (13), according to the mathematical toolbox (8)- (11).
where δ (1) denotes the first-order variation. Then, we used −δ E(φ) as the gradient descent flow for M(φ) and E(φ) to update φ 1 , as shown in the level set Equations (16) and (17).
where t is the virtual time.
The shrinking process of the DCTLSM controlled by Equation (16) is effective to retain the contour fragment of φ 0 that is located at the cellular boundary (contour fragment 1 in Figure 3f) and the morphological shape of the active contour φ 1 = 0 (contour fragment 2 in Figure 3f). Meanwhile, the contour fragment of φ 0 that is located away from the cellular boundary becomes distorted through the corrosion of the additional term δ 0 φ 1 in Equation (14), aiming to conform to the evolution of φ 1 , illustrated in Figure 4f. As a result, the fixed and active contours are naturally merged into a union contour φ = 0. Similar to the MSTF, the misallocation part of the initial contour (contour fragment 3) can be removed after the topological connectivity changes, as illustrated in Figure 3f. The topology filtering process is completely done after the expansion process of the DCTLSM, which is similar to the shrinking process. In Figure 3g, the transition region (contour fragment 2) between contour fragments 1 and 3 is naturally generated by the expanding process. Contour fragment 4, represented as a blob located at the inner region of the contour, will disappear naturally. Similar to the points missing problem in the MSTF, some points on active contour φ 1 are missing due to the morphological scaling effect of the DCTLSM. As a compensation for this case, we used the Canny edge detector to obtain the textured edge features of the cell image and applied the LTEDT-SDF on these features to generate a Canny mountain map, which is rich in edge information including the outline of the cellular boundary. As shown in Figure 4c, the contour points on the LTEDT-SDF slope will fall into the nearby Canny valley at unit velocity, which is realized by Equation (17) without considering the curvature term about φ 1 . In addition, the Canny valley is the textural edges that are detected by the Canny edge detector. Ultimately, we obtained the postprocessed contour as the initialization of the cell instances, shown as the blue contour in Figure 3h.

Overlapping Cell Segmentation
To perform overlapping cell segmentation, we took the union contour as the reinitialized contour for each overlapping cell and used the LSM to evolve it to the real boundary. Specifically, we first constructed a repellent, called a cutting line (see Figure 5d), based on which, the attractor is generated by applying LTEDT on the union contour. Then, we designed an evolution strategy, called contour scanning, to segment the overlapping region of the cytoplasm based on the evolution velocity field generated by merging the repellant and attractor into a distance map. The cutting line can ensure the contour does not seriously collapse inward into the cytoplasm. In addition, the contour scanning strategy can retain as much edge information as possible (see Figure 5e) to find the ground truth contour of the cytoplasm.

Cutting Line Detection
A union contour (red in Figure 5a) is generated by the combination of the cell contour processed by the DCTLSM (blue in Figure 5a) and the initial contour (green in Figure 5a) of its overlapping neighbor. Then, the MSTF is applied on this union contour to find its topological connectivity changes based on the LTEDT-SDF D u of the union contour (see Figure 6a). We constructed the LTEDT-SDF 1 and LTEDT-SDF 2 for the disconnected contours (green contours in Figure 5b), respectively, and obtained a criterion called the topology discriminant bridge (see Figure 6b) on which the value of the LTEDT-SDF 1 and LTEDT-SDF 2 are equal. The intersection of the topology discriminant bridge and the black contour (see Figure 5b), which is considered as the critical state before its topological connectivity changes, is represented as the discriminant point x d (see Figure 5c). Furthermore, we generated the LTEDT-SDF D d of this discriminant point, as shown in Figure 6c, to find the project points x p (square markers) on the union contour. Finally, the cutting line can be generated by the project points, shown as the blue line in Figure 5d. The project points can be easily found at the intersection of the union contour and the discriminant contour shown as the red circle in Figure 6c. Specifically, the discriminant contour is the isocontour with a certain distance D d ( x p ), where the project points x p can be calculated by.
The threshold T is used as the parameter of numerical smearing for the calculation error, which is set to 1 (pixel). This numerical recipe is based on the inspiration that D d ( x p ) and D u ( x d ) are considered to be theoretically equal. That is because the critical state contour (see the black contour in Figure 5b) will contain the mapping pixel of the project points if the union contour can be separated into different parts (see the green contour in Figure 5b) during the process of morphological shrinking.

Contour Scanning Strategy for Segmentation
Inspired by the watershed method, we used its reverse process in contour scanning. It is assumed that there is a lowest imaginary source inside the region of the cell to be segmented. Then, we considered this source as a drain (see Figure 6d), so that the process of watershed swelling became draining. The cutting line was considered as the judgment condition for closing the draining hole. When the water level dropped from the height of the union contour to the height of the dam constructed by the cutting line, the textured edges inside the neighbor cell were left, including the continuous contour associated with the cell to be segmented (see Figure 5e). To realize this draining process, there was no need to find this imaginary source mentioned above, because we used the level set method to directionally control the evolution of the contour. Specifically, the velocity field for the evolution of the union contour can be concluded by the following steps. First of all, the external region of the union contour in the velocity field was set at −1. Then, we set 0 to the internal region bounded by the cutting line and the cell to be segmented. Finally, we assigned the reversed LTEDT-SDF that was constructed by the cutting line to the region bounded by the neighbor cell and cutting line. After the construction of the velocity field, we controlled the evolution of the union contour by using Equation (2), where v is calculated by the velocity field illustrated in Figure 6d. The contour of the neighbor cell will be distorted and pushed towards the cutting line along the scanning vector (see Figure 6d). In this way, the contour scanning can be realized. The scanning vector with |∇v|= 1 is naturally produced by LTEDT-SDT. After contour scanning, the textural edges will be left as shown in Figure 5e. Furthermore, the internal part and the external part of the contour to be segmented can naturally form a connecting contour, and the inner textural edges of the neighbor cell are scattered into fragments. These fragments can be removed in the shrinking process of the DCTLSM (see Figure 5f). After the expanding process of the DCTLSM, we can obtain the coarse segmentation of the cytoplasm shown as Figure 5g. Then, the refined segmentation (see Figure 5h) of the cytoplasm will be obtained based on the Canny mountain map and evolution Equation (17). The segmentation based on this scanning strategy will be performed iteratively among the neighbor cells.

Image Datasets
To evaluate the performance, we used the ISBI-2014 dataset [71], which is provided as the first cervical cytology image segmentation challenge. This dataset was provided in the competition stage and divided into 45 synthetic training images (train-45 dataset) and 90 synthetic test images (test-90 dataset). The images in the datasets were synthesized from the 16 EDF images, each of which was generated by applying the discrete wavelet transform method [72] on the according 16 FOV image volumes that were acquired under an Olympus BX40 microscope system. The boundaries of the nuclei and cytoplasm in these 16 EDF images were manually annotated by an experienced cytotechnologist. The organizer divided the 16 annotated EDF images into 4 training EDF images and 12 test EDF images from which the 45 synthetic training and 90 synthetic test images (512 × 512 pixels, gray-scaled) were constructed separately. Specifically, after the manual selection of 12 isolated cells in the 4 training EDF images and 41 ones in the 12 test EDF images, an automatic synthesis procedure was applied to construct synthetic training and test image datasets with the following regulations: (a) the image patches, which contain the isolated cells, were clipped from the EDF images and then pasted into one 512 × 512 synthetic image after applying random rigid geometric transformation and random linear brightness transformation.

Evaluation Metrics
Evaluation metrics in the ISBI-2014 challenge [71] were introduced to assess the quantitative performance of our segmentation method for the "good" segmentation, including the pixel-based average Dice coefficient (DC), pixel-based average false positive rate (FP P ) and pixel-based average true positive rate (TP P ). In addition, we assessed not so "good" segmentation by using the object-based false negative rate (FN O ). DC is defined as where S SR denotes the segmentation result, S GT denotes the ground truth and |·| denotes the number of pixels in the object within this symbol. The "good" segmentation mentioned above is considered as those whose DC is greater than the thresholds ranging {0.6, 0.7, 0.8, 0.9}.
In addition, we introduced two more metrics to assess the effectiveness of the MSTF and DCTLSM, namely the pixel-based average false positive growth rate (∆FP) and pixelbased average false negative growth rate (∆FN). The FP P can not only represent the segmentation precision but also indicate the specificity of the morphology filtering effect, which occurs in the deletion of the fragments made up of the misallocation pixels for touching cells. In order to show the effect of touching cell splitting, we defined the growth rate of FP P as ∆FP. (18) where FP ini denotes the pixel-based average false positive rate of the initialization region described in Section 2.1 (see Figure 2g), and FP split denotes the pixel-based average false positive rate of the MSTF or DCTLSM in touching cell splitting. FN P was used to measure the loss of recall rate caused by the morphological scaling process. We defined the growth rate of FN P as ∆FN. (19) where FN ini denotes the pixel-based average false negative rate of the initialization region, and FN split denotes the pixel-based average false negative rate of the MSTF or DCTLSM in touching cell splitting.

The Determination of Morphological Scaling Threshold for MSTF and DCTLSM
To filter the false positive fragments of the touching cells caused by misallocation in contour initialization, it is essential to determine the suitable morphological scaling threshold TH. If TH is too small, the touching cells cannot be split. If TH is too large, overlapping cells will be segmented, which will stop the overlapping cells from being segmented based on the contour scanning strategy. We applied the MSTF on the cell images of the train-45 dataset with the overlap rate of [0,0.1], because touching cells exists in those whose overlap rate is between [0,0.1]. The results of ∆FP and ∆FN with different TH are shown in Table 1. The specificity of the morphology filtering effect manifests itself as the extremum of ∆FP, and we adopted the maximum principle to ensure the performance of touching cell splitting as much as possible. Therefore, we consider the transition of ∆FP from the extremum to an outlier value as the marker of threshold selection.
To supply a more visual description of the results, we show the changes of the ∆FP and ∆FN in Figure 7b. Compared with the MSTF, the ∆FN is obviously decreased by the DCTLSM, and there is almost no change in the ∆FP.
By performing outlier detection on the data in Table 1, we get a box diagram (Figure 7a) about the ∆FP. In the figure, the extremum occurs at TH = 5, and the transition from the extremum to the outlier value occurs at TH = 15. As the threshold increases, the outliers disappear at TH = 30, which indicates that the overlapping cells begin being segmented, and the ∆FP will not change after TH > 50, which indicates that the cells are all segmented. After the outlier detection and analysis, we set TH = 15 as the morphological scaling threshold. However, the FN P is increased by morphological scaling. To alleviate this problem, we apply the DCTLSM on the same images used by the MSTF to get the ∆FP and ∆FN at TH = 15, as shown in the last row of Table 1. The DCTLSM can reduce the ∆FN effectively compared with the MSTF.
To further verify the performance of the DCTLSM, we also performed comparative experiments on the test-90 dataset. As shown in Table 2, compared with the MSTF, the DCTLSM decreased the ∆FN obviously, with almost no change in the ∆FP. The ∆FP and ∆FN are equal to 0 when the number of cells is three or five. That is because there are no overlapping and touching cells in the test-90 dataset.
In summary, the experimental results demonstrate that the performance is improved by the DCTLSM compared with the MSTF in touching cell splitting when TH = 15.

Quantitative Comparison with Baseline Method
The quantitative evaluation of 810 cells is performed over the "good" cytoplasm segmentation with DC thresholds ranging {0.6,0.7,0.8,0.9}. The comparison between our proposed method (TH = 15) and the baseline method [33] is displayed in Table 3. To keep fair, they are based on the same initialization procedure. The experimental results show that, as the threshold increases, DC, TP P and FN O increase and FP P decreases, which indicates that a trade-off should be made between cell detection and the accuracy of segmentation for real-world missions. According to Table 3, we further visualized the advantages of the proposed method in Figure 6. Our segmentation results show a stable (i.e., lower variances) performance on both the train-45 dataset (see Figure 8a) and test-90 dataset (see Figure 8b) in the aspect of the improvement of the segmentation accuracy, namely DC, TP P and FP P . Only when DC > 0.8, the FP P of the proposed method on the test-90 dataset is slightly larger. As shown in Figure 6, our segmentation has a higher object-based average recall rate, because the FN O is obviously lower than [33] when the DC threshold is less than 0.9. In conclusion, our proposed method outperforms [33].
To further assess the performance of our proposed method, we graphically visualized DC, TP P , FP P and FN O as binary functions of the overlap ratio and the number of cells on the test-90 dataset with a "good" segmentation threshold of DC > 0.7 based on TH = 15, as shown in Figure 9. These functions show that the cellular clumps consisting of two to 10 cells under the overlap ratio less than 0.2 can be stably and accurately segmented by our proposed method, as illustrated by Figure 9a,d. Furthermore, the clumps consisting of two to four cells can be effectively segmented. FP P is small enough over the whole domain, even approaching zero as the overlap ration decreases (see Figure 9c). TP P is larger than 0.94 over most of the domain, as illustrated in Figure 9b. In addition, our proposed method has a more uniform performance, because the superior value of these functions is more invasive, distributing into the region of a larger number of cells and higher overlap ratio than the hierarchically exponential distribution in [33]. Therefore, our proposed method is more robust to the number of cells and overlap ratio.    Tables 4 and 5, including the methods of the ISBI-2014 challenge winner, Nosrati et al. [30], as well as their newly proposed method [31] and the baseline method [33]. In these two tables, the superior results are indicated by the bold numbers. Specifically, our proposed method achieves the best performance on the segmentation precision indicated by the DC value that achieves the highest value (90.4%) with almost up to 2% improvement over the average values (88.6%) of the other methods on the test-90 dataset and highest value (92.3%) with almost up to 3% improvement over the average value (89.7%) on the train-45 dataset. Moreover, our proposed method has the best object-based detection performance on the train-45 dataset, indicated by the lowest FN O , which achieved 10.4% with up to an almost 7% reduction over the average value of 17.3%. In addition, our FP P is the lowest on the test-90 dataset, which further shows that our proposed method has the highest segmentation precision over the other methods. Although the method [26] obtained a lower FP P on the train-45 dataset, it is not an accurate method for delineation of the cytoplasm boundaries due to the absence of outstanding DC and FP P . Similarly, in the absence of outstanding DC and FP P , a much higher TP P alone as obtained by Tareef [28,60] on both of the train-45 and test-90 datasets cannot indicate a promising performance on the segmentation precision. Huang [51] obtained the lowest FN O on the test-90 dataset, but their segmentation precision had no advantage, for the reason reported in [21].
In summary, our proposed method achieves the highest DC, lower FP P and reasonable TP P compared with the other methods. What is more, it can avoid extensive parameter tuning, such as [30,[32][33][34][35][36][37][50][51][52], which relied on manual selection in an arbitrary parameter set for multiple-level set functions. The methods used in this comparison, except for the other winner Ushizima [27] in the ISBI-2014 challenge, were all designed to delineate the cervical cells in Pap images, aiming to segment realistic cytoplasm boundaries. The authors of [27] obtained unrealistic geometric segmentation by the edges of a Voronoi diagram partitioning the overlapped pairs of cells with no respect of the realistic searching for the actual cell boundaries [73]. Although FP P in [27] was the lowest, we still made the second-lowest results (bold number in Table 5). This is not only because of the unrealistic geometric segmentation [27] described above but also that we consider the lower FP P alone cannot reflect the performance on segmentation accuracy, since the FP P would be zero if each nucleus was selected as the segmentation for the cytoplasm of each cell, supposing each nucleus was lying inside its cytoplasm.

Computational Complexity
We made an average running time comparison between our proposed algorithm and [33] on the test-90 dataset based on the hardware platform with a 2.30 GHz Intel Core i7-4712MQ processor and 12GB DDR3 RAM in the form of an unoptimized MATLAB (The MathWorks, Inc., Natick, MA, USA) code. Figure 10a shows the running time of the initialization described in Section 2.1, which takes 0.14-1.5 min, with an average of 0.79 (std: 0.37) minutes. Figure 10b,c show the running time of [33] and us, respectively, in segmenting cytoplasm contours from the background based on the initialization work. Our proposed algorithm is 2-26 times faster than [33], with an average of 5.4 (std: 4.0) times acceleration for the images on the test-90 dataset. Moreover, when the overlap ratio ranges

Qualitative Evaluation of Our Segmentation Results
Four typical segmentations compared with the ground truth and [33] are shown in Figure 11. The improvement can be listed in two aspects: the segmentation for the touching cells and the boundary convergence of the overlapping cells. The pseudo-fragments of the touching cells caused by the misallocation in the initialization of the cellular boundaries can be obviously filtered out by our proposed touching cell splitting algorithm, as shown in the first column of Figure 11. Over-segmentation caused by the absence of an effective guidance mechanism for LSF evolution in [33] is improved by our proposed contour scanning strategy, as shown in the second column of Figure 11. Similarly, the improvement for under-segmentation is shown in the third column of Figure 11. The orange contour and the dark blue contours in the fourth column of Figure 11 further illustrate the improvements by our method. Figure 11. Qualitative results of (a) Lu (reproduced from [33], Institute of Electrical and Electronics Engineers, 2015) and (b) our proposed method compared with (c) the ground truth.

Conclusions
In this paper, we proposed a novel algorithm based on the level set method (LSM) to automatically segment the touching and overlapping cells. The main goal of our proposed method is to address the nonconvergence problem of the baseline method in delineating the overlapping regions of cells, which is caused by ill-defined cytoplasm boundaries in the initialization. Our proposed method exhibits the following contents: (1) we proposed a morphological scaling-based topology filter (MSTF) and a new mathematical toolbox about vector calculus for the evolution of the level set function (LSF). (2) We combined the MSTF and the mathematical toolbox into a multifunctional filtering algorithm DCTLSM that can morphologically scale up and down the contour while keeping part of the contour points fixed. (3) We designed a novel contour scanning strategy for overlapping cell segmentation. The proposed algorithm is capable of fully segmenting cervical cells from the clumps of up to 10 cells with more accurate segmentation. Our segmentation provides an obviously lower false negative rate compared with the baseline method based on the same ill-defined cytoplasm boundaries. By quantitative and qualitive comparison on the ISBI-2014 dataset, our proposed method is superior to other state-of-the-art works.
In the future, our proposed method will be of increasing value in guiding automated slide analyses for cervical pathology after being intertwined with other promising techniques.
Author Contributions: G.L. contributed to this work by setting up the experimental environment, designing the algorithms, designing and performing the experiments, analyzing the data and writing the paper. Q.D., H.L., M.J., T.J., M.H. and G.D. contributed through research supervision and reviewer roles and by amending the paper. All authors have read and agreed to the published version of the manuscript.