A New Approach for Cell Detection and Tracking

In recent years, cell tracking methods by detection have become more and more popular because they outperformed cell tracking methods by contour evolution in most practical cell tracking applications. Yet, the most frequently used segmentation technique by cell detection methods is still threshold selection that is determined manually or by algorithms proposed in the 1970s. As a whole, these old threshold selection methods could not meet the accuracy requirement of cell detection adequately. In this paper, we propose a new approach of cell tracking by detection based on a multiple-threshold segmentation method that calculates multiple thresholds automatically and robustly. After cell detection, the proposed approach generates the timeline moving trajectory of a cell by connecting the cell positions along the time lapse image sequences based on morphological operations. We use four types of cells to verify the effectiveness of the proposed approach and the experimental results are favorable.


I. INTRODUCTION
Due to the advances in optical imaging and image processing, image cytometry has made great progress in recent years.As its central problem in many studies, cell segmentation faces the following world-widely recognized challenges.(1), the great variety of cells takes on various shapes, sizes, intensity distributions and edges.(2), the contrast between the cell and the background might be lower than the contrast between cells.The intensities of different cells in the same image might vary greatly and different parts of the same cell might also vary greatly.(3), different cells might overlap or connect each other.(4), the density of the cells changes greatly with the emergence of new cells and the disappearance of old cells.The first challenge requires the cell segmentation method to incorporate the priori knowledge.The second challenge requires robust image segmentation methods and threshold selection method has been recognized as the best choice.The third challenge requires robust morphological image processing methods to separate the connected or overlapped cells.The fourth challenge requires robust tracking methods to link the positions of the same cell correctly.It is obviously out of the capability of the existing methods to solve them as a whole.In this paper, we try to meet the generality of a great variety of cells and propose a new approach of cell detection and tracking.
Cell segmentation has achieved great advances in recent years that include the abundant emergence of new methods, the increased automation level and the increased accuracy.However, most state-of-the-art methods or software tools mainly rely on basic image processing algorithms developed many decades ago.For example, the watershed [1] algorithm developed in the 1970s is the one of the core image processing algorithms adopted by the Cellsegm method [2], the SMASH method [3], the ImageJ [4] and the CellProfiler [5].The Otsu's threshold selection algorithm [6] developed in the 1970s is one of the core image processing algorithms adopted by the SMASH method and the ImageJ.These basic image processing algorithms have been developed to solve the general computer vision problems and their abilities to address the challenging cell segmentation problem are limited by their inherent insufficiencies.For instance, the watershed algorithm could not overcome the severe over-segmentation problem and consequently it could not solve the third challenge well.The Otsu's algorithm is easily affected by low contrast and varied intensities and consequently it could not solve the second challenge well.In recent years, many new cell segmentation methods have been proposed and they have both advantages and disadvantages [7][8][9].For instance, the iterative erosion method proposed recently solves the third challenge well and it overcomes the oversegmentation problem by terminating the overlapping cell separation process with the area threshold [7].The slope difference distribution based threshold selection method proposed recently solves the problem caused by low contrast well [7].However, it could not solve the problem of varied intensities well with one single intensity threshold.Therefore, we propose a multiple-thresholds-based cell segmentation method in this paper.The number of multiple thresholds and the ways to determine the thresholds could be adjusted according to the type of the cell being segmented.
After cell segmentation, the trajectories of different cells should be generated for biological and medical analysis [10].There are two groups of cell tracking methods: (1) tracking by contour evolution [11][12] and (2) tracking by detection [13].Overall, tracking by detection methods perform better than tracking by contour evolution methods for cell tracking [14].Tracking by detection methods completely rely on the accuracy of cell segmentation.Thus, cell segmentation is most fundamental and important and it is also most challenging [15].However, the efforts put into cell segmentation by the researchers from all over the world make the methods of cell segmentation more and more divergent [16].Even for a given dataset, there is no simple way to point out the right algorithm after the objective comparisons [14], which indicates that neither of the evaluated algorithms could solve the presented problems completely from a biologist's view.In the point of view of an algorithm developer, it is not difficult to find the optimal solution for one specific type of cell dataset.What is much more challenging and makes the problem much more complex is to find the optimal solutions for so many types of cell datasets.In this paper, we propose a new approach to detect and track different types of cell datasets.Experimental results on four types of cells indicate that the proposed approach is promising.

