The identification of inliers is the key challenge in the kind of experimental data-set that we are discussing. When single-particle tracking methods require pre-registration, they may not be robust against lack of guaranteed inliers. On the other hand, when single-particle tracking does not need pre-registration, it is typically based on models of the objects’ motion, which are not optimal with transient objects and with low signal-to-noise ratio across extended image sequences.

We treated the position, appearance and disappearance of transient objects as a random spatial process. The frame-to-frame displacement of objects embedded in a tissue can be described as the effect of two random variables: (1) large-scale movements of the tissue within the field of view, due for example to its growth or to the inability of the microscope to focus on a fixed point of the field of view over time; (2) small-scale fluctuations of the objects of interest within the tissue.

The main purpose of the algorithm is to find the best large-scale transformation that explains most of the object correspondence, despite the noisy small-scale dynamics. The strategy adopted was to find correspondences between the random process at time t and the random process at time t-τ where τ is a lag variable. In this work, we only considered a τ=1 step, which is equivalent to approximating the tissue as effectively rigid at the given temporal scale. As mentioned above, objects that appear in both frames are called inliers and those that do not appear in both frames are called outliers and cannot be linked by a global frame-frame transformation. For example, outliers can be either debris, i.e. any object which is determined not to be an object of interest, or could be an object of interest that has just exited or entered the new frame.

There are two main steps in the execution of the algorithm: (i) candidate transformations are generated from the data; (ii) among these candidates, the optimal transformation is chosen with respect to some previously defined objective. Inliers are then mapped onto each other by the chosen global transformation and labelled (with identifiers). Finally, new identifiers are generated for the outliers within the region of interest. These new identifiers correspond to proposed objects of interest that may be filtered at other stages of the process as discussed later.

Algorithm 1 provides a terse, high-level overview of the method. The key algorithm stages, which amount to finding the best global transformation from the data, are discussed below.

Consider the sets of objects U(t),V(t-τ) at times offset by a lag τ. A subset of these corresponds to true objects of interest (i.e. not noise or debris) and a further subset of these exists in both frames. Given N objects at time t and M objects at time t-τ, a number of permutations P(N,k) are sampled from U and a number of permutations P(M,k) are sampled from V. We limited our method to k∈{1,3}, thus sampling either single objects (k=1) or triangular constellations (k=3) from each frame. Object identifiers take a natural ordering so that any set of the same k object identifiers are considered equivalent regardless of their ordering.

Let n be the number of k-constellations ci(t) sampled from U(t) and m be the number of k-constellations cj(t′=t-τ) sampled from V(t-τ). In general, for two tensors M1 and M2 related by M1=AM2, the affine transformation A can be uniquely determined through a least squares method (this is, in general, over-determined). Here, the transformation γij is the affine transformation, taking the set ci(t) to the set cj(t′): cj(t′)=∑γijci(t)

In lines 3–6 of the Algorithm 1, the constellation are sampled in frames U and V. Translations are then generated for each pair of objects v∈V and u∈U, and added to the list of candidate transformations.

Constellation sampling (line 9) is preferably seeded with linkages found by translations so that likely outliers are excluded from consideration. If a given translation maps n objects in U to n objects in V, then these object-object pairings can be considered candidate linkages. This not only reduces the candidate objects to those likely inliers, but constellation congruences can be identified as well, so the same set of objects in one frame can be paired to the same set of objects in the other frame when proposing transformations. Constellation pairs P(NL,3),P(ML,3) can then be sampled from the L inliers appearing in both frames. These will be used to generate proposal affine transformations, to provide a more optimal fit than translations alone.

Transformation ranking was carried out using a modified least-squares-loss function. The least-squares-loss method is generally described as the minimization of ∑|y-y^|2, where y is a proposal vector and y^ is the target vector. In our case, the proposal vectors are objects in the frame U(t) transformed back to t-τ by one of the proposed transformations. The target vectors are actual objects sampled in V(t-τ). The least-squares-loss objective was modified based on the following prior: (i) we expect a global rigid transformation to explain the movement of object centroids up to small fluctuations due to individual object movement which are deemed negligible. (ii) A distance >ε between objects’ centroids is assumed. The value of ε was determined from the data and is always slightly larger than the average object radius. When considering distances between proposal objects and target objects, only distances to the first nearest neighbour within a radial distance ε were considered. Each candidate transformation γ was applied to all objects ui(t)∈U(t),i∈{1,2,...,N}. Let u~i=γ(ui). We used a cost function

where nnV(U~i) is the nearest neighbour position in V(t-τ) to the projected object U~i. If the projected object had no nearest neighbour within a sphere of radius ε, the capped distance ε+1 was attributed. This is illustrated in Fig. 2.

The cost function will rank transformations by how well they explain the movement of the majority of objects. If there are small differences in scores, we may prefer transformations that have the smallest displacement. For this reason, we can regularise this cost function to add a penalty to larger displacements.

Noise and histogram levels were used on normalised grey-scale images taking pixel values in [0, 1]. The noise shown in the top panel of Additional file 5: Fig. S2 was estimated using a wavelet-based estimator of the Gaussian noise standard deviation [35, 36]. The middle panel in Additional file 5: Fig. S2 describes the signal. This is a simple quasi-signal metric defined as a ratio between the 99th and 95th histogram percentile boundary values, which are determined empirically for our data to separate background from signal. When this value approaches 0, the corresponding signal appears very low. For values above a calibrated threshold, the corresponding signal lies outside the 99th percentile of the data, which typically happens for isolated bright points in the image. A noise range was used to guide downstream strategies in the image processing pipeline. An empirically determined “normal” noise band [0.01, 0.045] was used in the sample parameters discussed here. When the noise was below 0.01 no denoising was necessary, otherwise wavelet denoising was applied. Noise exceeding 0.045, typically corresponded to light saturation and/or low signal-to-noise ratio, and was due to a low number of mitotic events. In such instances additional thresholding was applied to the data before wavelet denoising. When the noise exceeded 0.1, the signal-to-noise ratio was so low that the frame was marked as degenerate. Additional file 1: Fig. S1 contrasts different noise levels in sample images and Additional file 5: Fig. S2 plots image properties over time for a given experiment.

A 2D projection of the data was obtained by summing the 3D tensor along the z-axis. An adaptive threshold based on image histogram percentile ranges was used to construct a narrow-band filter for the 2D data. This range can vary erratically between frames as shown in the bottom time series in Additional file 5: Fig. S2. A thin slice of the image data range (shaded band) was selected and an aggressive “Gaussian smoothing” (averaging neighbouring pixel values using a Gaussian kernel). An empirically-determined sigma value (σ=8) was used to find a mesh-like connected component corresponding to the region of activity in the root tip. The 2D mask of this largest component was extended to a 3D mask by projecting the 2D region back into the z-plane. The largest connected component in the thresholded image was identified as the root tip. These stages are illustrated in Additional file 6: Fig. S6. The image processing pipeline continues to process data only within the ROI.

For cases where the region of activity is well-isolated (smaller region of interest), the processing is more efficient as the amount of volumetric data is reduced and debris beyond the root are filtered from the downstream pipeline. For our sample data, an ROI with area ≈105 pixels corresponded to a well-isolated root tip. Noise and light saturation can affect this part of the process making the ROI area a proxy for image quality. While the actual area will depend on the data, it should vary smoothly and be relatively small unless root activity genuinely extends to the entire field of view.

We performed a “pre-segmentation” procedure to detect blob centroids. Segmentation plays a central role in many image processing pipelines and typically involves (i) thresholding and identifying background, (ii) using distance/gradient transformations with peak detection to identify markers and (iii) routines such as watershed [37] to segment blob labels. For our data, we have found it appropriate not to carry out the final segmentation. Instead, we carried out the pre-segmentation steps from thresholding to peak detection in detecting object centroids. Given the variability in the data over a large frame sequence, we have found a simple “annealing thresholding” to be effective. This simply increases a threshold iteratively so as to remove large connected components beyond a maximum perimeter length (empirically determined from data to be in the range 1000–2000 pixels) thus removing the background from the image. Being iterative, this is not as efficient as simpler linear filters but it is simple and robust to noise and variability. In extreme cases where it was not possible to find a threshold level in this manner, we treated the 99th percentile as image background and removed it.

Having applied the threshold filter, we identified centroids by (i) performing a difference of Gaussians to emphasise blob-like objects, (ii) applying a maximum filter and (iii) returning the coordinates of the local maxima (peaks) in the image. To allow for a fully automated routine that can cope with arbitrary datasets (in the scope of our light-sheet microscopy datasets) the emphasis in the centroid detection stage has been to avoid spurious centroid detections at the risk of under-sampling, while optimising for objects to be identified for at least two frames somewhere during the peak of their light intensity arc over time. A video with the centroids identification (red = centroids identified in the previous frame; white = centroids identified in the current frame) in sync with the raw time-lapse data is available here [31].

Objects in the point cloud were expected to be separated by a minimum distance ε such that under global transformation single objects are matched.The parameter ε is of central importance at all stages of the image processing and tracking pipeline and can be chosen based on the data. For sparse point-cloud data there will typically be only one match within the ball of radius ε. An ε value can be chosen such that it is possible to find multiple objects matched to one object. The interpretation in this case is that all matched points in the catchment region correspond to the same object and that multiple points are the result of e.g. image processing anomalies. Merging lineages can be useful to avoid generating superfluous new identifiers. This adds robustness in the event of noisy data or image processing anomalies. Alternatives to using a ball of radius ε might consider anisotropic catchment regions that factor in the direction of the global transformation. The parameter ε value may result in the tracker merging objects that the object detection stage discriminates between if those objects are not sufficiently separated. In the case of accurately counting mitotic events, this reduces the introduction of false-positives i.e. the erroneous generation of new identifiers

Object position updates were treated as small Gaussian perturbations of a global transformation. If the position of an object is given as

where γ is the translation operator, τ is the lag and ξi∼N(0,σ0) is the random fluctuation for one object, taking the “consensus” effectively filters out the noise to reveal the global movement. For example, the top m translations could be averaged to produce a new translation.

Transient events were expected to persist for a minimum number of frames greater than or equal to 1. For our data the minimum number of frames is 2. Objects that have been identified for less than 2 frames, termed single-frame detections, were excluded from the tracking result.

In each frame, the objects of interest were expected to be found in the center of the field of view. Objects may be detected near the original frame boundaries depending on the size of the ROI. Objects can be treated as inliers/outliers based on standard outlier detection methods under the assumption of normality (for example within 1 or 2 standard deviations from the mean coordinate).

Object location and angle differences between object displacements and the global transformation may be used in a post-processing stage to remove misclassified outliers. In the PRIORS routine in Algorithm listing 1 additional transformations beyond the transformations generated from the data can be chosen based on modelling assumptions. Additionally, the filter routines on lines 16–18 can filter data from the final result rather than during the tracking process based on modelling assumptions.