TrackMate Algorithms
This page documents the current components of TrackMate. TrackMate has a modular design, meaning that it is made of different modules that each have a specific role. Developers can build their own module and re-used the other ones and the GUI to achieve a quick development. The module types are (in the order you meet them when executing the plugin):
- Spot detectors. Taking your image data, they detect spots in them.
- Spot analyzers. Each spot can receive a wide range of features, calculated from their location, radius and the image data. For instance: max intensity in spot, rough morphology, etc... They are then used to filter out spurious spots and retain only good ones for the subsequent tracking step.
- Views. Display the segmentation and tracking results overlaid on your image data.
- Spot trackers. Take the filtered spots and link them together to build tracks.
- Track analyzers. Like for spot analyzers, but operate on whole track. Can be used to report track mean velocity, displacement, etc... They are also used to filter spurious tracks.
- Actions. Miscellaneous actions you can take on the complete result of the tracking process. It can be used to copy the track overlay to another image, launch a 3D viewer, export the results to a simple format, generate a track stack, etc...
We describe here the best we can the current modules that are shipped with TrackMate.
Contents
Spot detectors
Spot features generated by the spot detectors
Behind this barbaric name stand the part responsible for spot detection in your image. The algorithm currently implemented are very generic and naive. They will most likely fail for complicated case such as touching objects, very weak SNR, etc... The 3 of them present are all based on Laplacian of Gaussian filtering, which we describe below.
All detectors must at least provide the following common spot features:
- X, Y, and Z: the spot coordinates in space. Coordinates are expressed in physical units (μm, ...).
- R the spot radius, also in physical units. The current detectors only set this radius value to be the one specified by the user. More advanced detectors - yet to be implemented - could retrieve each spot radius from the raw image.
- Quality: The implementation varies greatly from detector to detector, but this value reflects the quality of automated detection. It must be a positive real number, large values indicating good confidence in detection result for the target spot. This sole feature is then used in the initial filtering step, where spots with a quality lower that a specified threshold are purely and simply discarded.
The two other time features - T and Frame number - are set by TrackMate itself when performing detection on all the timepoints of the target raw data. T is expressed in physical units, and the Frame number - starting from 0 - is simply the frame the spot was found in.
Difference of Gaussian particle detection (DoG segmenter)
Given d an approximate expected particle diameter, determined upon inspection, two gaussian filters are produced with standard deviation σ₁ and σ₂:
- σ₁ = 1 / (1 + √2 ) × d
- σ₂ = √2 × σ₁
The image is filtered using these two gaussians, and the result of the second filter (largest sigma) is subtracted from the result of the first filter (smallest sigma). This yields a smoothed image with sharp local maximas at particle locations. A detection spot is then created for each of these maximas, and an arbitrary quality feature is assigned to the new spot by taking the smoothed image value at the maximum. If two spots are found to be closer than the expected radius d/2, the one with the lowest quality is discarded.
To improve the localization accuracy, and extra step is taken to yield a sub-pixel localization of the spots. The position of each spot is recalculated using a simple parabolic interpolation scheme, as in [Lowe]. The quality feature is also interpolated using this scheme.
A large number of spurious spots are created by finding local maximas. There spurious spots are discarded inn extra step, by applying a threshold on the quality feature computed during segmentation. The value of this threshold is set manually, to match the SNR of the input image. Thresholded spots are then retained for subsequent particle-linking.
Spot analyzers
Spot features, such as 'Max intensity', 'Estimated diameter', etc., are calculated for all spots just after the initial filtering step. They are then used to select spots, based on filters set to retain only spots with a given feature below or above a specified threshold. By the way initial filtering is a good way to limit spot feature calculation time on spots you know to be spurious. With the current features, this not such a reason to do it: most feature analyzer are cheap computationally.
Mean, Median, Min, Max, Total intensity and its Standard Deviation
The plain statistical estimates are simply calculated from all the values for pixels within the physical radius from the spot center.
Contrast & Signal/Noise ratio
This contrast followed Michelson contrast definition: C = (I_in - I_out) / (I_in + I_out)
where I_in
is the mean intensity inside the spot volume (using the physical radius), and I_out
is the mean intensity in a ring ranging from its radius to twice its radius.
The spots's SNR is computed a (I_in - I_out) / std_in
where std_in
is the standard deviation computed within the spot.
These two values depend greatly on the correctness of the radius of each spot. Negative values might be caused by incorrect radius.
Estimated diameter
This feature estimates an optimal diameter for each spot, based on contrast calculation.
The mean pixel intensity is calculated in 20 concentric, tangent rings, centered at the target spot location, and with radiuses ranging from a 10th of the spot radius to twice its radius. The contrast at a given radius is defined as the difference between the mean intensity of a ring with inner radius the radius sought, and the previous ring.
The estimated diameter is then defined as the radius that maximizes this contrast. A quadratic interpolation is performed to improve diameter accuracy.
Spot trackers
LAP trackers
The Linear Assignment Problem (LAP) trackers implemented here follow a stripped down version of the renowned method contributed by Jaqaman and colleagues^{[1]}. We repeat here the ideas found in the reference paper, then stresses the differences with the nominal implementation.
In TrackMate, the LAP framework backs up two instances of a tracker (four if you count the old ones, prior to v2.5.z):
- the Simple LAP tracker;
- the LAP tracker.
The first one is simply a simplified version of the second: it has less settings and only deal with particle that do not divide nor merge, and ignores any feature penalty (see below).
All the linking costs for these two trackers are based on the particle-to-particle square distance. If this tracker had to be summarized in one sentence, it would be the following: The SImple LAP tracker and the LAP tracker are well suited for particle undergoing Brownian motion. Of course, they will be fine for a non-Brownian motion as long as the particles are not too dense.
Particle-linking happens in two step: track segments creation from frame-to-frame particle linking, then track segments linking to achieve gap closing. The mathematical formulation used for both steps is linear assignment problem (LAP): a cost matrix is assembled contained all possible assignment costs. Actual assignments are retrieved by solving this matrix for minimal total cost. We describe first how cost matrices are arranged, then how individual costs are calculated.
Cost matrix for frame-to-frame linking
In the first step, two consecutive frames are inspected for linking. Each spot of the first frame is offered to link to any other spot in the next frame, or not to link. This takes the shape of a (n+m) x (n+m) matrix (n is the number of spots in the frame t, m is the number of spots in the frame t+1), that can be divided in 4 quadrants.
- The top-left quadrant (size n x m) contains the costs for linking a spot i in the frame t to any spot j in the frame 't+1.
- The top-right quadrant (size n x n) contains the costs for a spot i in the frame t not to create a link with next frame (yielding a segment stop).
- The bottom-left quadrant (size m x m) contains the costs for a spot j in the frame t+1 not to have any link with previous frame (yielding a segment start).
- The bottom-right quadrant (size m x n) is the auxiliary block mathematically required by the LAP formalism. A detailed explanation for its existence is given in the supplementary note 3 of ^{[1]}. This quadrant is built by taking the transpose of the top-left quadrant, and replacing all non-blocking costs by the minimal cost.
Solving LAP
To solve this LAP, we rely on the Munkres & Kuhn algorithm^{[2]}, that solves the problem in polynomial time (O(n^{3})). The algorithm returns the assignment list that minimizes the sum of their costs.
The frame-to-frame linking described above is repeated first for all frame pairs. This yields a series of non-branching track segments. A track segment may be start or stop because of a missing detection, or because of a merge or split event, which is not taken into account at this stage. A second step where track segments are offered to link between each other (and not only spots) is need, and described further down.
Calculating linking costs
In calculating costs, we deviate slightly from the original paper from Jaqaman et al.^{[1]}. In the paper, costs depend solely on the spot-to-spot distance, possibly weighted by the difference in spot intensity. Here, we offer to the user to tune costs by adding penalties on spot features, as explained below.
The user is asked for a maximal allowed linking distance (entered in physical units), and for a series of spot features, alongside with penalty weights. These parameters are used to tune the cost matrices. For two spots that may link, the linking cost is calculated as follow:
- The distance between the two spots D is calculated
- If the spots are separated by more than the max allowed distance, the link is forbidden, and the cost is set to infinity (i.e the blocking value). If not,
- For each feature in the map, a penalty p is calculated as
where W is the factor associated to the feature in the map. This expression is such that:- there is no penalty if the 2 feature values f1 and f2 are the same;
- with a factor of 1, the penalty is 1 is one value is the double of the other;
- the penalty is 2 if one is 5 times the other one.
- All penalties are summed, to form P = (1 + ∑ p )
- The cost is set to the square of the product: C = ( D × P )²
If the user feeds no penalty, the costs are simply the distances squared.
Calculating non-linking costs
The top-right and bottom-left quadrant of the frame-to-frame linking matrix contain costs associated with track segment termination or initiation (a spot is not linking to a spot in the next or previous frame). Each of these two blocks is a square matrix with blocking value everywhere, except along the diagonal for which an alternative cost is computed. Following Jaqaman^{[1]}, this cost is set to be
Calt = 1.05 × max( C )
where C is the costs of the top-left quadrant.
Cost calculation & Brownian motion
Without penalties and with a maximal linking allowed distance, the returned solution is the one that minimizes the sum of squared distances. This actually corresponds to the case where the motion of spots is governed by Brownian motion. See for instance Crocker and Grier^{[3]}.
By adding feature penalties, we aim at favoring linking particles that "resemble" each other. In brute single particle linking problems, spots are generally all the same, and they only differ by position. However, there is a variety of problems for which these feature penalties can add robustness to the tracking process.
For instance, we originally developed TrackMate for semi-automated lineaging of C.elegans embryos, using a strain fluorescent in the nucleus. Cells that are dividing have a fluorescence distribution which is very different from non-dividing cells, and this can be exploited for robust tracking.
Track segment linking
In a second step, the track segments built above are offered to link between each other. Jaqaman and colleagues proposes to exploit again the LAP framework for this step. A new cost matrix is generated, but this time the following events are considered:
- The end of a track segment is offered to link to any other track segment start. This corresponds to gap-closing events, where a link is created typically over two spots separated by a missed detection.
- The start of a track segment is offered to link to the spots in the central part (not start, not end) of any other track segment. This corresponds to splitting events, where a track branches in two sub-tracks.
- The end of a track segment is offered to link to the spots in the central part of any other track segment. This corresponds to merging events, where two tracks merges into one. I yet did not meet this case in Life-Sciences.
- A spot part of any track segment is offered not to create any link.
The second cost matrix has a shape that resembles the first cost matrix, calculated for frame-to-frame linking, and which is best described in the original article.
As before, we modified the way costs are calculated, and re-used the feature penalties framework described above. Also, the user must provide on top a maximal time-difference to link, over which linking will be provided. Careful: this maximal time is expressed in physical units and not in number of frames.
Main differences with the Jaqaman paper^{[1]}
The nominal implementation of the paper remains the one developed under Matlab by Khuloud Jaqaman et al. and published in Nature Methods. The software is called u-track and can be found on Khuloud Jaqaman homepage.
TrackMate was initially developed to simplify C.elegans lineaging. It therefore just bundles a stripped down version of this framework.
The notable differences are:
- The LAP framework is generic: Jaqaman and colleagues proposed a framework to approximate multiple-hypothesis tracking solutions using linear assignment problems. One just need to provide the link cost matrix. TrackMate properly implements the LAP framework, but the cost matrix calculation - which is specific to each problem - is much more simpler than in u-track.
For instance, in TrackMate all link costs are based on the square distance between two spots (weighted or not by feature differences, see above), which make it tailored for Brownian motion. In u-track, the user is proposed with different motion types, including a linear motion whose parameters are determined on the fly. See for instance CD36 tracking, in the supplementary note 7 of the paper^{[1]}. - In u-track, merging and splitting of tracks are used to describe two particles that temporally overlap spatially. These events' costs are weighted by the two particle intensities to properly catch the apparent increase in intensity due to the overlap. In TrackMate, we use splitting events to describe cell divisions, as we developed it initially to deal with C.elegans lineages. However is seems than Jaqaman and colleagues used it the same way to investigate CD36 dissociation and re-association.
- In TrackMate, distance and time cutoffs are specified manually by the user. In u-track they are derived for each particle automatically, providing self adaptation.
Linear motion tracker.
Starting from TrackMate v2.7.z, we shipped a new tracker that can deal specifically with linear motion, or particle moving with a roughly constant velocity. This velocity does not need to be the same for all particles. You can find it in TrackMate tracker selection under the name Linear motion LAP tracker.
Though it deals with a completely different motion model compared to the LAP trackers in TrackMate, it reuses the Jaqaman LAP framework, and it is similar to a tracker proposed in the Jaqaman paper as well: See the CD36 tracking, in the supplementary note 7 of the paper^{[1]}. But again, the version in TrackMate is simplified compared to what you can find in u-track.
Principle.
The linear motion tracker relies on the Kalman filter to predict the most probable position of a particle undergoing constant velocity movement.
Tracks are initiated from the first two frames, using the classical LAP framework with the Jaqaman cost matrix (see above), using the square distance as cost. The user can set what is the maximal distance allowed for the initial search with the Initial search radius setting.
Each track initiated from a pair of spots is used to create an instance of a Kalman filter. There are as many Kalman filters as tracks. In the next frames, each Kalman filter is used to generate a prediction of the most probable position of the particle. All these predictions are stored.
Then, all the predicted positions are linked against the actual spot positions in the frame, using again the Jaqaman LAP framework, with the square distance as costs. The user can set how far can be an actual position from a predicted position for linking with the Search radius setting.
Now of course, after linking, some Kalman filters might not get linked to a found spot. This event is called an occlusion: the predicted position did not correspond to an actual measurement (spot). The good thing with Kalman filters is that they are fine with this, and are still able to make a prediction for the next frame even with a missing detection. If the number of successive occlusions is too large, the track is considered terminated. The user can set the maximal number of successive occlusions allowed before a track is terminated with the Max frame gap setting.
Conversely, some spots might not get linked to a track. They will be used to initiate a new track in the next frame, as for the tracker initiation described above.
It is important to note here that the cost functions we use is the square distance, like for the Brownian motion, but from the predicted positions to the actual detections. Because the prediction positions are made assuming constant velocity, we are indeed dealing with an adequate cost function for linear motion. But since we are linking predicted positions to measured positions with the square distance cost function, we do as if the predicted positions deviate from actual particle position with an error that follows the gaussian distribution. This is a reasonable assumption and this is why this tracker will be robust.
Implementation
The code can be found on on GitHub. We now repeat the section above in pseudo-language. When you see the word link below, this means:
- Take all the source detections in frame t and the target detections in frame t+1.
- Compute the costs for all possible physical assignment (potential links) between source and target detections and store them in the cost matrix.
- Solve the LAP associated to this matrix.
- Create a link for each assignment found.
The particle linking algorithm would read as follow:
- Initialization:
- Link all the detections of frame 0 to the detections of frame 1, just based on the square distance costs (for instance).
- From each of the m links newly created, compute a velocity. This velocity is enough to initialize m Kalman filters.
- Initialize m tracks with the found detections and links, and store the associated Kalman filters.
- Track elongation:
- For each Kalman filter, run the prediction step. This will generate m predicted positions.
- Link the m predicted positions to the n detections in frame 2, based on square distance.
- Target detection that have been linked to a predicted position are added to the corresponding track.
- The accepted target detection is used to run the update step of the Kalman filter.
- Loop to next frame.
- Track termination:
- Some of the m predicted position might not find an actual detection to link to. In that case, we have an occlusion. The algorithm must decide whether it has to terminate the track or to bridge over a gap.
- If the number of successive occlusions for a Kalman filter is below a certain limit (typically 2 to 10), the track is not terminated, and the filter goes back to the track elongation step. Hopefully, from the new prediction a target particle will be found, and the detection in frame t will be linked to a detection in frame t+2 (or t+3 etc).
- Otherwise, the track is terminated and the Kalman filter object is dropped.
- Track initiation:
- Conversely, some detections in frame t+1 might not be linked to a predicted position. In this case, these orphan detections are stored to initiate a new track. But for this, other orphans detections are needed in frame t+2.
- This step is identical to the initiation step, but for subsequent frames. It requires to store orphan detections in current and previous frames.
- In frame t+2, priority must be given to detections linked to the predicted positions by the Kalman filters over orphan detections of frame t+1. So when you deal with frame t+2, you perform first the track elongation step, get a list of orphan detections in frame t+2, and then combine it to the orphan detections in frame t+1 to initiate new Kalman filters.
References
- ↑ ^{1.0} ^{1.1} ^{1.2} ^{1.3} ^{1.4} ^{1.5} ^{1.6} Jaqaman et al., "Robust single-particle tracking in live-cell time-lapse sequences", Nat Methods. 2008 Aug;5(8):695-702.
- ↑ J. Munkres, "Algorithms for the Assignment and Transportation Problems", Journal of the Society for Industrial and Applied Mathematics, 5(1):32–38, 1957 March
- ↑ Crocker and Grier. "Methods of Digital Video Microscopy for Colloidal Studies." J Colloid Interf Sci (1996) vol. 179 (1) pp. 298-310