II. The Proposed Cell Detection Method
Suppose the positions of cells in the image sequences determined automatically by the cell segmentation method is ( , ) where 1 T is the optimal intensity threshold calculated by slope difference distribution (SDD) [8]. 2 T is the low intensity threshold to segment the cells with lower intensities and 3 T is the high intensity threshold to segment the cells with higher intensities.4 T is the area threshold of the cell to remove noise blobs.5 T is the area threshold for the iterative erosion to avoid over-segmentation.6 T is the area threshold for the iterative erosion to remove noise blobs.Please note that Eq. ( 1) is only used to formulate cell segmentation.The ground truths are the actual cells' positions instead of the manually identified positions and the proposed method does not require any training sets to calculate the parameters.In addition, the proposed segmentation method aims to identify the positions of cells solely.The parameters 1 2 3 4 5 6 , , , , , T T T T T T are calculated from the current frame of the image sequences as follows.
After the slope difference distribution () Sxof the current frame is computed, the optimal intensity threshold 1 T of the cell image is calculated as follows.
where 1 g and 2 g are two positions in the slope difference distribution that correspond to the largest peak 1 P and the second largest peak 2 P .They are calculated as:   T is computed as: The low intensity threshold 2 T is calculated as: The max area distribution of the largest blobs in all the segmentations by the thresholds higher than 1 T is computed as: The area difference distribution of the largest blobs in all the segmentations by the thresholds higher than 1 T is computed as: The high intensity threshold 3 T is calculated as: With the optimal intensity threshold 1 T , the cell image is segmented as: T  , the areas of all the segmented blobs are computed as: The segmented blobs are determined as cells or not based on their areas.The blobs whose areas are smaller than the mean of the areas of all the segmented blobs are determined as cells.
{ }( ) where { }( ) A x i denotes the area of the th i blob in the segmentation result by the global threshold x .
x N denotes the total number of blobs in the segmentation result by the global threshold x .
x A denotes the mean of all the segmented blobs' areas computed by

 
Areas ( ) Sx and it is computed as: The final segmentation result with the low threshold 2 T is computed as: With the high threshold 3 T , the cell image is segmented as: 33 () In the segmented images, 2 S contains cells with lower intensities and 3 S contains cells with higher intensities.1 S contains most cells that might miss the cells with lower intensities.Different cells with higher intensities might be segmented as one cell in 1 S because the intensities surrounding these cells might be higher than the optimal threshold 1 T .Therefore, three thresholds are used to segment the cells with different intensities.Accordingly, the three segmentation results are fused as follows.
The cells excluding higher intensities are computed as: 4 S is then filtered by the noise blob removing filter [8]   to generate the segmentation ' 4

S
for the cells with middle intensities.The fusion segmentation is computed as: 4 T is the area threshold that distinguishes the cell and the noise blobs in the fusion segmentation and it is determined based on the priori knowledge of the specific type of cells to be segmented.In the fusion segmentation result, the blobs whose areas are greater than 4 T are determined as cells.The final segmentation result is computed as: () where is the total number of blobs in the fusion result.

5
T is the area threshold for the iterative erosion to avoid over-segmentation and it is computed as: denotes the total number of blobs in the final segmentation result.w is a weighting coefficient and its default value is 0.5.With the area threshold 5 T , the seeds seed S are computed with the iterative erosion [8].

6
T is the area threshold to remove fake seeds and it is determined based on the priori knowledge of the specific type of cells to be segmented.With the area threshold 6 T , the fake seeds are removed as follows: ()

 
Areas where s N denotes the total number of seeds.The positions auto P of cells are then computed from the true seeds true S .Fig. 1 (a) shows the flowchart of how to compute the low intensity threshold 2 T .Firstly, the optimal threshold 1 T of the current frame in the time-elapse image sequences is calculated by Eqs (2)(3)(4).Then, the current image is segmented by Eq. ( 5) In each segmentation result, the areas of all the segmented blobs are computed and the max area is computed by Eq. ( 6).In total, 1 T max areas are selected and form the max area distribution.The differences between adjacent max areas in the max area distribution are calculated and form the difference distribution denoted by Eq. ( 7).The value where the peak of the difference distribution is calculated and subtracted from 1 T as shown in Eq. ( 8) to get the low threshold T .Firstly, the optimal threshold 1 T of the current frame in the time-elapse image sequences is calculated by Eqs (2-4).Then, the current image is segmented by Eq. ( 5) with the thresholds 1 Tx  where 1, 2,..., 255 xT  .In each segmentation result, the areas of all the segmented blobs are computed and the max area is computed by Eq. ( 9).All max areas are selected and form the max area distribution.The differences between adjacent max areas in the max area distribution are calculated and form the difference distribution denoted by Eq. (10).The value where the peak of the difference distribution is calculated and added to 1 T as shown in Eq. ( 11) to get the high threshold 3 T .In Fig. 2 (a), we demonstrate how the low intensity threshold 2 T based on the optimal threshold 1 T .The max area distribution of the largest blobs is denoted by the blue line.The difference distribution is denoted by the red lines.As can be seen, there is a peak in the difference distribution.The x value corresponds to this peak is denoted as  T .When the intensities of the cells vary less significantly, equation 1 is shortened with fewer intensity thresholds.For example, two intensity thresholds are used to segment the pancreatic stem cells as formulated below.
Accordingly, the fusion segmentation is computed as: In addition, the high intensity threshold 3 T are calculated directly from the slope difference distribution instead of Eq. (11).Fig. 5 shows the computed cell positions at different times overlaying on the original cell images.As can be seen, the segmentation results by two intensity thresholds are accurate.
On the other hand, when the intensities of the cells are extremely low, the high intensity threshold becomes useless.Multiple low intensity thresholds are needed to segment the cells accurately.For example, the optimal intensity threshold and two low intensity thresholds are used to segment the mouse stem cells as formulated below.T is computed by Eq. ( 8) and 2 2 T is computed as: 1 TT  (31) Accordingly, the fusion segmentation is computed as: On the other hand, there must be cases where six parameters 1 2 3 4 5 6 , , , , , T T T T T T are too few to make Eq. ( 1 minimum enough and more thresholds are required.In such cases, slope difference distribution threshold selection method could be used to calculate multiple thresholds simultaneously and automatically [7].where B is the structuring element and B is its reflection.The structuring element is a unit sphere. The morphological dilation (Eq.( 35)) is repeated to connect the adjacent cell positions along the time axis and generate a connected cell trajectory.To make the cell tracking method tolerate of errors caused by imperfect cell segmentation, a trajectory filter is designed to remove small and invalid trajectories.
Immediately following the morphological dilation (Eq.( 35)), all the trajectories in the three-dimensional image The area threshold to separate the small trajectory and the large trajectory is calculated as: The coordinates of the small and invalid trajectories are computed as: The three-dimensional image ' 3d I is updated by the following equations to remove the small and invalid trajectories.

   
' 3 ( , , ) 0, , , , Because of the small and invalid trajectory removing process, the proposed cell tracking method could tolerate occasional errors caused by imperfect cell segmentation.After small and invalid trajectories are removed, the remaining large trajectories are the valid cell trajectories.The above adjacent cell position connecting process (Eqs.33-39) is repeated to connect adjacent cell trajectories and generate larger cell trajectories.

A. Results
We use 91 frames of Fluo-N3DH-CHO cells, 92 frames of Fluo-N2DL-Hela cells, 390 frames of pancreatic stem cells (PhC-C2DL-PSC) and 91 frames of mouse stem cells (Fluo-N2DH-GOWT1) to verify the effectiveness of the proposed cell detection and tracking approach.All the tested cell image sequences are available and could be downloaded at http://celltrackingchallenge.net/2ddatasets/.The identified positions of the cells are overlaid on the original images to evaluate the segmentation accuracy qualitatively.Fig. 7 shows some typical detection results for the Fluo-N3DH-CHO cells.Fig. 8 shows some typical detection results for the Fluo-N2DL-Hela cells.Fig. 9 shows some typical detection results for the PhC-C2DL-PSC cells.Fig. 10 shows some typical detection results for the Fluo-N2DH-GOWT1 cells.As can be seen, the proposed cell segmentation method could segment all these four types of cells robustly.Fig. 11 shows the generated cell trajectories for the 91 frames of Fluo-N3DH-CHO cells with different times of iterations of the proposed cell tracking method.Fig. 11 (a) shows the generated 12 cell trajectories by 3 iterations of the proposed cell tracking method.Fig. 11 (b) shows the generated 10 cell trajectories by 4 iterations.Fig. 11 (c) shows the generated 10 cell trajectories by 5 iterations and Fig. 11 (d) shows the generated 9 cell trajectories by 6 iterations.The generated trajectories with 6 iterations correlate with the practical cell trajectories most correctly, thus 6 iterations should be used for Fluo-N3DH-CHO cells tracking.The computation time in MATLAB for the proposed tracking methods with different iterations are shown in Table 1.Fig. 12 shows the generated cell trajectories for the 92 frames of Fluo-N2DL-Hela cells with different times of iterations of the proposed cell tracking method.Fig. 12 (a) shows the generated 174 cell trajectories by six iterations of the proposed cell tracking method.Fig. 12 (b) shows the generated 139 cell trajectories by eight iterations.Fig. 12 (c) shows the generated 138 cell trajectories by nine iterations and Fig. 12    Fig. 13 shows the generated cell trajectories for the 390 frames of pancreatic stem cells with different times of iterations.Fig. 13 (a) shows the generated 86 cell trajectories by one iteration of the proposed cell tracking algorithm.Fig. 13 (b) shows the generated 56 cell trajectories by two iterations of the proposed cell tracking algorithm.Fig. 13 (c) shows the generated 20 cell trajectories by three iterations of the proposed cell tracking algorithm.Fig. 13 (d) shows the generated 14 cell trajectories by four iterations of the proposed cell tracking algorithm.The generated trajectories with four iterations correlate with the practical cell trajectories well.From the green trajectory, the red trajectory and the blue trajectory, we could see cell division continuously gives rise to new cells where the daughter cells are linked to the mother cell and finally form a big tree.From the magenta trajectory, we see that a new cell moves in the view, divide into two cells and then moves out of the view.From the cyan trajectory, we see that a cell divide into two cells and one of the daughter cell moves out of the view.The computation time in MATLAB for the proposed tracking algorithm with different iterations are shown in Table 3.It is seen that it is very efficient for the proposed tracking method to track cells in low-resolution images (200×200).
Fig. 14 shows the generated cell trajectories for the 91 frames of mouse stem cells with different times of iterations.Fig. 14 (a) shows the generated 30 cell trajectories by eight iterations of the proposed cell tracking algorithm.Fig. 14 (b) shows the generated 30 cell trajectories by ten iterations of the proposed cell tracking algorithm.Fig. 14 (c) shows the generated 28 cell trajectories by twelve iterations of the proposed cell tracking algorithm.Fig. 14 (d) shows the generated 28 cell trajectories by sixteen iterations of the proposed cell tracking algorithm.For the iterations from 14 to 16, the generated cell trajectories match the ground truths of the cell trajectories with 100% accuracy.Although the true negative (TN) rate of the segmentation accuracy by the proposed method for the mouse stem cells is not zero, the cell tracking accuracy reaches 100% accuracy because of the error tolerance ability of the proposed cell tracking algorithm.The computation time in MATLAB for the proposed tracking algorithm with different iterations are shown in Table 4.It is seen that it is more time-consuming than tracking the pancreatic stem cells with higher density because of the relative high resolution of the processed images.For more qualitative results, please refer to the four supplementary videos with the tracked positions denoted by different markers overlaying on the cells.
We use two measures, SEG and TRA defined in [14] to evaluate the proposed cell detection and tracking approach quantitatively.The quantitative comparisons between the proposed framework and the best algorithms described in [14] are shown in Table 5.As can be seen, the proposed approach is more robust than the state-of-the-art approaches described in [14] in detecting and tracking some types of cells.

B. Discussion
The basic image processing algorithms in the proposed approach include slope difference distribution (SDD) based threshold selection method, the low intensity threshold calculation method, the high intensity threshold calculation method, the iterative erosion method, the morphological filtering method, the fusion method and the adjacency based cell tracking method.Among them, SDD based threshold selection method, the iterative erosion method and the adjacency based cell tracking method are fixed for any type of cell images while the other methods need to be adjusted to perform optimally for a specific type of cell image.The priori knowledge of the specific type of cells is very important for the accuracy of cell segmentation.The most accurate segmentation is usually achieved by making fully use of the priori knowledge of the cell type.Both the number of the intensity thresholds and the values of the area thresholds depend on the priori knowledge of the cell type.For some cell types, the methods to calculate the low threshold and the high threshold might also depend on the priori knowledge.For instance, the high threshold is calculated as the valley position with highest intensity for the pancreatic stem cell segmentation.For some types of cell images, e.g. the muscle cell images, only the gradient image is required for cell segmentation [7], which is also an application of priori knowledge.
In [14], most state of the art approaches including the approach with the best overall performance use threshold selection methods for cell detection.The threshold selection methods used in [14] include manually specifying method [17][18] and Otsu's method [19][20].However, Otsu's method could not segment cell images robustly when the contrast between the cells is higher than the contrast between the background and the cells [7].On the contrary, the SDD threshold selection method could calculate the optimal threshold robustly regardless of the contrast differences [7].Quantitative results show that SDD threshold selection method is more robust than all available state of the art image segmentation methods including the threshold selection methods [21][22].Because the proposed approach adopts the most robust threshold selection method for cell detection, it should be more robust than state-ofthe-art approaches based on other threshold selection methods theoretically.In addition, the low intensity and high intensity thresholds computation methods proposed in this paper make full use of the statistical information of the cells' areas and increase the segmentation robustness further.As a result, the proposed approach is more robust than state of the art approaches in detecting the cells tested in this paper.However, there are also some types of cells that could not be segmented by threshold selection methods with adequate accuracy.It is out of the capability of the proposed approach to segment and detect such types of cells.
Many researchers had claimed that deep learning performed the best for image segmentation [23].However, it is not widely validated and recognized.As a matter of fact, deep learning is only used for prediction and traditional image segmentation methods are also needed to segment the cells.Many experimental results showed that the ability of deep learning to improve the segmentation accuracy of traditional methods is limited.After combining deep learning with traditional methods, the accuracy of cell segmentation achieved little improvement [24].Research work of using deep learning for cell tracking indicates that deep learning is not a good option because it could only track one single cell while there are usually thousands of cells [25].In theory, deep learning is only good for classification problems when the number of outputs is significantly less than the number of inputs.To segment an image with resolution 200×200 into two classes, the total number of outputs is 40000  2 , which needs infinite training data.The slope difference distribution based threshold selection method performed better than other image segmentation methods for cell segmentation [7][8].However, it did not achieve the perfect performance yet because of the significantly varied intensities of cells in most cell images.Therefore, we propose a multiple-thresholds-based cell segmentation method in this paper to solve this problem.The iterative erosion could avoid the over-segmentation completely if the algorithms are designed with the priori-knowledge incorporated correctly.
Cell tracking is also very challenging task because of the great number of cells, the significantly varied moving speeds of cells and the unexpected emergence and disappearance of cells.In theory, cell tracking is a multipletargets non-linear and non-Gaussian tracking problem and its optimal solution is intractable for Kalman-filter or particle-filter [26].Multiple hypothesis method has been used frequently for cell or particle tracking [27][28].Unfortunately, the objective comparisons of these notable tracking methods revealed that the quest for better tracking methods remains [29].In this paper, we propose to track the cells by iteratively connecting adjacent cell positions by morphological operations along the time axis and removing small connected domains to eliminate the noise.Therefore, the proposed cell tracking method could tolerate occasionally introduced errors by cell segmentation.In addition, the proposed tracking method is validated with real cell datasets instead of the simulated data used by other methods [27][28][29], which might be more convincing.Once the perfect methods are developed, the automatically identified positions of the cells and the automatically generated cell trajectories should be more accurate than the manually identified positions and generated trajectories because the manual results are affected by random errors while the automatic results are uniform.Specifically, there are many cases where the naked eyes could not discern the cells well while the proposed methods could segment the cells robustly.Therefore, developing automatic methods to segment and track cells is an unavoidable task.Although it still has a long way to go before meeting the ultimate goal that the developed method can segment and identify different types of cells autonomously and the biologists can trust the segmentation and identification results blindly, the proposed approach has the potential to achieve this ultimate goal with the additional help of biologists and software engineers in this field.

V. Conclusion
In this paper, a new approach is proposed to detect and track different types of cells.The proposed multiplethresholds based cell segmentation method could segment the cells robustly after incorporating the previously proposed iterative erosion The proposed cell tracking method connects adjacent cell positions along the time axis to form connected domains and different cell trajectories are then clustered based on their distances to the labeled connected domains.The proposed approach is compared with state of the art approaches with four types of cells.Experimental results show that the proposed approach is more robust in segmenting and tracking some types of cells than state of the art approaches.

2 T
. Fig. 1 (b) shows the flowchart of how to compute the high intensity threshold 3

FIGURE 4 .optimal intensity threshold 1 T 3 T
The segmentation results of the demonstrated Hela cell image by single threshold segmentation method (a) The computed cell positions auto P with the optimal intensity threshold 1 T and iterative erosion overlaying on the original image.(b) The segmentation result by the and iterative erosion overlaying on the original image.(f) The segmentation result by the high intensity threshold 3 filtering.Fig. 6 (a) shows the computed mouse stem cell positions auto P and Fig. 6 (b) shows the fusion segmentation fusion S .As can be seen, some cells that are difficult to discern by naked eyes are identified accurately by the proposed method.

1 T , 1 2
typical result of the mouse stem cells by the proposed multiple thresholds based cell segmentation method (a) The computed cell positions auto P with the three intensity thresholds, three area thresholds and iterative erosion overlaying on the original image.(b) The fusion segmentation result by the three intensity thresholds Proposed Cell Tracking Method After the positions of the cells in all the frames of the image sequences are computed, they need to be linked to generate the trajectories of the cells.All the two-dimensional positions are stacked in the time direction to form the threedimensional dataset in the first dimension of the three-dimensional dataset are denoted by the set For each location with the value of 1, its adjacent locations are also assigned with the value of 1 by the following morphological dilation.
x y z x X y Y y Y I x y z j     (36) .e. the total number of pixels in the jth labeled track.
Fig.12shows the generated cell trajectories for the 92 frames of Fluo-N2DL-Hela cells with different times of iterations of the proposed cell tracking method.Fig.12 (a)shows the generated 174 cell trajectories by six iterations of the proposed cell tracking method.Fig.12 (b)shows the generated 139 cell trajectories by eight iterations.Fig.12 (c)shows the generated 138 cell trajectories by nine iterations and Fig.12 (d)shows the generated 133 cell trajectories by ten iterations.The generated trajectories with ten iterations correlate with the practical cell trajectories most correctly, thus ten iterations should be used for Fluo-N2DL-Hela cells tracking.The computation time in MATLAB for the proposed tracking methods with different iterations are shown in Table2.
Fig.12shows the generated cell trajectories for the 92 frames of Fluo-N2DL-Hela cells with different times of iterations of the proposed cell tracking method.Fig.12 (a)shows the generated 174 cell trajectories by six iterations of the proposed cell tracking method.Fig.12 (b)shows the generated 139 cell trajectories by eight iterations.Fig.12 (c)shows the generated 138 cell trajectories by nine iterations and Fig.12 (d)shows the generated 133 cell trajectories by ten iterations.The generated trajectories with ten iterations correlate with the practical cell trajectories most correctly, thus ten iterations should be used for Fluo-N2DL-Hela cells tracking.The computation time in MATLAB for the proposed tracking methods with different iterations are shown in Table2.

11 . 12 .FIGURE 13 . 14 .
FIGURE 13.Demonstration of generating different cell trajectories with different iterations by the proposed cell tracking method for the pancreatic stem cells in 390 frames (a) Generated cell trajectories with 1 iteration of the proposed cell tracking method.(b) Generated cell trajectories with 2 iterations of the proposed cell tracking method.(c) Generated cell trajectories with 3 iteration of the proposed cell tracking method.(d) Generated cell trajectories with 4 iterations of the proposed cell tracking method.

Table 1 .
Tracking results for the 91 frames of Fluo-N3DH-CHO cells

Table 2 .
Tracking results for the 92 frames of Fluo-N2DL-Hela cells

Table 3 .
Tracking results for the 390 frames of PhC-C2DL-PSC cells

Table 4 .
Tracking results for the 91 frames of Fluo-N2DH-GOWT1 cells

Table 5 .
[14]arison of the accuracy of the proposed framework with best results in[14].