As already mentioned previously, Karrenbauer and Wöll [78] follow the approach to minimize the squared distances between linked particles by choosing the optimal set of connections from all possible connec-tions. This is done at once for the whole data set of all positions at all times. The procedure can be divided into three parts:
• Create a list of all possible connections: For each particle (spot) in snapshot𝑡look for other particles (spots) that are closer than the search radius 𝛿𝑟maxin all subsequent snapshots from𝑡+ 1up to 𝑡+𝛿𝑡max. Do the search for all snapshots.
• Select those connections that are minimizing the squared distances between connected spots by solving a linear optimization problem (see below).
• Label connected spots with the same particle ID to save the connections.
Note that searching forward in time is enough, searching backward would give the same list of possible connections. This is also the reason why this tracking algorithm does not violate time-reversal symmetry.
How to construct the linear optimization problem is described in the next paragraph.
In the following the term “spot” or “point” is used for the position of one particle in a specific snapshot.
The task of tracking is thus to connect the spots to particle trajectories. The term “link” or “connection”
between two spots is used to describe that these are two spots in different snapshots that belong to the same particle trajectory. A “possible” link denotes a link that might connect two spots of the same particle, but the algorithm has not decided yet whether it considers this to be the case. A “selected” link denotes a link for which the algorithm has decided that it indeed correctly connects two spots of the same particle.
4.3.1 Linear optimization problem
Karrenbauer and Wöll [78] translated the tracking problem to a linear optimization problem (LP) and proofed the consistency of their solution. To keep it short and more simple the mathematical proofs are not given here, only the construction of the problem is presented.
Let us define𝑚as the number of possible links that have been found in the first step of the whole procedure (see above) and𝑛as the total number of spots in the data (including all snapshots of the measurement).
A general linear optimization problem can be written as:
optimize𝐱 for 𝐜⋅𝐱⟶min. while 𝐌 𝐱=𝐛 and 0≤𝐱≤1, (4.4)
4.3 Tracking according to Karrenbauer and Wöll Here𝐱is the(𝑚+2𝑛)-dimensional solution vector,𝐌is the2𝑛× (𝑚+2𝑛)matrix of constraints,𝐛is the 2𝑛-dimensional constraint vector and the product𝐜⋅𝐱describes the linear cost function, which is defined by the(𝑚+2𝑛)-dimensional cost vector𝐜. The general idea of an LP is to minimize the cost function by finding a globally optimal solution vector𝐱. The set of linear equations𝐌 𝐱=𝐛sets the rules for𝐱. Any possible solution must obey it.
The solution vector𝐱is defined to be a list of ones and zeros, that is divided into three parts: The first𝑚 entries tell us whether a possible link𝑘in the list of possible link is selected (then𝑥𝑘 = 1) or not (then 𝑥𝑘 = 0). The next𝑛elements provide a value for each spot in the data: A “1” indicates that the spot is the starting point of a trajectory and therefore has no predecessor; “0” means that it is not a starting point and therefore has a predecessor. Similarly, the last𝑛elements also provide a value for each spot: “1”
means that the spot is the endpoint of a trajectory and therefore has no successor; “0” means that the spot is not the endpoint of trajectory and has a successor. For instance, if spot number𝑙is a starting point (has no predecessor) we have𝑥𝑚+𝑙 = 1. The last2𝑛elements are not actually part of the solution (the set of selected connections would be enough), but they are necessary for the computation of the cost function.
We want to punish any point without successor and any point without predecessor with additional costs, so that connected points become more likely.
The cost vector𝐜has the same structure as the solution vector: It contains the costs for each of the possible links in the first𝑚entries, the costs for a starting point in the next𝑛entries and the costs for an endpoint in the last𝑛entries. With the above definition of𝐱 as a list of ones and zeros, the costs for any set of selected links, starting and endpoints is thus given by the scalar product𝐜⋅𝐱.
More complicated is the formation of the constraint matrix𝐌 and the constraint vector 𝐛. The rules that are enforced by the linear equation𝐌 𝐱 =𝐛must guarantee that always a valid set of connections, starting points and endpoints is marked as selected. In spoken words, the rules are:
1. Each point has at most one successor (may only be the left part of exactly one connection) 2. Each point has at most one predecessor (may only be the right part of exactly one connection) 3. Each point is either a starting point or it has a predecessor (is the right part of a connection) 4. Each point is either an endpoint or it has a successor (is the left part of a connection)
In the terminology used here a connection between two spots is a link that connects a spot in one snapshot 𝑡1(left part) with another spot in a later snapshot𝑡2(right part). Note that these rules also allow uncon-nected (single) spots; they are at the same time starting points and endpoints of a track and therefore have no predecessor and no successor.
All2𝑁 elements of the constraint vector𝐛on the right hand side of the linear equation (4.4) are set to
“1”. As a consequence the whole set of rules has to be encrypted in the constraint matrix𝐌in a way that the multiplication of a valid solution vector𝐱with any of𝐌’s rows must be equal to one. Like the cost vector and the solution vector, the matrix is also divided into three parts. This is seen in the definition of its columns in Equation4.5.
The first𝑚columns provide values for the each of the𝑚possible links. The constraint matrix𝐌connects each of the possible links (corresponding to the column number𝑘) to the left and right participating spots of the position data (corresponding to the row numbers𝑟and𝑛+𝑙): The𝑘th link connecting spot number 𝑟to spot number 𝑙 is encoded as two “1”s in column 𝑘, one in row𝑟 and one in row𝑛 +𝑙, so that 𝑀𝑟,𝑘 = 𝑀(𝑛+𝑙),𝑘 = 1. Accordingly, for each possible link there is a “1” in the upper part (first𝑛rows) and another “1” in the lower part (last𝑛rows) of𝐌, while the remaining matrix elements stay zero.
𝐌=
With this definition, rules 1 and 2 are already encoded in the matrix. To illustrate this, let us look at rule 1 which allows only one successor: If the optimizer would select the𝑘th and the𝑗th link in the list of possible links, we would have𝑥𝑘 = 1and𝑥𝑗 = 1in the solution vector. If both these links would have spot number𝑙as their left part, spot𝑙would have two different successors. But the multiplication with the(𝑛+𝑙)th row of𝐌would be equal to2, since connection𝑘leads to𝑀(𝑛+𝑙),𝑘 = 1and connection𝑗 leads to 𝑀(𝑛+𝑙),𝑗 = 1, so that𝑀(𝑛+𝑙),𝑘𝑥𝑘+𝑀(𝑛+𝑙),𝑗𝑥𝑗 = 2. This violates the constraints. Due to the definition𝐛 = 𝟏,1must be the result of the multiplication of𝐱with any of𝐌’s rows. Therefore, only one of the two connections may be selected to fulfill the constraints.
The next𝑛columns of𝐌(corresponding to column number𝑚+𝑟) are there to ensure that a spot that is already the right part of a connection (a spot with a predecessor) cannot be the starting point of a track (rule 3). This is done by setting𝑀𝑟,(𝑚+𝑟) = 1. The multiplication of row𝑟with a solution vector𝐱that has spot number𝑟enabled as a starting point (𝑥𝑚+𝑟 = 1) then already has a “1” in the sum, namely the product𝑀𝑟,(𝑚+𝑟)𝑥𝑚+𝑟 = 1. An additionally selected connection number𝑘with position𝑟as right part would yield another “1”, since then also𝑀𝑟,𝑘𝑥𝑘 = 1. The complete multiplication of row𝑟with𝐱would then be equal to 2, violating the constraints. Analogously, with the same arguments, the last𝑁 columns of𝐌ensure the validity of rule 4.
Another requirement is that the linear equation𝐌 𝐱=𝐛should only allow zeros and ones in the solution vector 𝐱, other values make no sense. As posed above, the linear optimization problem only sets the constraint ∀𝑖 ∶ 0 ≤ 𝑥𝑖 ≤ 1. The trick is that the minimization of the cost function 𝐜 ⋅𝐱 enforces binary values for all 𝑥𝑖. Let us look at an example: Due to the constraints set by the linear equation, selecting connection𝑘with left part𝑙 and right part𝑟leads to the equations𝑥𝑘 +𝑥𝑚+𝑟 = 1(row 𝑟in points and end points are more expensive than connections, so that setting𝑥𝑘 = 1would minimize the costs in this example.
4.3 Tracking according to Karrenbauer and Wöll
Figure 4.4: Without a time penalty the fragmentation of a track (blue and red) can be cheaper than the single track (green). This happens if the formation of one track involves paying for many longer spatial distances.
4.3.2 Cost function
In principle the definition of the linear optimization problem (Equ.4.4) does not put any restriction to the cost function𝐜⋅𝐱, there is a lot of flexibility. As discussed in the introduction to this chapter, the minimization of the squared distances between connected spots gives the most probable selection of links, at least for freely diffusing particles. This leads us to the following first guess for the cost vector:
𝑐𝑘=𝛿𝑥2𝑘+𝛿𝑦2𝑘+𝛿𝑧2𝑘 for 𝑘= 1,…, 𝑚 and 𝑐𝑘 =𝛿𝑟2max for 𝑘=𝑚+ 1,⋯, 𝑚+ 2𝑛. (4.6) Here the first𝑚entries of𝐜give the cost for all possible links and the last2𝑛entries are the punishment for a spot being a starting point or an endpoint. Note that Crocker and Grier [35] minimize the same function in their algorithm. The difference is that they only connect points from one snapshot with points from the next snapshot, instead of deciding on all possible connections in all snapshots at once. Since this cost function does not penalize the time interval between spots, the fragmentation of a track can sometimes be cheaper than a complete track. This is illustrated in Figure4.4.
To avoid fragmentation one includes the squared time difference𝛿𝑡2𝑘in the costs:
𝑐𝑘 =𝛿𝑥2𝑘+𝛿𝑦2𝑘+𝛿𝑧2𝑘+𝛿𝑡2𝑘 for 𝑘= 1,…, 𝑚 and 𝑐𝑘=𝛿𝑟2max+𝛿𝑡2max for 𝑘=𝑚+1,⋯, 𝑚+2𝑛.
(4.7) One can also adjust the time penalty by scaling𝛿𝑡2with an additional factor. This may become necessary if time units and space units differ a lot. For instance, if a particle on average moves20length units in one time step (𝛿𝑟𝑡 = 1), it can be a good idea to increase the time penalty by a factor of20. For freely diffusing particles, penalizing the time difference is contrary to the idea that expensive connections should be more unlikely. Actually the probability for a big spatial displacement increases with the time interval (see Equ.4.1). On the other hand it makes no sense to require lower costs for connections with bigger time intervals. The algorithm would produce tracks with as many gaps as possible and the fragmentation problem would take over completely. The concept of paying less for more likely connections only works in the ideal case, where a particle is constantly detected in each snapshot.
Another possibility is to weigh the costs for starting points and endpoints (points with not predecessor and/or no successor) with a factor proportional to the quality of the detection of that point, e.g. with its brightness. False positives (erroneously detected particles) would than have a lower probability to be included in a track. This possibility was not used within this work, since there were more ambiguities with false negatives (undetected particles) than with false positives.
There is no reason to restrict the dependence of the cost function to the displacements in space and time.
In addition one could use other particle features. For example the change in particle size𝛿𝜎 could be penalized to ensure that a particle does not change its size in the course of the trajectory. Again, this
was not used within this work since small and big particles can be separated well enough in the binary samples discussed here.
As pointed out above, the cost function presented here is designed for the tracking of freely diffusing particles. From the investigations in this work, one can say that glassy systems can also be tracked with good confidence with that function. For other types of motion like e.g. anomalous diffusion, drift or active transport, one could think of different cost functions leading to better results.
4.3.3 Joining
Karrenbauer and Wöll [78] also describe a method to deal with false positives. This is done by allowing connections in between a single snapshot (𝛿𝑡 = 0), so-called joins. In the LP, both directions of such a joining connection are added to the list of possible links. This means that a track can be continued with either of the two spots and the other spot is the predecessor for a possible continuation of the track.
In order to prohibit the selection of both directions at once, which would result in a very cheap circular track without ends, Karrenbauer and Wöll [78] proposed to raise the costs for a join by adding the mean penalty for being a starting point and an endpoint. This means that the costs for a circular track include paying for the circle being starting point and endpoint at the same time.
Obviously the search for joining neighbours in the snapshot should be limited to a smaller search radius 𝛿𝑟join< 𝛿𝑟max, so that spots being that close are supposed to come from a single real particle. Since one only knows something about detected spots, Karrenbauer and Wöll [78] argue that it is better to avoid false negatives at the expense of accepting more false positives. However, within this work it turned out that it is not easy to separate false positives from true detections. Furthermore, physical quantities calculated from space correlations suffer a lot more from false positives than from false negatives (cf. C.1). For instance, the pair distribution function does not change at all (only the statistics is worse), if one omits a random selection of particles in the computation.
4.3.4 Implementation
Originally, Karrenbauer and Wöll [78] implemented their tracking method only for two-dimensional data.
Their MATLAB code uses the IBM ILOG CPLEX Optimizer as linear programming solver, which is free of charge for academic use. Since it was necessary to extend the code to three dimensions anyway, it was also translated to Python/NumPy. This has the additional advantage that the code is easily embedded into the other software that was implemented within this work, for particle detection, computation of physical quantities, simulations and others. CPLEX smoothly integrates into Python as well, so that there is no significant difference in the computation speed comparing MATLAB and Python.
From computational point of view, a drawback of the method is clearly that the optimization has to be done for all questionable connections at once. A typical data from experiments with colloids, consist of4000particles in1000and more snapshots. The number of possible links easily exceeds10millions.
Although CPLEX works with sparse matrices and is tuned to use the smallest possible amount of memory, for large uncut data sets, the16GB machine used for the data analysis always runs out of RAM and is not able to run the code. This renders the method unusable for most of the measurements done here, but of course the code was tested for smaller records. For example a complex record with 500particles in 350snapshots is done in about 2 minutes on a 3 GHz machine, using∼ 5GB of RAM. The bottleneck is always the memory consumption not the CPU time. For the same data Crocker and Grier’s method only needs about 20 seconds using not much more RAM than necessary than to load the data from hard disk (∼ 200MB).