3D cell tracking with Ultrack¶
This document shows how Ultrack can track cells in large fluorescence microscopy data. We use the zebrafish tail data from zebrahub acquired with DaXi.
Here, we use the image processing functions provided with Ultrack, but any other segmentation method could be used. For example, refer to the stardist_2 example to see how to use Ultrack with segmentation labels.
Note: This is not the same methodology used in the zebrahub paper. There, we used a convolutional neural network to predict the nuclei boundaries and the initial version of Ultrack software.
First, we import the required packages. You can install them using the conda environment file on this folder.
[1]:
import napari
import numpy as np
import dask.array as da
from napari.utils.notebook_display import nbscreenshot
from rich.pretty import pprint
from ultrack import MainConfig, load_config, track, to_tracks_layer, tracks_to_zarr
from ultrack.imgproc import robust_invert, detect_foreground
from ultrack.utils.array import array_apply, create_zarr
We open the data remotely using the napari-ome-zarr plugin.
[2]:
viewer = napari.Viewer()
viewer.window.resize(1800, 1000)
im_layer, = viewer.open(
"http://public.czbiohub.org/royerlab/zebrahub/imaging/single-objective/ZSNS001_tail.ome.zarr/",
plugin="napari-ome-zarr",
rendering="attenuated_mip",
gamma=0.7,
contrast_limits=(0, 500),
)
nbscreenshot(viewer)
[2]:
We setup to process only a subset of frames from the dataset, feel free to change the values below.
[3]:
voxel_size = im_layer.scale[1:] # (z, y, x) pixel size
start_idx = 400 # starting frame
viewer.dims.set_point(0, start_idx + 5)
image = im_layer.data[0]
image = image[start_idx:(start_idx + 10)] # processing only a subset of time points
image.shape
[3]:
(10, 420, 1217, 1091)
We detect which pixels contains cells and which are background using the array_apply function to apply the detect_foreground to each individual time point while saving it to the detection.zarr path.
[4]:
detection = create_zarr(image.shape, bool, store_or_path="detection.zarr", overwrite=True)
array_apply(
image,
out_array=detection,
func=detect_foreground,
sigma=25.0,
voxel_size=voxel_size,
)
viewer.add_labels(
detection,
visible=True,
translate=(start_idx, 0, 0, 0),
scale=voxel_size,
).contour = 2
nbscreenshot(viewer)
Applying detect_foreground ...: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [01:38<00:00, 9.90s/it]
[4]:
Ultrack requires an estimate of the cell boundaries, we approximate this by inverting the image using the robust_invert function and saving it to boundaries.zarr.
[5]:
boundaries = create_zarr(image.shape, np.float16, store_or_path="boundaries.zarr", overwrite=True)
array_apply(
image,
out_array=boundaries,
func=robust_invert,
voxel_size=voxel_size,
)
viewer.add_image(boundaries, visible=False, translate=(start_idx, 0, 0, 0), scale=voxel_size)
Applying robust_invert ...: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [00:36<00:00, 3.61s/it]
[5]:
<Image layer 'boundaries' at 0x7f3f8bdf54e0>
Now we start the tracking stage. Ultrack parameters are provided in a configuration file which is equivalent to the MainConfig object.
Using a configuration file helps reproducing previous runs and it can be easily shared when using distributed processing.
In this example the configuration is already configured so we only need to load it. We change the segmentation threshold parameter as an example.
The configuration documentation can be found here.
[6]:
cfg = load_config("config.toml")
# cfg = MainConfig() # or load default config
cfg.segmentation_config.threshold = 0.5
pprint(cfg)
MainConfig( │ data_config=DataConfig(working_dir=PosixPath('.'), database='sqlite', address=None, n_workers=1), │ segmentation_config=SegmentationConfig( │ │ threshold=0.5, │ │ min_area=500, │ │ max_area=10000, │ │ min_frontier=0.0, │ │ anisotropy_penalization=0.0, │ │ max_noise=0.0, │ │ ws_hierarchy=<function watershed_hierarchy_by_area at 0x7f404268c670>, │ │ n_workers=1 │ ), │ linking_config=LinkingConfig( │ │ n_workers=1, │ │ max_neighbors=5, │ │ max_distance=5.0, │ │ distance_weight=0.0, │ │ z_score_threshold=5.0 │ ), │ tracking_config=TrackingConfig( │ │ appear_weight=-0.002, │ │ disappear_weight=-0.01, │ │ division_weight=-0.001, │ │ dismiss_weight_guess=None, │ │ include_weight_guess=None, │ │ window_size=100, │ │ overlap_size=5, │ │ solution_gap=0.0, │ │ time_limit=120000, │ │ method=0, │ │ n_threads=0, │ │ link_function='power', │ │ power=4.0, │ │ bias=0.0 │ ) )
Once the configuration is set, it can be used with the cell detection and boundaries to compute the tracks.
[7]:
track(
cfg,
detection=detection,
edges=boundaries,
scale=voxel_size,
overwrite=True,
)
Adding nodes to database: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [14:15<00:00, 85.53s/it]
Linking nodes.: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 9/9 [00:21<00:00, 2.34s/it]
Academic license - for non-commercial use only - expires 2024-08-06
Using GRB solver
Solving ILP batch 0
Constructing ILP ...
Solving ILP ...
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 694321 rows, 901725 columns and 2153390 nonzeros
Model fingerprint: 0xc491f3a0
Variable types: 0 continuous, 901725 integer (901725 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-16, 7e-01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 348598 rows and 358070 columns (presolve time = 8s) ...
Presolve removed 348598 rows and 358070 columns
Presolve time: 7.91s
Presolved: 345723 rows, 543655 columns, 1274833 nonzeros
Found heuristic solution: objective 5043.3524392
Variable types: 0 continuous, 543655 integer (543655 binary)
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Ordering time: 0.27s
Barrier statistics:
AA' NZ : 1.181e+06
Factor NZ : 5.885e+06 (roughly 400 MBytes of memory)
Factor Ops : 1.648e+08 (less than 1 second per iteration)
Threads : 6
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 4.21353902e+04 3.99028854e+05 2.47e+01 3.77e-01 1.78e+00 10s
1 2.12386917e+04 2.47108458e+05 1.03e+01 1.25e-15 6.30e-01 10s
2 7.87747235e+03 8.70116483e+04 6.09e-01 8.97e-16 7.84e-02 10s
3 8.08051876e+03 3.09158317e+04 5.61e-02 6.94e-16 1.79e-02 11s
4 9.84651687e+03 2.00505628e+04 5.16e-03 4.72e-16 7.71e-03 11s
5 1.08136667e+04 1.62059692e+04 2.42e-03 3.75e-16 4.07e-03 11s
6 1.16018041e+04 1.45098500e+04 1.01e-03 5.13e-16 2.19e-03 11s
7 1.20816818e+04 1.36345628e+04 4.26e-04 4.44e-16 1.17e-03 11s
8 1.24176057e+04 1.31090427e+04 1.26e-04 3.89e-16 5.20e-04 12s
9 1.25344957e+04 1.29002122e+04 5.41e-05 3.61e-16 2.75e-04 12s
10 1.25908505e+04 1.27840755e+04 2.44e-05 4.16e-16 1.45e-04 12s
11 1.26307457e+04 1.27115438e+04 6.41e-06 4.55e-16 6.07e-05 12s
12 1.26401826e+04 1.26726147e+04 2.85e-06 5.27e-16 2.44e-05 13s
13 1.26445408e+04 1.26610521e+04 1.37e-06 4.72e-16 1.24e-05 13s
14 1.26472290e+04 1.26553951e+04 4.52e-07 4.37e-16 6.14e-06 13s
15 1.26487004e+04 1.26520030e+04 7.88e-08 3.89e-16 2.48e-06 13s
Barrier performed 15 iterations in 13.02 seconds
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root relaxation: objective 1.264910e+04, 131331 iterations, 4.67 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 12649.0980 0 27003 5043.35244 12649.0980 151% - 15s
H 0 0 12459.677634 12649.0980 1.52% - 16s
H 0 0 12476.150138 12649.0980 1.39% - 17s
0 0 12580.7032 0 5199 12476.1501 12580.7032 0.84% - 21s
H 0 0 12535.283769 12580.7032 0.36% - 22s
H 0 0 12553.596884 12580.7032 0.22% - 22s
H 0 0 12562.441568 12580.7032 0.15% - 23s
H 0 0 12574.963919 12580.7032 0.05% - 24s
0 0 12580.7016 0 5198 12574.9639 12580.7016 0.05% - 25s
0 0 12575.9946 0 880 12574.9639 12575.9946 0.01% - 26s
H 0 0 12574.984242 12575.9946 0.01% - 27s
0 0 12575.9817 0 856 12574.9842 12575.9817 0.01% - 29s
0 0 12575.6498 0 126 12574.9842 12575.6498 0.01% - 32s
H 0 0 12575.620168 12575.6498 0.00% - 34s
0 0 12575.6412 0 104 12575.6202 12575.6412 0.00% - 35s
* 0 0 0 12575.635833 12575.6358 0.00% - 37s
Cutting planes:
Gomory: 142
Cover: 45
Clique: 2203
MIR: 38
Zero half: 4202
RLT: 5058
Explored 1 nodes (139255 simplex iterations) in 37.43 seconds
Thread count was 16 (of 16 available processors)
Solution count 10: 12575.6 12575.6 12575 ... 5043.35
Optimal solution found (tolerance 0.00e+00)
Best objective 1.257563583284e+04, best bound 1.257563583284e+04, gap 0.0000%
Saving solution ...
Done!
Finally, we export the resulting tracks into the napari tracks format and a zarr with the instance segmentation masks.
[8]:
tracks_df, graph = to_tracks_layer(cfg)
tracks_df.to_csv("tracks.csv", index=False)
segments = tracks_to_zarr(
cfg,
tracks_df,
store_or_path="segments.zarr",
overwrite=True,
)
viewer.add_tracks(
tracks_df[["track_id", "t", "z", "y", "x"]],
name="tracks",
graph=graph,
scale=voxel_size,
translate=(start_idx, 0, 0, 0),
visible=False,
)
viewer.add_labels(
da.from_zarr(segments),
name="segments",
scale=voxel_size,
translate=(start_idx, 0, 0, 0),
).contour = 2
nbscreenshot(viewer)
Exporting segmentation masks: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [00:19<00:00, 1.91s/it]
[8]:
