Rapid, topology-based particle tracking for high-resolution measurements of large complex 3D motion fields

Spatiotemporal tracking of tracer particles or objects of interest can reveal localized behaviors in biological and physical systems. However, existing tracking algorithms are most effective for relatively low numbers of particles that undergo displacements smaller than their typical interparticle separation distance. Here, we demonstrate a single particle tracking algorithm to reconstruct large complex motion fields with large particle numbers, orders of magnitude larger than previously tractably resolvable, thus opening the door for attaining very high Nyquist spatial frequency motion recovery in the images. Our key innovations are feature vectors that encode nearest neighbor positions, a rigorous outlier removal scheme, and an iterative deformation warping scheme. We test this technique for its accuracy and computational efficacy using synthetically and experimentally generated 3D particle images, including non-affine deformation fields in soft materials, complex fluid flows, and cell-generated deformations. We augment this algorithm with additional particle information (e.g., color, size, or shape) to further enhance tracking accuracy for high gradient and large displacement fields. These applications demonstrate that this versatile technique can rapidly track unprecedented numbers of particles to resolve large and complex motion fields in 2D and 3D images, particularly when spatial correlations exist.

Supplementary Note 6 Time-lapse particle motion with particles randomly disappearing and appearing 5 Supplementary Note 7 T-PT algorithm workflow 7 Supplementary Figure 1 Synthetically generated 3D volume of size 512x512x192 voxels consisting of 50,000 spherical particles. Supplementary Figure 4 Tracking performance of T-PT for an applied simple shear deformation with a varying level of non-affine deformation. 4 Supplementary Figure 5 Performance of various SPT algorithms in tracking random particle motion. 5 Supplementary Figure 6 Performance characterization of T-PT in tracking particles across time-lapse image series where particles randomly disappear and appear between image frames.

XZ Plane
Scaled Image Intensity 0 1 Figure 1: A simulated 3D volume of size 512x512x192 voxels consisting of 50,000 randomly distributed spherical beads.
x x Figure 2: A synthetic image of a single simulated spherical particle.

Nyquist resolution calculation
For a given particle seeding density in the image, ρ o , the mean particle distance in the images is, for 2D and 3D images respectively [1]. Here, ρ o,2D is the particle seeding density in the 2D image defined as the number of particles per unit image size and ρ 0,3D is the particle seeding density in the 3D image defined as the number of particles per unit 3D image volume.
According to the Nyquist-Shannon sampling theorem, the spatial resolution of the recovered displacement field is determined by half the sampling frequency. Therefore for a given mean particle distance, r, a first order estimate for the Nyquist spatial resolution of the recovered displacement field is 2r.

Mean displacement error versus signal-to-noise ratio
In Supplementary Figure 3, we plot the mean displacement error for our T-PT technique against synthetically generated images across a range of Poisson-distributed noise to account for variable experimental settings. The mean displacement error is defined as the mean of the difference between the imposed versus the measured particle displacement vector magnitude. The plot demonstrates that T-PT has low tracking error even for images with relatively high noise levels.

T-PT performance for non-affine motion in shear deformation
In the results section, we show the ability of T-PT to recover non-affine shear deformations. Here, we showcase tracking performance for the same shear deformation but for a range of different non-affinity. To define non-affinity in the images, we introduce a non-affinity parameter, N D, defined as the ratio of maximum non-affine displacement magnitude to the mean particle separation distance. As in the main text example, the non-affine shear deformation is generated by applying a homogeneous simple shear deformation of 10% nominal shearing strain along the x − y direction of a volumetric image of size 512x512x192 voxels. Superimposed on the affine displacement field is a non-affine displacement field of a truncated zero-mean normal distribution to two standard deviations in a range of standard deviations along the x, y and z directions. The displacement field was prescribed to the 50,000 particles randomly embedded in the synthetically generated images. T-PT is used to track this non-affine shear deformation from these synthetically generated images. In Supplementary Figure 4, we plot the recovery and mismatch ratio for T-PT for a range of non-affinity in the particle motion field, which shows that T-PT can accurately recover very large non-affine shear deformations.

Random particle motion
The tracking performance of T-PT for random particle motion was evaluated by simulating random particle motion in 3D images. We synthetically generated a 3D volume of 512x512x192 voxels seeded with 50,000 identically-sized spherical particles with an SNr of 25. The particles were prescribed with a random displacement vectors of a truncated zero-mean normal distribution to two standard deviations in a range of standard deviations along the x, y and z directions. We evaluated the performance of T-PT and other current SPT techniques as a function of increasing displacement parameter d for the prescribed displacement field (Supplementary Figure 5). For d < 0.25, T-PT and TrackMate exhibit very high recovery ratios η r ∼ 1. As d was increased further, the η r of TrackMate remains close to unity out to d ∼ 0.4, while T-PT decreases to η r ∼ 0.8. In comparison, Legant et al.'s method maintains a stable η r ∼ 0.97 until d ∼ 0.3, but with increasing d, recovery ratio decreases to η r ∼ 0.5. FVRM starts with a lower η r ∼ 0.9, which decreases with increasing d (Supplementary Figure 5). These results show that TrackMate, followed by T-PT, enables the most accurate tracking of particles undergoing random motion. However, the fast execution time of T-PT can make it preferable even for certain purely random motion applications, especially those involving small random particle motion at very high particle numbers.
6. Time-lapse particle motion with particles randomly disappearing and appearing We evaluate the performance of T-PT to track particle motion in a time-lapse image series in which particles randomly disappear and appear between image frames. By employing synthetic images we mimicked experimental particle tracking conditions in which particles can disappear, due to moving out of the image frame, detection failure or particle merging, or in which particles can appear due to moving into the image frame, failure of particle detection in the previous image frame or particle splitting. In the simulations, identical spherical particles were seeded in a space large enough such that the final image of 512x512x192 voxels region used as input for the algorithm was fully seeded with approximately 50,000 particles throughout the application of the motion field. We generated a time-lapse image series for six time points. The particles were prescribed Additionally, we randomly seeded and removed a varying fraction of particles in each frame to simulate particles appearing and disappearing in experimental multi-time point particle tracking. The number of particles seeded and removed in each image frame was picked from a Gaussian distribution having a mean equal to the number of particles intended to be seeded and removed with a standard deviation of 20% of the mean value.
To track particles across the time-lapse data, T-PT matches particles between two consecutive image frames. The particle trajectories across multiple frames are later merely computed from joining the frame-to-frame particle tracking results. In Supplementary Figure 6, we plot the recovery and mismatch ratios of the particles tracked between the initial image frame (t = 0) to the current image frame (t = t o ). The recovery ratio of the algorithm decreases with increasing image frame. This is expected because T-PT computes particle trajectories from consecutive frame-to-frame particle matches without employing any gap closing schemes. As a result, the recovery ratio is compounded over image frames. Moreover, the recovery ratio of the algorithm decreases with an increase in the fraction of seeded/removed particles. The algorithm maintains a very high η r ∼ 0.98 until the 5th image frame for 10% of particles seeded/removed. The recovery ratio at the 5th image frame drops to η r ∼ 0.91 for 20% of seeded/removed particles for each image frame. The mismatch ratio for all the cases is very low, η m < 1.4 × 10 −4 . Overall, these results demonstrate that T-PT can accurately track particles across time-lapse data even when a significant fraction of particles disappear and appear between image frames. These results could be further improved using gapclosing schemes.

T-PT algorithm Workflow
In this section we describe the workflow of the T-PT algorithm (Supplementary Figure 7). In particular, we explain each step in the algorithm and the effect of various parameters on the algorithm performance.

Input image pair:
T-PT loads a consecutive image pair for every time point in the time-lapse image series. Each time point image is stored within a cell container 'vol' inside a binary .mat file (the standard file extension to store MATLAB formatted data).

Particle detection and localization:
First, T-PT normalizes the input images such that a voxel has a maximum intensity value of 1. To detect particles in the images, T-PT converts the input images into binary images using image thresholding operations. Image thresholding operations are used to segment the central particle region as bright voxels and background as dark voxels. The default image thresholding value is 0.5. However, the user should set the image thresholding value such that adjacent particles are separated by at least one dark voxel.
From the binary images, the algorithm computes simply-connected regions of bright voxels. To eliminate simply-connected regions corresponding to noise or multiply connected particles, the user needs to define the minimum and maximum voxel size for the simply-connected region of the particle based on its size within the image. The default value for the minimum and the maximum number of voxels is 4 and infinity, respectively. The algorithm computes the centroid of each filtered simply-connected region to estimate a voxel-level particle center.
The radial symmetry method [7,8] computes each particle center with sub-voxel accuracy from the original image subset surrounding each estimated particle center. The user needs to set the image subset size such that each particle just fits inside the image subset. The default image subset size is 7x7x7 voxels. 3. Generating the particle descriptor: As described in the Results section, T-PT computes the particle descriptor for each particle in the reference and deformed image. By default, the algorithm uses 16 nearest neighboring particles, n, and two concentric shells, k, to create the particle descriptor. The reported default values of n and k have been found empirically to work well for a variety of motion fields and particle seeding densities in 3D images. While changing values of n and k do not significantly affect the algorithm's performance, the users can optimize these parameters for the motion fields within their own images.
For small n, the particle descriptor bins get sparsely filled and thus reduce the uniqueness of particle descriptor. For large n, each particle descriptor bin gets filled with many particles. Since the particles are randomly distributed, each bin in the particle descriptor gets a similar number of particles, and thus reduces the uniqueness of the particle descriptor.
Tests on simulated motion fields showed that k = 1 only allowed for eight bins in the particle descriptor. The length of this particle descriptor is too short in most cases to create a unique particle descriptor signature for a large number of particles. For large k, the particle descriptor bin size reduces to a small volume, which does not allow particles enough room to move without significantly altering the particle descriptor.

Particle linking:
The reference and deformed volumes are divided into cubic volume subsets for particle linking. The particles from the reference image subset are linked to the corresponding deformed image subset. To allow particles near the edge of the reference image subset to undergo large displacement while still being in the same corresponding deformed image subset, the deformed image subset size is always set to be twice the size of the reference image subset. However, the reference and deformed image have the same centers. The initial subset size should be larger than the maximum particle motion in the images. The default value of the initial subset size is 256 voxels. In the iterative particle matching process, the size of the image subset remains the same or decreases by a factor of two, until the convergence criteria are met.
The particles are linked between the reference image subsets and the corresponding deformed image subsets such that a chosen particle pair has the smallest L 2 particle descriptor distance. Assuming, r and s to be a particle descriptor for a particle in the reference and deformed images respectively, the L 2 particle descriptor distance between r and s is defined as, Here, r and s are expressed as a vector having 8k (number of bins) components. Then, the algorithm eliminates any particle links which are not bijective to enforce a one-to-one particle linking condition.

Particle link verification:
T-PT uses the similarity of neighborhood test and the universal median filter test as outlier removal schemes to remove spurious particle links. For the similarity of neighborhood test, the default values of p = 2 and q = 5 are empirically found to work well for a variety of motion fields and particle seeding densities in our 3D images. A small p/q ratio relaxes the outlier removal criterion for the similarity of neighborhood test, which is suited for high spatial frequency motion as the local relative particle positions change vastly. A high p/q ratio is appropriate for low spatial frequency motion as the the local relative particle position does not change significantly.
In the universal median test [9], 27 neighbors are used to compute the residual from the displacement for each particle link. The residual is normalized by the median residual of these 27 neighboring particles. Any particle link exceeding the normalized residual value above a user-defined threshold value is removed. By default, the residual threshold value is 4. Increasing the threshold value relaxes the outlier removal criteria and is suited for very high spatial frequency motion fields while a low threshold value is appropriate for a low spatial frequency motion field.
Overall, relaxing the outlier removal criteria increases the number of particle matches found by the algorithm, but it can also lead to a higher number of spurious particle matches. Whereas setting stricter outlier removal criteria will decrease the number of particle matches found by the algorithm, but it also reduces the number of false particle matches. The user can optimize the extent of the outlier removal process to level appropriate to their application and motion field.

Iterative deformation warping:
Iterative deformation warping is used to warp particle positions in reference and deformed image volumes based on cumulative particle matches. IDM helps in the recovery of large particle displacements, and improves the similarity of the particle descriptor in the reference and deformed volumes, which further enhances the accuracy of the particle linking process. Assume x k i and y k i to be the particle positions in the reference and deformed images in the k th iteration of deformation warping. From matched particle pairs, the displacement field at matched particle positions in the reference image is calculated as u m i = y m i − x m i . The displacement field is linearly interpolated to all particle positions x i and y i as u k i for the k th iteration of deformation warping. The particle positions in the reference and deformed frame are updated in the k th + 1 iteration step as, 7. Convergence criteria: In each iteration step of the iterative particle matching and IDM process, T-PT computes the total fraction of particles matched between the image pair. After each iteration, the fraction of particle matches increases. The algorithm meets the convergence criteria when the fraction of matched particles after each iteration reaches a steady value. Specifically, the convergence criteria is reached if the current subset size is 16 voxels and if any of the following conditions are met: • The algorithm has run through a maximum of 5 iterations for a subset size of 16 voxels.
• The fraction of matched particles has improved from the previous iteration by an amount smaller than 0.25%.
• The ratio of the fraction of matched particles between the current and previous iterations is less than 0.1.
If the convergence criteria are not met, the algorithm undergoes another iteration of particle matching and IDM. In the next iteration of particle matching, the subset size reduces by a factor of 2 if any of the following conditions are achieved. Otherwise, the subset size for the next iteration remains the same.
• The algorithm has run through a maximum of 5 iterations for the same current subset size.
• The fraction of matched particles has improved from the previous iteration by an amount smaller than 4%.
• The ratio of the fraction of matched particles between the current iteration and the previous iteration is less than 0.6.
The minimum subset size is limited to 16 voxels. If the algorithm tries to set the subset size below 16 voxels, it is reset back to 16 voxels, unless the user specifies a smaller subset size.
8. Particle search in predicted regions: As a final step to improve the recovery ratio of the algorithm, displacements from the matched particles are interpolated to predict the position of unmatched particles in the deformed volume with a default search radius of one-third of the nearest neighbor particle separation distance to the unmatched particle in the reference image. If the algorithm finds a single particle in the predicted region in the deformed image, it saves the particle pair as a temporary particle link. 9. Particle link verification: The temporary particle links formed in the previous step are subjected to the earlier described outlier removal scheme comprised of the similarity of neighborhood test and the universal median test [9]. The links verified by the outlier removal schemes are stored in the list of successful particle matches between the image pair.
10. Output: The particle centers from the reference and deformed images and the particle matches are saved as the output of the algorithm.