# A Survey of Miss-Ratio Curve Construction Techniques

Daniel Byrne  Michigan Technological UniversityHoughtonMichigan
<EMAIL_ADDRESS>
###### Abstract.

Miss-ratio curve (MRC), or equivalently hit-ratio curve (HRC), construction
techniques have recently gathered the attention of many researchers. Recent
advancements have allowed for approximating these curves in constant time,
allowing for online working-set-size (WSS) measurement. Techniques span the
algorithmic design paradigm from classic dynamic programming to artificial
intelligence inspired techniques. Our survey produces broad classification of
the current techniques primarily based on _what_ locality metric is being
recorded and _how_ that metric is stored for processing.

Applications of theses curves span from dynamic cache partitioning in the
processor, to improving block allocation at the operating system level. Our
survey will give an overview of the historical, exact MRC construction
methods, and compare them with the state-of-the-art methods present in today’s
literature. In addition, we will show where there are still open areas of
research and remain excited to see what this domain can produce with a strong
theoretical background.

##  1\. Introduction

From the working set theory proposed by Denning in 1968 (Denning, 1968) to
recent advancements in cache locality theory (Ding and Xiang, 2012; Hu et al.,
2016) modeling data locality has been essential to cache design and OS
resource allocation choices. In the past decade there has been a large
emphasis on using these models to construct plots relating miss rate to cache
size, called miss-ratio curves (MRCs) or equivalently hit-ratio curves (HRCs).

Miss ratio curves have proven to be extremely useful in estimating _how much_
data is being used by a particular workload, known as the working set size.
Knowing the working set size of a workload helps estimate the _utility_ of
adding more memory or increasing the cache size. An example miss-ratio curve
is shown in Figure 1. The knee-point in this particular MRC is just over cache
size 3000 for this workload. Notice how that an increase cache size to 4000
yields no significant decrease in miss-rate.

Figure 1. Hit Ratio Curve for two common synthetic web application workloads.
ETC models the Facebook workload from (Atikoglu et al., 2012). It contains 50
million accesses (GETs) to 7 million unique objects. PSA is the D’Carra
workload from (Carra and Michiardi, 2014), it contains 50 million accesses to
7 million unique objects following an Pareto distribution.

Given the recent theoretic and implementation advancements in miss-ratio curve
construction, we find a survey on these methods to be useful update to the
community and present some identified open lines of research yet to be fully
explored in this domain.

##  2\. Taxonomy

In this section we present a taxonomy of the current miss-ratio curve
construction techniques based on _which_ metric is recorded, _how_ the data is
stored, and the _derivation method_ of the miss-ratio curve.

###  2.1. Definitions

In the cache modeling literature there are variations on meaning of terms such
as _reuse distance_. To avoid this, we give definitions for the following
terms used in our taxonomy:

  * •

Stack Distance - The number of _unique_ accesses between two accesses to the
same piece of data. This has also been referred to as reuse distance.

  * •

Reuse Time - The _total_ number of accesses between two accesses to the same
piece of data. This is consistent with (Hu et al., 2016; Hu et al., 2015;
Zhong et al., 2009; Xiang et al., 2011).

  * •

Buckets - A set of variable sized bins that contain access to a given piece of
data.

  * •

Access Counters - Hardware or software monitors that count unique or total
accesses to a piece of data.

###  2.2. Taxonomy of MRC Constructions

[.Metric [blank]; [.Data Structure [blank]; [.MRC Derivation [blank]; [.Rep.
Work ]]]][.Stack Distance [.Buckets [.Stack Distance Distr. MIMIR (Saemundsson
et al., 2014) ] [.Hill Climbing Cliffhanger (Cidon et al., 2016) ] ] [.Access
Counters [.Stack Processing Counter Stacks (Wires et al., 2014) UMON (Qureshi
and Patt, 2006) ] ] [.Tree [.Search Tree Olken (Olken, 1981) ] [.Interval Tree
SHARDS (Waldspurger et al., 2015) ]]] Figure 2. Methods based on Stack
Distance [.Metric [blank]; [.Data Structure [blank]; [.MRC Derivation [blank];
[.Rep. Work ]]]][.Reuse Time [.Histogram [.Average Eviction AET (Hu et al.,
2016) ] [.Footprint LAMA (Hu et al., 2015) ] [.Hit-Evict Distribution Beyond
LRU (Beckmann and Sanchez, 2016) ]] [.Access Counters [.Expected Reuse
Statcache (Berg and Hagersten, 2004) Statstack (Eklov and Hagersten, 2010) ]]]
Figure 3. Methods based on Reuse Time

Figure 2 and 3 classify the techniques that use stack distance and reuse time
respectively, by how they store and process the data into a miss-ratio curve.
We also give a representative work from each class.

##  3\. Stack Distance Algorithms

This class of algorithms for calculating miss-ratio curves takes advantage of
the least recently used (LRU) replacement policy which maintains the stack
order. Table 1, demonstrates this for a set of requests.

In order to generate a miss-ratio curve for a given cache size, we create a
_distribution_ of stack distances. Once we have a distribution of the stack
distances, we can calculate the number of hits in a cache of size $d$ as the
sum of all the references with stack distances $\leq d$ and equally the number
of misses as the sum of all references with stack distance $\geq d$.

request | time | stack | stack dist.
---|---|---|---
a | 1 | a | $\infty$
b | 2 | b,a | $\infty$
c | 3 | c,b,a | $\infty$
d | 4 | d,c,b,a | $\infty$
a | 5 | a,d,c,b | 3 (d,c,b)
a | 6 | a,d,c,b | 0 (none)
d | 7 | d,a,c,b | 1 (a)
b | 8 | b,d,a,c | 3 (d,a,c)
Table 1. Example stack distance calculation for a set of requests

Formally, let cache size be $d$, $N$ is the total number of requests, and $M$
is the total number of unique requests, then the $missrate(d)$ is express as:

$(1-\frac{\sum\limits_{i=1}^{d}freq(i)}{N})=\frac{\sum\limits_{i=d}^{M}freq(i)}{N}$

It is important to note that some represent this as an integration of the
probability density function given by the distribution. The reason is to
highlight the order relationship between stack distance and miss-rate, that
is, miss-rate is an order higher than stack distance. This result was proved
in Ding’s Higher Order Theory of Locality (HOTL) (Ding and Xiang, 2012).

###  3.1. Mattson’s Stack Distance Algorithm

Mattson et. al. (Mattson et al., 1970) developed the first algorithm to use
stack distance as a method for calculating the miss-rate for a given cache
size. It has five main steps for each reference $x$:

  1. (1)

Search stack to find the current location of $x$ (if any)

  2. (2)

Calculate the distance $d$ from the top of the stack

  3. (3)

Update cache miss counters for all sizes $\geq d$

  4. (4)

Update cache hit counters for all sizes $\leq d$

  5. (5)

Push $x$ onto the stack and remove the last location of $x$ (if any)

The algorithm runs in $O(N*M)$ time, where $N$ is total number of references
and $M$ is unique number of references. Many works improve on this algorithm
by using more efficient data structures, sampling, and estimated stack
distances.

Fang et. al. (Fang et al., 2005) showed that the stack distance distribution
of a program’s typical workload could be predicted accurately given a
significantly smaller version of the workload. They predicted the stack
distance distribution for the _reference_ input set of the SPEC2K benchmark
using only the _train_ and _test_ input sets. These results encouraged
researchers to move towards sampling the references in order to build the
distribution while still achieving accurate stack distance distributions.

###  3.2. Tree

In tree-based stack distance calculations we exploit the logarithmic search
time in balanced tree structure. Methods in this class have developed towards
approximation of distances by getting to a leaf node that represents a close
enough stack distance.

####  Partial Sums Tree by Bennet et. al. (Bennett and Kruskal, 1975)

Bennet introduces a partial sums tree based on binary vectors representing the
stack distance of a reference. But first they introduce the concept of using a
hash table to store the _previous_ location $p$ of a reference. Second, they
maintain a binary vector, $B$, at time $t$ for a given reference $x$ with the
following property:

$B_{t}=\begin{cases}1&\text{if $x_{p}\neq x_{i}(i=p,p+1,...,t)$}\\\
0&\text{otherwise.}\end{cases}$

Now the number of 1s in $B_{t}$ is the stack distance of reference $x$. As $t$
increases the idea is to create a hierarchy of partial sums out of $B$ so that
we don’t have to scan $B$ each time we calculate the stack distance for $x$.
For example, imagine that every 3 time intervals we add up all the 1s counted
so far and store them as a partial sum $B^{i}$, then for every 3 partial sums
$B^{i}$ we make a new final sum. To calculate the stack distance now, we do a
tree traversal from our current time $t$ to the previous use time $p$.

Since each lookup in the hash table is $O(1)$ and each stack distance
calculation is $O(log(n))$, then for a trace of $N$ references the running
time is now $O(N*log(n))$.

####  Interval tree by Almasi et. al. (Almási et al., 2002)

Almasi interprets Bennet’s binary vectors in a different manner. They use only
the 0s from the vector, calling them holes. The calculation for stack distance
now becomes:

$stackdist(x)=t-p-holes(x_{p})$

Where $t$ is current time, $p$ is previous time, and $holes(x_{p})$ represents
the number of holes (0s) between $t$ and $p$ in $B$. Now we represent a hole’s
index as $h_{i}$ and let the interval $[h_{i},h_{j}]$ represent the time
between two consecutive holes.

Since these intervals of holes are non-overlapping, they can be represented as
AVL or red-black trees. For $N$ references this takes $O(N*log(n))$, similar
to Bennet’s asymptotic bound, but in practice it performs significantly better
because the tree is well balanced. Although we are still at the point of using
a large amount of memory, $O(M)$, where $M$ is unique number of accesses.

####  Search Tree (Olken, 1981)

Olken proposed ordering the LRU stack as a binary tree based on access times
so that an in-order traversal gives the LRU stack position. For a given
reference at the root, the right side of the tree represents references used
more recently and the left side represents older references. At each node in
the tree, we store a count of the number of nodes in the right subtree (more
recent) and the number of nodes in the left subtree (older).

In order to calculate a stack distance, on a reference to $x$, we have to
search the tree for its location and traverse our way to the top of the tree,
recording the number of nodes in the right subtree of each node that we
visited. In order to keep this tree balanced we should use an AVL or Splay
tree as suggested by Olken. This algorithm runs in $O(N*log(n))$ time and
needs to store $M$ number of nodes, where $M$ is the number of unique
accesses.

####  Scale Tree (Zhong et al., 2009)

In Zhong et. al. they produce a modified version of Olken’s search tree by
modifying it to support a _range_ of times at each node. This maps several
references to one node and allows for much faster searches. The amount of
error incurred by using a range of times at each node is bounded by the range.
The overall run time is reduced to $O(N*log(log(M))$, but we still require
$O(M)$ space.

####  SHARDS (Waldspurger et al., 2015)

Spatially hashed reuse distances (stack distances), or SHARDS, is a recent
advancement to constant space use. They achieve this by sampling accesses
based on the hashed value of their location, hence spatially hashed. These
samples are then put into an interval tree to compute their stack distance.
After computing stack distance, it is recorded in the stack distance histogram
so that we can solve for the miss-ratio curve. The achievement in this work is
running in $O(1)$ through sampling and linear run time since it only requires
one scan through the set of accesses.

###  3.3. Buckets

These algorithms employ some variant on reducing the number of stack positions
that exist in the LRU stack. The key insight is to reduce $n$ stack positions
into $m$ buckets while still maintaining the LRU stack order.

####  MIMIR by Saemundsson et. al. (Saemundsson et al., 2014)

The MIMIR algorithm created fixed number of buckets to store references. Each
bucket $B$ maintained the LRU stack order property, that is $B_{i}$ has lower
stack distance than $B_{j}$ since $i<j$. While buckets can vary in size, they
introduce an _aging_ procedure to ensure that buckets stay balanced. A common
metric used to resize the buckets is the _average_ stack distance.

The miss-ratio curve is generated from the estimated stack distance
distribution. For any bucket we know that the actual stack distance lies
somewhere between the sum of the size of the buckets before $B_{i}$, call it
$n$ and $n+sizeof(B_{i})$. Now we create the normal stack distance histogram
from the intervals and perform the integration to get the miss rate at for a
cache size (stack distance) $d$.

####  Cliffhanger by Cidon et. al. (Cidon et al., 2016)

The Cliffhanger approach uses the same MIMIR variable bucket algorithm, but
their insight is to use _shadow queues_. Shadow queues can be thought of
victim caches but without the value of the data, just the _key_. For example,
if we were to a miss in a cache with size 10, but then hit in the shadow queue
of size 20 then that would correspond to a hit with cache size 30. Therefore,
in the MIMIR algorithm each bucket corresponds to a shadow queue of size $n$.

In order to solve for the miss-ratio curve, they use _hill climbing_ to
incrementally build a miss-ratio curve. Hill climbing is a technique to find
the local slope of a curve. Recall that a hit in the shadow queue corresponds
to a miss in the cache, hence and increase in the miss-rate for that given
cache size. Therefore, we can increase the size of this cache a small amount
in order to combat future misses. In resource partitioning, this also means
decreasing the size of a randomly chosen _different_ cache. This process
continues until there is no overall improvement in the total hit rate of the
system.

###  3.4. Access Counters

These methods either use: (1) a set of counters to represent stack distances
at each cache size or (2) a set of counters to represent _total_ accesses and
_unique_ accesses.

####  Powers of 2 by (Kim et al., 1991)

First purposed by Kim et. al. in 1991, they use Mattson’s algorithm with two
key differences:

  * •

Counters are used representing cache sizes of $2^{k}$

  * •

Use a hash table to record previous stack location and stack distance for each
reference

Therefore on reference $x$, we lookup its position and stack distance $i$.
Then increment the counter corresponding to cache size $i$, this means that
$x$ would hit in all caches size $>i$. Now we have to push $x$ on the top of
the stack and update stack distances for all references in the stack, the
worst case on this is $M$. Therefore the running time is in the same
asymptotic bound, $O(N*M)$ as Mattson’s, but in practice it runs faster since
it is insensitive to a program’s average stack distance.

This algorithm has been applied to build miss-ratio curves in Zhou et. al.
(Zhou et al., 2004). They present a hardware and software implementation for
calculating page miss-ratio curves to guide memory management in a multi-
program environment, a known NP-hard problem.

At the hardware level, Tam et. al. (Tam et al., 2009) have used this algorithm
to manage L2 cache partitioning. Here they create an access trace to L2 by
recording the misses of L1, then they send the access trace log into the
algorithm given by Kim. Their miss-ratio curve is then used to guide how much
of the L2 cache should an application receive.

####  CounterStacks (Wires et al., 2014)

This recent development by Wires et. al. uses probabilistic counters
(HyperLogLogs) and Bloom filters (for set inclusion) in order to maintain a
list of unique accesses over a period of time. Specifically, they build a list
for time $t$ that contains the number of unique accesses. This allows us to
keep track of how much a counter increased (stack distance) over a number of
non-distinct accesses. An example is given in Table 2.

In order to derive the stack distance, we need to look at the _intra-counter_
change, that is to find the last reference, we must find the newest counter.
In Table 2 we find the last reference to $a$ is at position 1, hence the stack
distance of $a$ lies at position (1,4) = 3.

Table 2. Counter Stacks a | b | c | a
---|---|---|---
1 | 2 | 3 | 3
| 1 | 2 | 3
|  | 1 | 2
|  |  | 1

CounterStacks uses $O(N*log(M))$ to process a set of references, since we
still need to find the last reference to each element. But CounterStacks
achieves a great reduction in space complexity, $O(log(M))$ through its use of
probabilistic counters which set a small bound on the overall error in stack
distance calculation.

####  UMONs (Qureshi and Patt, 2006)

Qureshi et. al. developed utility monitor circuits (UMONs) that approximate
hit rate vs. number of sets in an LRU cache. They use access counters for a
fixed number of stack positions for each set. Since the cache obeys the LRU
policy, we can get the hit rate for a given number of sets by reading the
counters for $n$ sets.

This process is used in partitioning caches by set per application. If
application A has a higher utility for a number of sets (i.e. its hit rate
curve is not as steep) then it will be assigned that number of sets.
Application B will be left with the remaining number of sets.

####  Fractals (He et al., 2012)

L. He et. al. use a fractal equation to calculate the miss-rate of a given
cache size from the time between two consecutive misses. This is known as the
inter-miss gap described by Denning and related to LRU miss-rate by Ding in
HOTL (Ding and Xiang, 2012). From the distribution of these inter-miss gaps we
can get the LRU miss rate for a given cache size, by the following:

$intermiss(c)=\frac{1}{missrate(c)}$

Their collection inter-miss gaps used hardware counters to record a cache
miss, the number of references $n$, and a new cache miss. $n$ is then saved
into a buffer. This model allows for online processing at the cost of a 2%
slowdown.

Unfortunately, their model ignores cold misses which leads to inaccuracy (the
cold miss equation is not fractal). On average in the SPEC2006 benchmark they
report 76% percent miss-ratio curve accuracy.

####  PARDA (Niu et al., 2012)

Q. Niu et. al introduce the first explicit parallel stack distance algorithm
based on Mattson’s stack processing. PARDA exploits the independences of the
_define-use_ chains in references. That is, since we only care about the
reference to the last location, all previous references are independent.
Therefore, we can create chunks that contain all references between the last
use and the current use for a given access. The algorithm is outlined below:

  1. (1)

Start at the end of the trace with reference $x$

  2. (2)

Scan until you see another $x$

  3. (3)

Send that interval off for processing via Mattson’s algorithm

  4. (4)

Repeat with the next reference after $x$

##  4\. Reuse Time Algorithms

More recently, reuse time has emerged as a metric used to approximate miss-
ratio curves. In 1968, Denning proposed that the reuse time distribution can
yield the working set size at time $t$, but until recently it was hard to
achieve accurate results with the analytical models available.

###  4.1. Reuse Time

Table 3 gives an example of the reuse time calculation. Note we do not need to
maintain the LRU stack, but as a result the miss-rate for a given cache size
is no longer an exact integration over the distribution of reuse times.

request | time | last use | reuse time
---|---|---|---
a | 1 | $\infty$ | $\infty$
b | 2 | $\infty$ | $\infty$
c | 3 | $\infty$ | $\infty$
d | 4 | $\infty$ | $\infty$
a | 5 | 1 | 4 (5-1)
a | 6 | 5 | 1 (6-5)
d | 7 | 4 | 3 (7-4)
b | 8 | 2 | 6 (8-2)
Table 3. Example reuse time calculation for a set of requests

###  4.2. Histogram

These methods bin the data recorded into a histogram, and build a miss-ratio
curve from that histogram by integrating the probability density function or
summing the total frequencies over a given range.

####  Average Eviction Time (Hu et al., 2016)

Hu et. al. presents as series of kinetic equations related to average data
eviction in the cache. Based on the reuse time histogram, they estimate the
probability that reference $x$ has reuse time greater than $t$, which is then
related to a stack movement. From here we can solve for the _average_ eviction
time for a given cache size. The average eviction time for a cache now relates
miss-rate to time, so we have the following equation:

$\int_{0}^{AET(c)}P(t)dt=c$

Where $P(x)$ is the probability that a reference has reuse time greater than
$t$.

In order to gain constant space complexity, they perform sampling over random
random intervals, keeping the amount of data in the reuse time histogram
constant with respect to the number of requests.

####  Footprint (Xiang et al., 2011)

Xiang et. al. introduces the _footprint_ function as a mapping from an
execution window to the volume of references during that time. For a given
number of references $n$ there are $\binom{n}{2}$ number of execution windows.
Therefore, the _average_ footprint is measured for a set of references. This
can be done by measuring the following:

  * •

Reuse Time Histogram

  * •

The first use of every reference

  * •

The last use of every reference

The miss-ratio curve under the footprint metric is constructed for a given
cache size $d$ as the fraction of reuse times that have an average footprint
smaller than $d$. In HOTL, it was shown that the footprint metric is the cache
locality theory equivalent of Denning’s working set size.

The footprint algorithm runs in $O(N)$ and $O(M)$ where $N$ is the total
number of accesses and $M$ is the total unique accesses. The algorithm was
recently applied in LAMA (Locality Aware Memory Allocation) (Hu et al., 2015),
to guide `memcached`’s slab reallocation problem.

####  Hit-Evict Distribution (Beckmann and Sanchez, 2016)

Beckmann et. al. introduces two new histograms of information, hit and evict.
The hit distribution is the frequency of hits for given reuse times. Likewise,
the evict distribution is the frequency of evicts for given reuse times. The
motivation behind this work is to reduce the dependence on the LRU stack
property that many analytic models assume. In their work they create an _age_
distribution based off of the reuse times. This age distribution when combined
with the hit and evict distribution can give an accurate LRU model for cache
miss rates.

###  4.3. Access Counters

These methods use counters to record reuse times for a particular set of
references, these counters are then aggregated to form a reuse time
distribution.

####  Statcache (Berg and Hagersten, 2004)

Statcache by Berg et. al. uses random sampling to select some memory
references. For each memory reference it does the following:

  1. (1)

Set a (hardware) watchpoint on the address and record current number of memory
accesses at that time $n$.

  2. (2)

On the watchpoint trap read the new current number of memory access $N$. The
reuse time of this address is now $N-n$.

  3. (3)

Update the reuse time distribution with the new reuse time.

The Statcache model assumes a cache with random replacement, because of this
we can assume that the probability that a cache line is still in a cache after
a cache miss is uniform. Since the model calculates the probability that a
reference will cause a cache miss, by assuming miss-ratio does not change over
time, we can fix miss-rate and derive the probability that a number of
references cause a miss from its reuse distance from their reuse time.

####  Statstack (Eklov and Hagersten, 2010)

The Stackstack model builds off of the Statcache work, this time modeling an
LRU cache. They define the _expected_ stack distance of a reference with reuse
time $t$ to be the _average_ stack distance of all references with reuse time
$t$. The miss-ratio curve is then constructed by computing the expected stack
distance of each reuse time weighted by their frequency, this gives us a stack
distance distribution.

If one were to interpret the Statstack model’s relationship between average
stack distance vs. time, then we would arrive at the kinetic equation present
in the average eviction model (AET). Stackstack performs in the same time
bound as the AET model, $O(N)$ time and $O(1)$ space due to their sampling
techniques.

##  5\. Open Problems

We have described the history of the miss-ratio curve construction techniques
that rely on reuse time and stack distance measurements. In a broader scope,
we have the enormous amount of research done in the cache behavior modeling
domain. Fixing our limit on constructing practical miss-ratio curves and their
applications, we find that results in the following areas would excite many
researchers in this field.

  * •

To what extent can we relate logical access time to physical clock time?
Currently, many models assume logical access time, that is each reference is
considered a unique point in time and there are no inactive periods. However,
in real systems it is often the case where the number of references per unit
of time vary significantly. An investigation into how much error is produced
when using physical clock time vs. the overhead caused by measuring logical
time would give insights on how to better formulate these models for use on
real systems.

  * •

How to extend these models to non-LRU caches?
While some models could be used for alternative cache policies, there is a
large reliance on the nice LRU stack property, this is especially painful
since Intel’s switch to RRIP policy has made cache modeling significantly more
difficult. In latest work by Beckmann et. al. (Beckmann and Sanchez, 2016),
they claim a general model, and that policy specific models will be the
subject of future work.

  * •

Despite the above remarks, the theory remains mature and well understood. How
can we leverage these new models in day-to-day systems?
The implementation of these models leaves much to be desired, rather high
overhead has stopped many of these miss-ratio curve guided techniques from
making it to practice. What improvements can we expect with miss-ratio curve
constructions being done in $O(1)$ space in an online and dynamic workload
environment?

##  6\. Conclusion

Our survey has covered the two main metrics, stack distance and reuse time,
and their associated cache models from a miss-ratio curve construction
standpoint. In Denning’s 1968 work (Denning, 1968), he defined page residency
as how long a page will stay in main memory. In the cache modeling field, we
have called it the stack distance. More recently, it has been called the
average eviction time of a cache. You could say that modeling locality has
been reused over and over again.

## References

  * (1)
  * Almási et al. (2002) George Almási, Cǎlin Caşcaval, and David A. Padua. 2002.  Calculating Stack Distances Efficiently.  SIGPLAN Not. 38, 2 supplement (June 2002), 37–43.  https://doi.org/10.1145/773039.773043
  * Atikoglu et al. (2012) Berk Atikoglu, Yuehai Xu, Eitan Frachtenberg, Song Jiang, and Mike Paleczny. 2012.  Workload analysis of a large-scale key-value store. In ACM SIGMETRICS Performance Evaluation Review, Vol. 40. ACM, 53–64.
  * Beckmann and Sanchez (2016) Nathan Beckmann and Daniel Sanchez. 2016.  Modeling cache performance beyond lru. In High Performance Computer Architecture (HPCA), 2016 IEEE International Symposium on. IEEE, 225–236.
  * Bennett and Kruskal (1975) Bryan T Bennett and Vincent J. Kruskal. 1975.  LRU stack processing.  IBM Journal of Research and Development 19, 4 (1975), 353–357.
  * Berg and Hagersten (2004) Erik Berg and Erik Hagersten. 2004.  StatCache: a probabilistic approach to efficient and accurate data locality analysis. In Performance Analysis of Systems and Software, 2004 IEEE International Symposium on-ISPASS. IEEE, 20–27.
  * Carra and Michiardi (2014) Damiano Carra and Pietro Michiardi. 2014.  Memory partitioning in memcached: An experimental performance analysis. In Communications (ICC), 2014 IEEE International Conference on. IEEE, 1154–1159.
  * Cidon et al. (2016) Asaf Cidon, Assaf Eisenman, Mohammad Alizadeh, and Sachin Katti. 2016.  Cliffhanger: scaling performance cliffs in web memory caches. In 13th USENIX Symposium on Networked Systems Design and Implementation (NSDI 16). 379–392.
  * Denning (1968) Peter J Denning. 1968\.  The working set model for program behavior.  Commun. ACM 11, 5 (1968), 323–333.
  * Ding and Xiang (2012) Chen Ding and Xiaoya Xiang. 2012.  A higher order theory of locality. In Proceedings of the 2012 ACM SIGPLAN Workshop on Memory Systems Performance and Correctness. ACM, 68–69.
  * Eklov and Hagersten (2010) David Eklov and Erik Hagersten. 2010.  StatStack: Efficient modeling of LRU caches. In Performance Analysis of Systems & Software (ISPASS), 2010 IEEE International Symposium on. IEEE, 55–65.
  * Fang et al. (2005) Changpeng Fang, S Can, Soner Onder, and Zhenlin Wang. 2005\.  Instruction based memory distance analysis and its application to optimization. In Parallel Architectures and Compilation Techniques, 2005. PACT 2005. 14th International Conference on. IEEE, 27–37.
  * He et al. (2012) L. He, Z. Yu, and H. Jin. 2012.  FractalMRC: Online Cache Miss Rate Curve Prediction on Commodity Systems. In 2012 IEEE 26th International Parallel and Distributed Processing Symposium. 1341–1351.  https://doi.org/10.1109/IPDPS.2012.121
  * Hu et al. (2015) Xiameng Hu, Xiaolin Wang, Yechen Li, Lan Zhou, Yingwei Luo, Chen Ding, Song Jiang, and Zhenlin Wang. 2015\.  LAMA: Optimized locality-aware memory allocation for key-value cache. In 2015 USENIX Annual Technical Conference (USENIX ATC 15). 57–69.
  * Hu et al. (2016) Xiameng Hu, Xiaolin Wang, Lan Zhou, Yingwei Luo, Chen Ding, and Zhenlin Wang. 2016\.  Kinetic Modeling of Data Eviction in Cache. In 2016 USENIX Annual Technical Conference (USENIX ATC 16).
  * Kim et al. (1991) Yul H. Kim, Mark D. Hill, and David A. Wood. 1991.  Implementing Stack Simulation for Highly-associative Memories. In Proceedings of the 1991 ACM SIGMETRICS Conference on Measurement and Modeling of Computer Systems (SIGMETRICS ’91). ACM, New York, NY, USA, 212–213.  https://doi.org/10.1145/107971.107995
  * Mattson et al. (1970) Richard L. Mattson, Jan Gecsei, Donald R. Slutz, and Irving L. Traiger. 1970.  Evaluation techniques for storage hierarchies.  IBM Systems journal 9, 2 (1970), 78–117.
  * Niu et al. (2012) Q. Niu, J. Dinan, Q. Lu, and P. Sadayappan. 2012\.  PARDA: A Fast Parallel Reuse Distance Analysis Algorithm. In 2012 IEEE 26th International Parallel and Distributed Processing Symposium. 1284–1294.  https://doi.org/10.1109/IPDPS.2012.117
  * Olken (1981) Frank Olken. 1981\.  Efficient methods for calculating the success function of fixed-space replacement policies.  Technical Report. Lawrence Berkeley Lab., CA (USA).
  * Qureshi and Patt (2006) Moinuddin K. Qureshi and Yale N. Patt. 2006.  Utility-Based Cache Partitioning: A Low-Overhead, High-Performance, Runtime Mechanism to Partition Shared Caches. In Proceedings of the 39th Annual IEEE/ACM International Symposium on Microarchitecture (MICRO 39). IEEE Computer Society, Washington, DC, USA, 423–432.  https://doi.org/10.1109/MICRO.2006.49
  * Saemundsson et al. (2014) Trausti Saemundsson, Hjortur Bjornsson, Gregory Chockler, and Ymir Vigfusson. 2014.  Dynamic Performance Profiling of Cloud Caches. In Proceedings of the ACM Symposium on Cloud Computing (SOCC ’14). ACM, New York, NY, USA, Article 28, 14 pages.  https://doi.org/10.1145/2670979.2671007
  * Shi et al. (2009) Xudong Shi, Feiqi Su, Jih-Kwon Peir, Ye Xia, and Zhen Yang. 2009.  Modeling and stack simulation of CMP cache capacity and accessibility.  IEEE Transactions on Parallel and Distributed Systems 20, 12 (2009), 1752–1763.
  * Tam et al. (2009) David K. Tam, Reza Azimi, Livio B. Soares, and Michael Stumm. 2009\.  RapidMRC: Approximating L2 Miss Rate Curves on Commodity Systems for Online Optimizations. In Proceedings of the 14th International Conference on Architectural Support for Programming Languages and Operating Systems (ASPLOS XIV). ACM, New York, NY, USA, 121–132.  https://doi.org/10.1145/1508244.1508259
  * Waldspurger et al. (2015) Carl A Waldspurger, Nohhyun Park, Alexander T Garthwaite, and Irfan Ahmad. 2015.  Efficient MRC Construction with SHARDS.. In FAST. 95–110.
  * Wires et al. (2014) Jake Wires, Stephen Ingram, Zachary Drudi, Nicholas JA Harvey, Andrew Warfield, and Coho Data. 2014.  Characterizing Storage Workloads with Counter Stacks.. In OSDI. 335–349.
  * Xiang et al. (2011) X. Xiang, B. Bao, C. Ding, and Y. Gao. 2011\.  Linear-time Modeling of Program Working Set in Shared Cache. In 2011 International Conference on Parallel Architectures and Compilation Techniques. 350–360.  https://doi.org/10.1109/PACT.2011.66
  * Zhong et al. (2009) Yutao Zhong, Xipeng Shen, and Chen Ding. 2009.  Program Locality Analysis Using Reuse Distance.  ACM Trans. Program. Lang. Syst. 31, 6, Article 20 (Aug. 2009), 39 pages.  https://doi.org/10.1145/1552309.1552310
  * Zhou et al. (2004) Pin Zhou, Vivek Pandey, Jagadeesan Sundaresan, Anand Raghuraman, Yuanyuan Zhou, and Sanjeev Kumar. 2004\.  Dynamic tracking of page miss ratio curve for memory management. In ACM SIGOPS Operating Systems Review, Vol. 38. ACM, 177–188.


11institutetext: Chan Zuckerberg Biohub, San Francisco, USA 11email:
<EMAIL_ADDRESS>
# Large-Scale Multi-Hypotheses Cell Tracking
Using Ultrametric Contours Maps

Jordão Bragantini  0000-0001-7652-2735 Merlin Lange  0000-0003-0534-4374 Loïc
Royer  0000-0002-9991-9724

###### Abstract

In this work, we describe a method for large-scale 3D cell-tracking through a
segmentation selection approach. The proposed method is effective at tracking
cells across large microscopy datasets on two fronts: (i) It can solve
problems containing millions of segmentation instances in terabyte-scale 3D+t
datasets; (ii) It achieves competitive results with or without deep learning,
bypassing the requirement of 3D annotated data, that is scarce in the
fluorescence microscopy field. The proposed method computes cell tracks and
segments using a hierarchy of segmentation hypotheses and selects disjoint
segments by maximizing the overlap between adjacent frames. We show that this
method is the first to achieve state-of-the-art in both nuclei- and membrane-
based cell tracking by evaluating it on the 2D epithelial cell benchmark and
3D images from the cell tracking challenge. Furthermore, it has a faster
integer linear programming formulation, and the framework is flexible,
supporting segmentations from individual off-the-shelf cell segmentation
models or their combination as an ensemble. The code is available as
supplementary material.

##  1 Introduction

Tracking cells in terabyte-scale time-lapse volumetric image datasets of
developing organisms acquired with the latest generation of light-sheet
microscopes [20, 24, 55] is still an open problem. This conceptual and
practical challenge involves segmenting individual cells and reconstructing
their trajectories and lineages across multiple rounds of cell divisions. The
scale of the data hinders the application of computationally intensive
algorithms as these datasets easily contain thousands to millions of
segmentation instances. Moreover, the lack of 3D ground-truth annotations also
limits the application and performance of deep learning (DL) and other machine
learning methodologies.

Methods usually tackle the problem in two different steps: computing the
segmentation (or detection) and linking the previously detected instances
given an association score. These methods can employ off-the-shelf
segmentation models [52, 56, 42], and might use association scores from hand-
crafted features [29]. However, segmentation mistakes undermine the
association score, harming the final tracking result. Because state-of-the-art
segmentation methods are based on deep learning and require hard-to-obtain
annotated data, their accuracy is greatly impaired when annotations are not
available.

More recent methods based on flow-field estimation [17, 18, 57, 30, 19]
circumvent these limitations by avoiding segmentation altogether and
estimating the position and motion of cells between frames, relying only on a
marker in the center of each cell and a link between the same cell and their
divisions in adjacent frames. The association score is computed through the
motion estimation. This kind of cell-centroid supervision requires much less
effort than annotating volumetric segmentation instances.

Rather than avoiding segmentation, other approaches leverage the spatial
information and compute the cell segments and their tracks jointly by modeling
the cell nuclei as a mixture of Gaussians [3, 51], computing multiple
segmentation hypotheses using ellipsoid fitting on their contours [60], or
formulating the splitting and linking of regions (_i.e_. cells) in the image
as a graph-cut problem in space and time [22, 58, 15]. While these approaches
provide an elegant framework for joint segmentation and tracking, they are
much more computationally intensive when compared to flow-based methods.

In practice, cell segmentations are invaluable because it assists subsequent
downstream image analysis steps such as track validation and quantitative
analysis (_e.g_. morphological features, image intensity). Thus, we propose a
novel method with the following goals in mind:

  * •

The segmentation instances should be computed jointly with tracking to avoid
dependence on the accuracy of the segmentation method and to leverage both
spatial and temporal information.

  * •

The method should be computationally efficient and support out-of-memory
processing of multi-terabyte datasets.

  * •

It should leverage existing pre-trained models or traditional image processing
for segmentation while achieving reasonable performance when ground-truth
labels are unavailable.

To accomplish this, we propose a new algorithm for joint segmentation and
tracking. It builds upon the concept of hierarchical segmentation and a novel
method for extracting disjoint segments from a hierarchy and, at the same
time, formulating it as an integer linear program (ILP) that respects the
biological constraints (_e.g_., cell divides into two).

We perform multiple experiments evaluating different use cases of the proposed
method and validate its performance against state-of-the-art-methods using the
cell tracking challenge [61] and the epithelial cell benchmark [15]. Moreover,
we show our method can scale to large multi-terabyte scale real-world
microscopy datasets [66] with millions of cell instances.

##  2 Background and Related Works

###  2.1 Hierarchical Segmentation

Hierarchical segmentation has previously been used in cell tracking [60] but
in a very limited scenario. Here, we review a more general representation of
hierarchies and contours and their application to image segmentation.

Flat image segmentation (_i.e_. non-hierarchical) concerns partitioning an
image into disjoint segments, while hierarchical segmentation produces nested
regions. Let $V$ be our basic units, the image voxels, a hierarchy
$\mathcal{H}=\\{S_{1},S_{2},...,S_{n}\\}$ contains a sequence of
segmentations, where every $S_{i}$ is a collection of disjoint non-empty sets
and their union is $V$. The sequence is ordered, such that $S_{i-1}$ is a
refinement of $S_{i}$, so any element in $S_{i-1}$ is a subset of an
individual region of $S_{i}$. Therefore, at the finest level, $S_{1}$, each
region contains a single unit, and at the coarsest, there is a single region
$S_{n}=\\{V\\}$, the whole image. By assigning an increasing value
$\lambda_{i}$ for each partition $S_{i}$, the hierarchy can be represented as
a dendrogram.

The hierarchical segmentation cannot be easily represented with a set of
discrete labels as done for flat segmentation. Instead, they are commonly
represented using an Ultrametric Contour Maps (UCM) [38, 36, 4] where the
boundary between different regions receives the value $\lambda_{i}$ of where
their merging happened.

Several methodologies explore this duality between fuzzy contour maps (_i.e_.
UCM, edge detection, affinity map) and hierarchies. Notably, Pont-Tuset  _et
al_. [48] accelerated the computation of normalized cut [53] of affinity maps
by starting the clustering from a lower-resolution image and increasing its
resolution, and finally aligning and merging the results into a UCM. Maninis
_et al_. [31] used a convolutional neural network (CNN) to estimate boundaries
of different orientations and applied the oriented watershed transform [5] to
obtain the hierarchical segmentation.

In the more general context, these nested segmentations are a hierarchical
clustering of an edge-weighted graph $G=(V,E)$, where $V$ are the image voxels
and $E$ the edges connecting nearby elements in the image planar lattice. For
example, the watershed framework from [38] was extended to edge-weighted
_image graphs_ [11, 36].

In the microscopy domain, segmentation from edge weights is extremely popular
[10, 16, 64]. Because most applications expect a single label per instance,
methodologies based on correlation clustering [23, 25, 67, 43] are preferred
over hierarchical [41] because they avoid the extremely challenging step of
selecting the segment from the hierarchy without using a naive horizontal cut
in the dendrogram. Kiran and Serra [26] proposed an algorithm that finds the
optimum flat partition in a hierarchy for a restricted set of energies. More
recently, GASP [6] was proposed to unite correlation and hierarchical
clustering.

Our approach is not constrained to any hierarchical segmentation method; we
chose the watershed framework due to its simplicity and computational
performance.

###  2.2 Cell Tracking

Cell tracking has been a central problem for the bioimage community for more
than two decades [69]. Several different strategies have been explored, a non-
exhaustive list being: contour evolution [13]; sequential fitting of mixtures
of Gaussians [3]; dynamic programming using the Viterbi algorithm [29]; and
probabilistic graphical models that account for under/over-segmentation errors
[51].

Most related to our approach are joint segmentation and tracking methods,
where the final cell segments are found during the tracking step [9]. In [60],
the authors compute a hierarchical segmentation by fitting ellipsoids to their
cells detection binary map contours, then, they estimate their association
score using gradient boosting [14]; and obtain the tracks by solving a network
flow ILP.

Jug  _et al_. [22] proposed (and later improved [49]) the moral lineage
tracing (MLT) problem by formulating tracking as a minimum cost multi-cut
subject to biological constraints: cells can divide but not merge, and that
they can appear or disappear from the field of view. This approach extended
the wide literature of multi-cut problems in the spatial domain (_i.e_., image
segmentation) [23, 25, 67, 43] to include the time axis and additional
constraints. Despite the elegant framework, the NP-hardness of the multi-cut
problem limits its application to medium-size datasets, even with specialized
algorithms [49].

More recently, flow-based methods [17, 40, 57] were proposed and have shown
the most competitive performance. In [18], they estimate both the centroid
detection and the motions of cells using a CNN. The tracking is computed by
advecting backward in time using motion estimation and starting from the
cells’ detections. Malin-Mayor  _et al_. [30, 19] extends this method by
proposing a loss function that allows sparse annotation of the cells’
centroids and their links, additionally, the lineages are obtained with an ILP
solver rather than a heuristic as in earlier work [17, 18, 57].

Other notable DL-learning-based tracking methods employ different strategies:
[44] tracks cells by clustering pixels from adjacent frames using their
coordinates and embeddings from a recurrent hourglass network. Starting from a
set of segments, Ben  _et al_. [7] construct a directed graph and solves the
tracking using a graph neural network to classify the edges that belong to the
cell lineages.

##  3 Method

Our method was developed for fluorescence microscopy images, where cells are
engineered to express localized fluorescence, making them or their boundaries
easily identifiable from the background. The challenge is to split the cells’
instances, which cannot be easily resolved because of optical artifacts, and
link these instances through time.

To avoid committing to a single potentially erroneous candidate segment per
cell, we construct a hierarchy of segments. Ideally, the true cell instance
will be one of these candidate segments. This differs from methods that
aggregate superpixels during tracking [22, 49] because the set of possible
segments is fixed during the ILP optimization.

While it is rather optimistic to assume that the hierarchy will always closely
fit the true segments, this approach can always find a single region at a
finer level that belongs to a unique true segment at the cost of not covering
it completely.

We formulate the problem of finding the instances in the hierarchy as an
optimization problem to maximize the overlap between the unknown true
instances and a selection of disjoint segments from the hierarchy, which we
will call maximum total intersection over union (MTIoU).

Given a hierarchy $\mathcal{H}$ and a set of ground-truth segments
$\mathcal{G}=\\{g_{1},g_{2},...,g_{n}\\}$, the set of disjoint segmentation
with MTIoU can be found with the following ILP:

$\displaystyle\operatorname*{arg\,max}_{\mathbf{x}}$
$\displaystyle=\sum_{p\in\mathcal{H}}\sum_{q\in\mathcal{G}}w_{pq}x_{pq}$ (1)
$\displaystyle\text{s.t.}\quad y_{p}$
$\displaystyle=\sum_{q\in\mathcal{G}}x_{pq}\quad\forall\ p\in\mathcal{H}$ (2)
$\displaystyle\sum_{p\in\mathcal{H}}x_{pq}$ $\displaystyle\leq 1\quad\forall\
q\in\mathcal{G}$ (3) $\displaystyle y_{p}+y_{q}$ $\displaystyle\leq
1\quad\forall\ p,q\in\mathcal{H}|p\subset q$ (4) $\displaystyle y_{p}$
$\displaystyle\in\\{0,1\\}\quad\forall\ p\in\mathcal{H}$ (5) $\displaystyle
x_{pq}$ $\displaystyle\in\\{0,1\\}\quad\forall\
p,q\in\mathcal{H}\times\mathcal{G}$ (6)

where $w_{pq}$ is the intersection over union (IoU) between segment $p$ and
$q$; $x_{pq}$ indicates if $p$ is assigned to $q$; $y_{p}$ indicates if
segment $p$ was selected ($y_{p}=1$) with Eq. 2 forcing a single assignment
per segment; Eq. 3 forces a single assignment per ground-truth segment; Eq. 4
restricts the solution to only disjoint segments (_i.e_. they belong to
different branches of the hierarchy).

In practice, information about the true instances is not available, so we
approximate this problem by solving the segmentation selection between
hierarchies of adjacent frames under the assumption that they are similar.
Thus, the cell segments will overlap in some nodes in the hierarchies.

By formulating the ILP considering every time point, not just a pair of
frames, and including the biological constraints, such as cell division,
disappearing and appearing cells [60], the resulting segmentation and their
association between hierarchies can reconstruct the cell lineages.

Therefore, the ILP formulation of the MTIoU between hierarchies for cell
tracking becomes:

$\displaystyle\operatorname*{arg\,max}_{\mathbf{x}}$
$\displaystyle=\sum_{t\in[2..T-1]}\sum_{p\in\mathcal{H}_{t-1}}\sum_{q\in\mathcal{H}_{t}}w_{pq}x_{pq}+$
$\displaystyle\sum_{p\in\bm{\mathcal{H}}}w_{\alpha}x_{\alpha
p}+w_{\beta}x_{p\beta}+w_{\delta}x_{\delta p}$ (7)
$\displaystyle\text{s.t.}\quad y_{q}$ $\displaystyle=x_{\alpha q}+\sum_{p\in
H_{t-1}}x_{pq}\ \ \forall\ q\in\mathcal{H}_{t}\ \forall\ t\in[1..T]$ (8)
$\displaystyle y_{p}+x_{\delta p}$ $\displaystyle=x_{p\beta}+\sum_{q\in
H_{t+1}}x_{pq}\ \ \forall\ p\in\mathcal{H}_{t}\ \forall\ t\in[1..T]$ (9)
$\displaystyle y_{p}$ $\displaystyle\geq x_{\delta p}\ \ \forall\ \
p\in\bm{\mathcal{H}}$ (10) $\displaystyle y_{p}+y_{q}$ $\displaystyle\leq
1\quad\forall\ p,q\in\mathcal{H}_{t}\ |\ p\subset q\quad\forall\ t\in[1..T]$
(11) $\displaystyle y_{p}$ $\displaystyle\in\\{0,1\\}\ \ \forall\
p\in\bm{\mathcal{H}}$ (12) $\displaystyle x_{pq}$ $\displaystyle\in\\{0,1\\}\
\ \forall\ p,q\in\mathcal{H}_{t-1}\times\mathcal{H}_{t}\ \forall\ t\in[2..T]$
(13)

where $\mathcal{H}_{t}$ is the candidate segments (_i.e_. hierarchy) of frame
$t$ and $\bm{\mathcal{H}}=\cup_{t\in[1..T]}\mathcal{H}_{t}$. We consider
$\mathcal{H}_{t}=\emptyset$ when $t$ is out of range $[1..T]$; $x_{\alpha
p},x_{p\beta},x_{\delta p}$ are the segment $p$ appearance, disappearance, and
division binary variables, and their non-positive penalization parameters
$w_{\alpha},w_{\beta},w_{\delta}$, respectively. For segments from $t=1$, we
set $w_{\alpha}=0$, and $w_{\beta}=0$ for $t=T$ as boundary conditions;
$y_{p}$ is a slack binary variable indicating if segment $p$ is part of the
solution. Equation 8 and 9 are the constant flow constraints [60] such that
for the selected segments, there is an incoming segment from the previous
frame or it is an appearing segment, and the outflow allows up to 2 out-going
segments (_i.e_. division) or its disappearance; Eq. 10 asserts there are only
divisions from an existing cell; Eq. 11 restricts the solution to disjoint
segments as in Eq. 4. In Section 4.3, we compare the run time of our ILP
formulation with previous work.

###  3.1 Hierarchy construction

While not limited to any hierarchical segmentation method, we chose to compute
the hierarchies using the hierarchical watershed framework [46, 45], which
allows constructing hierarchies in log-linear time (sorting complexity) with
respect to the number of edges in a graph with the minimum-spanning tree
algorithm [37] and Tarjan union-find [59].

For each frame $t$, we estimate a set of candidate segments,
$\mathcal{H}_{t}$, from a binary foreground, $F_{t}$, and a fuzzy contour map,
$C_{t}$, provided by the user. The foreground $F_{t}$ indicates the location
of cells against the background. The contour map indicates where the cell
boundaries (_i.e_. contours) may lie. Since the connected components in
$F_{t}$ are disconnected by definition, a hierarchy is computed for each one
and integrated into a single $\mathcal{H}_{t}$ while ignoring the background.

Each hierarchy is constructed using the watershed hierarchy by area [34, 45]
on the undirected weighted graph from their respective connected component in
$F_{t}$. The edge weight between a pair of pixels is their average intensity.

Using the full hierarchies on the tracking ILP (Eq. 7) can easily make the
problem too large to solve. Therefore, three criteria were used to filter
irrelevant segments [47] and decrease the number of candidate solutions:

  1. 1.

Minimum size: the hierarchy is pruned below a minimum size threshold;

  2. 2.

Maximum size: the hierarchy is cut above a maximum size threshold;

  3. 3.

Non-relevant contours: segments with average contour strength below a
threshold are fused.

The minimum and maximum cell sizes are usually known for each organism, and
the non-relevant contours filtering can be tuned by observing the under (over)
segmentation of the solution. Figure 1(a) shows a hierarchy filtering example.

(a) Example of hierarchy filtering by minimum size, maximum size, and contour
strength; dashed lines show deleted elements. Leaves of the hierarchy are not
shown

(b) Example of linkage between two hierarchies and their weights matrix
$w_{..}$ which includes all $w_{p,q}$. Solution considering $w_{..}$ in blue
and in red solution of $w^{p}_{..}$ when $p>=2$

The MTIoU can suffer from over-segmentation when spurious segments with high
overlap between frames are present in the hierarchies. For example, in Figure
1(b), despite $AB$ and $CD$ having perfect overlap ($w_{AB,CD}=1.0$) when
solving Eq. 7, the solution selects segments $A,B,C,D$ (blue solution) because
their sum of weights is larger. By taking a power $p$ of $w_{..}$, we reduce
the over-segmentation as it favors fewer segments with larger $w_{..}$ values
(red solution) rather than many with smaller $w_{..}$ values.

###  3.2 Scaling to terabyte-scale datasets

As described previously, one major challenge of joint segmentation and
tracking is the intense computational requirements. This section will discuss
how our approach can effectively scale to terabyte-scale datasets.

The hierarchical segmentation is computed per connected component of each
foreground map $F_{t}$, greatly reducing the number of pixels considered. For
each hierarchy, the segments outside of the minimum and maximum size range are
removed. Next, the disjoint segmentation constraint, Eq. 11, is computed by
going through each segment to its last ancestor. For each segment, we store
their center coordinates and bounding boxes for computing their IoU. Since the
segmentation of each volume (plane) is independent, it can be easily run in
parallel.

The association weight between segments is only computed for a subset of
$k\times N_{t}$ pairs, where $N_{t}$ is the total number of segments on frame
$t$. For $(t-1,t)$-pair of frames, we start by computing a KDTree [62] with
the segments’ coordinates from $t-1$ and querying the $2k$-nearest neighbors
within a certain radius for each segment coordinate in $t$. From these $2k$
samples, we select the $k$ with the highest IoU, thus, we avoid computing the
$w_{..}$ for every pair in $(t-1,t)$, as $w_{..}$ is significantly sparse.
Additionally, since the IoU only considers the non-zero values of the
segmentation masks, it can be computed using a crop around the segments’
bounding box. When the bounding boxes between a segment and its neighbor do
not intersect, its IoU is zero and does not require further computation. The
weights between different $(t-1,t)$ are independent, so each pair of frames
can be run in parallel.

The ILP can easily reach millions of variables due to multiple candidate
segments per cell. When that happens, it cannot be solved in a reasonable
amount of memory or time. However, its global solution can be approximated by
solving the problem in multiple overlapping time windows as in [30]. The
regions between overlapping windows are not used in the final solution except
for a single frame, where the segments from already computed solutions are
used as constraints, enforcing continuity of tracks, essentially stitching
adjacent windows.

Malin-Mayor  _et al_. [30] employed their own distributed computing scheduler
to process these windows in parallel while avoiding concurrent processing of
overlapping windows. We propose a simpler two-pass interleaved scheduling
scheme: i) Split the timelapse only along the time-axis into $N$ overlapping
windows, where the overlap is smaller than a single window; ii) Solve the ILPs
of the even-numbered window; iii) Add the solutions from the even-numbered
ILPs as constraints; and iv) solve the odd-numbered ILPs. This simple strategy
guarantees that no overlapping window is processed simultaneously, and it can
be trivially executed in parallel at the cost of waiting for every even-
numbered window to finish before starting the odd-numbered executions. Further
improvements can be made by waiting for only adjacent pairs to finish (_e.g_.
window 3 depends on 2 and 4).

##  4 Experiments

We evaluate our claims and the proposed method in four different experiments:

  1. 1.

Multiple quantitative comparisons of the proposed algorithm against state-of-
the-art methods for cell tracking.

  2. 2.

Analysis of an ensemble of segmentations with traditional segmentation
algorithms and off-the-shelf segmentation CNNs.

  3. 3.

Runtime comparison of different ILP formulations.

  4. 4.

A multi-terabyte dataset cell tracking use case.

When available, the experiments were evaluated using the metrics proposed by
the Cell Tracking Challenge (CTC) [61]: TRA, a score of operations required to
transform the resulting lineage graph to the ground-truth lineage [32]; SEG,
the average IoU between the predicted and the ground-truth segmentation
instances (a predicted instance is assigned to a ground-truth instance when it
overlaps with more than half of its pixels); and CTB, the average between the
TRA and SEG scores, this is the primary ranking score in the competition and
used to characterize state-of-the-art method.

###  4.1 Quantitative comparison

Cell Tracking Challenge: This first quantitative experiment was run on five
datasets from CTC: Fluo-C3DL-MDA231: Infected human breast carcinoma cells;
Fluo-N3DH-CE: C. Elegans developing embryo [35]; Fluo-N3DL-DRO: D.
Melanogaster developing embryo [3]; Fluo-N3DL-TRIF (TRIC): Developing
Tribolium Castaneum embryo (and its cartographical projection). In each case,
two timelapses are available for training and validation purposes, and two
other timelapses are used for testing. Test set results were evaluated by the
competition organizers, and their labels are not publicly available.

The MDA231 and CE datasets have ground-truth tracks and silver-standard
segmentation annotation automatically generated by a consensus [2] of the top-
ranking algorithms from previous iterations of CTC. The DRO, TRIF, and TRIC
datasets only have ground-truth tracks for a specific group of cells and a few
2D labels. In these cases where not every lineage are labeled, the initial
positions of these cells’ tracks are provided during submission and must be
used to select only the cells that belong to these lineages. These sparse
annotations make the tracking much more challenging because less training data
are available, and during evaluation, other lineages must be removed, so
missing links on the predicted lineages will prune the rest of the tracks and
greatly decrease the final score.

Table 1:  Results on the hidden test dataset of the cell tracking challenge at
the time of article submission; subscripts indicate the method used according
to enumeration; our submission is in bold (anonymous submission on the
competition website until paper acceptance)

Dataset | Rank |  CTB  |  TRA  |  SEG
---|---|---|---|---
MDA | 1st |  0.797 2 |  0.884 2 |  0.710 2
2nd |  0.761 5 |  0.882 1 |  0.642 5
3rd |  0.757 1 |  0.880 5 |  0.632 1
7th |  0.683 |  0.818 |  0.549
CE | 1st |  0.844 |  0.987 3 |  0.759 6
2nd |  0.829 1 |  0.979 4 |  0.729 2
3rd |  0.808 2 |  0.975 7 |  0.722 1
4th | - |  0.967 |  0.722
DRO | 1st |  0.708 |  0.802 |  0.613
2nd |  0.617 1 |  0.785 4 |  0.613 1
3rd |  0.591 4 |  0.785 3 |  0.397 4
TRIF | 1st |  0.841 |  0.955 3 |  0.746
2nd |  0.804 3 |  0.936 |  0.684 8
3rd |  0.785 8 |  0.886 8 |  0.654 3
TRIC | 1st |  0.867 1 |  0.952 3 |  0.791 1
2nd |  0.816 3 |  0.942 1 |  0.766 2
3rd |  0.787 2 |  0.854 |  0.680 3
4th |  0.754 | - |  0.654

Table 2: Epithetial Cell Benchmark segmentation (SEG metric) results, top
scores in bold, higher is better. Our tracking algorithm is indicated by ✝

Method |  | Variation |  | PER |  | PRO
---|---|---|---|---|---|---
FFN |  | Default |  | $0.879$ |  | $0.796$
MTL |  | Default |  | $0.904$ |  | $\mathbf{0.818}$
TA |  | Default |  | - |  | $0.758$
PlantSeg |  | Default model |  | $0.787$ |  | $0.761$
| Fine tuned |  | $0.885$ |  | $0.800$
MALA |  | Original |  | $\mathbf{0.907}$ |  | $\mathbf{0.817}$
| Our repeat |  | $0.902$ |  | $0.810$
|  Our tracking${}^{\text{\char 61}}$ |  | $\mathbf{0.909}$ |  | $\mathbf{0.817}$
Blur |  |  Our tracking${}^{\text{\char 61}}$ |  | $0.853$ |  | $0.804$

Table 3: Tracking and segmentation results for individual segmentation methods
and their ensemble

File | Method |  | TRA | SEG | CTB
---|---|---|---|---|---
1 | Watershed |  | 0.905 | 0.593 | 0.749
Cellpose |  | 0.863 | 0.671 | 0.767
Stardist |  | 0.910 | 0.724 | 0.817
Ensemble |  | 0.956 | 0.740 | 0.848
2 | Watershed |  | 0.898 | 0.695 | 0.796
Cellpose |  | 0.908 | 0.743 | 0.825
Stardist |  | 0.941 | 0.818 | 0.880
Ensemble |  | 0.956 | 0.826 | 0.891

We used a U-Net CNN [50] for predicting the cell and edge detection maps, and
both maps were optimized using their own dice loss [12] with additional dice
losses on side edge-detection predictions [65] during training.

The DRO, TRIF, and TRIC do not have any 3D segmentation annotations.
Therefore, we tried two approaches. i) Using our U-Net and training it using
labels computed from a traditional segmentation algorithm. We denoted this as
pseudo-labels (PL). ii) Using filters to detect the cells and compute the
contour topology, denoted by image processing (IM).

In both cases, the foreground (_i.e_. cells) was detected using a difference
of Gaussians and thresholding it with Otsu. The contour maps were simply the
inverse of the Gaussian blurred image. For the pseudo-labels, we used the
contour maps and foregrounds to run the hierarchical watershed and performed a
horizontal cut in the hierarchy to obtain the flat segments used for training
our CNN. In IM, the foreground and contour maps were used as the tracking
input.

We optimized the hyper-parameters using grid search and cross-validation on
the training set to avoid a biased estimate of our final accuracy, so models
trained on dataset 01 were applied to dataset 02 and the other way around.
More details about these procedures, including their parameters, can be found
in the supplementary materials.

The non-deep-learning solutions achieved better results (see supplementary) on
these datasets than deep learning combined with pseudo-labels, therefore, IM
was used in our final submission. Next, the proposed method was evaluated on
the hidden test set.

The competition organizers report the top three methods for each metric and
provide access to our ranking regardless of our position. The final results
are reported in Table 4.1.

The leaderboard’s methods at the time of submission were:

  1. 1.

KTH-SE: [29] image processing solution with tracking using Viterbi algorithm;

  2. 2.

KIT-GE: EmbedTrack [28], a joint flow-based tracking and pixel embedding
segmentation method;

  3. 3.

MPI-GE-CBG: Description is not publicly available, source code indicates it is
a DL-based method with nearest neighbor linking;

  4. 4.

JAN-US: [19] it is an improved version of [30], where they classify the cell
states to detect division and assist the ILP formulation and use a structured
SVM to find the optimal weights of the ILP automatically;

  5. 5.

LEID-NL: Model-evolution [13] with routines to split dividing cells and track
cells entering the field of view;

  6. 6.

MU-CZ: H-minima watershed [33] plus U-Net for cell contour and foreground
prediction, tracking uses an intersection-based heuristic;

  7. 7.

IGFL-FR: ELEPHANT [57] uses the cell detection and advection of DL flow
prediction linking as item 3;

  8. 8.

RWTH-GE: [54] Laplacian of Gaussian cell detection, with nearest neighbor
linking with manual curation for missing detections.

The large variety of top-ranking methods shows that there is not a single
best-tracking algorithm. Despite this, our approach obtained the best
performance in the combined score, CTB, on three out of the five datasets: CE
with supervised DL and TRIF and DRO with image processing only. In the TRIC
dataset, most of our mistakes were on the periphery of the projection, where
the cells’ movements are inflated, making our IoU weights significantly less
effective. In the MDA231, we achieved our worst ranking. We suspect our edge-
detection network did not generalize to the test dataset.

Epithelial Cell Benchmark [15] evaluates a different tracking modality, where
the cells are densely packed, and the cells’ membranes are visible rather than
their nucleus. It contains eight 2D timelapses, divided into two cell types,
three peripodial (PRE) and five proper discs (PRO) cells’ timelapses [1]. In
this context, graph-cut methodologies extended from the electron microscopy
segmentation domain are commonly used [16, 64].

We compared with the original baselines MTL [49], FFN [21], TA [1], MALA [16]
and PlantSeg [64] which uses Multicut from ilastik [8] and GASP [6].
Additionally, we reproduced MALA’s experiment
111https://github.com/funkey/flywing, achieving a score lower than the
original run but still competitive, Table 4.1. Refer to [15] for a more in-
depth description of baselines.

We evaluated our algorithm in two scenarios, marked with a ✝ in Table 4.1: (i)
the original frames filtered with a Gaussian kernel ($\sigma=1.0$) as our
contour map, (ii) using MALA’s edge predictions as input.

The minimal solution without deep learning (i) is comparable to some of the DL
approaches, especially in the PRO datasets, which have diverse cell sizes. Our
solution from MALA’s boundaries (ii) meets or surpasses other methods despite
having inferior predictions than the original MALA work. The performance could
be further improved by accessing MALA’s original network weights. This
displays the effectiveness of our approach at performing the correct
segmentation against state-of-the-art DL and graph segmentation schemes.
Videos are attached as supplementary.

###  4.2 Segmentation ensemble experiment

Most biologists cannot train their segmentation models or design complex
image-processing pipelines. Therefore, leveraging broadly available
segmentation tools is essential. We show that our method can use segmentations
from different algorithms and combine them to achieve superior performance
using fuzzy contour map representations on the training set of Fluo-N2DL-HeLa
[39], a 2D dataset of HeLa cells expressing H2b-GFP from CTC. We chose this
dataset since it contains high-quality images and annotations. Because it is a
2D dataset, we can leverage state-of-the-art off-the-shelf 2D segmentation
methods such as Cellpose [56] and Stardist [52].

(a) GT

(b) ensemble

(c) watershed

(d) cellpose

(e) stardist

Figure 2: Results from individual methods versus their ensemble, mistakes are
indicated by the red arrows. (a) ground-truth cell detection; (b) result;
(c-e) segmentation results of different methods

We evaluate three different segmentation methods: a classical watershed from
h-minima using a distance transform as the watershed basins [27, 63]; Cellpose
[56] using their nuclei pre-trained model; and Stardist [52] using their ”2d
versatile fluorescence” pre-trained model.

For each segmentation map, we extract a binary map indicating foreground
(_i.e_. cell) versus background and the distinct labels’ binary contours. The
ensemble maps were created by combining the foregrounds using the logical OR
operator and the contours by averaging. Table 4.1 presents the scores from
running our proposed algorithm on the foreground and contour maps from the
labels of each method. The top score for the ensemble model shows that the
contour map representation is an effective way to combine multiple
segmentation, and our ILP picks segments that are better than any individual
segmentation result, consequently improving cell tracking. Figure 2 shows how
this approach benefits segmentation and tracking.

###  4.3 ILP run time comparison

Figure 3: Runtime comparison between ILP formulations

Cells moving in and out of the field of view and cell division constraints are
critical aspects of the cell-tracking ILP and can be formulated differently.
In this experiment, we compare the performance of our formula against other
recent ILP-based methods for optimizing Eq. 7, Türetken  _et al_. [60] and
Malin-Mayor  _et al_. [30], with the same experimental setup and parameters.
We included Eq. 11 to Malin-Mayor formulation to allow multiple segmentation
hypotheses and set $w_{\beta}=0$ for all methods because Malin-Mayor does not
penalize disappearance events.

We compared the average run time of 3 runs with increasingly longer timelapses
(_i.e_. window) of the Fluo-N3DH-CE datasets. The window started from the
middle of the timelapse and increased at both ends. A linear fit for each
method reported slope coefficients of $1.75\pm 0.12$ for Türetken,
$\mathbf{1.37\pm 0.06}$ for Ours, and $1.84\pm 0.06$ for Malin-Mayor, the
smaller, the better. Figure 3 display the measurements and their linear
regression. Therefore, our formulation is faster and scales better with
increasingly larger problems. All sizes except the initial size of 10 produce
reasonable tracking results.

###  4.4 Terabyte-scale tracking use case

We ran our algorithm on a 3D time-lapse dataset of a developing zebrafish
embryo [66] recorded with a light-sheet microscope. The recording starts at 12
hrs post-fertilization and lasts 6.6 hrs with an acquisition interval of 30
secs. After cropping and stabilizing the embryo region, the dataset’s
uncompressed size is about 3.4 TB, and its shape is $791\times 448\times
2174\times 2423$ voxels (t, z, y, x).

The algorithm was executed in a compute cluster using SLURM [68], following
the strategies described in Section 3.2. For the tracking step, a window of
size 50 with an overlap of 5 frames on each side was used. For a fair
comparison, we limited our resources to 100 CPU cores and 20 GPUs as used in
[30].

Processing took approximately a total of 9.5 hours of distributed computing,
including: 1 hour and 20 min. to detect the foreground and predict the
contours using a U-Net as described in the other experiments, and trained on
pseudo-labels from a watershed segmentation; 2 hours and 7 min. to compute the
hierarchical segmentation for all frames; 40 minutes for the association
between segments; 3 hours for the tracking first pass; 2 hours and 8 minutes
for the second pass; and 18 minutes to export the data. The first step used
GPUs, the remaining only CPU cores, processing in total 25 hours of GPU time
and 31 days 22 hours of CPU time. Resulting in more than 21.5 million cell
instances. The ILP size is partially proportional to the number of resulting
instances.

Results are shown in Figure 4 and in a video in supplementary materials.
Annotations are not available for a quantitative analysis. However, we expect
an accuracy similar to embryonic datasets from the cell tracking challenge,
Tab. 4.1.

In comparison, with similar resources, Malin-Mayor  _et al_. [30] tracked a
mouse dataset in an earlier developmental stage (therefore, fewer and less
densely distributed nuclei) in 44 hours, segmenting about 7 million cell
instances. Information about the total GPU and CPU time is not available.

(a)

(b)

(c)

(d)

Figure 4: Whole embryo segmentation results: (a) Maximum intensity projection;
(b) 3D rendering of segmentation labels; (c) YX slice; and (d) ZY slice.

##  5 Limitations & Social Impact

While the methodology is more robust to segmentation mistakes than non-joint
segmentation and tracking methods, the space of possible segmentations is
constrained from the input contour map. Therefore, contour maps that poorly
corresponds to the cell boundaries greatly impar our methods performance.

The negative social impact of this work would be the employment of this method
to assist unethical/immoral biological research, as with any other
bioinformatics tool.

##  6 Conclusion

In this paper, we introduced a novel approach to cell tracking in microscopy
data that selects cell segments from hierarchical segmentations. Our method
demonstrated state-of-the-art performance on three out of five datasets from
the Cell Tracking Challenge [61] and on the Epithelial Cell Benchmark [15]. To
date, it is the only method to excel in both nuclei and membrane-based cell
tracking while also scaling to terabyte-scale datasets.

However, the benefits of using contour maps for cell tracking extend beyond
surpassing benchmark metrics and improving segmentation performance. The
contour (_i.e_. affinity) space offers more flexibility than the traditional
domain of discrete labels, enabling the use of an ensemble of segmentation
models and differential models for edge detection.

Therefore, the outlined goal of developing an effective joint segmentation and
tracking method that supports a variety of inputs was achieved.

Overall, we believe our approach has significant potential for advancing cell
tracking in complex microscopy data and opening new avenues for exploration in
bioimaging research. In future work, a model could be used to predict the ILP
weights.

## Acknowledgements

We thank Ahmed Abbas, Ahmet Can Solak, Alexandre Falcão, Guillaume Le Treut,
Ilan Francisco, Laurent Najman, Paul Swoboda, and Sandy Schmid for discussions
regarding this work and Chan Zuckerberg Biohub for the support.

