Tracking highly motile cells using flow field estimation and Ultrack¶
Tracking highly motile cells can be challenging. Here, we show how you can estimate the cells’ flow to assist in tracking motile cells.
For this purpose, we analyze the Tribolium Castaneum embryo 3D cartographic projection from the cell tracking challenge (CTC); the cells in this dataset show a clear migration pattern, plus their movement are further amplified on the left and right edges of the image due to projection artifacts.
First, we download the data from the CTC website.
[1]:
!wget -nc http://data.celltrackingchallenge.net/training-datasets/Fluo-N3DL-TRIC.zip
!unzip -n Fluo-N3DL-TRIC.zip
File ‘Fluo-N3DL-TRIC.zip’ already there; not retrieving.
Archive: Fluo-N3DL-TRIC.zip
We import the required packages. You can install them using the conda environment file in this folder.
[2]:
from pathlib import Path
import napari
import numpy as np
import pandas as pd
import dask.array as da
import scipy.ndimage as ndi
from napari.utils.notebook_display import nbscreenshot
from rich.pretty import pprint
from tifffile import imread
from ultrack import MainConfig, add_flow, segment, link, solve, to_tracks_layer, tracks_to_zarr, to_ctc
from ultrack.utils.array import array_apply, create_zarr
from ultrack.imgproc import robust_invert, detect_foreground
from ultrack.imgproc.flow import timelapse_flow, advenct_from_quasi_random, trajectories_to_tracks
from ultrack.utils.cuda import on_gpu
from ultrack.tracks.stats import tracks_df_movement
We set the download dataset path, open a napari viewer, and load the images using napari’s default reader.
If preferred, you can restrict the analysis to the last 50 frames where the interesting part is. However, you won’t be able to compute the accuracy since the ground-truth dimension won’t match.
[3]:
dataset_path = Path("Fluo-N3DL-TRIC/02")
# im_files = sorted(dataset_path.glob("*.tif"))[-50:] # uncomment to run the last 50 frames
im_files = sorted(dataset_path.glob("*.tif")) # uncomment to run the whole dataset
viewer = napari.Viewer()
viewer.window.resize(1800, 1000)
im_layer = viewer.open(im_files, stack=True)
image = viewer.layers[0].data
nbscreenshot(viewer)
[3]:
We detect which voxels 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, "detection.zarr", overwrite=True)
array_apply(
image,
out_array=detection,
func=on_gpu(detect_foreground),
)
viewer.add_labels(detection, visible=True).contour = 2
nbscreenshot(viewer)
/hpc/mydata/ilan.silva/miniforge3/envs/ultrack-flow-field/lib/python3.11/site-packages/ultrack/utils/array.py:190: FutureWarning: The NestedDirectoryStore is deprecated and will be removed in a Zarr-Python version 3, see https://github.com/zarr-developers/zarr-python/issues/1274 for more information.
store = zarr.NestedDirectoryStore(str(store_or_path))
Applying detect_foreground ...: 100%|███████████████████████████| 210/210 [04:05<00:00, 1.17s/it]
[4]:
Ultrack also 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, "boundaries.zarr", overwrite=True)
array_apply(
image,
out_array=boundaries,
func=on_gpu(robust_invert),
sigma=3.0,
)
viewer.add_image(boundaries, visible=False)
/hpc/mydata/ilan.silva/miniforge3/envs/ultrack-flow-field/lib/python3.11/site-packages/ultrack/utils/array.py:190: FutureWarning: The NestedDirectoryStore is deprecated and will be removed in a Zarr-Python version 3, see https://github.com/zarr-developers/zarr-python/issues/1274 for more information.
store = zarr.NestedDirectoryStore(str(store_or_path))
Applying robust_invert ...: 100%|███████████████████████████████| 210/210 [03:55<00:00, 1.12s/it]
[5]:
<Image layer 'boundaries' at 0x14e659b06650>
Using the original images we compute the movement flow. We use a lower number of n_scales because the z dimension is very short and using the default n_scales=3 results in a zero-length dimension.
[6]:
!rm -r flow.zarr # removing previous flow
flow = timelapse_flow(image, store_or_path="flow.zarr", n_scales=2, lr=1e-2, num_iterations=2_000)
viewer.add_image(
flow,
contrast_limits=(-0.001, 0.001),
colormap="turbo",
visible=False,
scale=(4,) * 3,
channel_axis=1,
name="flow field",
)
/hpc/mydata/ilan.silva/miniforge3/envs/ultrack-flow-field/lib/python3.11/site-packages/ultrack/utils/array.py:190: FutureWarning: The NestedDirectoryStore is deprecated and will be removed in a Zarr-Python version 3, see https://github.com/zarr-developers/zarr-python/issues/1274 for more information.
store = zarr.NestedDirectoryStore(str(store_or_path))
Computing flow: 0%| | 0/209 [00:00<?, ?it/s]/hpc/mydata/ilan.silva/miniforge3/envs/ultrack-flow-field/lib/python3.11/site-packages/torch/nn/functional.py:4436: UserWarning: Since version 1.3.0, affine_grid behavior has changed for unit-size grids when align_corners=True. This is not an intended use case of affine_grid. See the documentation of affine_grid for details.
warnings.warn(
Computing flow: 100%|███████████████████████████████████████████| 209/209 [11:34<00:00, 3.32s/it]
[6]:
[<Image layer 'flow field' at 0x14e878241ad0>,
<Image layer 'flow field [1]' at 0x14e87811bcd0>,
<Image layer 'flow field [2]' at 0x14e873a1c190>]
We advect random samples through the flow field and visualize them using the napari tracks layer to evaluate the flow. Plus, we color the tracks according to their angle on the xy-plane.
[7]:
trajectory = advenct_from_quasi_random(flow, detection.shape[-3:], n_samples=1000)
flow_tracklets = pd.DataFrame(
trajectories_to_tracks(trajectory),
columns=["track_id", "t", "z", "y", "x"],
)
flow_tracklets[["z", "y", "x"]] += 0.5 # napari was crashing otherwise, might be an openGL issue
flow_tracklets[["dz", "dy", "dx"]] = tracks_df_movement(flow_tracklets)
flow_tracklets["angles"] = np.arctan2(flow_tracklets["dy"], flow_tracklets["dx"])
flow_tracklets.to_csv("flow_tracklets.csv", index=False)
viewer.add_tracks(
flow_tracklets[["track_id", "t", "z", "y", "x"]],
name="flow vectors",
visible=True,
tail_length=25,
features=flow_tracklets[["angles", "dy", "dx"]],
colormap="hsv",
).color_by="angles"
nbscreenshot(viewer)
[7]:
Now that we have our detection, boundaries and flow, and have checked that they look ok, we will start the tracking step. Ultrack’s tracking with flow relies on 4 steps and their respective python functions: - segment: Computes the segmentation hypothesis for tracking; - add_flow: Adds the flow to each segmentation hypothesis, must be called before link; - link: Links and assign edge weights to the segmentation hypothesis, taking into consideration the previously
added flow; - solve: Solves the tracking problem by selecting the strongly associated segmentation hypothesis.
These steps use our configuration object, MainConfig, which we’ll set up below. Its documentation can be found here.
The parameters were chosen manually by inspection.
NOTE: If you’re running out of memory, you should decrease the n_workers parameters. If you want to speed up the processing and have spare memory, you can increase it.
We will evaluate the tracking accuracy in two scenarios, with and without flow.
[8]:
cfg = MainConfig()
cfg.data_config.n_workers = 8
cfg.segmentation_config.n_workers = 8
cfg.segmentation_config.min_area = 250
cfg.segmentation_config.max_area = 15_000
cfg.linking_config.n_workers = 12
cfg.linking_config.max_neighbors = 5
cfg.linking_config.max_distance = 50.0
cfg.linking_config.distance_weight = 0.0001
cfg.tracking_config.window_size = 70
cfg.tracking_config.overlap_size = 5
cfg.tracking_config.appear_weight = -0.01
cfg.tracking_config.disappear_weight = -0.001
cfg.tracking_config.division_weight = 0
pprint(cfg)
MainConfig( │ data_config=DataConfig(working_dir=PosixPath('.'), database='sqlite', address=None, n_workers=8), │ segmentation_config=SegmentationConfig( │ │ threshold=0.5, │ │ min_area=250, │ │ max_area=15000, │ │ min_frontier=0.0, │ │ anisotropy_penalization=0.0, │ │ max_noise=0.0, │ │ ws_hierarchy=<function watershed_hierarchy_by_area at 0x14e8dc2fb060>, │ │ n_workers=8 │ ), │ linking_config=LinkingConfig( │ │ n_workers=12, │ │ max_neighbors=5, │ │ max_distance=50.0, │ │ distance_weight=0.0001, │ │ z_score_threshold=5.0 │ ), │ tracking_config=TrackingConfig( │ │ appear_weight=-0.01, │ │ disappear_weight=-0.001, │ │ division_weight=0, │ │ dismiss_weight_guess=None, │ │ include_weight_guess=None, │ │ window_size=70, │ │ overlap_size=5, │ │ solution_gap=0.001, │ │ time_limit=36000, │ │ method=0, │ │ n_threads=-1, │ │ link_function='power', │ │ power=4, │ │ bias=-0.0 │ ) )
Tracking without flow¶
First, we evaluate the tracking without the flow information.
Independently of using the flow or not, we compute the set of candidate segmentation hypotheses from the detection and boundaries maps.
[9]:
segment(detection, boundaries, cfg, overwrite=True)
Adding nodes to database: 100%|█████████████████████████████████| 210/210 [07:24<00:00, 2.12s/it]
Next, we compute the link (association) between frames and solve the tracking problem.
[10]:
link(cfg, overwrite=True)
solve(cfg)
Linking nodes.: 100%|███████████████████████████████████████████| 209/209 [05:02<00:00, 1.45s/it]
Set parameter TokenServer to value "license01"
Using Gurobi solver
Solving ILP batch 0
Constructing ILP ...
Set parameter TimeLimit to value 36000
Solving ILP ...
Set parameter NodeLimit to value 1073741824
Set parameter SolutionLimit to value 1073741824
Set parameter IntFeasTol to value 1e-06
Set parameter Method to value 3
Set parameter MIPGap to value 0.001
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Rocky Linux 8.9 (Green Obsidian)")
CPU model: Intel(R) Xeon(R) Platinum 8480CL, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 112 physical cores, 112 logical processors, using up to 32 threads
Optimize a model with 1746085 rows, 3989594 columns and 8324193 nonzeros
Model fingerprint: 0xe23f4ed2
Variable types: 0 continuous, 3989594 integer (3989594 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-31, 1e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 247123 rows and 462622 columns (presolve time = 5s) ...
Presolve removed 248562 rows and 466232 columns (presolve time = 10s) ...
Presolve removed 248562 rows and 466263 columns (presolve time = 15s) ...
Presolve removed 248878 rows and 466680 columns
Presolve time: 19.22s
Presolved: 1497207 rows, 3522914 columns, 7515316 nonzeros
Found heuristic solution: objective 2389.0099480
Variable types: 0 continuous, 3522914 integer (3522914 binary)
Found heuristic solution: objective 5999.1902205
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Elapsed ordering time = 6s
Elapsed ordering time = 14s
Elapsed ordering time = 15s
Elapsed ordering time = 20s
Ordering time: 22.49s
Barrier statistics:
AA' NZ : 5.157e+06
Factor NZ : 3.532e+08 (roughly 5.0 GB of memory)
Factor Ops : 1.162e+12 (roughly 1 second per iteration)
Threads : 30
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 5.31657476e+05 5.05434771e+06 4.04e+01 6.91e-01 4.17e+00 54s
1 2.20617073e+05 3.27805942e+06 1.57e+01 2.78e-15 1.44e+00 56s
2 5.32418704e+04 1.22061627e+06 1.28e+00 1.78e-15 1.97e-01 58s
3 5.14601935e+04 5.44316238e+05 1.56e-01 1.40e-15 6.46e-02 59s
Barrier performed 3 iterations in 59.50 seconds (44.27 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
991798 1.3147298e+05 0.000000e+00 0.000000e+00 60s
Root relaxation: objective 1.314730e+05, 991798 iterations, 37.85 seconds (28.71 work units)
Total elapsed time = 77.37s (DegenMoves)
Total elapsed time = 82.32s (DegenMoves)
Total elapsed time = 93.57s (DegenMoves)
Total elapsed time = 103.16s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 131472.984 0 41039 5999.19022 131472.984 2092% - 121s
H 0 0 130968.94953 131472.984 0.38% - 122s
H 0 0 131031.14354 131472.984 0.34% - 129s
H 0 0 131031.95393 131472.984 0.34% - 133s
0 0 131322.537 0 6698 131031.954 131322.537 0.22% - 142s
H 0 0 131197.38611 131322.537 0.10% - 149s
Cutting planes:
Gomory: 24
Clique: 8504
Zero half: 1375
RLT: 756
BQP: 1
PSD: 609
Explored 1 nodes (1364600 simplex iterations) in 149.55 seconds (130.91 work units)
Thread count was 32 (of 112 available processors)
Solution count 7: 131197 131032 131031 ... -0
Optimal solution found (tolerance 1.00e-03)
Best objective 1.311973861129e+05, best bound 1.313225369358e+05, gap 0.0954%
Saving solution ...
Done!
Set parameter TokenServer to value "license01"
Using Gurobi solver
Solving ILP batch 2
Constructing ILP ...
Set parameter TimeLimit to value 36000
Solving ILP ...
Set parameter NodeLimit to value 1073741824
Set parameter SolutionLimit to value 1073741824
Set parameter IntFeasTol to value 1e-06
Set parameter Method to value 3
Set parameter MIPGap to value 0.001
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Rocky Linux 8.9 (Green Obsidian)")
CPU model: Intel(R) Xeon(R) Platinum 8480CL, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 112 physical cores, 112 logical processors, using up to 32 threads
Optimize a model with 932575 rows, 1935048 columns and 4095825 nonzeros
Model fingerprint: 0x176063d1
Variable types: 0 continuous, 1935048 integer (1935048 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-32, 9e-01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 154514 rows and 243805 columns (presolve time = 5s) ...
Presolve removed 158979 rows and 250066 columns (presolve time = 10s) ...
Presolve removed 158979 rows and 250066 columns
Presolve time: 10.02s
Presolved: 773596 rows, 1684982 columns, 3626318 nonzeros
Found heuristic solution: objective 1851.7615723
Variable types: 0 continuous, 1684982 integer (1684982 binary)
Found heuristic solution: objective 3345.5559972
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Ordering time: 5.72s
Barrier statistics:
AA' NZ : 2.589e+06
Factor NZ : 5.047e+07 (roughly 1.4 GB of memory)
Factor Ops : 3.380e+10 (less than 1 second per iteration)
Threads : 30
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 2.10927373e+05 2.34924116e+06 3.44e+01 6.77e-01 3.55e+00 21s
1 8.87703784e+04 1.44012212e+06 1.29e+01 2.22e-15 1.16e+00 22s
2 2.44900519e+04 4.81334787e+05 8.52e-01 1.78e-15 1.48e-01 22s
3 2.96797528e+04 2.44047922e+05 5.58e-02 1.12e-15 5.60e-02 23s
4 4.29087850e+04 1.16620185e+05 4.39e-03 7.91e-16 1.89e-02 23s
5 4.95027846e+04 9.02778898e+04 1.64e-03 6.66e-16 1.04e-02 24s
6 5.41376718e+04 7.54593383e+04 5.84e-04 6.04e-16 5.45e-03 24s
7 5.70187427e+04 6.85596615e+04 1.99e-04 4.98e-16 2.95e-03 25s
8 5.85484914e+04 6.40872588e+04 6.14e-05 4.86e-16 1.42e-03 25s
9 5.90355155e+04 6.18934498e+04 2.89e-05 5.67e-16 7.31e-04 26s
10 5.94298225e+04 6.07408188e+04 7.14e-06 4.86e-16 3.35e-04 26s
11 5.95606501e+04 6.00229993e+04 2.15e-06 5.79e-16 1.18e-04 27s
12 5.96191797e+04 5.98430232e+04 7.67e-07 9.78e-16 5.72e-05 28s
13 5.96449321e+04 5.97633740e+04 3.46e-07 5.55e-16 3.03e-05 28s
Barrier performed 13 iterations in 28.10 seconds (32.30 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
497994 5.9674107e+04 0.000000e+00 0.000000e+00 29s
Root relaxation: objective 5.967411e+04, 497994 iterations, 17.18 seconds (15.47 work units)
Total elapsed time = 31.90s (DegenMoves)
Total elapsed time = 36.55s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 59674.1068 0 45215 3345.55600 59674.1068 1684% - 47s
H 0 0 59192.118737 59674.1068 0.81% - 47s
H 0 0 59247.786071 59674.1068 0.72% - 51s
H 0 0 59248.360872 59674.1068 0.72% - 52s
0 0 59529.7254 0 8810 59248.3609 59529.7254 0.47% - 58s
H 0 0 59409.093216 59529.7254 0.20% - 61s
H 0 0 59482.560766 59529.7254 0.08% - 61s
Cutting planes:
Gomory: 19
Clique: 7399
Zero half: 2014
RLT: 1423
PSD: 386
Explored 1 nodes (595005 simplex iterations) in 62.05 seconds (57.82 work units)
Thread count was 32 (of 112 available processors)
Solution count 8: 59482.6 59409.1 59248.4 ... -0
Optimal solution found (tolerance 1.00e-03)
Best objective 5.948256076595e+04, best bound 5.952972542605e+04, gap 0.0793%
Saving solution ...
Done!
Set parameter TokenServer to value "license01"
Using Gurobi solver
Solving ILP batch 1
Constructing ILP ...
Set parameter TimeLimit to value 36000
Solving ILP ...
Set parameter NodeLimit to value 1073741824
Set parameter SolutionLimit to value 1073741824
Set parameter IntFeasTol to value 1e-06
Set parameter Method to value 3
Set parameter MIPGap to value 0.001
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Rocky Linux 8.9 (Green Obsidian)")
CPU model: Intel(R) Xeon(R) Platinum 8480CL, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 112 physical cores, 112 logical processors, using up to 32 threads
Optimize a model with 1632247 rows, 3495711 columns and 7439607 nonzeros
Model fingerprint: 0x9e9ec629
Variable types: 0 continuous, 3495711 integer (3495711 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-40, 9e-01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -31.9460000
Presolve removed 300654 rows and 440013 columns (presolve time = 5s) ...
Presolve removed 301011 rows and 440394 columns (presolve time = 10s) ...
Presolve removed 301032 rows and 442789 columns (presolve time = 18s) ...
Presolve removed 302780 rows and 442789 columns
Presolve time: 18.59s
Presolved: 1329467 rows, 3052922 columns, 6559606 nonzeros
Found heuristic solution: objective 3006.5426638
Variable types: 0 continuous, 3052922 integer (3052922 binary)
Found heuristic solution: objective 4027.9245508
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Elapsed ordering time = 6s
Ordering time: 10.97s
Barrier statistics:
AA' NZ : 4.677e+06
Factor NZ : 2.118e+08 (roughly 3.5 GB of memory)
Factor Ops : 4.288e+11 (roughly 1 second per iteration)
Threads : 30
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 3.37107194e+05 4.22064829e+06 4.25e+01 6.70e-01 3.66e+00 40s
1 1.36247779e+05 2.60326334e+06 1.60e+01 2.22e-15 1.21e+00 41s
2 3.38292788e+04 8.70244824e+05 1.16e+00 1.76e-15 1.55e-01 42s
3 4.16258376e+04 3.30255435e+05 7.42e-02 1.33e-15 4.22e-02 44s
4 6.21952013e+04 1.77579395e+05 1.57e-02 1.05e-15 1.65e-02 45s
5 7.19083691e+04 1.30314198e+05 6.78e-03 6.31e-16 8.35e-03 46s
6 8.14985464e+04 1.09123125e+05 7.56e-04 6.94e-16 3.93e-03 47s
7 8.38834929e+04 9.81137282e+04 3.22e-04 5.76e-16 2.02e-03 48s
8 8.55530039e+04 9.11308740e+04 9.48e-05 6.11e-16 7.93e-04 50s
9 8.62445813e+04 8.81564974e+04 2.94e-05 6.78e-16 2.72e-04 51s
10 8.65331434e+04 8.74403987e+04 1.05e-05 5.76e-16 1.29e-04 53s
11 8.66519447e+04 8.72081484e+04 4.73e-06 5.55e-16 7.90e-05 54s
12 8.66986052e+04 8.70102087e+04 2.91e-06 5.55e-16 4.43e-05 55s
13 8.67487547e+04 8.68946900e+04 1.17e-06 6.11e-16 2.07e-05 57s
Barrier performed 13 iterations in 56.67 seconds (68.01 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
879608 8.6791712e+04 0.000000e+00 0.000000e+00 57s
Root relaxation: objective 8.679171e+04, 879608 iterations, 36.37 seconds (29.60 work units)
Total elapsed time = 66.61s (DegenMoves)
Total elapsed time = 70.62s (DegenMoves)
Total elapsed time = 78.39s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 86791.7125 0 65853 4027.92455 86791.7125 2055% - 93s
H 0 0 86277.199874 86791.7125 0.60% - 94s
H 0 0 86349.501117 86791.7125 0.51% - 101s
H 0 0 86350.110141 86791.7125 0.51% - 104s
0 0 86631.9526 0 12921 86350.1101 86631.9526 0.33% - 116s
H 0 0 86508.427651 86631.9526 0.14% - 121s
H 0 0 86580.436126 86631.9526 0.06% - 122s
Cutting planes:
Gomory: 21
Clique: 9451
Zero half: 2951
RLT: 2205
PSD: 438
Explored 1 nodes (1075388 simplex iterations) in 122.53 seconds (110.55 work units)
Thread count was 32 (of 112 available processors)
Solution count 8: 86580.4 86508.4 86350.1 ... -31.946
Optimal solution found (tolerance 1.00e-03)
Best objective 8.658043612555e+04, best bound 8.663195257059e+04, gap 0.0595%
Saving solution ...
Done!
The cell tracking challenge evaluates only a subset of cells or this dataset t, so we supply the first_frame parameter to the to_ctc function.
[11]:
gt_path = dataset_path.parent / "02_GT"
reference_frame = imread(gt_path / "TRA/man_track000.tif")
no_flow_res_path = Path("02_NO_FLOW_RES/TRA")
to_ctc(no_flow_res_path, cfg, first_frame=reference_frame, overwrite=True)
Selecting tracks from first trame: 100%|█████████████████████| 196/196 [00:00<00:00, 85163.53it/s]
Exporting segmentation masks: 100%|█████████████████████████████| 210/210 [00:49<00:00, 4.22it/s]
Tracking with flow¶
The flow step should be executed after the segment function. Since it has already been computed, we don’t need to start from scratch; we can apply the flow and recompute from the link step with overwrite=True.
[12]:
add_flow(cfg, flow)
100%|███████████████████████████████████████████████████████████| 210/210 [09:41<00:00, 2.77s/it]
As with the previous run, we solve the tracking and export using the first_frame parameter, but to a different directory.
[13]:
link(cfg, overwrite=True)
solve(cfg)
flow_res_path = Path("02_FLOW_RES/TRA")
to_ctc(flow_res_path, cfg, first_frame=reference_frame, overwrite=True)
Linking nodes.: 100%|███████████████████████████████████████████| 209/209 [04:57<00:00, 1.42s/it]
Set parameter TokenServer to value "license01"
Using Gurobi solver
Solving ILP batch 0
Constructing ILP ...
Set parameter TimeLimit to value 36000
Solving ILP ...
Set parameter NodeLimit to value 1073741824
Set parameter SolutionLimit to value 1073741824
Set parameter IntFeasTol to value 1e-06
Set parameter Method to value 3
Set parameter MIPGap to value 0.001
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Rocky Linux 8.9 (Green Obsidian)")
CPU model: Intel(R) Xeon(R) Platinum 8480CL, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 112 physical cores, 112 logical processors, using up to 32 threads
Optimize a model with 1746085 rows, 3990144 columns and 8325293 nonzeros
Model fingerprint: 0xd1ed2924
Variable types: 0 continuous, 3990144 integer (3990144 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e-30, 9e-01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 249114 rows and 464614 columns (presolve time = 5s) ...
Presolve removed 250805 rows and 468772 columns (presolve time = 10s) ...
Presolve removed 250805 rows and 468801 columns (presolve time = 15s) ...
Presolve removed 251073 rows and 469151 columns
Presolve time: 18.63s
Presolved: 1495012 rows, 3520993 columns, 7510847 nonzeros
Found heuristic solution: objective 1221.0507973
Variable types: 0 continuous, 3520993 integer (3520993 binary)
Found heuristic solution: objective 3752.4989648
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Elapsed ordering time = 8s
Elapsed ordering time = 10s
Elapsed ordering time = 14s
Elapsed ordering time = 15s
Elapsed ordering time = 20s
Ordering time: 22.68s
Barrier statistics:
AA' NZ : 5.153e+06
Factor NZ : 3.772e+08 (roughly 5.0 GB of memory)
Factor Ops : 1.327e+12 (roughly 2 seconds per iteration)
Threads : 30
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 4.09053039e+05 5.11412500e+06 3.87e+01 7.02e-01 3.88e+00 54s
1 1.67717044e+05 3.18387371e+06 1.48e+01 2.33e-15 1.30e+00 55s
2 4.06976318e+04 1.12030426e+06 1.20e+00 1.77e-15 1.78e-01 57s
3 4.22705882e+04 5.33991450e+05 9.25e-02 1.17e-15 6.28e-02 59s
Barrier performed 3 iterations in 59.43 seconds (47.47 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
1019045 1.0772673e+05 0.000000e+00 0.000000e+00 60s
Root relaxation: objective 1.077267e+05, 1019045 iterations, 38.68 seconds (29.32 work units)
Total elapsed time = 76.92s (DegenMoves)
Total elapsed time = 83.00s (DegenMoves)
Total elapsed time = 96.04s (DegenMoves)
Total elapsed time = 106.92s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 107726.725 0 59212 3752.49896 107726.725 2771% - 125s
H 0 0 107202.76717 107726.725 0.49% - 126s
H 0 0 107259.30859 107726.725 0.44% - 133s
H 0 0 107260.08201 107726.725 0.44% - 137s
0 0 107580.715 0 9950 107260.082 107580.715 0.30% - 149s
H 0 0 107452.43109 107580.715 0.12% - 155s
H 0 0 107530.66414 107580.715 0.05% - 156s
Cutting planes:
Gomory: 21
Clique: 8557
Zero half: 2899
RLT: 1874
BQP: 3
PSD: 582
Explored 1 nodes (1387096 simplex iterations) in 156.52 seconds (133.98 work units)
Thread count was 32 (of 112 available processors)
Solution count 8: 107531 107452 107260 ... -0
Optimal solution found (tolerance 1.00e-03)
Best objective 1.075306641382e+05, best bound 1.075807148838e+05, gap 0.0465%
Saving solution ...
Done!
Set parameter TokenServer to value "license01"
Using Gurobi solver
Solving ILP batch 2
Constructing ILP ...
Set parameter TimeLimit to value 36000
Solving ILP ...
Set parameter NodeLimit to value 1073741824
Set parameter SolutionLimit to value 1073741824
Set parameter IntFeasTol to value 1e-06
Set parameter Method to value 3
Set parameter MIPGap to value 0.001
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Rocky Linux 8.9 (Green Obsidian)")
CPU model: Intel(R) Xeon(R) Platinum 8480CL, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 112 physical cores, 112 logical processors, using up to 32 threads
Optimize a model with 932575 rows, 1934341 columns and 4094411 nonzeros
Model fingerprint: 0xbd6fa3c2
Variable types: 0 continuous, 1934341 integer (1934341 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [9e-31, 9e-01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -0.0000000
Presolve removed 156876 rows and 246320 columns (presolve time = 5s) ...
Presolve removed 159157 rows and 252929 columns (presolve time = 10s) ...
Presolve removed 161565 rows and 252929 columns
Presolve time: 10.42s
Presolved: 771010 rows, 1681412 columns, 3618374 nonzeros
Found heuristic solution: objective 759.8521681
Variable types: 0 continuous, 1681412 integer (1681412 binary)
Found heuristic solution: objective 2210.0673266
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Ordering time: 5.90s
Barrier statistics:
AA' NZ : 2.582e+06
Factor NZ : 5.335e+07 (roughly 1.4 GB of memory)
Factor Ops : 3.997e+10 (less than 1 second per iteration)
Threads : 30
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 2.13938719e+05 2.28938875e+06 4.61e+01 6.64e-01 4.46e+00 22s
1 8.94406369e+04 1.58874774e+06 1.86e+01 2.11e-15 1.62e+00 22s
2 2.10474703e+04 6.48039701e+05 1.83e+00 1.99e-15 2.37e-01 23s
3 1.80920419e+04 2.16715288e+05 1.20e-01 1.25e-15 5.33e-02 23s
4 2.89398907e+04 9.66911263e+04 1.87e-02 8.97e-16 1.75e-02 24s
5 3.39301768e+04 6.91445769e+04 7.39e-03 6.52e-16 9.08e-03 24s
6 3.80313176e+04 5.81616533e+04 2.83e-03 5.41e-16 5.18e-03 25s
7 4.04854600e+04 5.14013140e+04 1.16e-03 4.72e-16 2.80e-03 25s
8 4.23318631e+04 4.69232994e+04 3.03e-04 5.74e-16 1.18e-03 26s
9 4.28489795e+04 4.54276442e+04 1.52e-04 4.58e-16 6.61e-04 26s
10 4.32878775e+04 4.45269127e+04 4.21e-05 4.82e-16 3.18e-04 27s
11 4.34648088e+04 4.38909547e+04 9.95e-06 6.11e-16 1.09e-04 28s
Barrier performed 11 iterations in 27.75 seconds (31.22 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
501448 4.3573138e+04 0.000000e+00 0.000000e+00 28s
Root relaxation: objective 4.357314e+04, 501448 iterations, 16.32 seconds (13.57 work units)
Total elapsed time = 31.55s (DegenMoves)
Total elapsed time = 36.39s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 43573.1377 0 50605 2210.06733 43573.1377 1872% - 47s
H 0 0 43147.122423 43573.1377 0.99% - 47s
H 0 0 43209.968678 43573.1377 0.84% - 50s
H 0 0 43210.549853 43573.1377 0.84% - 52s
0 0 43449.3762 0 8302 43210.5499 43449.3762 0.55% - 58s
H 0 0 43360.629407 43449.3762 0.20% - 61s
H 0 0 43402.712704 43449.3762 0.11% - 61s
H 0 0 43404.433168 43449.3762 0.10% - 62s
H 0 0 43409.647127 43449.3762 0.09% - 62s
Cutting planes:
Gomory: 26
Clique: 7915
Zero half: 2764
RLT: 1593
BQP: 3
PSD: 411
Explored 1 nodes (596457 simplex iterations) in 63.00 seconds (57.34 work units)
Thread count was 32 (of 112 available processors)
Solution count 10: 43409.6 43404.4 43402.7 ... -0
Optimal solution found (tolerance 1.00e-03)
Best objective 4.340964712700e+04, best bound 4.344937616001e+04, gap 0.0915%
Saving solution ...
Done!
Set parameter TokenServer to value "license01"
Using Gurobi solver
Solving ILP batch 1
Constructing ILP ...
Set parameter TimeLimit to value 36000
Solving ILP ...
Set parameter NodeLimit to value 1073741824
Set parameter SolutionLimit to value 1073741824
Set parameter IntFeasTol to value 1e-06
Set parameter Method to value 3
Set parameter MIPGap to value 0.001
Gurobi Optimizer version 11.0.2 build v11.0.2rc0 (linux64 - "Rocky Linux 8.9 (Green Obsidian)")
CPU model: Intel(R) Xeon(R) Platinum 8480CL, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 112 physical cores, 112 logical processors, using up to 32 threads
Optimize a model with 1632240 rows, 3496624 columns and 7441426 nonzeros
Model fingerprint: 0xb0820919
Variable types: 0 continuous, 3496624 integer (3496624 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [2e-29, 9e-01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+00]
Found heuristic solution: objective -31.9300000
Presolve removed 305528 rows and 445262 columns (presolve time = 5s) ...
Presolve removed 306838 rows and 446601 columns (presolve time = 10s) ...
Presolve removed 306860 rows and 449107 columns (presolve time = 18s) ...
Presolve removed 308688 rows and 449107 columns
Presolve time: 18.58s
Presolved: 1323552 rows, 3047517 columns, 6546969 nonzeros
Found heuristic solution: objective 986.0253768
Variable types: 0 continuous, 3047517 integer (3047516 binary)
Found heuristic solution: objective 1881.0688555
Concurrent LP optimizer: primal simplex, dual simplex, and barrier
Showing barrier log only...
Root barrier log...
Elapsed ordering time = 5s
Ordering time: 11.14s
Barrier statistics:
AA' NZ : 4.665e+06
Factor NZ : 2.170e+08 (roughly 3.5 GB of memory)
Factor Ops : 4.463e+11 (roughly 1 second per iteration)
Threads : 30
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 3.34766804e+05 4.05877123e+06 4.80e+01 6.47e-01 3.76e+00 40s
1 1.45304635e+05 2.60108649e+06 1.98e+01 2.88e-15 1.35e+00 41s
2 3.55900369e+04 9.91809666e+05 2.09e+00 3.77e-15 2.02e-01 42s
3 3.09423575e+04 2.93647967e+05 9.48e-02 3.10e-15 3.88e-02 44s
4 5.09533528e+04 1.59270727e+05 2.46e-02 1.77e-15 1.56e-02 45s
5 6.18389976e+04 1.16613933e+05 1.08e-02 9.91e-16 7.87e-03 46s
6 6.86470739e+04 9.91739007e+04 5.91e-03 9.35e-16 4.38e-03 47s
7 7.50884279e+04 8.61372312e+04 1.47e-03 7.13e-16 1.58e-03 49s
8 7.68713914e+04 8.18747305e+04 6.06e-04 6.57e-16 7.16e-04 50s
9 7.76784163e+04 8.03371292e+04 2.77e-04 7.68e-16 3.80e-04 51s
10 7.80329190e+04 7.95079810e+04 1.49e-04 7.68e-16 2.11e-04 52s
11 7.82694020e+04 7.90042341e+04 7.53e-05 8.80e-16 1.05e-04 54s
Barrier performed 11 iterations in 53.75 seconds (64.45 work units)
Barrier solve interrupted - model solved by another algorithm
Solved with dual simplex
Root simplex log...
Iteration Objective Primal Inf. Dual Inf. Time
886607 7.8578038e+04 0.000000e+00 0.000000e+00 54s
Root relaxation: objective 7.857804e+04, 886607 iterations, 33.40 seconds (26.51 work units)
Total elapsed time = 63.77s (DegenMoves)
Total elapsed time = 67.99s (DegenMoves)
Total elapsed time = 76.04s (DegenMoves)
Total elapsed time = 83.33s (DegenMoves)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 78578.0376 0 69583 1881.06886 78578.0376 4077% - 119s
H 0 0 77992.814417 78578.0376 0.75% - 120s
H 0 0 78071.827000 78578.0376 0.65% - 126s
H 0 0 78073.009638 78578.0376 0.65% - 130s
0 0 78414.4541 0 10318 78073.0096 78414.4541 0.44% - 141s
H 0 0 78301.124238 78414.4354 0.14% - 147s
H 0 0 78360.637798 78414.4354 0.07% - 148s
Cutting planes:
Gomory: 22
Clique: 10728
Zero half: 3472
RLT: 2132
BQP: 3
PSD: 650
Explored 1 nodes (1107150 simplex iterations) in 148.47 seconds (113.77 work units)
Thread count was 32 (of 112 available processors)
Solution count 8: 78360.6 78301.1 78073 ... -31.93
Optimal solution found (tolerance 1.00e-03)
Best objective 7.836063779777e+04, best bound 7.841443542375e+04, gap 0.0687%
Saving solution ...
Done!
Selecting tracks from first trame: 100%|█| 196/
Exporting segmentation masks: 100%|█| 210/210 [
Evaluation¶
We use the traccuracy package to replicate the cell tracking challenge metrics and compare the results.
To finalize, we export all tracks and segmentation masks to napari.
[17]:
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.layers["flow vectors"].visible = False
viewer.layers["detection"].visible = False
viewer.add_tracks(
tracks_df[["track_id", "t", "z", "y", "x"]],
name="tracks",
graph=graph,
visible=True,
)
viewer.add_labels(da.from_zarr(segments), name="segments").contour = 2
nbscreenshot(viewer)
/hpc/mydata/ilan.silva/miniforge3/envs/ultrack-flow-field/lib/python3.11/site-packages/ultrack/utils/array.py:190: FutureWarning: The NestedDirectoryStore is deprecated and will be removed in a Zarr-Python version 3, see https://github.com/zarr-developers/zarr-python/issues/1274 for more information.
store = zarr.NestedDirectoryStore(str(store_or_path))
Exporting segmentation masks: 100%|█| 210/210 [
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
[17]:
