Multi-task ant system for multi-object parameter estimation and its application in cell tracking

Applied Soft Computing 35 (2015) 449–469

Contents lists available at ScienceDirect

Applied Soft Computing journal homepage: www.elsevier.com/locate/asoc

Multi-task ant system for multi-object parameter estimation and its application in cell tracking Benlian Xu a,∗ , Mingli Lu a,b , Yayun Ren a , Peiyi Zhu a , Jian Shi a , Dahai Cheng a a b

School of Electrical & Automatic Engineering, Changshu Institute of Technology, 215500 Changshu, China School of Automation, Nanjing University of Science & Technology, 210094 Nanjing, China

a r t i c l e

i n f o

Article history: Received 7 December 2012 Received in revised form 4 June 2015 Accepted 29 June 2015 Available online 6 July 2015 Keywords: Ant system Cell tracking Parameter estimation

a b s t r a c t Inspired by ant’s stochastic behavior in search for multiple food sources, we propose a cooperating multitask ant system for tracking multiple synthetic objects as well as multiple real cells in a bio-medical field. In our framework, each ant colony is assumed and assigned to fulfill a given task to estimate the state of an object. Furthermore, two ant levels are used, i.e., ant individual level and ant cooperation level. In the ant individual level, ants within one colony perform independently, and the motion of each individual is probabilistically determined by both its intended motion modes and the likelihood function score. In the ant cooperation level, each ant adjusts individual state within its influence region according to heuristic information of all other ants within the same colony, while the global best template at current iteration is found among all ant colonies and utilized to update ant model probability, influence region, and probability of fulfilling task. Our algorithm is validated by comparing it to the-state-of-art algorithms, and specifically the improved tracking performance in terms of false negative rate (up to 10.0%) and false negative rate (up to 2.1%) is achieved based on the studied three real cell image sequences. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The study of analysis of cellular behavior is rapidly becoming a requisite for describing biological processes and diagnosing diseases, and some promising research results on cellular image segmentation or tracking have been reported [1–10], ranging from single cell to multiple cells [3,4], from deterministic analysis method to Bayesian estimate technique [5–7], from phase-contrast microscopy cells to fluorescent image cells [8–10], and so on. Since the study of cellular behavior analysis involves many challenges, such as model uncertainty, morphological variance, and overlapping and colliding between cells, both conventional and manual techniques become tedious and time consuming processes for computing qualitative and quantitative features of a large number of cells. For efficiency and accuracy, the automated tracking methods that eliminate the bias and variability to a certain degree are of great importance. Many efforts have been made in automated cell tracking over the past decades, and could be divided into three groups, i.e., model propagation based method, detection based method, and

∗ Corresponding author. Tel.: +86 51252251170. E-mail address: xu [email protected] (B. Xu). http://dx.doi.org/10.1016/j.asoc.2015.06.045 1568-4946/© 2015 Elsevier B.V. All rights reserved.

Random Finite Sets (RFS) based method. As for the model propagation based method, snakes or active contours are parametric or nonparametric, closed or open curves that can capture object boundary and further track the object of interest [11–14], but it requires cells to be partially overlapping in adjacent frames; level set [15] is able to tackle the changes in object topology, but it will merge two contacting contours into a single one and it also requires cells to be partially overlapping in different frames; meanshift algorithms [16] give a fast solution for object tracking in video sequences (e.g., vehicle tracking, closed loop video), but usually do not give object contours. In terms of the detection based method, also called “detect-before-track” technique [7,17–19], particle filter is usually adopted since it is able to cope with heavy clutter, and is very easy to implement. However, various accurate models, such as state transition model, observation model, object intensity model and shape model, are required to be formulated besides the detection segmentation module. Furthermore, this algorithm usually encounters some problems during the temporal data association stage especially in the case of cell division and cell high density, which lead to tracking failures. For instance, as reported in [19], although the segmentation accuracy is 99.5%, the tracking accuracy is only 90.0%. The last category of tracking algorithm is the RFS-based method, which has demonstrated that it is suitable for tackling the multi-object state estimate problem of spawn events, objects entering and leaving image region

450

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

[20,21]. Therefore, RFS-based filters are naturally suitable for tracking multiple cells in various cases, such as E. coli bacteria tracking and T cell tracking [22–24]. Specifically, the multi-Bernoulli filtering algorithm used in [23] is very competitive and also called a “track-before-detect” method, which by-passes the detection module and exploits the spatiotemporal information directly from the image sequence, and yields better results in low contrast cell image sequences. Ant colony optimization (ACO) is a population-based and metaheuristic approach to model the behavior of almost blind ants in establishing the shortest path from their colony to their feeding sources and back by laying a chemical substance called pheromone. By virtue of the learning and searching power of an autocatalytic self-organization of an ant colony, ACO-based algorithms are receiving increasing attention and have enjoyed great success in the solution of traditionally difficult optimization problems [25,26]. To date, most of these algorithms are improved and extended to applications such as vehicle routing [27], scheduling [28], clustering [29–31], classification task [32], and target tracking [33–36], but the application to biological cell tracking is seldom reported. As we know, in addition to the aforementioned learning and searching power of ant colony, there are other salient features for the track of multiple biological cells, and they could be characterized as follows. First, there usually exist dynamic differences between tiny cells (such as cells 1 and 2 in Fig. 1), and we have to resort to well-designed motion models for acquiring better performance in the general tracking framework. However, through the use of all-orientation searching strategy of ant colony, the constraint of requiring accurate models for different tracking modes can be relaxed, i.e., we utilize the proposed ant decision of simple form to accommodate the differences in cell dynamics. Second, the number of cells varies over time, i.e., tiny cells leaving or entering the image (cell 3 in Fig. 1), and this requires the proposed approach could adaptively and jointly estimate the number of cells and its individual states. In the framework of ant colony, we define the mapping between colonies and cells, namely, some ant colonies are modeled to correspond to existing or new entering cells with a big value of colony existing probability, whereas some to disappearing cells with a small existing probability, which will be removed in the next frame. Finally, the variation in cell distribution density often happens (cells 1 and 2 in Fig. 1), and this easily leads to the “hijacking” problem for adjacent or interacting cells. However, for a group of ant colonies, the independent and cooperating mechanism is designed and expected to be suitable for solving the problem. The main contributions of our method can be summarized as follows: (1) introducing an idea of bioinspired estimator, i.e., multi-task ant system, for estimating the states of multiple cells (objects), which opens another interesting research branch for both swarm intelligence and object tracking; (2) introducing a series of ant-related behavior variables, such as moving mode, influence region and ant colony existing probability, and the corresponding recursive forms are formulated as well; (3) considering two ant levels, namely, the ant individual level and the ant cooperation level. The ant individual level propels ants with the same task to move toward the location of the same interested cell, whereas in the ant cooperation level the cooperation mechanism among colonies is utilized to update and regulate ant-related behavior parameters for accurately locating the cells of interest in the following iterations as a whole. Preliminary results have been announced in the conference paper [37]. This paper presents a more complete theoretic study and abundant experimental analysis. We first address the background on multi-Bernoulli filter for multi-object image tracking in Section 2. In Section 3, a multi-task ant system algorithm is then proposed to track multiple cells and determine the number of cells,

which utilizes the heuristic information and cooperation mechanism among ant colonies to find potentials in a recursive way. In Section 4, the experiment results on various cases are presented to demonstrate the effectiveness of our algorithm. Finally, conclusions are summarized in Section 5.

2. Background on multi-Bernoulli filter for multi-object image tracking In the framework of multi-object tracking, our focus is to determine the posterior multi-object distribution given all history observations. With the assumption that each object follows Markov process, the Bayesian filtering algorithm offers a concise way to describe the multi-object tracking problem. Since both the states and measurements are represented as sets of random vectors, the concept of random finite set (RFS) is introduced to solve the multiobject estimation problem in the Bayesian framework. For tracking n objects at time k, the multi-object state is denoted by a RFS X k = { x1,k , x2,k , . . ., xn,k }, where xi,k is the state vector of i-th object. Let z k = [z1,k , z2,k , . . ., zm,k ¯ ] denote the image observation com¯ pixel values at time k, and z 1:k is defined prising an array of m as the cumulative image observations up to time k. Therefore, if the ˆ k|k (·|z 1:k ), the recurposterior multi-object density is denoted by  sive form of multi-object state estimation problem can be defined as

 ˆ k|k−1 (X k |z 1:k−1 ) = ˘ ˆ k|k (X k |z 1:k ) = ˘



ˆ k−1|k−1 (X|z 1:k−1 )ıX f (X k |X)˘

ˆ k|k−1 (X k |z 1:k−1 ) h(z k |X k )˘

(1)

ˆ k|k−1 (X|z 1:k−1 )ıX h(z k |X)˘

where f(·) is the multi-object transition density function, h(·) denotes the observation likelihood function, and ı is an appropriate reference measure on some state space. Observe that, due to the fact that there are multiple integrals implied in the above recursion, the estimate of posterior multiobject density is in essence an intractable problem in most cases. However, an approximated solution to this problem is the probability hypothesis density (PHD) recursion [20], in which the first order statistical moment of the posterior multiple-object state is propagated in time in the Bayesian recursion form, instead of propagating the multiple-object posterior density as a whole. More recently, the multi-Bernoulli filter developed from the RFS framework is proposed and utilized as a more effective visual tracking technique [22,23], and the salient features of the filter are that it operates in the single-target state space and in parallel. To give an intuitive understanding of the multi-Bernoulli filter, which will be discussed in Section 4 for performance comparison with our approach, we briefly introduce the principle and its notions. With the assumption of multi-object density being characterized by multi-Bernoulli parameters, the multi-Bernoulli filter for image sequences can be represented, respectively, by a prediction step and an update step. In the prediction step, if the multi-Bernoulli parameters of X k Mk−1 (i) (i) ˆ k−1|k−1 = {(rk−1 , pk−1 (·))} is denoted by  , where the paramei=1

(i)

(i)

ter rk−1 is the existence probability of the i-th object and pk−1 represents the probability density of the state conditioned on its existence, the predicted multi-object state is also a multi-Bernoulli RFS (i)

(i)

ˆ k|k−1 = {(rk|k−1 , pk|k−1 (·))} 

Mk−1 i=1

(i)

(i)

∪ {(r,k , p,k (·))}

M,k i=1

(2)

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

451

Fig. 1. The possible coincidence between the evolvement of working ant colony and cell motion.

(i)

(i)

where {(r,k , p,k (·))}

M,k i=1

describe the parameters of multi-

Bernoulli RFS of births at time k, and the predicted parameters Mk−1 (i) (i) {(rk|k−1 , pk|k−1 (·))} are given by (i) rk|k−1

=

(i) rk−1



i=1



(i) pk−1 (·), ps,k (·)

 (i) pk|k−1 (x)



(i)

(3)

f (x|·), pk−1 (·)ps,k (·)



=



(i)

pk−1 (·), ps,k (·)

with the object existence probability ps,k (·) at time k and the    standard inner product notion v, h = v(x)h(x)dx. (i)

(i)

Given the predicted parameters {(rk|k−1 , pk|k−1 (·))}

Mk|k−1 i=1

of

multi-Bernoulli RFS, the updated parameters of multi-Bernoulli RFS are formulated as (i)

(i)

Mk

ˆ k|k = {(rk , pk (·))}i=1  where (i)



(i)

rk|k−1 pk|k−1 (·), h(z k )

(i)

rk =



(4)

(i)

(i)



(i)



1 − rk|k−1 + rk|k−1 pk|k−1 (·), h(z k ) (i)

(i) pk (·)

=



pk|k−1 (·)h(z k ) (i)

(5)



pk|k−1 (·), h(z k )

With the above recursive form and a problem-dependent likelihood function h(·), the generic multi-Bernoulli filter can be easily extended to track multiple objects in a sequence of raw image frames [23].

and some of ant colonies are merged as one colony due to the same task, which means the old colonies disappear or the existing prob(i) ability tends to zero (corresponding to rk−1 → 0); whereas some ant colonies work independently all the time with different tasks, which means that these ant colonies have been working with big (i) values of existing probability (corresponding to 0.5 < rk−1 ≤ 1). On the other side, for each ant colony, its distribution can be considered as an approximation to the state density function of an object, (i) (·) (·) and we further have pk−1 ∼{xi,k−1 }, where xi,k−1 denotes the state of an ant in colony i at time k − 1. Therefore, due to the underlying similarities between the definition of working mechanism of real ant colonies and the multi-Bernoulli filter, it is natural for us to develop a multi-task ant system to track multiple objects in a series of low contrast image sequences. 3.1. Ant individual and colony behavior description An ant colony can do various jobs in a random but wholly regular way, such as gathering corpses to clean up their nest, transporting larvae to place them by size, and moving between nest and food back and forth by chemical substance called “pheromone”. These tasks ranging from simplicity to complexity attribute to the power of ant individual and colony cooperative searching behaviors. The ant individual behavior shows the independence and randomness of an ant, while the ant colony behavior defines the communication and cooperation among ants. In our proposed algorithm, to fulfill different tasks, we define two ant behaviors as below. 3.1.1. Ant individual behavior description Assume that each ant has three modes, i.e., moving forward, turning right, and turning left, and the mode probability of ant i in colony s at time k is denoted by a triple vector (i),3 T

3. Multi-task ant system for multi-cell tracking

s,k = [s,k , s,k , s,k ]

A real ant system usually consists of several ant colonies working together or independently for the same goal or different tasks,

where elements in the vector represent, the probabilities of moving forward, turning right and turning left respectively, and we have

(i)

(i),1

(i),2

(6)

452

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

different variances or influence regions. For a big variance, an ant can move within a larger circle region centered on the mean, but as the variance decreases, the movement of the ant is restricted to a relatively smaller region in the end.

3.1.2. Ant colony behavior description In our ant tracking system, we usually assign m ant colonies to search for n(m > n) potential objects, and all ants within the same colony have the same goal, i.e., locating the same object of interest. If two different colonies, whose distance is less than a given threshold, try to find two objects or the overlap ratio between two colonies is in excess of a prior threshold, the two colonies are then merged together, and these two colonies are treated as one in the following iterations. The main considerations of ant colony behavior are as below. To exploit more useful information, we introduce a communication mechanism among colonies to help ants locate more objects (dash line in Fig. 2). Therefore, a histogram template pool1 T{k+1} is built in a recursive way T{k+1} = [T{k} , Tkbest ]

Fig. 2. Ant behavior in our ant system.

if h(Tkbest ) > t0

(8)

where Tkbest denotes the best histogram template at the current iteration with the highest likelihood score, h(·) denotes the likelihood function, t0 is a predefined threshold, and we set the cardinality of histogram template pool to be |T{k+1} | ≤ Nmax . Also, to characterize how a given ant colony s evolves as a whole, we introduce a variable qs,k (qs,k ∈ [0, 1]) to evaluate the existing probability of ant colony s. If qs,k ≥ q¯ 0 (where q¯ 0 is a predefined threshold), the ant colony s is supposed to have found the object of interest and could be used to extract the corresponding cell state. In our proposed algorithm, the variable propagates over time as below qs,k+1|k = qs,k

(9)

qs,k+1 = f¯ (qs,k+1|k , h(·)) where the update function f¯ (·) will be discussed in Section 3.3.

3.2. Ant individual decision level Fig. 3. Ant moving ability corresponding to different variances.

 j=1,2,3

(i),j

s,k = 1. At the initial stage of iteration, for instance, each (i)

element in s,k is assigned the same value due to lack of prior information. As iteration evolves, the prediction of mode probability of ant i can be written as

⎡

⎤

P1,1

P1,2

P1,3

s,k+1|k = PM × s,k = ⎣ P2,1

P2,2

P2,3 ⎦ × s,k

P3,1

P3,2

P3,3

(i)

(i)

⎢

⎥

(i)

(7)

where PM is the mode jump matrix, and Pi1,j1 denotes the probability of mode j1 jumping to mode i1. With the evolvement of mode probability, as shown in Fig. 2, an ant initially located at position (1) will probably turn left, then move forward to position (2), and finally to position (3) where an object is located. Once the ant’s current location is close to the object, it will keep track of the object in the following frames through adjusting its moving ability. Therefore, we define each ant has a varying moving ability over time, which is characterized by its influence (i) region s (k). Fig. 3 illustrates an ant is assumed to have different moving abilities represented by the same mean or position but with

In the ant individual decision level, we assume each colony work independently, and multi-object tracking system can be treated as several single-object tracking sub-problems in which one ant colony corresponds to one object. Next, without lossof generality,  we only describe the decision behavior of ant i(i ∈ 1, ..., N ) in





colony s(s ∈ 1, ...m ). Indeed, since our ant system belongs to the class of swarm intelligence algorithm, it still follows some generic definitions: (1) an appropriate description of the problem, so that it lends itself to incrementally building a solution to the problem, such as selecting state in sequence from time k to k + 1, and so on; (2) a probabilistic transition rule based on the value of the heuristic function and/or on the pheromone information, such as moving mode, state transition probability, likelihood score, etc., as shown in Fig. 4. (i) At time k, the state of ant i, xs,k , is known and usually represented (i)

(i)

(i)

(i)

by its position (xs,k , ys,k ) and velocity (x˙ s,k , y˙ s,k ), thus the decision of ant i at next step can be defined with a probability as

1 Note that there are only a limited number of samples in the initial template poolT{k+1} , which are obtained and extracted from the first frame.

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

453

Table 1 The pseudo-code of local searching strategy of each ant. (j)

Initialization: given an ant j with position xs,k+1|k and influence region (j) s,k ;

threshold 0 Generate a random number generated in [0,1] If < 0 (j) Generate new candidate x¯ s,k+1 according to Gaussian distribution (j)

(j)

N(x; xs,k+1|k , s,k ) (j)

(j)

if h(z k+1 |¯xs,k+1 ) > h(z k+1 |xs,k+1|k ) (j) xs,k+1|k

←

(j) x¯ s,k+1

End if Else (j) Keep xs,k+1|k End if

Fig. 4. The topology of one step of ant decision.

where , ,  and  are the adjustment coefficients, and gk is defined as below |T{k} | p

 (i,j)

ps,k+1|k =

(i),j

(j)

(i)

˛

s,k+1|k f (j) (xs,k+1|k |xs,k )

 (i) l∈˝s



(i),l

(l)

(i)

(j)

[h(zk+1 |xs,k+1|k )]

˛

s,k+1|k f (l) (xs,k+1|k |xs,k )

(l)

1  gk = min(uk (j), u˜ i (j)) |T{k} |

ˇ

[h(zk+1 |xs,k+1|k )]

where u˜ i (j) denotes the value of j-th element of i-th histogram u˜ i in a template pool, p is the total number of elements in a histogram and takes the value of 256 × 3 for an RGB image, and |T{k} | denotes the number of samples in histogram template pool T{k} .

ˇ

(10) (i),j

(i)

where mode probability s,k denotes the j-th element of s,k ; (j)

(i)

xs,k+1|k is the one step prediction of state xs,k following a Markov process with a transition density f(j) (·), observation image z (·) is (j)

usually modeled by the likelihood function h(z (·) |xs,(·) ) conditioned (j)

on the current state xs,(·) , j corresponds to one of predicted positions from three modes, ˛ and ˇ are the importance adjustment parami denotes the space of three predicted positions of ant eters, and ˝(s) i. Remarks

• Mode probability describes the inclination degree to which an ant moves toward one of three predicted positions. A big value indicates that the ant gradually moves toward the location of an object of interest, and even though it might make several undesirable decisions during movement, it still keeps on adjusting itself and eventually moves closely to an object. • State transition density encompasses all possible moving modes of an object, and, for simplicity, all modes are assumed to follow Gaussian distributions but with different means and variances. From the perspective of cell tracking, since each cell moves randomly in the observation image, the combination of three modes is naturally suitable for capturing the dynamic uncertainty of a cell, which may perform a strong or weak maneuver at any time. • Measurement likelihood function plays an important role in ant’s decision, and it is usually selected according to empirical analysis in computer vision. In terms of track-before-detect tracking, literature [23] proposed a likelihood function based on color space conversion using kernel density method, and such technique resorts to a large number of training data. However, our algorithm employs a novel likelihood function to discriminate potential regions in a convenient way, which operates directly on RGB space and requires a small number of templates. Therefore, for a given histogram uk of an RGB image, the corresponding likelihood score is computed as



h(uk ) =  e−(1−gk )



(12)

i=1 j=1

 (11)

To view the effectiveness of our proposed likelihood function in a direct manner, Fig. 5 illustrates an example of comparison results on various regions of a real cell image. In terms of two regions occupied by two cells (i.e. case 1 in Fig. 5(a)), the similarity measure is about 0.99974, as shown in Fig. 5(b). If one region is fully occupied by cell and the other is covered by cell nearly in half (i.e. case 2 in Fig. 5(a)), the similarity decreases to 0.73539, as illustrated in Fig. 5(c). However, if one cell leaving the interested area, we find that the similarity is close to zero, as shown in Fig. 5(d). It is observed that the above function can discriminate the difference of any two regions of the same size, and it is, therefore, utilized to act as a similarity measure between the candidates and template samples. As mentioned before, some ants might make undesirable choices due to the probabilistic selection. To remedy such consequence and guide more ants moving toward the region of interest, a local adjustment strategy is proposed within the influence region of each ant. Note that the influence region is represented by the variance of each ant, which means that an ant will search for an object in a large scale of area with a big variance and vice versa. Without loss of generality, we consider the search behavior in one dimension, and other dimensions follow the same rule. For a given random number generated in [0,1], if is smaller than a predefined threshold 0 , the local search step is done, as illustrated in Table 1. 3.3. Ant cooperation decision level The salient feature of swarm intelligence, such as the ACO algorithm, is the information exchange among agents. In our algorithm, we use multiple ant colonies to track multiple objects and each colony in principle when implementation corresponds to one object, thus information exchange happens both   within each colony and among colonies. For ant i(i ∈ 1, ..., N ) in colony s(s ∈





1, ..., m ) at time k, we first use a c1 × c2 rectangle region tem(i)

(i)

plate whose center is closest to (xs,k , ys,k ), then we calculate its corresponding histogram likelihood score according to Eq. (11), and the best template with the highest likelihood score is finally best . Following this way, the best template found and denoted by Ts,k

454

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

Fig. 5. Comparison results on various regions using our proposed likelihood function.

Tkbest at the current iteration can be found and exchanged among

m colonies. Moreover, if h(Tkbest ) > t0 , a new augmented template is declared until |T{k+1} | = Nmax , which means that a complete template pool has already been built. In contrast to the strategy proposed in [23], our template pool is built in a dynamic way through

learning rule with less computation burden, and Table 2 gives the generation procedure of template pool T{k} . In the generic ACO algorithm, update procedure using the global best information is requisite to further regulate ant behavior in the following iterations. Since this step constitutes a positive feedback

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

455

Table 2 The generation procedure of histogram template pool in a dynamic manner. Initialization: Get the RGB histograms using a c1 × c2 template block, and each histogram corresponds to one cell in the first frame. Therefore, the initial histogram template pool is built and denoted by T{0} , and the number of T{0} is equal to that of cells in the first frame. For time k = 1 : t For ant colony s = 1 : m  

best within ant colony s according to h(uk ) = (e−(1−gk ) ) . Find the best histogram template Ts,k End for Find the histogram template among m ant colonies, and denoted by Tkbest ;

 

Calculate the similarity, denoted by gaver , between Tkbest and T{k−1} according to h(uk ) = (e−(1−gk ) ) ; If gaver > t0 and |T{k−1} | < Nmax T{k} = [T{k−1} , Tkbest ] Else T{k} = T{k−1} End if End for

mechanism, a well-established updating strategy not only guarantees the quality of solution, but also accelerates convergence speed. Without loss of generality, for ant colony s, once the searching behavior of each ant is finished, the likelihood score corresponding to each predicted position is calculated for re-evaluating the importance of each ant. Therefore, we have (i)

(i) temp,k+1

ωs

=

(i)

h(z k+1 |xs,k+1|k )ωs,k+1|k N 

(13)

(i),1

(j) (j) h(z k+1 |xs,k+1|k )ωs,k+1|k

(i) s,k+1

j=1 (i)

velo(i) = k

(x

(j)



if (i)

− xs,k )/T

s,k+1|k ⎪ ⎪ ⎪ ⎩ −velomax

(j)

(i)



xs,k+1|k − xs,k /T > velomax

(j)

(i)

if |xs,k+1|k − xs,k |/T ≤ velomax



if

(j)

(i)



xs,k+1|k − xs,k /T < −velomax (14)

where velomax denotes the maximum speed of ant i, and T is the sampling interval. Meanwhile, one step distance of ant i is also restricted within (i) (i) the range of [min , s,k ], where s,k ≤ max ,  min and  max denotes the minimum and maximum influence regions, respectively. Once all ants finish their individual state updates, denoted by



(i)

xs,k+1

N

(i)

(i),3

3 

T

(i),j

s,k+1

j=1

(16)

(j)

where

(i),j s,k+1

=

h(z k+1 |xs,k+1 ) 3 

(i),j s,k+1|k

(l)

h(z k+1 |xs,k+1 )

l=1

Note that the above likelihood score is computed on the histogram template pool T{k+1} , and such a definition encourages more ants to select the mode with a higher mode probability in the following decisions. In our proposed multi-task ant system, each colony in principle corresponds to one object, and our focus, therefore, is how to describe the degree to which an ant colony has found the corresponding object, which is also represented by the existing probability of ant colony qs,k+1 = f¯ (qs,k+1|k , h(·))



=



qs,k+1|k f (xs,k+1 ), h(z k+1 |xs,k+1 )





1 − qs,k+1|k + qs,k+1|k f (xs,k+1 ), h(z k+1 |xs,k+1 )

(17)

where 1 − qs,k+1|k describes the degree to which ant colony s has not found the corresponding object, or the disappearing probability of ant colony s at time k + 1. 3.4. Remarks and implementation issues

, we further use the information to update the influence i=1

region represented by s,k+1 =

=

(i),2

[s,k+1 , s,k+1 , s,k+1 ]

(i)

where the prediction weight ωs,k+1|k is equal to the previous ωs,k , and the calculation of likelihood function is based on the template pool T{k} . In order to fully characterize the individual difference of each ant, two parameters, i.e., ant state and its influence region, are jointly updated to achieve better tracking performance based on (i) the above ωs temp,k+1 . In terms of the state update, we follow the same rule as Eqs. (5)–(9) in [35], but both the velocity and influence region restrictions are introduced to make ants move toward area (i) (j) of interest step by step. For instance, xs,k and xs,k+1|k , respectively, represent the previous state of ant i and the state to be selected by ant i, and the following rule is applied to ant i

⎧ ⎪ velomax ⎪ ⎪ ⎨

score over all ants in colony s, and K is the adjustment coefficient. It can be seen that the bigger the likelihood score is, the smaller the influence region is, which means that if an ant is close to an object of interest, it will continue searching for a better candidate within a small area around its current position. Because the new position of each ant is finally obtained, its corresponding deterministic moving mode is evident. Therefore, in terms of updating mode probability, we use three different mode candidates as a basis to normalize weights

K ·ω ¯ s,k+1 (i)

ωs,k+1

(i)

s,k+1|k

(15)

(i) where s,k+1|k denotes the prediction of influence region with (i) (i) (i) s,k+1|k = s,k , ωs,k+1 is the update likelihood score using the same (i) ¯ s,k+1 denotes the average likelihood form as Eq. (13) with xs,k+1 , ω

Although we have addressed the main attributes of our algorithm in previous sections, some contributions, which do not appear in [23,34,35], are summarized and listed below. • Candidates. In terms of the ant individual selection level, our algorithm assumes that each ant decision is based on three candidates corresponding to three modes, however, in [34], the stochastic selection is evaluated on all ants’ predicted N positions where N 3;

456

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

including synthetic image data and real cell image sequences, and all examples are conducted on a DELL 6 GHz processor with 1.99 GB RAM. 4.1. Application to synthetic image objects tracking 4.1.1. Parameter settings This section aims to test theoretically the performance of our algorithm on a series of known synthetic images. To obtain an objective comparison with the multi-Bernoulli filter proposed in [22], we employ the observation model as Eq. (24) in [22] instead of Eq. (11) in our algorithm. In the first example, we only track one object moving in an observation image with the size of 200 pix × 200 pix (one pixel is equivalent to 1 m, as shown in Fig. 7), and one ant colony of birth is generated at each iteration with a mean state x0 = [160 m, 1 m/s, 160 m, − 1 m/s]T . The simulation uses a maximum of Lmax = 500 and minimum of Lmin = 200 ants per hypothesized track, three ant modes are, respectively, defined as a uniform linear motion and two nearly constant turn models with ω1 = /5 and ω2 = − /5, and the maximum and minimum influence regions are taken as  max = 20pix and  min = 2pix, respectively. Finally, a fixed model  0.6 0.2 0.2 jump matrix is introduced as PM = 0.2 0.6 0.2 , parameters 0.2 0.2 0.6 ˛ and ˇ are set to be 1 and 0.01, respectively, and the remaining parameters are selected as in [22].

Fig. 6. The main framework of our proposed multi-task ant system.

• Movement. In our algorithm, each ant moves gradually toward region of interest, but in [35] some ants are allowed to move at a big step, which results in big errors in velocity element of state; • Efficiency. Our algorithm develops a novel likelihood function and builds the template sample pool in a dynamic way through learning rule. However, in [23], a large number of samples (850 or more) are first trained off-line, and the traditional kernel density estimate with Bhattacharyya distance is then used to calculate the likelihood score, which leads to an inaccurate tracking and less efficiency. In terms of the issues on algorithm implementation, ant colony at each step consists of two parts, i.e., the remaining ant colonies (except for the first iteration) and a fix number of ant colonies added at each step. We also assume that the total number of ant colonies is larger than that of objects, therefore, as the technique in [21–23], the merge process should be considered at the stage of state estimate for the sake of computation burden and prior constraints, and if the overlap ratio of two blobs corresponding two ants, respectively, is greater than a given threshold, the merge process is declared. To visualize our proposed algorithm in a full view, we summarize the main blocks as in Fig. 6, and corresponding implementation procedure is illustrated in Table 3. 4. Experiments In this section, the tracking ability of our proposed algorithm is investigated through a number of tracking examples

4.1.2. Results Fig. 8 illustrates the tracking results in x and y directions using our algorithm and the multi-Bernoulli filter respectively, and it can be observed that our algorithm is competitive and works well when it is utilized to track image object. Fig. 9 presents the 2-D plane tracking results, and the above conclusion is further verified, but it seems that the multi-Bernoulli filter performs better in this case. The underlying reason is that the multi-Bernoulli filter in the settings overestimates the number of real targets (as shown in Fig. 10), and the locations of these estimated targets are in essence close to each other. Therefore, in order to fully characterize the tracking performance, we introduce the metric of OSPA distance [36,38], and it can be seen, as shown in Fig. 10, that the OSPA distance of our algorithm is significantly smaller than that of the multi-Bernoulli filter. In our second example, we investigate a relatively challenging tracking scenario with a time varying but a maximum number of four targets, and simulation results are presented as well, as shown in Fig. 11 and Figs. 12–14, respectively. All parameters are set as the above except the birth ant colony term. In this example, we use four birth ant colonies at each iteration with mean state x1,0 = [120 m, 1 m/s, 120 m, −1 m/s]T , x2,0 = [70 m, 1 m/s, 70 m, −1 m/s]T , x3,0 = [50 m, 1 m/s, 50 m, −1 m/s]T , and x4,0 = [20 m, 1 m/s, 20 m, −1 m/s]T . It is noted that both our algorithm and the multi-Bernoulli filter can simultaneously give the state estimates of multiple objects. However, as the results of the first example, our algorithm performs better in terms of the cardinality estimate and the tracking accuracy, as shown in Fig. 14 and Figs. 12 and 13, respectively. So FPR, we have evaluated the performance of multi-target tracking on a series of synthetic image sequences, and it is found that our algorithm is working well and has achieved satisfactory results as we expected. Next, we will employ our algorithm to the problem of real cell tracking on various scenarios in the following section.

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

457

Fig. 7. The object moving trajectory. 200 x-coordinate (m)

x-coordinate (m)

200 Our algorithm True tracks

150 100 50 0

10

15

20

25 Time

30

35

40

45

50

5

50

10

15

20

25 Time

30

35

40

45

50

200 Our algorithm True tracks

150

y-coordinate (m)

y-coordinate (m)

100

0 5

200

100 50 0

Multi-Bernoulli filter True t racks

150

Multi-Bernoulli filter True t racks

150 100 50 0

5

10

15

20

25 Time

30

35

40

45

5

50

a) The tracking result using our algorithm

10

15

20

25 Time

30

35

40

45

50

b) The tracking result using the multi-Bernoulli filter

Fig. 8. The tracking results in x and y directions using different methods. 180

180

Multi-Bernoulli filter True track s

Our algorithm True tracks

160

160

140 y-coordinate (m)

y-coordinate (m)

140

120

120

100

100

80

80

60

0

20

40

60

80 100 x-coordinate (m)

120

140

a) The tracking result using our algorithm

160

180

60

0

20

40

60

80 100 x-coordinate (m)

120

140

160

180

b) The tracking result using the multi-Bernoulli filter

Fig. 9. The tracking results in 2-D image using different methods.

5 Estimated value True value

4 3 2 1 0

5

10

15

20

25

30

35

40

45

50

Number of targets

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

Number of targets

458

5 Estimated value True value

4 3 2 1 0

5

10

15

20

1

10

0

10

-1

5

10

15

20

25

30

35

40

45

50

OSPA distance (in m)

OSPA distance (in m)

2

10

10

25

30

35

40

45

50

35

40

45

50

time step

time step 2

10

1

10

0

10

-1

10

5

10

15

20

time step

25

30

time step

a) The tracking result using our algorithm

b) The tracking result using the multi-Bernoulli filter

Fig. 10. Target number estimate and OSPA distance using different methods (p = 2, c = 100 m).

Fig. 11. The objects moving trajectory.

200

Our algorithm True tracks

100 50 0

5

10

15

20

25 Time

30

35

40

45

50

100 50

5

10

15

20

25 Time

30

35

40

45

50

200 y-coordinate (m)

Our algorithm True tracks

150 100 50 0

150

0

200 y-coordinate (m)

Multi-Bernoulli filter True t racks

150 x-coordinate (m)

x-coordinate (m)

200

5

10

15

20

25 Time

30

35

40

a) The tracking result using our algorithm

45

50

Multi-Bernoulli filter True tracks

150 100 50 0

5

10

15

20

25 Time

30

35

40

45

b) The tracking result using the multi-Bernoulli filter

Fig. 12. The tracking results in x and y directions using different methods.

50

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

459

Table 3 The implementation procedure of our algorithm of multi-task ant system. 0. Initialization Given a frame image, we first uniformly divide it into several sub-regions (4,9, or more), and assign a number of ants to each sub-region, which constitutes an initial ant colony. From the second frame image, the existing ant colonies also include the predicted ant colonies from the previous frame. 1. In the ant individual decision level For ant colony s = 1 : m For each ant i in colony s (i = 1 : L)



Select potential location according to

(i,j) ps,k+1|k

=



l∈˝

(j)

˛

(i),j (j) (i) f (j) (x |x ) s,k+1|k s,k+1|k s,k





(i) s

[h(z k+1 |x

˛

(i),l (l) (i) f (l) (x |x ) s,k+1|k s,k+1|k s,k



ˇ (j) )] s,k+1|k

[h(z k+1 |x

ˇ (l) )] s,k+1|k

;

In order to find better location around xs,k+1|k , a local searching strategy is implemented if < 0 holds, as illustrated in Table 1; End Once all ants in the given colony finish the above two decisions, find the highest likelihood score in the colony according to h(uk ) =  e−(1−gk )



!

, and

the corresponding template Tsbest as well; End 2. In the ant cooperation decision level Find the global highest likelihood score among m colonies and update the template sample pool, as illustrated in Table 2; For each ant in any given colony, the corresponding state is updated using the same rule as in [35], and Eq. (14) as well; (i) (i) (i) Update ant’s influence region s , mode probability s , importance weight ωs , and ant colony existing probability qs,k+1 according to Eqs. (13), (15)–(17); If qs > q¯ 0 (a threshold, we set it to be 0.01), ant colony s, we say, is retained for the next frame image, otherwise delete it. 3. Output According to the centroid distance measured between colonies, we merge two or more spatially neighboring ant colonies that they might be interested in the same object. After merging process, one colony corresponds to one object, and the state of an object is easily obtained through averaging all states of ants within the same colony.

4.2. Application to real multiple cells tracking 4.2.1. Parameter settings In this section, we will test the performance of our algorithm on several low contrast cell image sequences, and these sequences include various challenging scenarios, such as the dynamic differences between cells, the variations in cell morphology, and the changes in cell density, etc. In order to track all possible cells, we introduce a given number of birth ant colonies (four ant colonies in the first three scenarios and nine in the last scenario) in each frame, and these birth ant colonies are initially assigned to different parts in current image. We use a maximum of Lmax = 200 and minimum of Lmin = 100 ants per hypothesized track, and the maximum and minimum influence regions are  max = 80 pix and  min = 15 pix, respectively. Furthermore, we assume that the initial existing probability of each colony is qs,0 = 0.1, and the maximum of histogram template pool, i.e., Nmax , is equal to 40.

In our examples, each cell is represented by a square blob, and (i) the state of each ant xs,k is denoted by position and velocity in x and y directions, respectively. Three ant modes are defined as uniform linear motion and two nearly constant turn models with ω1 = /3 and ω2 = − /3. Therefore, the corresponding dynamics (i),j (i) of each ant is assumed to follow xs,k+1 = Aj xs,k + wk (j = 1, 2, 3), where Aj is the state transition matrix of mode j, and wk is the process noise with zero mean and variance Q = diag (x2 , x2˙ , y2 , y2˙ ) = diag (1, 0.1, 1, 0.1). Other related parameters are listed in Table 4. 4.2.2. Results In scenario 1, as shown in Fig. 15(a), cell 1 with nearly round shape moves slowly, while cell 2 undergoes relatively fast speed and its shape keeps on changing. Since the state is only represented by position and velocity, and other features such as size and shape are not considered, this enhances the tracking difficulties in low contrast image sequences. However, using our proposed multi-task

140

140 Our algorithm True track s

100

100

80

80

60 40

60 40

20

20

0

0

-20

0

20

40

Multi-Bernoulli filter True tracks

120

y-coordinate (m)

y-coordinate (m)

120

60 80 100 x-coordinate (m)

120

140

a) The tracking result using our algorithm

160

-20

0

20

40

60 80 100 x-coordinate (m)

120

140

160

b) The tracking result using the multi-Bernoulli filter

Fig. 13. The tracking results in 2-D image using different methods.

460

Number of targets

Number of targets

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469 6 4 2 Estimated value 0

5

10

15

20

25

30

True value

35

40

45

10

5

Estimated value 0

50

5

10

15

20

2

10

1

10

0

10

-1

10

5

10

15

20

25

30

25

30

35

40

45

50

35

40

45

50

time step

35

40

45

50

time step

OSPA distance (in m)

OSPA distance (in m)

time step

True value

2

10

1

10

0

10

-1

10

5

10

15

20

25

30

time step

a) The tracking result using our algorithm

b) The tracking result using the multi-Bernoulli filter

Fig. 14. Target number estimate and OSPA distance using different methods (p = 2, c = 100 m).

ant system, dynamics, shape and size differences between cells are overcome and satisfactory tracking results have been achieved, as illustrated in Fig. 15(b). Fig. 16 presents the above position estimate and the number of cells in each frame, and it further verifies the ability of tracking multiple cells in low contrast image sequences. Without loss of generality, two measures are introduced, namely, false negative rate (FNR) and false positive rate (FPR), and defined as FNR =

t  k=1

t 

ˆ (k) / H miss

(k)

Htrue and FPR =

k=1

t 

t 

ˆ (k) / H extra

k=1

(k)

Htrue , respec-

k=1

ˆ (k) denotes the number of cells that are missed at tively, where H miss (k) ˆ time k, H the number of non-existing cells that are tracked, and (k)

extra

Htrue the number of true cells at time k. We first record all false negative reports and false positive reports in each frame over 20 Monte-Carlo simulations, then the averaged FNR and FPR are computed as around 1.33% and 0.67%, respectively, using our algorithm. It is noted that the false positive mainly comes from a small shade region at the left side of cell 1 in Fig. 15(a), since with iteration evolving more ants gather around cell 1 and are subject to mistaking this small region as a real cell. We also find that the missing cells (false negative reports) take place at the initial tracking phase mainly due to lack of prior information such as the limited elements in T{k} which is built dynamically over time. To compare the performance with the state-of-the-art multi-Bernoulli filter, we use the same likelihood function as Eq. (11) but total 40 histogram templates are extracted off-line in advance, and the averaged FNR and FPR are calculated as 4.50% and 0.83%, respectively. In terms of both the false negative reports and the false positive reports, or both the averaged FNR and FPR, our algorithm is competitive with the multi-Bernoulli filter as a whole. In scenario 2, our aim is to further test the ability of tracking spatially adjacent moving cells with our proposed algorithm. As shown in Fig. 17(a), it is observed that cell 2 undergoes shape change and sudden variation in dynamics, and it moves around cell 1 from frames 20 to 29. An example of successful tracking is shown in

Fig. 17(b). Our algorithm can track spatially adjacent cells and it can discriminate the clutter and cells of interest. Similarly, over 20 Monte-Carlo runs, the averaged FNR and FPR are 8.44% and 2.03%, respectively, using our algorithm, and 18.43% and 2.65%, respectively, using the multi-Bernoulli filter. In particularly, these false negative reports in our algorithm mostly happen between frame 22 and frame 29, but the multi-Bernoulli filter remains higher reports during all tracking frames. In scenario 3, our focus is to evaluate the tracking ability of cells entering and/or leaving observation image. In terms of the ant-individual-related parameters, we take the same values as in scenarios 1 and 2, while other parameters are set as in Table 4. Altogether three cases are considered for a given image sequences, namely, new cell entering, existing cell leaving and both new cell entering and existing cell leaving. In case 1, as shown in Fig. 18, the new cell 2 enters first at the lower rim of image and our algorithm captures the cell instantly when it is partially in the observation area in frame 40. Afterwards, cell 2, as well as the original cell 1, is kept on being tracked with our algorithm in the following frames. In case 2, as shown in Fig. 19, cell 2 keeps on moving left, and partially leaves observation region in frame 47 and fully in frame 48. It is noted that our algorithm can track all cells even though one cell is partially leaving the image. In case 3, as shown in Fig. 20, cell 1 moves right, partially leaves the image in frames 52 and 53, and fully leaves the image in frame 54. Meanwhile, cell 3 partially enters from left upper part of image in frames 53 and 54. It can be observed that our algorithm can track existing cells and capture new entering one instantly in frame 53. It is noted that the averaged FNR and FPR are 10.20% and 0.20%, respectively, using our algorithm, and 20.00% and 0.40%, respectively, using the multi-Bernoulli filter. It can be seen that our algorithm shows better performance than the multi-Bernoulli filter, especially between frame 42 and frame 48. In frame 53, the two algorithms record relatively higher false negative reports, and this is because both are insensitive to cell 3 that only a small part

Table 4 Parameter settings on various scenarios. Scenarios

Parameter values

1 and 2 3 4

|velomax | = 5pix/frame; = 100 ;  = 15 ;  =  = 10; threshold to merge tracks is set to be 0.6; cell size is 7 pix × 7 pix square blob; |velomax | = 2 pix/frame; = 100 ;  = 10 ;  =  = 20; threshold to merge tracks is set to be 0.1; cell size is 10 pix × 10 pix square blob; |velomax | = 2 pix/frame; = 120 ;  = 15 ;  =  = 30; threshold to merge tracks is set to be 0.1; cell size is 12 pix × 12 pix square blob;

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

461

Fig. 15. Two cells tracking with different dynamics (128 pix × 128 pix).

appears in the left upper region. To full view how ant colonies work, Fig. 21 illustrates the corresponding ant colony distributions of each extracted cell state in selected frames discussed above. It can be observed that each ant colony moves around their individual interested cell, and each is bounded within the image. Finally, we will test the tracking performance of multiple cells (up to 7 cells) using the proposed algorithm in scenario 4. Fig. 22 presents an example of successfully tracking results. It is seen that our algorithm can track simultaneously 5 or more cells in this case, and partially entering and leaving cells, such as cells 5, 6, 7 and 3, are tracked as well. According to the statistic results over 20 Monte-Carlo runs, the averaged FNR and FPR are 3.11% and 2.56%, respectively, using our algorithm, and 3.56% and 4.67%, respectively, using the multi-Bernoulli filter, which means that our algorithm achieves a robust tracking and remains competitive in our all studied cases.

4.2.3. Discussions 4.2.3.1. Issues on ant-related parameters. It can be observed that the above results are acquired based on a series of ant-related parameters, which include parameters ˛ and ˇ related to ant decision in Eq. (10), ,  and  related to ant discrimination as Eq. (11), velomax , max and  min related to ant motion in ant state update process as [35], and Lmax , Lmin related to ant colony scale. To achieve better tracking performance, our algorithm resorts to appropriate selections of these ant parameters, thereby, a multi-object multiparameter optimization problem is then raised. Although analytical determination of these values is out of the scope of this research, the corresponding empirical values in scenario 1 are determined through the following analysis. (a) Ant decision parameters ˛ and ˇ In order to investigate the effect of sensitive parameters ˛ and ˇ in Eq. (10) on tracking performances, we test different

462

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469 100

3 x-coordinate in pixel

80

2.6

40

cell1 cell2

2.4

4

6

8 Frame

10

12

No. of targets

20

2

14

140 y-coordinate in pixel

Target number The estimated target number

2.8 60

cell1 cell2

120

2.2 2 1.8 1.6

100 80

1.4

60 40

1.2

20 2

4

6

8

10

12

14

1

16

Frame

2

4

6

8

10

12

14

time step

Fig. 16. Position and the number of cells estimate using our algorithm (left and right, respectively).

combinations of the two parameters, and Table 5 illustrates the statistic results on the averaged FNR and FPR if other parameters are set as scenario 1 in Table 4 (b) Ant discrimination parameters , and  We will investigate the influence of another set of parameter combinations on tracking results, i.e., , ,  in the likelihood function in Eq. (11). As indicated in Fig. 5, the ability of ant discriminating the difference or similarity between any two regions mainly depends on the likelihood function form, as well as the corresponding values of , , . As a matter of fact, a better choice of these parameters is beneficial to the improvement of tracking performance. For simplicity, we only present the statistic results of various combinations of  and  =  with ˛ = 2 and ˇ = 5, and other parameters are set as scenario 1 in Table 4. Table 6 shows that the tracking performance is significantly affected by the above parameters’ different combinations, and inappropriate definition of likelihood function will definitely lead to large values of FNR and FPR. Furthermore, the FNR is almost inversely proportion to the FPR, and partly because detecting all possible cells of region means more false positives are generated,

which can be circumvented by the introduction of other features. We find that the combination of  = 15 and  =  = 10 is preferred for our cell image sequences. (c) Ant motion parameters velomax ,  max and  min Parameters, such as velomax , max and  min , are used to describe the ant motion range and moving ability, which can be acquired based on the prior information according to the studied cell image sequences. The parameter velomax describes the maximum speed of each cell in image, which is restricted by the maximum cell moving distance between frames, and it takes a large value (5 pix/frame) in scenario 1. Meanwhile,  max and  min determine the maximum and minimum movement abilities of each ant, respectively, and both of them are determined according to the statistic results in Table 7. It is implicated in Table 7 that both the FNR and FPR are affected by the ant moving ability, and the combination of  max = 80 and  min = 15 could achieve better tracking results. (d) Ant colony scale parameters Lmax and Lmin In various ant-related algorithms, a moderate number of ants are required and essential for achieving better performance. On the one hand, to reduce the computational burden (o((Lm)2 )), Lmax ,

Table 5 Statistic results using various ˛ and ˇ with  = 15 and  =  = 10 over 20 runs. ˛ = 0.5

ˇ=1 ˇ=2 ˇ=3 ˇ=4 ˇ=5

˛=1

˛ = 1.5

˛=2

FNR

FPR

FNR

FPR

FNR

FPR

FNR

FPR

0.0267 0.0100 0.0217 0.0467 0.0283

0.1067 0.0850 0.0867 0.1067 0.1050

0.0383 0.0383 0.0433 0.0383 0.0217

0.1067 0.0850 0.0883 0.0683 0.0933

0.0533 0.0317 0.0317 0.0483 0.0417

0.0983 0.1050 0.1183 0.0783 0.1083

0.0400 0.0633 0.0333 0.0283 0.0131

0.0983 0.1083 0.0817 0.0967 0.0062

Note that the maximum values of FNR and FPR could be reached, with lower than 7% and 12%, respectively, for various combinations of ˛ and ˇ. Besides, we find that the combination of ˛ = 2 and ˇ = 5 is a better choice due to smaller values of FNR and FPR. The numbers in bold and italic form only means that we find an acceptable and favorable FNR and FPR results for various combinatios. Table 6 Statistic results using various  and  =  with ˛ = 2 and ˇ = 5 over 20 runs.  =  = 10

=5  = 10  = 15  = 20  = 25

 =  = 20

 =  = 30

 =  = 40

FNR

FPR

FNR

FPR

FNR

FPR

FNR

FPR

0.0133 0.0133 0.0135 0.1333 0.2000

0.4667 0.3000 0.0064 0.0333 0.0333

0.0000 0.0000 0.0000 0.0000 0.0000

1.3333 1.1667 0.7000 0.6667 0.5000

0.0000 0.0000 0.0000 0.0333 0.0000

4.5000 2.6333 1.6333 1.5000 1.3000

0.0000 0.0000 0.0000 0.0000 0.0000

5.9333 5.3667 4.5000 3.0000 2.0000

The numbers in bold and italic form only means that we find an acceptable and favorable FNR and FPR results for various combinatios.

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

463

Fig. 17. Close moving cells tracking (128 pix × 128 pix). Table 7 Statistic results of tracking performance using different  max and  min over 20 runs.  max = 20

 min = 5  min = 10  min = 15  min = 20

 max = 50

 max = 80

 max = 110

FNR

FPR

FNR

FPR

FNR

FPR

FNR

FPR

0.0267 0.0283 0.0483 0.0317

0.0233 0.0183 0.0217 0.0233

0.0400 0.0233 0.0233 0.0300

0.0217 0.0333 0.0200 0.0300

0.0283 0.0283 0.0133 0.0183

0.0283 0.0100 0.0069 0.0107

0.0333 0.0283 0.0267 0.0400

0.0233 0.0217 0.0250 0.0133

The numbers in bold and italic form only means that we find an acceptable and favorable FNR and FPR results for various combinatios. Table 8 Statistic results of tracking performance using different Lmax and Lmin over 20 runs. Lmax = 150

Lmin = 50 Lmin = 100 Lmin = 150

Lmax = 200

Lmax = 250

Lmax = 300

FNR

FPR

FNR

FPR

FNR

FPR

FNR

FPR

0.0467 0.0183 0.0135

0.0056 0.0060 0.0071

0.0483 0.0131 0.0133

0.0059 0.0064 0.0068

0.0487 0.0136 0.0130

0.0084 0.0075 0.0069

0.0492 0.0132 0.0130

0.0079 0.0074 0.0067

The numbers in bold and italic form only means that we find an acceptable and favorable FNR and FPR results for various combinatios.

464

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

Fig. 18. Cells tracking when new cell enters images (100 pix × 100 pix).

a maximum of ants per hypothesized task, is introduced to ensure the required estimation accuracy; on the other hand, to search for more candidates, Lmin , a minimum of ants per hypothesized task, is defined and remained when implementation. Table 8 illustrates the statistic results on the averaged FNR and FPR if other parameters take the proposed values according to the aforementioned discussions. As we observe that both the FNR and FPR are relatively stable for the combinations of Lmax ≥ 200 and Lmin ≥ 100, while the combination of Lmax = 200 and Lmin = 100 is preferred due to less execution time required. Remark: We have demonstrated the influence of various parameter combinations on tracking performance, and some parameters, such as the ant discrimination parameters, play a significant role in

the resulting FNRs and FPRs. In our studied cell image sequences, even for a general cell tracking problem, the corresponding parameters are recommended as follows: (1) Parameters ˛ and ˇ indicate the importance between ant dynamics and ant discrimination ability. If cells move in an erratical way, we usually have ˛ > ˇ, otherwise ˛ < ˇ, which means that cells move smoothly and each ant tries to find potential cells via heuristic function. (2) Parameters , ,  in the likelihood function play a significant role in discriminating the difference of any two regions of the same size, especially for a region containing a cell and a region including clutter or being background. Therefore, an empirical method for adjusting parameters , ,  is that the value of similarity decreases linearly in proportion to cell area within a fixed region. (3) Parameter velomax

Fig. 19. Cells tracking when new cell leaves images (100 pix × 100 pix).

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

465

Fig. 20. Cells tracking when cell both leaves and enters images (100 pix × 100 pix).

indicates the maximum speed of a cell in a sequence, which is restricted by the maximum cell moving distance between frames; parameters  max and  min describe the ant influence regions at the current location. In general,  min is acquired approximately by the average diameter of the studied cells, and  max is determined roughly by the distance between the two closest neighboring cells in a sequence. (4) Parameters Lmax and Lmin denotes the number of ants allocated to each cell or task, and Lmin is approximately equal to the number of occupied pixels of a cell, and Lmax generally has

twice the number of occupied pixels of a cell, i.e., Lmax = 2Lmin , for the sake of execution efficiency. 4.2.3.2. Performance. Besides the aforementioned multi-Bernoulli filter, the PFS-based GM-PHD filter [24] and the particle filter (PF) [39], are further introduced for position tracking accuracy comparison. Since the two approaches belong to the category of “detect-before-track” technique, an extra cell detection module is required to generate detections using our previously proposed

Fig. 21. Ant colony distributions of selected frames in scenario 3.

466

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

Fig. 22. Multi-cell tracking in low contrast image sequences (200 pix × 200 pix).

hybrid cell detection algorithm [40]. As for the ground truth, we first draw the boundary of each cell within each frame, and the values of the marked pixels on cell boundaries are then set to be ones and zeros for the remaining pixels. In this way, a binary image is acquired. We continue to fill the image regions and holes, remove small objects, and produce another binary image with cell boundaries. In this way, the centroid corresponding to each cell boundary can be easily calculated as the ground truth. As for the label management, the traditional nearest neighboring data association is applied to all mentioned algorithms and performed frame by frame to obtain states of all detected cells. Without loss of generality, we use the sum of absolute position errors of all cells per frame in scenario 4 to compare tracking

accuracy. Fig. 23 gives the statistic results of the mean of position errors and the corresponding deviations over 20 runs, and it is observed that the proposed method can generate more accurate estimates in most cases compared with other three approaches. Furthermore, a ROC-like methodology, i.e., F-Measure in [41] is introduced in our work, which is capable of generating a meaningful scalar measure of performance. The F-Measure looks at the performance on object-level in precision-recall space without calculating the true negatives, and is defined as the weighted harmonic mean of the precision (P) and recall (R) metrics F=

1

v P1 + (1 − v) R1

(18)

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

467

Fig. 23. The statistic results of the mean of position errors and the corresponding deviations over 20 runs.

where R=

NTP NTP ,P = , NTP + NFN NTP + NFP

(19)

the parameter v controls the relative importance between P and R, and NTP , NFN and NFP are the number of true positives, false negatives and false positives, respectively. Fig. 24 shows the F-Measure values corresponding to various values of control parameter v. Results illustrate that our algorithm

outperforms the other three approaches when we mainly focus on the precision, i.e., corresponding to a larger v. Although both the PF and GM-PHD filter come very closely to our algorithm if the recall is mainly taken into account, the two filters heavily depends on an extra and problem-specific detection module. As a consequence, in terms of skipping detection module, our algorithm and the multi-Bernoulli filter still remain competitive for exploiting the spatiotemporal information directly from the image sequence.

468

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469

1 0.99

The value of F-Measure

0.98 0.97 0.96 0.95 PF GM-PHD filter multi-Bernoulli filter Our method

0.94 0.93 0.92 0.91 0.9

0

0.1

0.2

0.3

0.4 0.5 0.6 Parameter v

0.7

0.8

0.9

1

Fig. 24. Performance comparison based on F-Measure metric using various approaches.

5. Conclusions In this paper, we propose a novel parameter estimate algorithm in the framework of track-before-detect, i.e., a multi-task ant system, to track multiple cells in various challenging scenarios. By introducing two different levels, namely, the ant individual level and the ant cooperation level, the ant individual motion is well regulated and the best information is exchanged both within colony itself and among ant colonies. The proposed algorithm has been tested on synthetic image data as well as on several real cell image sequences. Results based on the synthetic data experiments clearly show that the superiority of the proposed algorithm over the stateof-the-art multi-Bernoulli filter in terms of the cardinality estimate and the OSPA distance. The real cell image experiments have also confirmed the validity of multi-cell tracking using the proposed algorithm in terms of tracking errors and F-Measure. In addition, the guidelines of ant-related parameter setting are given for a general cell tracking task.

Acknowledgments The authors are grateful to Dr. Dahai Cheng from NICTA in Australia for providing them with cell images for tests and also for his many helpful comments on the draft of this paper. The authors would like to thank the anonymous reviewers for their valuable comments that greatly improved the content of this paper. The authors would further thank the support by National Natural Science Foundation of China (No. 61273312), Natural Science Foundation of Higher Education Colleges in Jiangsu (14KJB510001), and the Project of Talent Peak of Six Industries (DZXX-013) as well.

References [1] M. Chicurel, Cell migration research is on the move, Science 295 (2002) 606–609. [2] E. Meijering, O. Dzyubachyk, I. Smal, W.A. van Cappellen, Tracking in cell and developmental biology, Semin. Cell Dev. Biol. 20 (2009) 894–902. [3] N.H. Nguyen, S. Keller, E. Norris, T.T. Huynh, M.G. Clemens, M.C. Shin, Tracking colliding cells in vivo microscopy, IEEE Trans. Biomed. Eng. 58 (2011) 2391–2400.

[4] M.A.A. Dewan, M.O. Ahmad, M.N.S. Swamy, Tracking biological cells in timelapse microscopy: an adaptive technique combining motion and topological features, IEEE Trans. Biomed. Eng. 58 (2011) 1637–1647. [5] X. Chen, X. Zhou, S.T.C. Wong, Automated segmentation, classification, and tracking of cancer cell nuclei in time-lapse microscopy, IEEE Trans. Biomed. Eng. 53 (2006) 762–766. [6] X. Yang, H. Li, X. Zhou, Nuclei segmentation using marker-controlled watershed, tracking using mean-shift, and Kalman filter in time-lapse microscopy, IEEE Trans. Circuits Syst. 53 (2006). [7] H. Shen, G. Nelson, S. Kennedy, D. Nelson, J. Johnson, D. Spiller, M.R.H. White, D.B. Kell, Automatic tracking of biological cells and compartments using particle filters and active contours, Chemometr. Intell. Lab. Syst. 82 (2006) 276–282. [8] A.J. Hand, T. Sun, D.C. Barber, Automated tracking of migrating cells in phasecontrast video microscopy sequences using image registration, J. Microsc. 234 (2009) 62–79. [9] J. Mansfield, K. Gossage, C. Hoyt, R. Levenson, Auto fluorescence removal, multiplexing, and automated analysis methods for in-vivo fluorescence imaging, J. Biomed. Opt. 10 (2005) 41207. [10] O. Dzyubachyk, W.A.v. Cappellen, J. Essers, W.J. Niessen, E. Meijering, Advanced level-set-based cell tracking in time-lapse fluorescence microscopy, IEEE Trans. Med. Imaging 29 (2010) 852–867. [11] T.F. Chan, L.A. Vese, Active contours without edges, IEEE Trans. Image Process. 10 (2001) 266–277. [12] C. Zimmer, J.-C. Olivo-Marin, Coupled parametric active contours, IEEE Trans. Pattern Anal. Mach. Intell. 27 (2005) 1838–1842. [13] C. Xu, J.L. Prince, Snakes, shapes, and gradient vector flow, IEEE Trans. Image Process. 7 (1998) 359–369. [14] N. Ray, S.T. Acton, K. Ley, Tracking leukocytes in vivo with shape and size constrained active contours, IEEE Trans. Med. Imaging 21 (1998) 1222–1235. [15] D.P. Mukherjee, N. Ray, S.T. Acton, Level set analysis for leukocyte detection and tracking, IEEE Trans. Image Process. 13 (2004) 562–572. [16] O. Debeir, P.V. Ham, C. Decaestecker, Tracking of migrating cells under phasecontrast video microscopy with combined mean-shift processes, IEEE Trans. Med. Imaging 24 (2005) 697–711. [17] L.L. Ong, J. Dauwels, M.H. Ang Jr., H.H. Asada, A Bayesian filtering approach to incorporate 2D/3D time-lapse confocal images for tracking angiogenic sprouting cells interacting with the gel matrix, Med. Image Anal. 18 (2014) 211–227. [18] O. Demirel, I. Smal, W.J. Niessen, E. Meijering, I.F. Sbalzarini, IET Conference on on Piecewise Constant Sequential Importance Sampling for Fast Particle Filtering, Data Fusion & Target Tracking 2014: Algorithms and Applications (DF&TT 2014), 2014, pp. 1–8. [19] F. Li, X. Zhou, S. Wong, Multiple nuclei tracking using integer programming for quantitative cancer cell cycle analysis, IEEE Trans. Med. Imaging 29 (2010) 96–105. [20] B.-N. Vo, S. Singh, A. Doucet, Sequential Monte Carlo methods for multi-target filtering with random finite sets, IEEE Trans. Aerosp. Electron. Syst. 41 (2005) 1224–1245. [21] B.-N. Vo, W.-K. Ma, The Gaussian mixture probability hypothesis density filter, IEEE Trans. Signal Process. 54 (2006) 4091–4104. [22] B.-N. Vo, B.-T. Vo, N.-T. Pham, D. Suter, Joint detection and estimation of multiple objects from image observations, IEEE Trans. Signal Process. 58 (2010) 5129–5141.

B. Xu et al. / Applied Soft Computing 35 (2015) 449–469 [23] R. Hoseinnezhad, B.-N. Vo, B.-T. Vo, D. Suter, Visual tracking of numerous targets via multi-Bernoulli filtering of image data, Pattern Recognit. 45 (2012) 3625–3635. [24] R.R. Juang, A. Levchenko, P. Burlina, Tracking cell motion using GM-PHD, in: Sixth IEEE International Conference on Symposium on Biomedical Imaging, 2009, pp. 1154–1157. [25] M. Dorigo, V. Maniezzo, A. Colorni, The ant system: optimization by a colony of cooperating agents, IEEE Trans. Syst. Man Cybern. B 26 (1996) 29–42. [26] M. Dorigo, L.M. Gambardella, Ant colony system: a cooperative learning approach to the traveling salesman problem, IEEE Trans. Evol. Comput. 1 (1997) 53–66. [27] J.E. Bell, P.R. McMullen, Ant colony optimization techniques for the vehicle routing problem, Adv. Eng. Inf. 18 (2001) 41–48. [28] C. Blum, Beam-ACO-hybridizing ant colony optimization with beam search: an application to open shop scheduling, Comput. Oper. Res. 32 (2005) 1565–1591. [29] J.-L. Deneubourg, S. Goss, N. Franks, The dynamics of collective sorting robotlike ants and ant-like robots, First Int. Conf. Simul. Adapt. Behav. Anim. Anim. 10 (1991) 356–365. [30] P.M. Kanade, L.O. Hall, Fuzzy ants and clustering, IEEE Trans. Syst. Man Cybern. A 37 (2007) 758–769. [31] J. Handl, J. Knowles, M. Dorigo, Strategies for the increased robustness of antbased clustering, LNCS 2977 (2004) 90–104.

469

[32] M.D. Backer, R. Haesen, Classification with ant colony optimization, IEEE Trans. Evol. Comput. 11 (2007) 651–665. [33] L. Nolle, On a novel ACO-Estimator and its application to the target motion analysis problem, Knowl. Based Syst. 21 (2008) 225–231. [34] B. Xu, Q. Chen, Z. Wang, Ant estimator with application to target tracking, Signal Process. 90 (2010) 1496–1509. [35] B. Xu, Q. Chen, Z. Wang, A novel estimator with moving ants, Simul. Model. Pract. Theory 17 (2009) 1663–1677. [36] B. Xu, H. Xu, J. Zhu, Ant clustering PHD filter for multiple-target tracking, Appl. Soft Comput. 11 (2011) 1074–1086. [37] B. Xu, Q. Chen, M. Lu, P. Zhu, Two ant decision levels and its application to multi-cell tracking, LNCS 7928 (2013) 288–296. [38] D. Schuhmacher, B.-T. Vo, B.-N. Vo, A consistent metric for performance evaluation of multi-object filters, IEEE Trans. Signal Process. 56 (2008) 3447–3457. [39] I. Smal, K. Draegestein, N. Galjart, W. Niessen, E. Meijering, Particle filtering for multiple object tracking in dynamic fluorescence microscopy images application to microtubule growth analysis, Med. Imaging 27 (2008). [40] M. Lu, B. Xu, A. Sheng, Cell automatic tracking technique with particle filter, in: Adv. Swarm Intell., 2012, pp. 589–595. [41] N. Lazarevic-McManus, J.R. Renno, D. Makris, G.A. Jones, An object-based comparative methodology for motion detection based on the F-Measure, Comput. Vision Image Underst. 111 (2008) 74–85.