## References

  * [1] Aigouy, B., Umetsu, D., Eaton, S.: Segmentation and quantitative analysis of epithelial tissues. Drosophila: Methods and Protocols pp. 227–239 (2016)
  * [2] Akbaş, C.E., Ulman, V., Maška, M., Jug, F., Kozubek, M.: Automatic fusion of segmentation and tracking labels. In: Computer Vision–ECCV 2018 Workshops: Munich, Germany, September 8-14, 2018, Proceedings, Part VI 15. pp. 446–454. Springer (2019)
  * [3] Amat, F., Lemon, W., Mossing, D.P., McDole, K., Wan, Y., Branson, K., Myers, E.W., Keller, P.J.: Fast, accurate reconstruction of cell lineages from large-scale fluorescence microscopy data. Nature Methods 11(9), 951–958 (2014)
  * [4] Arbelaez, P.: Boundary extraction in natural images using ultrametric contour maps. In: Conference on Computer Vision and Pattern Recognition Workshop. pp. 182–182. IEEE (2006)
  * [5] Arbelaez, P., Maire, M., Fowlkes, C., Malik, J.: Contour detection and hierarchical image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence 33(5), 898–916 (2010)
  * [6] Bailoni, A., Pape, C., Hütsch, N., Wolf, S., Beier, T., Kreshuk, A., Hamprecht, F.A.: GASP, a generalized framework for agglomerative clustering of signed graphs and its application to instance segmentation. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 11645–11655 (2022)
  * [7] Ben-Haim, T., Raviv, T.R.: Graph neural network for cell tracking in microscopy videos. In: European Conference on Computer Vision. pp. 610–626. Springer (2022)
  * [8] Berg, S., Kutra, D., Kroeger, T., Straehle, C.N., Kausler, B.X., Haubold, C., Schiegg, M., Ales, J., Beier, T., Rudy, M., et al.: Ilastik: interactive machine learning for (bio) image analysis. Nature Methods 16(12), 1226–1232 (2019)
  * [9] Bise, R., Maeda, Y., Kim, M.h., Kino-Oka, M.: Cell tracking under high confluency conditions by candidate cell region detection-based-association approach. Biomedical Engineering pp. 1004–1010 (2013)
  * [10] Briggman, K., Denk, W., Seung, S., Helmstaedter, M., Turaga, S.C.: Maximin affinity learning of image segmentation. vol. 22 (2009)
  * [11] Cousty, J., Bertrand, G., Najman, L., Couprie, M.: Watershed cuts: Minimum spanning forests and the drop of water principle. IEEE Transactions on Pattern Analysis and Machine Intelligence 31(8), 1362–1374 (2009)
  * [12] Deng, R., Shen, C., Liu, S., Wang, H., Liu, X.: Learning to predict crisp boundaries. In: European Conference on Computer Vision. pp. 562–578 (2018)
  * [13] Dzyubachyk, O., Van Cappellen, W.A., Essers, J., Niessen, W.J., Meijering, E.: Advanced level-set-based cell tracking in time-lapse fluorescence microscopy. IEEE Transactions on Medical Imaging 29(3), 852–867 (2010)
  * [14] Friedman, J.H.: Stochastic gradient boosting. Computational statistics & data analysis 38(4), 367–378 (2002)
  * [15] Funke, J., Mais, L., Champion, A., Dye, N., Kainmueller, D.: A benchmark for epithelial cell tracking. In: Proceedings of The European Conference on Computer Vision (ECCV) Workshops. pp. 0–0 (2018)
  * [16] Funke, J., Tschopp, F., Grisaitis, W., Sheridan, A., Singh, C., Saalfeld, S., Turaga, S.C.: Large scale image segmentation with structured loss based deep learning for connectome reconstruction. IEEE Transactions on Pattern Analysis and Machine Intelligence 41(7), 1669–1680 (2018)
  * [17] Hayashida, J., Bise, R.: Cell tracking with deep learning for cell detection and motion estimation in low-frame-rate. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 397–405. Springer (2019)
  * [18] Hayashida, J., Nishimura, K., Bise, R.: Mpm: Joint representation of motion and position map for cell tracking. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 3823–3832 (2020)
  * [19] Hirsch, P., Malin-Mayor, C., Santella, A., Preibisch, S., Kainmueller, D., Funke, J.: Tracking by weakly-supervised learning and graph optimization for whole-embryo c. elegans lineages. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 25–35. Springer (2022)
  * [20] Huisken, J., Swoger, J., Del Bene, F., Wittbrodt, J., Stelzer, E.H.: Optical sectioning deep inside live embryos by selective plane illumination microscopy. Science 305(5686), 1007–1009 (2004)
  * [21] Januszewski, M., Kornfeld, J., Li, P.H., Pope, A., Blakely, T., Lindsey, L., Maitin-Shepard, J., Tyka, M., Denk, W., Jain, V.: High-precision automated reconstruction of neurons with flood-filling networks. Nature Methods 15(8), 605–610 (2018)
  * [22] Jug, F., Levinkov, E., Blasse, C., Myers, E.W., Andres, B.: Moral lineage tracing. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 5926–5935 (2016)
  * [23] Kappes, J.H., Speth, M., Andres, B., Reinelt, G., Schn, C.: Globally optimal image partitioning by multicuts. In: Energy Minimization Methods in Computer Vision and Pattern Recognition. pp. 31–44. Springer (2011)
  * [24] Keller, P.J., Schmidt, A.D., Wittbrodt, J., Stelzer, E.H.: Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy. Science 322(5904), 1065–1069 (2008)
  * [25] Keuper, M., Levinkov, E., Bonneel, N., Lavoué, G., Brox, T., Andres, B.: Efficient decomposition of image and mesh graphs by lifted multicuts. In: IEEE International Conference on Computer Vision. pp. 1751–1759 (2015)
  * [26] Kiran, B.R., Serra, J.: Global–local optimizations by hierarchical cuts and climbing energies. Pattern Recognition 47(1), 12–24 (2014)
  * [27] Lotufo, R.d.A., Falcão, A.X., Zampirolli, F.A.: IFT-watershed from gray-scale marker. In: Conf. on Graphics, Patterns and Images (SIBGRAPI). pp. 146–152. IEEE (2002)
  * [28] Löffler, K., Mikut, R.: Embedtrack—simultaneous cell segmentation and tracking through learning offsets and clustering bandwidths. IEEE Access 10, 77147–77157 (2022). https://doi.org/10.1109/ACCESS.2022.3192880
  * [29] Magnusson, K.E.: Segmentation and tracking of cells and particles in time-lapse microscopy. Ph.D. thesis, KTH Royal Institute of Technology (2016)
  * [30] Malin-Mayor, C., Hirsch, P., Guignard, L., McDole, K., Wan, Y., Lemon, W.C., Kainmueller, D., Keller, P.J., Preibisch, S., Funke, J.: Automated reconstruction of whole-embryo cell lineages by learning from sparse annotations. Nature Biotechnology pp. 1–6 (2022)
  * [31] Maninis, K.K., Pont-Tuset, J., Arbeláez, P., Van Gool, L.: Convolutional oriented boundaries: From image segmentation to high-level tasks. IEEE Transactions on Pattern Analysis and Machine Intelligence 40(4), 819–833 (2017)
  * [32] Matula, P., Maška, M., Sorokin, D.V., Matula, P., Ortiz-de Solórzano, C., Kozubek, M.: Cell tracking accuracy measurement based on comparison of acyclic oriented graphs. PloS one 10(12), e0144959 (2015)
  * [33] Meyer, F.: Morphological segmentation revisited. Space, Structure and Randomness: Contributions in Honor of Georges Matheron in the Field of Geostatistics, Random Sets and Mathematical Morphology pp. 315–347 (2005)
  * [34] Meyer, F., Oliveras Vergés, A., Salembier Clairon, P.J., Vachier, C.: Morphological tools for segmentation: connected operators and watersheds. Annales des télecommunications. Annals of telecommunications 52(7-8), 366–379 (1997)
  * [35] Murray, J.I., Bao, Z., Boyle, T.J., Boeck, M.E., Mericle, B.L., Nicholas, T.J., Zhao, Z., Sandel, M.J., Waterston, R.H.: Automated analysis of embryonic gene expression with cellular resolution in c. elegans. Nature Methods 5(8), 703–709 (2008)
  * [36] Najman, L.: On the equivalence between hierarchical segmentations and ultrametric watersheds. Journal of Mathematical Imaging and Vision 40(3), 231–247 (2011)
  * [37] Najman, L., Cousty, J., Perret, B.: Playing with kruskal: algorithms for morphological trees in edge-weighted graphs. In: Mathematical Morphology and Its Applications to Signal and Image Processing. pp. 135–146. Springer (2013)
  * [38] Najman, L., Schmitt, M.: Geodesic saliency of watershed contours and hierarchical segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence 18(12), 1163–1173 (1996)
  * [39] Neumann, B., Walter, T., Hériché, J.K., Bulkescher, J., Erfle, H., Conrad, C., Rogers, P., Poser, I., Held, M., Liebel, U., et al.: Phenotypic profiling of the human genome by time-lapse microscopy reveals cell division genes. Nature 464(7289), 721–727 (2010)
  * [40] Nishimura, K., Hayashida, J., Wang, C., Ker, D.F.E., Bise, R.: Weakly-supervised cell tracking via backward-and-forward propagation. In: European Conference on Computer Vision. pp. 104–121. Springer (2020)
  * [41] Nunez-Iglesias, J., Kennedy, R., Parag, T., Shi, J., Chklovskii, D.B.: Machine learning of hierarchical clustering to segment 2d and 3d images. PloS one 8(8), e71715 (2013)
  * [42] Ouyang, W., Beuttenmueller, F., Gómez-de Mariscal, E., Pape, C., Burke, T., Garcia-López-de Haro, C., Russell, C., Moya-Sans, L., de-la Torre-Gutiérrez, C., Schmidt, D., et al.: Bioimage model zoo: A community-driven resource for accessible deep learning in bioimage analysis. bioRxiv (2022)
  * [43] Pape, C., Beier, T., Li, P., Jain, V., Bock, D.D., Kreshuk, A.: Solving large multicut problems for connectomics via domain decomposition. In: Proceedings of the IEEE International conference on computer vision workshops. pp. 1–10 (2017)
  * [44] Payer, C., Štern, D., Neff, T., Bischof, H., Urschler, M.: Instance segmentation and tracking with cosine embeddings and recurrent hourglass networks. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 3–11. Springer (2018)
  * [45] Perret, B., Chierchia, G., Cousty, J., Guimaraes, S.J.F., Kenmochi, Y., Najman, L.: Higra: Hierarchical graph analysis. SoftwareX 10, 100335 (2019)
  * [46] Perret, B., Cousty, J., Guimaraes, S.J.F., Maia, D.S.: Evaluation of hierarchical watersheds. IEEE Transactions on Image Processing 27(4), 1676–1688 (2017)
  * [47] Perret, B., Cousty, J., Guimarães, S.J.F., Kenmochi, Y., Najman, L.: Removing non-significant regions in hierarchical clustering and segmentation. Pattern Recognition Letters 128, 433–439 (2019)
  * [48] Pont-Tuset, J., Arbelaez, P., Barron, J.T., Marques, F., Malik, J.: Multiscale combinatorial grouping for image segmentation and object proposal generation. IEEE Transactions on Pattern Analysis and Machine Intelligence 39(1), 128–140 (2016)
  * [49] Rempfler, M., Lange, J.H., Jug, F., Blasse, C., Myers, E.W., Menze, B.H., Andres, B.: Efficient algorithms for moral lineage tracing. In: IEEE International Conference on Computer Vision. pp. 4695–4704 (2017)
  * [50] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [51] Schiegg, M., Hanslovsky, P., Kausler, B.X., Hufnagel, L., Hamprecht, F.A.: Conservation tracking. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 2928–2935 (2013)
  * [52] Schmidt, U., Weigert, M., Broaddus, C., Myers, G.: Cell detection with star-convex polygons. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 265–273. Springer (2018)
  * [53] Shi, J., Malik, J.: Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence 22(8), 888–905 (2000)
  * [54] Stegmaier, J., Otte, J.C., Kobitski, A., Bartschat, A., Garcia, A., Nienhaus, G.U., Strähle, U., Mikut, R.: Fast segmentation of stained nuclei in terabyte-scale, time resolved 3d microscopy image stacks. PloS one 9(2), e90036 (2014)
  * [55] Stelzer, E.H., Strobl, F., Chang, B.J., Preusser, F., Preibisch, S., McDole, K., Fiolka, R.: Light sheet fluorescence microscopy. Nature Reviews Methods Primers 1(1), 73 (2021)
  * [56] Stringer, C., Wang, T., Michaelos, M., Pachitariu, M.: Cellpose: a generalist algorithm for cellular segmentation. Nature Methods 18(1), 100–106 (2021)
  * [57] Sugawara, K., Çevrim, Ç., Averof, M.: Tracking cell lineages in 3d by incremental deep learning. eLife 11, e69380 (2022)
  * [58] Tan, M., Le, Q.: Efficientnet: Rethinking model scaling for convolutional neural networks. In: International Conference on Machine Learning. pp. 6105–6114 (2019)
  * [59] Tarjan, R.E.: Efficiency of a good but not linear set union algorithm. Journal of the ACM (JACM) 22(2), 215–225 (1975)
  * [60] Türetken, E., Wang, X., Becker, C.J., Haubold, C., Fua, P.: Network flow integer programming to track elliptical cells in time-lapse sequences. IEEE Transactions on Medical Imaging 36(4), 942–951 (2016)
  * [61] Ulman, V., Maška, M., Magnusson, K.E., Ronneberger, O., Haubold, C., Harder, N., Matula, P., Matula, P., Svoboda, D., Radojevic, M., et al.: An objective comparison of cell-tracking algorithms. Nature Methods 14(12), 1141–1152 (2017)
  * [62] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al.: Scipy 1.0: fundamental algorithms for scientific computing in python. Nature Methods 17(3), 261–272 (2020)
  * [63] Van der Walt, S., Schönberger, J.L., Nunez-Iglesias, J., Boulogne, F., Warner, J.D., Yager, N., Gouillart, E., Yu, T.: scikit-image: image processing in python. PeerJ 2, e453 (2014)
  * [64] Wolny, A., Cerrone, L., Vijayan, A., Tofanelli, R., Barro, A.V., Louveaux, M., Wenzl, C., Strauss, S., Wilson-Sánchez, D., Lymbouridou, R., et al.: Accurate and versatile 3d segmentation of plant tissues at cellular resolution. eLife 9, e57613 (2020)
  * [65] Xie, S., Tu, Z.: Holistically-nested edge detection. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 1395–1403 (2015)
  * [66] Yang, B., Lange, M., Millett-Sikking, A., Zhao, X., Bragantini, J., Vijay Kumar, S., Kamb, M., Gómez-Sjöberg, R., Solak, A.C., Wang, W., et al.: Daxi—high-resolution, large imaging volume and multi-view single-objective light-sheet microscopy. Nature Methods 19(4), 461–469 (2022)
  * [67] Yarkony, J., Ihler, A., Fowlkes, C.C.: Fast planar correlation clustering for image segmentation. In: European Conference on Computer Vision. pp. 568–581. Springer (2012)
  * [68] Yoo, A.B., Jette, M.A., Grondona, M.: SLURM: Simple linux utility for resource management. In: Job Scheduling Strategies for Parallel Processing: 9th International Workshop. pp. 44–60. Springer (2003)
  * [69] Zimmer, C., Labruyere, E., Meas-Yedid, V., Guillén, N., Olivo-Marin, J.C.: Segmentation and tracking of migrating cells in videomicroscopy with parametric active contours: A tool for cell-based drug testing. IEEE Transactions on Medical Imaging 21(10), 1212–1221 (2002)