Stable excitatory-inhibitory synapse balance despite dynamic turnover.

Code ↔ Paper

11 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 11 matches

[1] § Methods › Experimental models and subject details › Live image analysis › Drift correction ↔ ConsensusDrift/registration.py, lines 1626–1691 · score 0.77 · tiled phase correlation, composite image, consensus cluster, DBSCAN, Otsu, shift
[2] § Methods › Experimental models and subject details › Live image analysis › Puncta tracking analysis ↔ unisynapse/track/mincostflow.py, lines 2234–2286 · score 0.76 · sink node, min cost flow, source node, solver, bypass, MCF
[3] § Methods › Experimental models and subject details › Live image analysis › Puncta tracking analysis ↔ unisynapse/track/mincostflow.py, lines 70–164 · score 0.69 · nearest neighbor, Candidate links, missed frame, gap, displacement, detections
[4] § Methods › Experimental models and subject details › Live image analysis › Puncta track categorization ↔ scripts/analyze_tracks.py, lines 560–682 · score 0.65 · detection persistence, paired fraction, analyzed tracks, unpaired, duration, frame
[5] § Methods › Experimental models and subject details › Live image analysis › Puncta segmentation ↔ unisynapse/image/masks.py, lines 1187–1262 · score 0.59 · local background statistics, seed, confidence, scores, weighted, segmented
[6] § Methods › Experimental models and subject details › Live image analysis ↔ scripts/imagej_puncta_preprocess.py, lines 1–61 · score 0.59 · ImageJ, synaptic puncta, FIJI, subtracted, background
[7] § Methods › Experimental models and subject details › Live image analysis › Pairing analysis ↔ unisynapse/pair/hungarian.py, lines 18–71 · score 0.58 · overlap score, segmented masks, intersection, union, distance, pre
[8] § Methods › Experimental models and subject details › Live image analysis › Puncta tracking analysis ↔ unisynapse/config/schema.py, lines 468–528 · score 0.56 · motion smoothness weight, penalized, velocity, transit, displacement, intensity
[9] § Methods › Experimental models and subject details › Live image analysis › Puncta tracking analysis ↔ unisynapse/track/mincostflow.py, lines 70–164 · score 0.56 · motion smoothness weight, penalized, velocity, transit, displacement, intensity
[10] § Methods › Experimental models and subject details › Live image analysis › Puncta segmentation ↔ unisynapse/detect/segmentation.py, lines 209–282 · score 0.53 · low thresholds, seed, Touching, peak, blob, segmented
[11] § Methods › Experimental models and subject details › Live image analysis › Puncta tracking analysis ↔ unisynapse/track/mincostflow.py, lines 1940–2003 · score 0.51 · death costs, Google, transit, Birth, MCF, flow

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 2,591 lines · 111 KB · MIT · 4 matches

"""
Min-cost flow tracking using Google OR-Tools.
Also includes implementations of jaqaman et al. 2008 style tracking algorithm
and a hungarian tracking style algorithm for comparison purposes.
Implements multi-object tracking with:
- Distance-based linking
- Gap closing
- Birth/death handling
"""
from __future__ import annotations
from pathlib import Path
from typing import Dict, List, Optional, Tuple, Union, Set
import numpy as np
import pandas as pd
try:
from scipy.spatial import cKDTree
SCIPY_AVAILABLE = True
except ImportError:
SCIPY_AVAILABLE = False
# Try to import OR-Tools
try:
from ortools.graph.python import min_cost_flow
ORTOOLS_AVAILABLE = True
except ImportError:
ORTOOLS_AVAILABLE = False
# Try to import tqdm for progress bars
try:
from tqdm import tqdm
TQDM_AVAILABLE = True
except ImportError:
TQDM_AVAILABLE = False
# Import visualization from separate module
from .visualization import save_debug_visualizations
class MinCostFlowTracker:
"""
Multi-object tracker using min-cost flow optimization.
Parameters
----------
max_displacement_um : float
Maximum displacement per frame in µm
max_gap_frames : int
Maximum gap to close (frames)
cost_distance : float
Weight for distance cost
cost_intensity_diff : float
Weight for intensity difference cost
cost_size_diff : float
Weight for size difference cost
birth_cost : float
Cost for track initiation
death_cost : float
Cost for track termination
pixel_size_xy : float
XY pixel size in µm
pixel_size_z : float or None
Z pixel size in µm
"""
def __init__(
self,
greedy_max_displacement_um: float = 0.5,
greedy_link_sigma_um: Optional[float] = None,
greedy_max_gap_frames: int = 1,
greedy_cost_distance: float = 1.0,
greedy_cost_intensity_diff: float = 0.1,
greedy_cost_size_diff: float = 0.05,
greedy_motion_smoothness_weight: float = 0.0,
greedy_max_neighbors: Optional[int] = 5,
gap_close_penalty_per_missed_frame: float = 0.0,
use_simple_displacement_extraction: bool = False,
birth_cost: float = 10.0,
death_cost: float = 10.0,
pixel_size_xy: float = 1.0,
pixel_size_z: Optional[float] = None,
# MCF-specific parameters (optional overrides)
mcf_birth_cost: Optional[float] = None,
mcf_death_cost: Optional[float] = None,
mcf_cost_distance: Optional[float] = None,
mcf_cost_intensity_diff: Optional[float] = None,
mcf_cost_size_diff: Optional[float] = None,
mcf_max_displacement_um: Optional[float] = None,
mcf_switch_penalty: Optional[float] = None,
mcf_crowding_multiplier: Optional[float] = None,
mcf_max_neighbors: Optional[int] = None,
mcf_slack_cost_component: Optional[float] = None,
mcf_quality_bonus_weight: Optional[float] = None,
mcf_detection_bonus: Optional[float] = None,
mcf_motion_smoothness_weight: Optional[float] = None,
mcf_max_gap_frames: Optional[int] = None,
mcf_use_gaussian_link_cost: bool = False,
mcf_link_sigma_um: Optional[float] = None,
# Verbosity
verbosity_level: int = 0
):
if not ORTOOLS_AVAILABLE:
raise ImportError(
"Google OR-Tools is required for tracking. "
"Install with: pip install ortools"
)
# Greedy tracking parameters
self.max_displacement_um = greedy_max_displacement_um # deprecated name kept for internal reuse
self.max_gap_frames = greedy_max_gap_frames # deprecated name kept for internal reuse
self.cost_distance = greedy_cost_distance # deprecated name kept for internal reuse
self.cost_intensity_diff = greedy_cost_intensity_diff
self.cost_size_diff = greedy_cost_size_diff
self.greedy_link_sigma_um = greedy_link_sigma_um
self.greedy_motion_smoothness_weight = greedy_motion_smoothness_weight
self.greedy_max_neighbors = greedy_max_neighbors
self.gap_close_penalty_per_missed_frame = gap_close_penalty_per_missed_frame
self.use_simple_displacement_extraction = use_simple_displacement_extraction
self.birth_cost = birth_cost
self.death_cost = death_cost
# Default pixel size to 1.0 µm/px if unspecified to avoid NaNs
self.pixel_size_xy = pixel_size_xy if pixel_size_xy is not None else 1.0
self.pixel_size_z = pixel_size_z
# MCF-specific parameters (with intelligent defaults)
# Birth/death costs, usually set higher than greedy tracking to discourage singletons
# Otherwise flow can easily go from source->detecton->sink
# Default MCF birth/death to greedy values when unspecified
self.mcf_birth_cost = mcf_birth_cost if mcf_birth_cost is not None else self.birth_cost
self.mcf_death_cost = mcf_death_cost if mcf_death_cost is not None else self.death_cost
# Other MCF params default to same as greedy
self.mcf_cost_distance = mcf_cost_distance if mcf_cost_distance is not None else greedy_cost_distance
self.mcf_cost_intensity_diff = mcf_cost_intensity_diff if mcf_cost_intensity_diff is not None else greedy_cost_intensity_diff
self.mcf_cost_size_diff = mcf_cost_size_diff if mcf_cost_size_diff is not None else greedy_cost_size_diff
self.mcf_max_displacement_um = mcf_max_displacement_um if mcf_max_displacement_um is not None else greedy_max_displacement_um
# Optional: use Gaussian negative log-likelihood distance costs for MCF transition edges
# (off by default to preserve backwards-compatible behavior)
self.mcf_use_gaussian_link_cost = bool(mcf_use_gaussian_link_cost)
self.mcf_link_sigma_um = mcf_link_sigma_um
self.mcf_switch_penalty = mcf_switch_penalty if mcf_switch_penalty is not None else 0.0
self.mcf_crowding_multiplier = mcf_crowding_multiplier if mcf_crowding_multiplier is not None else 0.2
# Limit candidate links to k nearest neighbors (per gap). None = unlimited
self.mcf_max_neighbors = mcf_max_neighbors if mcf_max_neighbors is not None else 5
self.mcf_slack_cost_component = mcf_slack_cost_component if mcf_slack_cost_component is not None else 0.0
# Quality bonus for including detections (negative cost on transit edges)
self.mcf_quality_bonus_weight = mcf_quality_bonus_weight if mcf_quality_bonus_weight is not None else 0.5
# Flat detection bonus - incentivizes including any detection in tracks (negative cost)
self.mcf_detection_bonus = mcf_detection_bonus if mcf_detection_bonus is not None else 0.0
# Motion smoothness penalty - penalizes sudden changes in velocity (acceleration)
self.mcf_motion_smoothness_weight = (
mcf_motion_smoothness_weight
if mcf_motion_smoothness_weight is not None
else greedy_motion_smoothness_weight
)
# MCF-specific gap closing (independent from greedy tracking)
self.mcf_max_gap_frames = mcf_max_gap_frames if mcf_max_gap_frames is not None else greedy_max_gap_frames
# Verbosity level (propagated from config; 0=quiet, 1=normal, 2=verbose)
self.verbosity_level: int = verbosity_level
def _compute_link_cost(
self,
det1: pd.Series,
det2: pd.Series,
is_3d: bool,
frame_gap: int = 1,
use_mcf_costs: bool = False,
prev_displacement: Optional[np.ndarray] = None,
return_breakdown: bool = False
) -> Union[float, Tuple[float, dict]]:
"""
Compute cost to link two detections.
Parameters
----------
det1, det2 : pd.Series
Detections to link
is_3d : bool
Whether data is 3D
frame_gap : int
Number of frames between detections
use_mcf_costs : bool
If True, use MCF-specific cost parameters instead of greedy costs
prev_displacement : np.ndarray, optional
Previous displacement vector (µm) for motion smoothness calculation.
If provided, penalizes sudden changes in velocity (acceleration).
return_breakdown : bool
If True, return tuple of (cost, breakdown_dict) with cost components
Returns
-------
float or tuple
Cost value, or (cost, breakdown_dict) if return_breakdown=True
"""
# Basic finite checks
if not (np.isfinite(det1.get('x')) and np.isfinite(det1.get('y')) and np.isfinite(det2.get('x')) and np.isfinite(det2.get('y'))):
return np.inf
use_z = (
is_3d and
self.pixel_size_z is not None and
'z' in det1 and 'z' in det2 and
np.isfinite(det1.get('z')) and np.isfinite(det2.get('z'))
)
# Distance cost (in µm)
if use_z:
dist = np.sqrt(
((det1['z'] - det2['z']) * self.pixel_size_z)**2 +
((det1['y'] - det2['y']) * self.pixel_size_xy)**2 +
((det1['x'] - det2['x']) * self.pixel_size_xy)**2
)
else:
dist = np.sqrt(
((det1['y'] - det2['y']) * self.pixel_size_xy)**2 +
((det1['x'] - det2['x']) * self.pixel_size_xy)**2
)
# If distance is non-finite, treat as invalid link
if not np.isfinite(dist):
return np.inf
# Select cost parameters based on algorithm
if use_mcf_costs:
max_disp = self.mcf_max_displacement_um
cost_dist = self.mcf_cost_distance
cost_int = self.mcf_cost_intensity_diff
cost_size = self.mcf_cost_size_diff
switch_penalty = self.mcf_switch_penalty
else:
max_disp = self.max_displacement_um
cost_dist = self.cost_distance
cost_int = self.cost_intensity_diff
cost_size = self.cost_size_diff
switch_penalty = 0.0
# Scale max displacement by sqrt(frame_gap) for multi-frame gaps
# This reflects random walk scaling: displacement ~ sqrt(time)
max_disp_scaled = max_disp * np.sqrt(frame_gap)
if dist > max_disp_scaled:
return np.inf
# Distance cost
# Greedy/Hungarian (Jaqaman-style): Gaussian negative log-likelihood with variance scaling ~ frame_gap
# cost = (d^2) / (2*sigma^2), where sigma^2 = sigma_1frame^2 * frame_gap
# MCF default: keep legacy linear distance cost unless you explicitly change it.
distance_cost = 0.0
if not use_mcf_costs:
sigma_1 = getattr(self, 'greedy_link_sigma_um', None)
if sigma_1 is None:
# Interpret the max displacement gate as ~3-sigma by default
sigma_1 = max_disp / 3.0 if max_disp > 0 else 1.0
sigma2 = (sigma_1 * sigma_1) * float(frame_gap)
# Avoid division by zero
if sigma2 <= 0:
return np.inf
distance_cost = (dist * dist) / (2.0 * sigma2)
cost = distance_cost * cost_dist
else:
if getattr(self, 'mcf_use_gaussian_link_cost', False):
# MCF (optional): Gaussian negative log-likelihood with variance scaling ~ frame_gap
# cost = (d^2) / (2*sigma^2), where sigma^2 = sigma_1frame^2 * frame_gap
sigma_1 = getattr(self, 'mcf_link_sigma_um', None)
if sigma_1 is None:
# Interpret the max displacement gate as ~3-sigma by default
sigma_1 = max_disp / 3.0 if max_disp > 0 else 1.0
sigma2 = (sigma_1 * sigma_1) * float(frame_gap)
if sigma2 <= 0:
return np.inf
distance_cost = (dist * dist) / (2.0 * sigma2)
else:
# Legacy linear scaling for min-cost flow (backwards compatible)
distance_cost = dist * float(frame_gap)
cost = distance_cost * cost_dist
# Intensity difference cost
intensity_cost = 0.0
if 'signal_mean' in det1 and 'signal_mean' in det2:
int1 = det1.get('signal_mean', 0)
int2 = det2.get('signal_mean', 0)
if int1 > 0 and int2 > 0:
int_diff = abs(int1 - int2) / max(int1, int2)
intensity_cost = int_diff * cost_int
cost += intensity_cost
# Size/area difference cost
# Prefer physical size, fallback to pixel size, then sigma
size_cost = 0.0
if 'region_area_um2' in det1 and 'region_area_um2' in det2:
area1 = det1.get('region_area_um2', 0)
area2 = det2.get('region_area_um2', 0)
if area1 > 0 and area2 > 0:
area_diff = abs(area1 - area2) / max(area1, area2)
size_cost = area_diff * cost_size
cost += size_cost
elif 'region_volume_um3' in det1 and 'region_volume_um3' in det2:
vol1 = det1.get('region_volume_um3', 0)
vol2 = det2.get('region_volume_um3', 0)
if vol1 > 0 and vol2 > 0:
vol_diff = abs(vol1 - vol2) / max(vol1, vol2)
size_cost = vol_diff * cost_size
cost += size_cost
elif 'region_size_px' in det1 and 'region_size_px' in det2:
area1 = det1.get('region_size_px', 0)
area2 = det2.get('region_size_px', 0)
if area1 > 0 and area2 > 0:
area_diff = abs(area1 - area2) / max(area1, area2)
size_cost = area_diff * cost_size
cost += size_cost
elif 'sigma' in det1 and 'sigma' in det2:
s1 = det1.get('sigma', 1)
s2 = det2.get('sigma', 1)
if s1 > 0 and s2 > 0:
size_diff = abs(s1 - s2) / max(s1, s2)
size_cost = size_diff * cost_size
cost += size_cost
# Amplify cost based on ambiguity: more link candidates = more ambiguity = higher cost
# This directly addresses identity switching in dense synapse fields
ambig_col = f'ambiguity_gap{frame_gap}'
ambig_factor = det1.get(ambig_col, 1.0)
if ambig_factor < 0 or not np.isfinite(ambig_factor):
ambig_factor = 1.0
cost *= ambig_factor
# Constant penalty to discourage switching between unrelated puncta
cost += switch_penalty
# Motion smoothness penalty: penalize sudden changes in velocity (acceleration)
# Applies to both greedy (Hungarian) and MCF paths
smoothness_penalty = 0.0
smoothness_weight = self.mcf_motion_smoothness_weight if use_mcf_costs else getattr(self, 'greedy_motion_smoothness_weight', 0.0)
if prev_displacement is not None and smoothness_weight and smoothness_weight > 0:
# Current displacement vector (µm)
if use_z:
curr_displacement = np.array([
(det2['x'] - det1['x']) * self.pixel_size_xy,
(det2['y'] - det1['y']) * self.pixel_size_xy,
(det2['z'] - det1['z']) * self.pixel_size_z
])
else:
curr_displacement = np.array([
(det2['x'] - det1['x']) * self.pixel_size_xy,
(det2['y'] - det1['y']) * self.pixel_size_xy
])
# Compute change in velocity (acceleration vector)
# For multi-frame gaps, normalize by frame_gap to get velocity change per frame
acceleration = (curr_displacement - prev_displacement) / frame_gap
acceleration_magnitude = np.linalg.norm(acceleration)
# Penalty proportional to acceleration magnitude
# Higher acceleration = less smooth motion = higher cost
smoothness_penalty = smoothness_weight * acceleration_magnitude
cost += smoothness_penalty
if return_breakdown:
# Return cost along with breakdown of components
breakdown = {
'distance': dist,
'distance_cost': float(distance_cost) * float(cost_dist),
'intensity_cost': intensity_cost,
'size_cost': size_cost,
'total_cost': cost,
'ambiguity_factor': ambig_factor,
'smoothness_penalty': smoothness_penalty,
'switch_penalty': switch_penalty,
'max_disp_scaled': max_disp_scaled
}
return cost, breakdown
return cost
def _add_displacement_from_greedy_tracks(
self,
detections_by_frame: Dict[int, pd.DataFrame],
greedy_tracks: pd.DataFrame,
is_3d: bool
) -> Dict[int, pd.DataFrame]:
"""
Add previous displacement vectors to detections based on greedy tracking results.
This enables motion smoothness penalties in MCF by providing the expected
velocity from Stage 1 greedy tracking.
Notes
-----
We key greedy-track rows back to detections using (frame, det_id). If det_id is
missing in greedy_tracks (older outputs), we fall back to using the greedy_tracks
row index as det_id (best-effort), but the preferred/accurate path is to include
det_id in greedy outputs.
"""
if greedy_tracks is None or len(greedy_tracks) == 0:
return detections_by_frame
gt = greedy_tracks.copy()
# Ensure we have a stable per-detection id
if 'det_id' not in gt.columns:
# Best-effort fallback: use row index
gt = gt.reset_index(drop=False).rename(columns={'index': 'det_id'})
# Compute per-row displacement from previous frame within the same greedy track
gt = gt.sort_values(['track_id', 'frame']).reset_index(drop=True)
gt['prev_frame'] = gt.groupby('track_id')['frame'].shift(1)
gt['prev_x'] = gt.groupby('track_id')['x'].shift(1)
gt['prev_y'] = gt.groupby('track_id')['y'].shift(1)
if is_3d and self.pixel_size_z is not None and 'z' in gt.columns:
gt['prev_z'] = gt.groupby('track_id')['z'].shift(1)
valid_prev = gt['prev_frame'].notna() & (gt['prev_frame'] == (gt['frame'] - 1))
gt['prev_disp_x_um'] = np.nan
gt['prev_disp_y_um'] = np.nan
if is_3d and self.pixel_size_z is not None and 'z' in gt.columns:
gt['prev_disp_z_um'] = np.nan
gt.loc[valid_prev, 'prev_disp_x_um'] = (gt.loc[valid_prev, 'x'] - gt.loc[valid_prev, 'prev_x']) * self.pixel_size_xy
gt.loc[valid_prev, 'prev_disp_y_um'] = (gt.loc[valid_prev, 'y'] - gt.loc[valid_prev, 'prev_y']) * self.pixel_size_xy
if is_3d and self.pixel_size_z is not None and 'z' in gt.columns:
# Only assign Z displacement when both z and prev_z are finite
z_ok = valid_prev & gt['z'].notna() & gt['prev_z'].notna()
gt.loc[z_ok, 'prev_disp_z_um'] = (gt.loc[z_ok, 'z'] - gt.loc[z_ok, 'prev_z']) * self.pixel_size_z
disp_cols = ['frame', 'det_id', 'prev_disp_x_um', 'prev_disp_y_um']
if is_3d and self.pixel_size_z is not None and 'prev_disp_z_um' in gt.columns:
disp_cols.append('prev_disp_z_um')
disp = gt[disp_cols].copy()
# Apply to detections_by_frame
result: Dict[int, pd.DataFrame] = {}
for frame, df in detections_by_frame.items():
df_copy = df.copy()
# Make sure detections have det_id (use original df index if absent)
if 'det_id' not in df_copy.columns:
df_copy['det_id'] = df_copy.index.astype(int)
# Initialize columns
df_copy['prev_disp_x_um'] = np.nan
df_copy['prev_disp_y_um'] = np.nan
if is_3d and self.pixel_size_z is not None:
df_copy['prev_disp_z_um'] = np.nan
dframe = disp[disp['frame'] == frame].drop(columns=['frame'], errors='ignore')
if not dframe.empty:
df_copy = df_copy.merge(dframe, on=['det_id'], how='left', suffixes=('', '_gt'))
# Prefer merged values (from greedy) where available
df_copy['prev_disp_x_um'] = df_copy['prev_disp_x_um_gt'].combine_first(df_copy['prev_disp_x_um'])
df_copy['prev_disp_y_um'] = df_copy['prev_disp_y_um_gt'].combine_first(df_copy['prev_disp_y_um'])
if is_3d and self.pixel_size_z is not None and 'prev_disp_z_um_gt' in df_copy.columns:
df_copy['prev_disp_z_um'] = df_copy['prev_disp_z_um_gt'].combine_first(df_copy.get('prev_disp_z_um'))
# Drop temp columns
drop_cols = [c for c in df_copy.columns if c.endswith('_gt')]
df_copy = df_copy.drop(columns=drop_cols)
# Keep the original index ordering
result[frame] = df_copy
return result
def _add_ambiguity_factors(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool,
max_gap: int
) -> Dict[int, pd.DataFrame]:
"""
Compute per-detection ambiguity factors based on number of link candidates.
For each detection at frame t, counts how many potential link targets exist
at frame t+dt within max_disp_scaled = mcf_max_displacement_um * sqrt(dt).
This directly addresses identity switching: more candidates = more ambiguity = higher cost.
Parameters
----------
detections_by_frame : dict
Detections by frame
is_3d : bool
Whether 3D data
max_gap : int
Maximum gap to consider (max_gap_frames + 1)
Returns
-------
dict
Detections with ambiguity_gap{dt} columns added
"""
if not SCIPY_AVAILABLE:
return detections_by_frame
max_disp = self.mcf_max_displacement_um if self.mcf_max_displacement_um is not None else self.max_displacement_um
if max_disp <= 0:
return detections_by_frame
frames = sorted(detections_by_frame.keys())
result = {}
for frame in frames:
df = detections_by_frame[frame]
if len(df) == 0:
result[frame] = df
continue
df_copy = df.copy()
# Get coords for current frame
if is_3d and self.pixel_size_z is not None and 'z' in df.columns:
curr_coords = np.stack([
df['x'].to_numpy() * self.pixel_size_xy,
df['y'].to_numpy() * self.pixel_size_xy,
df['z'].to_numpy() * self.pixel_size_z
], axis=1)
else:
curr_coords = np.stack([
df['x'].to_numpy() * self.pixel_size_xy,
df['y'].to_numpy() * self.pixel_size_xy
], axis=1)
# For each gap (1 to max_gap), count candidates at next_frame
for gap in range(1, max_gap + 1):
next_frame = frame + gap
if next_frame not in detections_by_frame:
# No next frame, set ambiguity to 0
df_copy[f'ambiguity_gap{gap}'] = 0.0
continue
df_next = detections_by_frame[next_frame]
if len(df_next) == 0:
df_copy[f'ambiguity_gap{gap}'] = 0.0
continue
# Get coords for next frame
if is_3d and self.pixel_size_z is not None and 'z' in df_next.columns:
next_coords = np.stack([
df_next['x'].to_numpy() * self.pixel_size_xy,
df_next['y'].to_numpy() * self.pixel_size_xy,
df_next['z'].to_numpy() * self.pixel_size_z
], axis=1)
else:
next_coords = np.stack([
df_next['x'].to_numpy() * self.pixel_size_xy,
df_next['y'].to_numpy() * self.pixel_size_xy
], axis=1)
# Build tree for next frame
next_tree = cKDTree(next_coords)
# Scale max displacement by sqrt(gap) for multi-frame gaps
max_disp_scaled = max_disp * np.sqrt(gap)
# Count candidates for each detection in current frame
candidate_counts = []
for coord in curr_coords:
candidates = next_tree.query_ball_point(coord, r=max_disp_scaled)
candidate_counts.append(len(candidates))
candidate_counts = np.array(candidate_counts, dtype=float)
# Ambiguity factor: 1 + multiplier * min(max(candidates-1, 0), 10)
# candidates-1 because 1 candidate = no ambiguity
ambiguity = np.maximum(candidate_counts - 1, 0)
ambig_factor = 1.0 + self.mcf_crowding_multiplier * np.minimum(ambiguity, 10.0)
df_copy[f'ambiguity_gap{gap}'] = ambig_factor
result[frame] = df_copy
return result
def track(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool = False,
verbose: bool = False,
use_min_cost_flow: bool = False,
mcf_crowding_multiplier: Optional[float] = None,
mcf_max_neighbors: Optional[int] = None,
mcf_slack_cost_component: Optional[float] = None,
output_dir: Optional[Path] = None,
track_type: str = "",
type_id: int = 0,
debug_images: bool = False,
precomputed_displacements: Optional[Dict[int, pd.DataFrame]] = None,
debug_dir: Optional[Path] = None
) -> pd.DataFrame:
"""
Track detections across frames.
When use_min_cost_flow=True:
1. Extract frame-pair displacements
2. Run min-cost flow on all detections with motion-aware costs
When use_min_cost_flow=False:
1. Extract frame-pair displacements
2. Run greedy tracking with Hungarian assignment and motion-aware costs
Parameters
----------
detections_by_frame : dict
Mapping from frame number to DataFrame of detections
is_3d : bool
Whether data is 3D
verbose : bool
Print progress information
use_min_cost_flow : bool
Use min-cost flow optimization (with quality bonus filtering)
output_dir : Path or None
Directory to save track statistics CSV (for debugging)
debug_dir : Path or None
Override directory for tracking debug images
Returns
-------
tracks : DataFrame
Tracking results with track_id assigned
"""
frames = sorted(detections_by_frame.keys())
if len(frames) == 0:
return pd.DataFrame()
debug_output_dir = debug_dir if debug_dir is not None else output_dir
has_debug_output = debug_output_dir is not None
if use_min_cost_flow:
if verbose:
print(" [MCF] Using min-cost flow tracking")
# Check if we need motion calculations (motion smoothness)
motion_smoothness_enabled = (
self.mcf_motion_smoothness_weight is not None and self.mcf_motion_smoothness_weight != 0
)
need_motion = motion_smoothness_enabled
if need_motion:
# Use pre-computed displacements if provided, otherwise compute them
if precomputed_displacements is not None:
if verbose:
print(f" [MCF] Using pre-computed displacement vectors")
detections_with_displacement = precomputed_displacements
else:
if verbose:
print(f" [MCF] Stage 1: Frame-pair extraction for motion estimation")
detections_with_displacement = self._extract_displacement_parallel(
detections_by_frame,
is_3d,
verbose=verbose
)
# Add ambiguity factors to penalize high-ambiguity links (discourage identity switches)
detections_with_ambiguity = self._add_ambiguity_factors(
detections_with_displacement, is_3d, max_gap=self.mcf_max_gap_frames + 1
)
if verbose:
total_dets = sum(len(df) for df in detections_with_ambiguity.values())
print(f" [MCF] Stage 2: Running min-cost flow on all {total_dets} detections")
print(f" [MCF] Using mcf_max_gap_frames={self.mcf_max_gap_frames} for gap closing")
print(f" [MCF] Quality bonus will filter low-confidence detections")
else:
if verbose:
print(f" [MCF] Skipping Stage 1 (motion calculations disabled)")
print(f" [MCF] mcf_motion_smoothness_weight={self.mcf_motion_smoothness_weight}")
# Skip Stage 1 entirely - only add ambiguity factors for crowding penalty
detections_with_ambiguity = self._add_ambiguity_factors(
detections_by_frame, is_3d, max_gap=self.mcf_max_gap_frames + 1
)
if verbose:
total_dets = sum(len(df) for df in detections_with_ambiguity.values())
print(f" [MCF] Running min-cost flow on all {total_dets} detections")
print(f" [MCF] Using mcf_max_gap_frames={self.mcf_max_gap_frames} for gap closing")
print(f" [MCF] Quality bonus will filter low-confidence detections")
# Stage 2: Run min-cost flow on ALL detections
mcf_tracks = self._min_cost_flow_track(
detections_with_ambiguity,
is_3d,
verbose=verbose,
mcf_max_neighbors=mcf_max_neighbors,
mcf_slack_cost_component=mcf_slack_cost_component,
output_dir=output_dir,
debug_dir=debug_dir,
track_type=track_type,
type_id=type_id,
debug_images=debug_images
)
if len(mcf_tracks) > 0:
n_mcf_tracks = mcf_tracks['track_id'].nunique()
if verbose:
mcf_lengths = mcf_tracks.groupby('track_id').size()
print(f" [MCF] Output: {n_mcf_tracks} tracks, lengths: min={mcf_lengths.min()}, mean={mcf_lengths.mean():.1f}, max={mcf_lengths.max()}")
else:
print(f" [MCF] WARNING: Min-cost flow returned NO tracks!")
print(f" [MCF] Falling back to greedy tracking (full 2-pass)...")
# Run full greedy tracking (motion_only=False to get actual tracks)
return self._greedy_track(
detections_by_frame,
is_3d,
verbose=verbose,
motion_only=False,
output_dir=output_dir,
track_type=track_type,
type_id=type_id,
debug_images=debug_images
)
return mcf_tracks
else:
if verbose:
print(" Using greedy tracking (fast, locally optimal)")
return self._greedy_track(
detections_by_frame,
is_3d,
verbose=verbose,
motion_only=False,
output_dir=output_dir,
debug_dir=debug_dir,
track_type=track_type,
type_id=type_id,
debug_images=debug_images
)
def _greedy_track(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool,
verbose: bool = False,
motion_only: bool = False,
output_dir: Optional[Path] = None,
debug_dir: Optional[Path] = None,
track_type: str = "",
type_id: int = 0,
debug_images: bool = False
) -> pd.DataFrame:
"""
Simple greedy tracking with gap closing.
Two-stage flow (used for both motion-only pass and full greedy tracking):
1. Stage 1: Batch frame-pair displacement extraction.
2. Stage 2: Hungarian assignment with optional smoothness penalty.
Parameters
----------
motion_only : bool
If True, only compute motion vectors (Pass 1) and return empty DataFrame.
Used when MCF will do the actual tracking.
"""
frames = sorted(detections_by_frame.keys())
if len(frames) == 0:
return pd.DataFrame()
# Stage 1: displacement extraction (if motion smoothness enabled)
need_motion = getattr(self, 'greedy_motion_smoothness_weight', 0.0) > 0
detections_stage1 = detections_by_frame
if need_motion:
if verbose:
print(f" [GREEDY] Stage 1: Extracting displacements")
detections_stage1 = self._extract_displacement_parallel(
detections_by_frame,
is_3d,
verbose=verbose
)
if motion_only:
return pd.DataFrame()
# Stage 2: Hungarian assignment with motion-aware costs
try:
from scipy.optimize import linear_sum_assignment # noqa: F401
return self._greedy_track_hungarian(
detections_stage1,
is_3d,
verbose=verbose,
output_dir=output_dir,
debug_dir=debug_dir,
track_type=track_type,
type_id=type_id,
debug_images=debug_images
)
except Exception as e:
if verbose and not getattr(self, "_hungarian_warned", False):
print(f" [GREEDY] Warning: Hungarian assignment unavailable ({e}); falling back to simple greedy pass")
print(f" (If this repeats, try increasing greedy_max_displacement_um or reducing smoothness)")
self._hungarian_warned = True
return self._greedy_track_single_pass(detections_by_frame, is_3d, verbose)
def _greedy_track_hungarian(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool,
verbose: bool = False,
output_dir: Optional[Path] = None,
debug_dir: Optional[Path] = None,
track_type: str = "",
type_id: int = 0,
debug_images: bool = False
) -> pd.DataFrame:
"""
Jaqaman-style 2-pass greedy tracking:
Pass 1: frame-to-frame LAP (Hungarian), ONLY frame_gap==1 -> produces tracklets
Pass 2: post hoc gap closing LAP between tracklet ends and starts
"""
try:
from scipy.optimize import linear_sum_assignment
except ImportError:
raise
frames = sorted(detections_by_frame.keys())
if not frames:
return pd.DataFrame()
debug_output_dir = debug_dir if debug_dir is not None else output_dir
has_debug_output = debug_output_dir is not None
all_results: List[dict] = []
next_track_id = 0
# Active track heads for frame-to-frame linking ONLY
# list of (track_id, last_frame, last_det_series)
active_tracks: List[Tuple[int, int, pd.Series]] = []
# Debug bookkeeping (kept consistent with your current interface)
start_detections: List[Tuple[int, int]] = []
link_distances: List[float] = []
death_to_nextframe_dists: List[float] = []
track_costs: Dict[int, float] = {}
cost_breakdowns: List[dict] = []
ambiguity_factors: List[float] = []
smoothness_penalties: List[float] = []
BIG_M = 1e9
# Setup progress bar if verbose
if verbose and TQDM_AVAILABLE:
frame_iter = tqdm(frames, desc=" [GREEDY-P1] Frame-to-frame linking", unit="frame")
else:
frame_iter = frames
for frame in frame_iter:
df = detections_by_frame[frame]
if len(df) == 0:
# In tracklet mode, if a frame has no detections, all active tracks terminate.
active_tracks = []
continue
# Preserve original detection id for later mapping (e.g., displacement extraction)
if 'det_id' not in df.columns:
df = df.reset_index(drop=False).rename(columns={'index': 'det_id'})
else:
df = df.reset_index(drop=True)
# If no active tracks, every detection starts a new tracklet
if not active_tracks:
for det_idx, det in df.iterrows():
tid = next_track_id
next_track_id += 1
rec = det.to_dict()
rec["track_id"] = tid
rec["frame"] = frame
all_results.append(rec)
active_tracks.append((tid, frame, det))
start_detections.append((frame, det_idx))
continue
num_tracks = len(active_tracks)
num_dets = len(df)
# Cost matrix tracks x detections (only for frame_gap==1)
cost_matrix = np.full((num_tracks, num_dets), np.inf, dtype=float)
# Build k-NN structure over active tracks from previous frame
k_neighbors = self.greedy_max_neighbors
use_knn = SCIPY_AVAILABLE and k_neighbors is not None and k_neighbors > 0
track_coords = []
track_indices = []
if use_knn:
for t_idx, (_, last_frame, last_det) in enumerate(active_tracks):
if frame - last_frame != 1:
continue
if not (np.isfinite(last_det.get("x")) and np.isfinite(last_det.get("y"))):
continue
if (
is_3d and self.pixel_size_z is not None
and "z" in last_det.index and np.isfinite(last_det.get("z"))
):
coord = [
last_det["x"] * self.pixel_size_xy,
last_det["y"] * self.pixel_size_xy,
last_det["z"] * self.pixel_size_z
]
else:
coord = [
last_det["x"] * self.pixel_size_xy,
last_det["y"] * self.pixel_size_xy
]
track_coords.append(coord)
track_indices.append(t_idx)
if not track_coords:
use_knn = False
else:
tree = cKDTree(np.array(track_coords))
# Fill costs: compute candidates per detection
for det_idx, det in df.iterrows():
if not (np.isfinite(det.get("x")) and np.isfinite(det.get("y"))):
continue
# Candidate tracks (all or k nearest)
if k_neighbors is None or k_neighbors == 0:
candidate_order = range(num_tracks)
elif use_knn:
if (
is_3d and self.pixel_size_z is not None
and "z" in det.index and np.isfinite(det.get("z"))
):
det_coord = [
det["x"] * self.pixel_size_xy,
det["y"] * self.pixel_size_xy,
det["z"] * self.pixel_size_z
]
else:
det_coord = [
det["x"] * self.pixel_size_xy,
det["y"] * self.pixel_size_xy
]
k = min(int(k_neighbors), len(track_coords))
dists, idxs = tree.query(det_coord, k=k)
if np.isscalar(idxs):
idxs = [int(idxs)]
else:
idxs = [int(i) for i in np.atleast_1d(idxs)]
candidate_order = [track_indices[i] for i in idxs]
else:
dists = []
for t_idx, (_, last_frame, last_det) in enumerate(active_tracks):
if frame - last_frame != 1:
continue
if not (np.isfinite(last_det.get("x")) and np.isfinite(last_det.get("y"))):
dist_val = np.inf
else:
if (
is_3d and self.pixel_size_z is not None
and "z" in det.index and "z" in last_det.index
and np.isfinite(det.get("z")) and np.isfinite(last_det.get("z"))
):
dx = (det["x"] - last_det["x"]) * self.pixel_size_xy
dy = (det["y"] - last_det["y"]) * self.pixel_size_xy
dz = (det["z"] - last_det["z"]) * self.pixel_size_z
dist_val = float(np.sqrt(dx * dx + dy * dy + dz * dz))
else:
dx = (det["x"] - last_det["x"]) * self.pixel_size_xy
dy = (det["y"] - last_det["y"]) * self.pixel_size_xy
dist_val = float(np.sqrt(dx * dx + dy * dy))
dists.append((dist_val, t_idx))
dists.sort(key=lambda x: x[0])
k = max(int(k_neighbors), 0)
candidate_order = [idx for _, idx in dists[:k]] if k > 0 else []
for t_idx in candidate_order:
track_id, last_frame, last_det = active_tracks[t_idx]
# Enforce frame-to-frame linking only
if frame - last_frame != 1:
continue
prev_disp = None
if "prev_disp_x_um" in last_det.index and pd.notna(last_det.get("prev_disp_x_um")):
if is_3d and self.pixel_size_z is not None and "prev_disp_z_um" in last_det.index:
prev_disp = np.array([
last_det["prev_disp_x_um"],
last_det["prev_disp_y_um"],
last_det["prev_disp_z_um"],
])
else:
prev_disp = np.array([
last_det["prev_disp_x_um"],
last_det["prev_disp_y_um"],
])
base_cost = self._compute_link_cost(
last_det, det, is_3d, frame_gap=1,
use_mcf_costs=False,
prev_displacement=prev_disp,
return_breakdown=bool(debug_images)
)
if isinstance(base_cost, tuple):
cost_val, breakdown = base_cost
else:
cost_val, breakdown = base_cost, None
if breakdown is not None:
bd = breakdown.copy()
bd['total_cost'] = cost_val
cost_breakdowns.append(bd)
if 'ambiguity_factor' in bd and np.isfinite(bd['ambiguity_factor']):
ambiguity_factors.append(float(bd['ambiguity_factor']))
if bd.get('smoothness_penalty', 0) > 0:
smoothness_penalties.append(float(bd['smoothness_penalty']))
cost = cost_val
if np.isfinite(cost):
cost_matrix[t_idx, det_idx] = cost
# --- Pass 1 Hungarian with explicit birth/death costs (square assignment) ---
# base is cost_matrix (T x D) with BIG_M where invalid
base = cost_matrix.copy()
base[~np.isfinite(base)] = BIG_M
T, D = base.shape
# Square matrix:
# rows: 0..T-1 are tracks, T..T+D-1 are "birth rows" (one per detection)
# cols: 0..D-1 are detections, D..D+T-1 are "termination cols" (one per track)
M = np.full((T + D, D + T), BIG_M, dtype=float)
# Track -> detection link costs
M[:T, :D] = base
# Track -> terminate (death) costs on diagonal
# Each track row can choose its own termination column at cost death_cost
term_block = M[:T, D:D+T]
np.fill_diagonal(term_block, float(self.death_cost))
# Birth -> detection costs on diagonal
# Each detection can be "born" by its own birth row at cost birth_cost
birth_block = M[T:T+D, :D]
np.fill_diagonal(birth_block, float(self.birth_cost))
# Birth rows -> termination cols are dummy-to-dummy; set 0 so unused births/terms don't matter
M[T:T+D, D:D+T] = 0.0
row_ind, col_ind = linear_sum_assignment(M)
matched_dets = set()
birth_dets = set()
new_active_tracks = []
# Apply assignments
for r, c in zip(row_ind, col_ind):
# Track rows
if r < T:
# Track -> detection (extend tracklet)
if c < D and M[r, c] < BIG_M * 0.5:
track_id, last_frame, prev_det = active_tracks[r]
det = df.iloc[c]
rec = det.to_dict()
rec["track_id"] = track_id
rec["frame"] = frame
all_results.append(rec)
new_active_tracks.append((track_id, frame, det))
matched_dets.add(int(c))
# Debug link distance bookkeeping (unchanged from your code)
if (
np.isfinite(det.get("x")) and np.isfinite(det.get("y"))
and np.isfinite(prev_det.get("x")) and np.isfinite(prev_det.get("y"))
):
if (
is_3d and self.pixel_size_z is not None
and "z" in det.index and "z" in prev_det.index
and np.isfinite(det.get("z")) and np.isfinite(prev_det.get("z"))
):
dx = (det["x"] - prev_det["x"]) * self.pixel_size_xy
dy = (det["y"] - prev_det["y"]) * self.pixel_size_xy
dz = (det["z"] - prev_det["z"]) * self.pixel_size_z
link_distances.append(float(np.sqrt(dx * dx + dy * dy + dz * dz)))
else:
dx = (det["x"] - prev_det["x"]) * self.pixel_size_xy
dy = (det["y"] - prev_det["y"]) * self.pixel_size_xy
link_distances.append(float(np.sqrt(dx * dx + dy * dy)))
# else: track terminated (death) -> do not carry forward
# Birth rows
else:
# Birth row for detection j is row r = T + j, and only diag to col j is finite
j = int(r - T)
if c == j and c < D and M[r, c] < BIG_M * 0.5:
birth_dets.add(j)
# Create new tracks for detections chosen as births
for det_idx in sorted(birth_dets):
det = df.iloc[det_idx]
tid = next_track_id
next_track_id += 1
rec = det.to_dict()
rec["track_id"] = tid
rec["frame"] = frame
all_results.append(rec)
new_active_tracks.append((tid, frame, det))
start_detections.append((frame, det_idx))
# IMPORTANT: detections assigned to tracks are already handled above.
# Any detection not matched AND not explicitly birthed should not happen here
# (it would imply the solver paid BIG_M), but keep a safety fallback:
for det_idx, det in df.iterrows():
if det_idx in matched_dets or det_idx in birth_dets:
continue
if verbose:
print(f"[GREEDY] WARNING: detection {det_idx} forced into fallback birth (likely BIG_M assignment)")
tid = next_track_id
next_track_id += 1
rec = det.to_dict()
rec["track_id"] = tid
rec["frame"] = frame
all_results.append(rec)
new_active_tracks.append((tid, frame, det))
start_detections.append((frame, det_idx))
active_tracks = new_active_tracks
if not all_results:
return pd.DataFrame()
tracks_df = pd.DataFrame(all_results)
# ---- Pass 2: post hoc gap closing (Jaqaman-style 2nd LAP) ----
tracks_df = self._gap_close_jaqaman_lap(
tracks_df=tracks_df,
is_3d=is_3d,
verbose=verbose
)
# Debug visualizations (greedy path)
if debug_images and has_debug_output:
try:
save_debug_visualizations(
detections_by_frame=detections_by_frame,
all_results=tracks_df.to_dict("records"),
link_distances=link_distances,
start_detections=start_detections,
death_to_nextframe_dists=death_to_nextframe_dists,
track_costs=track_costs,
is_3d=is_3d,
output_dir=debug_output_dir,
track_type=track_type,
type_id=type_id,
verbose=verbose,
pixel_size_xy=self.pixel_size_xy,
max_displacement_um=self.max_displacement_um,
mcf_max_displacement_um=self.mcf_max_displacement_um,
mcf_cost_distance=self.mcf_cost_distance,
mcf_crowding_multiplier=self.mcf_crowding_multiplier,
cost_breakdowns=cost_breakdowns if cost_breakdowns else None,
ambiguity_factors=ambiguity_factors if ambiguity_factors else None,
quality_bonuses=None,
smoothness_penalties=smoothness_penalties if smoothness_penalties else None,
dpi=150,
debug_dir=debug_dir
)
except Exception as e:
print(f" [GREEDY] Warning: could not save debug visualizations: {e}")
return tracks_df
def _gap_close_jaqaman_lap(
self,
tracks_df: pd.DataFrame,
is_3d: bool,
verbose: bool = False
) -> pd.DataFrame:
"""
Pass 2 (Jaqaman-like): gap closing LAP between tracklet ends and starts.
We:
- compute each tracklet start (first detection) and end (last detection)
- allow end->start links for gaps in [2, max_gap_frames+1]
- cost = _compute_link_cost(end_det, start_det, frame_gap=gap) + optional gap penalty
- compare against terminate+initiate by using an effective cost:
eff_cost = link_cost - (death_cost + birth_cost)
and dummy no-link cost = 0.
This means we only link when it is cheaper than ending+starting.
Notes:
- This implements gap closing only (no merge/split LAP here).
- It produces one-in/one-out per tracklet -> chains, which we then relabel.
"""
if tracks_df.empty:
return tracks_df
try:
from scipy.optimize import linear_sum_assignment
except ImportError:
return tracks_df
# If max_gap_frames==0, there is no post hoc gap closing to do.
max_gap = int(getattr(self, "max_gap_frames", 0)) + 1 # frame_gap allowed: 2..max_gap
if max_gap < 2:
return tracks_df
BIG_M = 1e9
# Optional per-missed-frame penalty (tuneable; default 0)
# This discourages long gaps even if distance is small.
gap_penalty_per_missed_frame = float(getattr(self, "gap_close_penalty_per_missed_frame", 0.0))
# Sort and get start/end rows per tracklet
tracks_sorted = tracks_df.sort_values(["track_id", "frame"]).reset_index(drop=True)
grp = tracks_sorted.groupby("track_id", sort=False)
starts = grp.head(1).copy()
ends = grp.tail(1).copy()
# Map from row index in starts/ends -> track_id
start_tids = starts["track_id"].to_numpy()
end_tids = ends["track_id"].to_numpy()
start_frames = starts["frame"].to_numpy().astype(int)
end_frames = ends["frame"].to_numpy().astype(int)
# Build lookup of starts by frame for fast candidate search
starts_by_frame: Dict[int, List[int]] = {}
for j, f in enumerate(start_frames):
starts_by_frame.setdefault(int(f), []).append(j)
n_ends = len(ends)
n_starts = len(starts)
if n_ends == 0 or n_starts == 0:
return tracks_df
baseline = float(self.birth_cost) + float(self.death_cost)
# Extract NumPy arrays before loop for vectorized operations
ends_x = ends['x'].to_numpy()
ends_y = ends['y'].to_numpy()
ends_frame = ends['frame'].to_numpy().astype(int)
starts_x = starts['x'].to_numpy()
starts_y = starts['y'].to_numpy()
# Handle 3D coordinates
use_z = is_3d and self.pixel_size_z is not None
if use_z and 'z' in ends.columns and 'z' in starts.columns:
ends_z = ends['z'].to_numpy()
starts_z = starts['z'].to_numpy()
else:
ends_z = None
starts_z = None
use_z = False
# Optional properties
has_signal = 'signal_mean' in ends.columns and 'signal_mean' in starts.columns
has_size = 'region_size_px' in ends.columns and 'region_size_px' in starts.columns
if has_signal:
ends_signal = ends['signal_mean'].fillna(0).to_numpy()
starts_signal = starts['signal_mean'].fillna(0).to_numpy()
if has_size:
ends_size = ends['region_size_px'].fillna(0).to_numpy()
starts_size = starts['region_size_px'].fillna(0).to_numpy()
# Cost parameters (greedy, not MCF)
max_disp = self.max_displacement_um
cost_dist = self.cost_distance
cost_int = self.cost_intensity_diff
cost_size = self.cost_size_diff
# Optional k-NN pruning for gap closing (defaults to greedy_max_neighbors)
gap_knn = getattr(self, "gap_close_max_neighbors", getattr(self, "greedy_max_neighbors", None))
# Collect sparse candidate costs to avoid huge dense matrices
candidates: List[Tuple[int, int, float]] = []
used_rows: set = set()
used_cols: set = set()
# Setup progress bar for gap closing
if verbose and TQDM_AVAILABLE:
end_iter = tqdm(range(n_ends), desc=" [GREEDY-P2] Gap closing", unit="end")
else:
end_iter = range(n_ends)
# Fill candidate costs with gating (vectorized inner loop)
for i in end_iter:
f_end = ends_frame[i]
# Consider gap frames 2..max_gap (i.e., missing 1..max_gap-1 frames)
for gap in range(2, max_gap + 1):
f_start = f_end + gap
js = starts_by_frame.get(f_start, [])
if not js:
continue
# Convert js list to numpy array for vectorized indexing
js_arr = np.array(js, dtype=int)
# Vectorized distance computation
dx = (starts_x[js_arr] - ends_x[i]) * self.pixel_size_xy
dy = (starts_y[js_arr] - ends_y[i]) * self.pixel_size_xy
if use_z:
dz = (starts_z[js_arr] - ends_z[i]) * self.pixel_size_z
distances = np.sqrt(dx*dx + dy*dy + dz*dz)
else:
distances = np.sqrt(dx*dx + dy*dy)
# Scale max displacement by sqrt(frame_gap) for multi-frame gaps
max_disp_scaled = max_disp * np.sqrt(gap)
# Filter by distance gating
valid_mask = distances <= max_disp_scaled
if not np.any(valid_mask):
continue
# Base costs (scale distance cost by frame gap)
link_costs = distances * cost_dist * gap
# Add intensity differences (vectorized)
if has_signal:
end_sig = ends_signal[i]
start_sigs = starts_signal[js_arr]
# Only add cost if both signals are positive
sig_valid = (end_sig > 0) & (start_sigs > 0)
if np.any(sig_valid):
int_diffs = np.abs(start_sigs - end_sig) / np.maximum(start_sigs, end_sig)
link_costs = np.where(sig_valid, link_costs + int_diffs * cost_int, link_costs)
# Add size differences (vectorized)
if has_size:
end_sz = ends_size[i]
start_szs = starts_size[js_arr]
# Only add cost if both sizes are positive
size_valid = (end_sz > 0) & (start_szs > 0)
if np.any(size_valid):
size_diffs = np.abs(start_szs - end_sz) / np.maximum(start_szs, end_sz)
link_costs = np.where(size_valid, link_costs + size_diffs * cost_size, link_costs)
# Apply gap penalty per missed frame
if gap_penalty_per_missed_frame != 0.0:
link_costs += gap_penalty_per_missed_frame * (gap - 1)
# Compute effective costs
eff_costs = link_costs - baseline
# Apply distance gating
eff_costs = np.where(valid_mask, eff_costs, np.inf)
# Optional k-NN pruning within this gap for this end
if gap_knn is not None and gap_knn > 0:
# Select top-k lowest finite costs
finite_idxs = np.where(np.isfinite(eff_costs))[0]
if len(finite_idxs) == 0:
continue
k = min(gap_knn, len(finite_idxs))
topk_local = finite_idxs[np.argpartition(eff_costs[finite_idxs], k-1)[:k]]
js_candidates = js_arr[topk_local]
eff_candidates = eff_costs[topk_local]
else:
js_candidates = js_arr
eff_candidates = eff_costs
# Store sparse candidates
for j_idx, eff in zip(js_candidates, eff_candidates):
if not np.isfinite(eff):
continue
candidates.append((i, int(j_idx), float(eff)))
used_rows.add(i)
used_cols.add(int(j_idx))
if not candidates:
return tracks_df
# Build compact cost matrix using only involved rows/cols
row_list = sorted(used_rows)
col_list = sorted(used_cols)
row_index = {r: idx for idx, r in enumerate(row_list)}
col_index = {c: idx for idx, c in enumerate(col_list)}
R = len(row_list)
C = len(col_list)
# Add dummy terminate columns (one per row)
aug = np.full((R, C + R), BIG_M, dtype=float)
# Initialize terminate columns to 0 (no-link baseline)
aug[:, C:] = 0.0
for r_orig, c_orig, eff in candidates:
r = row_index[r_orig]
c = col_index[c_orig]
if eff < aug[r, c]:
aug[r, c] = eff
# Solve on compact augmented matrix (only involved rows/cols)
row_ind, col_ind = linear_sum_assignment(aug)
# Build end->start stitching map (tracklet-level)
outgoing: Dict[int, int] = {}
incoming: Dict[int, int] = {}
n_links = 0
for r, c in zip(row_ind, col_ind):
if c < C:
eff = aug[r, c]
if not np.isfinite(eff) or eff >= BIG_M * 0.5:
continue
# Only accept links that actually beat baseline (negative eff_cost)
if eff >= 0:
continue
tid_end = int(end_tids[row_list[r]])
tid_start = int(start_tids[col_list[c]])
# Safety: enforce time ordering
if int(start_frames[col_list[c]]) <= int(end_frames[row_list[r]]):
continue
# Enforce one-in/one-out (Hungarian gives this, but keep safe)
if tid_end in outgoing:
continue
if tid_start in incoming:
continue
outgoing[tid_end] = tid_start
incoming[tid_start] = tid_end
n_links += 1
if verbose:
print(f" [GREEDY-P2] Gap closing: candidate ends={n_ends}, starts={n_starts}, stitched_links={n_links}")
if not outgoing:
return tracks_df
# Rebuild merged tracks by concatenating stitched tracklets
tracklets = {
int(tid): grp.sort_values("frame").reset_index(drop=True)
for tid, grp in tracks_sorted.groupby("track_id")
}
new_tracks: List[pd.DataFrame] = []
visited: Set[int] = set()
new_id = 0
def walk_chain(start_tid: int) -> None:
nonlocal new_id
chain: List[int] = []
cur = start_tid
while cur is not None and cur not in visited:
visited.add(cur)
chain.append(cur)
cur = outgoing.get(cur, None)
if not chain:
return
# Concatenate tracklets in temporal order
merged = pd.concat([tracklets[tid] for tid in chain], ignore_index=True)
merged = merged.sort_values("frame").reset_index(drop=True)
merged = merged.copy()
merged["track_id"] = new_id
new_tracks.append(merged)
new_id += 1
# Start with chain heads (no incoming links)
for tid in tracklets.keys():
if tid in incoming:
continue
walk_chain(tid)
# Handle any remaining tracklets (cycles or isolated)
for tid in tracklets.keys():
if tid not in visited:
walk_chain(tid)
if not new_tracks:
return tracks_df
out = pd.concat(new_tracks, ignore_index=True)
out = out.sort_values(["track_id", "frame"]).reset_index(drop=True)
return out
def _greedy_track_single_pass(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool,
verbose: bool = False
) -> pd.DataFrame:
"""Single pass of greedy tracking."""
frames = sorted(detections_by_frame.keys())
if len(frames) == 0:
return pd.DataFrame()
# Estimate complexity
total_dets = sum(len(df) for df in detections_by_frame.values())
avg_dets_per_frame = total_dets / len(frames) if frames else 0
if verbose:
print(f" Total detections: {total_dets}")
print(f" Frames: {len(frames)}")
print(f" Avg detections/frame: {avg_dets_per_frame:.1f}")
# Complexity estimate
max_active = int(avg_dets_per_frame * (self.max_gap_frames + 1))
links_per_frame = avg_dets_per_frame * max_active
total_links = links_per_frame * len(frames)
print(f" Estimated link evaluations: {total_links:.0f}")
if total_links > 10_000_000:
print(f" WARNING: This may take a while (estimated {total_links/1e6:.1f}M link evaluations)")
all_results = []
next_track_id = 0
# Active tracks: list of (track_id, last_frame, last_detection)
active_tracks: List[Tuple[int, int, pd.Series]] = []
# Initialize progress bar
pbar = None
if TQDM_AVAILABLE and verbose:
pbar = tqdm(total=len(frames), desc=" Greedy tracking", unit="frame",
bar_format='{l_bar}{bar}| {n_fmt}/{total_fmt} [{elapsed}<{remaining}, {postfix}]')
for frame in frames:
df = detections_by_frame[frame]
if len(df) == 0:
continue
# Preserve original detection id for later mapping (e.g., displacement extraction)
if 'det_id' not in df.columns:
df = df.reset_index(drop=False).rename(columns={'index': 'det_id'})
else:
df = df.reset_index(drop=True)
# Try to link each detection to existing tracks
linked_det_indices = set()
linked_track_indices = set()
# Compute all possible links
possible_links = []
for det_idx, det in df.iterrows():
for track_idx, (track_id, last_frame, last_det) in enumerate(active_tracks):
frame_gap = frame - last_frame
# Allow gaps up to max_gap_frames
# frame_gap=1 -> consecutive frames (0 skipped)
# frame_gap=2 -> 1 frame skipped (allowed if max_gap_frames >= 1)
# frame_gap=3 -> 2 frames skipped (allowed if max_gap_frames >= 2)
if frame_gap > self.max_gap_frames + 1:
continue
cost = self._compute_link_cost(last_det, det, is_3d, frame_gap)
if cost < np.inf:
possible_links.append((cost, det_idx, track_idx))
# Sort by cost and assign greedily
possible_links.sort(key=lambda x: x[0])
for cost, det_idx, track_idx in possible_links:
if det_idx in linked_det_indices or track_idx in linked_track_indices:
continue
# Link detection to track
track_id = active_tracks[track_idx][0]
det = df.iloc[det_idx]
result = det.to_dict()
result['track_id'] = track_id
result['frame'] = frame
all_results.append(result)
# Update active track
active_tracks[track_idx] = (track_id, frame, det)
linked_det_indices.add(det_idx)
linked_track_indices.add(track_idx)
# Start new tracks for unlinked detections
for det_idx, det in df.iterrows():
if det_idx not in linked_det_indices:
track_id = next_track_id
next_track_id += 1
result = det.to_dict()
result['track_id'] = track_id
result['frame'] = frame
all_results.append(result)
active_tracks.append((track_id, frame, det))
# Remove dead tracks (keep for max_gap_frames + 1 to allow proper gap closing)
active_tracks = [
t for t in active_tracks
if frame - t[1] <= self.max_gap_frames + 1
]
# Update progress bar if available
if TQDM_AVAILABLE and pbar is not None:
pbar.update(1)
pbar.set_postfix({'active_tracks': len(active_tracks), 'assignments': len(all_results)})
# Close progress bar
if TQDM_AVAILABLE and pbar is not None:
pbar.close()
if not all_results:
return pd.DataFrame()
return pd.DataFrame(all_results)
@staticmethod
def _process_single_frame_pair(
t_curr: int,
t_next: int,
curr_df: pd.DataFrame,
next_df: pd.DataFrame,
tracker_params: dict,
is_3d: bool,
return_pairs: bool = False,
use_simple: bool = False,
) -> Tuple[int, int, List[Tuple[int, int, float, float, float]], List[Tuple[pd.Series, pd.Series]]]:
"""
Process a single frame pair for displacement extraction.
"""
# Reset index to ensure we can use integer indexing
curr_df = curr_df.reset_index(drop=True)
next_df = next_df.reset_index(drop=True)
matches: List[Tuple[int, int, float, float, float]] = []
det_pairs: List[Tuple[pd.Series, pd.Series]] = []
# Extract parameters
pixel_size_xy = tracker_params["pixel_size_xy"]
pixel_size_z = tracker_params.get("pixel_size_z")
max_displacement_um = float(tracker_params["max_displacement_um"])
cost_distance = float(tracker_params["cost_distance"])
cost_intensity_diff = float(tracker_params["cost_intensity_diff"])
cost_size_diff = float(tracker_params["cost_size_diff"])
if "use_simple_displacement_extraction" in tracker_params:
use_simple = bool(tracker_params["use_simple_displacement_extraction"])
# Precompute coordinate arrays in um
has_z = bool(is_3d and pixel_size_z is not None and "z" in curr_df.columns and "z" in next_df.columns)
if has_z:
curr_coords = np.stack(
[
curr_df["x"].to_numpy() * pixel_size_xy,
curr_df["y"].to_numpy() * pixel_size_xy,
curr_df["z"].to_numpy() * pixel_size_z,
],
axis=1,
) if len(curr_df) else np.zeros((0, 3))
next_coords = np.stack(
[
next_df["x"].to_numpy() * pixel_size_xy,
next_df["y"].to_numpy() * pixel_size_xy,
next_df["z"].to_numpy() * pixel_size_z,
],
axis=1,
) if len(next_df) else np.zeros((0, 3))
else:
curr_coords = np.stack(
[
curr_df["x"].to_numpy() * pixel_size_xy,
curr_df["y"].to_numpy() * pixel_size_xy,
],
axis=1,
) if len(curr_df) else np.zeros((0, 2))
next_coords = np.stack(
[
next_df["x"].to_numpy() * pixel_size_xy,
next_df["y"].to_numpy() * pixel_size_xy,
],
axis=1,
) if len(next_df) else np.zeros((0, 2))
use_kdtree = SCIPY_AVAILABLE and len(curr_coords) > 0 and len(next_coords) > 0
tree = cKDTree(curr_coords) if use_kdtree else None
# Optional extra costs (only if columns exist)
has_signal = "signal" in curr_df.columns and "signal" in next_df.columns
has_size = "size" in curr_df.columns and "size" in next_df.columns
curr_signal = curr_df["signal"].to_numpy() if has_signal else None
curr_size = curr_df["size"].to_numpy() if has_size else None
if use_simple:
# Simple greedy matching (kept as you had it)
for next_idx in range(len(next_df)):
next_det = next_df.iloc[next_idx]
if use_kdtree:
candidate_idxs = tree.query_ball_point(next_coords[next_idx], r=max_displacement_um)
else:
candidate_idxs = np.arange(len(curr_df))
if len(candidate_idxs) == 0:
continue
candidate_idxs = np.asarray(candidate_idxs, dtype=int)
disp_vectors = next_coords[next_idx] - curr_coords[candidate_idxs]
distances = np.linalg.norm(disp_vectors, axis=1)
valid_mask = np.isfinite(distances) & (distances <= max_displacement_um)
if not np.any(valid_mask):
continue
disp_vectors = disp_vectors[valid_mask]
distances = distances[valid_mask]
candidate_idxs = candidate_idxs[valid_mask]
costs = distances * cost_distance
if has_signal:
next_signal = next_det["signal"]
costs += np.abs(curr_signal[candidate_idxs] - next_signal) * cost_intensity_diff
if has_size:
next_size = next_det["size"]
costs += np.abs(curr_size[candidate_idxs] - next_size) * cost_size_diff
best_rel_idx = int(np.argmin(costs))
best_cost = float(costs[best_rel_idx])
if not np.isfinite(best_cost):
continue
best_curr_idx = int(candidate_idxs[best_rel_idx])
disp_x, disp_y = float(disp_vectors[best_rel_idx][0]), float(disp_vectors[best_rel_idx][1])
disp_z = float(disp_vectors[best_rel_idx][2]) if has_z else 0.0
matches.append((best_curr_idx, next_idx, disp_x, disp_y, disp_z))
if return_pairs:
det_pairs.append((curr_df.iloc[best_curr_idx], next_det))
return (t_curr, t_next, matches, det_pairs if return_pairs else [])
# LAP-based matching
try:
from scipy.optimize import linear_sum_assignment
except ImportError:
return MinCostFlowTracker._process_single_frame_pair(
t_curr, t_next, curr_df, next_df, tracker_params, is_3d, return_pairs, use_simple=True
)
n_curr = len(curr_df)
n_next = len(next_df)
if n_curr == 0 or n_next == 0:
return (t_curr, t_next, matches, det_pairs if return_pairs else [])
BIG_M = 1e9
# Dummy cost: must be > any "good" gated link cost, so matching is preferred when possible.
# Gated links have distance <= max_displacement_um, so base distance cost <= max_displacement_um * cost_distance.
# Add a small margin plus some allowance for signal/size terms.
DUMMY_COST = (max_displacement_um * cost_distance) + (cost_intensity_diff + cost_size_diff) + 1e-3
cost_matrix = np.full((n_next, n_curr), BIG_M, dtype=float)
any_candidate = False
for next_idx in range(n_next):
if use_kdtree:
candidate_idxs = tree.query_ball_point(next_coords[next_idx], r=max_displacement_um)
else:
disp_all = next_coords[next_idx] - curr_coords
dist_all = np.linalg.norm(disp_all, axis=1)
candidate_idxs = np.where((np.isfinite(dist_all)) & (dist_all <= max_displacement_um))[0]
if len(candidate_idxs) == 0:
continue
candidate_idxs = np.asarray(candidate_idxs, dtype=int)
next_det = next_df.iloc[next_idx]
disp_vectors = next_coords[next_idx] - curr_coords[candidate_idxs]
distances = np.linalg.norm(disp_vectors, axis=1)
valid_mask = np.isfinite(distances) & (distances <= max_displacement_um)
if not np.any(valid_mask):
continue
candidate_idxs = candidate_idxs[valid_mask]
distances = distances[valid_mask]
costs = distances * cost_distance
if has_signal:
next_signal = next_det["signal"]
costs += np.abs(curr_signal[candidate_idxs] - next_signal) * cost_intensity_diff
if has_size:
next_size = next_det["size"]
costs += np.abs(curr_size[candidate_idxs] - next_size) * cost_size_diff
valid2 = np.isfinite(costs)
if not np.any(valid2):
continue
any_candidate = True
cost_matrix[next_idx, candidate_idxs[valid2]] = costs[valid2]
# If nothing was gated-in, bail early (no displacement to extract)
if not any_candidate:
return (t_curr, t_next, matches, det_pairs if return_pairs else [])
# Augment with dummy columns to allow unmatched next detections, but NOT for free.
aug_cost = np.full((n_next, n_curr + n_next), DUMMY_COST, dtype=float)
aug_cost[:, :n_curr] = cost_matrix
row_ind, col_ind = linear_sum_assignment(aug_cost)
for r, c in zip(row_ind, col_ind):
if c < n_curr:
if aug_cost[r, c] >= BIG_M * 0.5:
continue
next_idx = int(r)
curr_idx = int(c)
disp_vector = next_coords[next_idx] - curr_coords[curr_idx]
disp_x, disp_y = float(disp_vector[0]), float(disp_vector[1])
disp_z = float(disp_vector[2]) if has_z else 0.0
matches.append((curr_idx, next_idx, disp_x, disp_y, disp_z))
if return_pairs:
det_pairs.append((curr_df.iloc[curr_idx], next_df.iloc[next_idx]))
return (t_curr, t_next, matches, det_pairs if return_pairs else [])
def _extract_displacement_parallel(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool,
verbose: bool = False
) -> Dict[int, pd.DataFrame]:
"""
Extract displacement vectors from consecutive frame pairs in parallel.
Parameters
----------
detections_by_frame : dict
Dict mapping frame number to detection DataFrame
is_3d : bool
Whether tracking is 3D
verbose : bool
Print progress messages
Returns
-------
dict
Dict with same structure, but detections have prev_disp_x_um, prev_disp_y_um, prev_disp_z_um
"""
# Get sorted frame numbers
frame_nums = sorted(detections_by_frame.keys())
if len(frame_nums) < 2:
return detections_by_frame
# Initialize displacement storage
# Key: (frame_t, det_index) -> (disp_x, disp_y, disp_z)
displacement_map = {}
if verbose:
print(f" [Stage 1] Extracting displacement from {len(frame_nums)-1} frame pairs")
# Create frame pairs: [(t0, t1), (t1, t2), ...]
frame_pairs = [(frame_nums[i], frame_nums[i+1]) for i in range(len(frame_nums)-1)]
tracker_params = {
'pixel_size_xy': self.pixel_size_xy,
'pixel_size_z': self.pixel_size_z if is_3d else None,
'max_displacement_um': self.max_displacement_um,
'cost_distance': self.cost_distance,
'cost_intensity_diff': self.cost_intensity_diff,
'cost_size_diff': self.cost_size_diff
}
# Execute frame-pair matching sequentially (channel-level parallelization is upstream)
for t_curr, t_next in frame_pairs:
_, _, matches, _ = self._process_single_frame_pair(
t_curr,
t_next,
detections_by_frame[t_curr],
detections_by_frame[t_next],
tracker_params,
is_3d,
return_pairs=False,
use_simple=self.use_simple_displacement_extraction
)
# Store displacements for frame t_next
for _, next_idx, disp_x, disp_y, disp_z in matches:
# Use DataFrame original index (before reset_index)
next_original_idx = detections_by_frame[t_next].index[next_idx]
displacement_map[(t_next, next_original_idx)] = (disp_x, disp_y, disp_z)
# Apply displacements to detection DataFrames
detections_with_disp = {}
for frame, df in detections_by_frame.items():
df_copy = df.copy()
# Initialize displacement columns
df_copy['prev_disp_x_um'] = np.nan
df_copy['prev_disp_y_um'] = np.nan
if is_3d:
df_copy['prev_disp_z_um'] = np.nan
# Apply stored displacements
for idx in df_copy.index:
key = (frame, idx)
if key in displacement_map:
disp_x, disp_y, disp_z = displacement_map[key]
df_copy.at[idx, 'prev_disp_x_um'] = disp_x
df_copy.at[idx, 'prev_disp_y_um'] = disp_y
if is_3d:
df_copy.at[idx, 'prev_disp_z_um'] = disp_z
detections_with_disp[frame] = df_copy
if verbose:
total_displacements = len(displacement_map)
print(f" [Stage 1] Extracted {total_displacements} displacement vectors")
return detections_with_disp
def _min_cost_flow_track(
self,
detections_by_frame: Dict[int, pd.DataFrame],
is_3d: bool,
verbose: bool = False,
mcf_max_neighbors: Optional[int] = None,
mcf_slack_cost_component: Optional[float] = None,
output_dir: Optional[Path] = None,
debug_dir: Optional[Path] = None,
track_type: str = "",
type_id: int = 0,
debug_images: bool = False
) -> pd.DataFrame:
"""
Min-cost flow tracking using Google OR-Tools.
Constructs a flow network where:
- Source connects to all detections (birth cost)
- Detections connect to detections in next frames (link cost)
- All detections connect to sink (death cost)
- Bypass arc allows unused flow to go source->sink at cost 0
Flow budget is set to total number of detections. The solver only creates
tracks when their total cost (birth + transits + links + death) is negative,
beating the bypass baseline of 0.
Parameters
----------
detections_by_frame : dict
Detections organized by frame
is_3d : bool
Whether data is 3D
verbose : bool
Print progress information
mcf_max_neighbors : int or None
Maximum neighbors to consider for linking (k-NN constraint)
mcf_slack_cost_component : float or None
Deprecated (bypass cost is now always 0)
output_dir : Path or None
Directory for debug output
debug_dir : Path or None
Override directory for debug output
track_type : str
Label for debug output
type_id : int
ID for debug output
debug_images : bool
Generate debug visualizations
"""
# Check if OR-Tools is available
if not ORTOOLS_AVAILABLE:
print(f" ERROR: OR-Tools is not installed!")
print(f" Min-cost flow tracking requires OR-Tools.")
print(f" Install with: pip install ortools")
print(f" Returning empty result (will trigger greedy fallback)")
return pd.DataFrame()
frames = sorted(detections_by_frame.keys())
if len(frames) == 0:
return pd.DataFrame()
debug_output_dir = debug_dir if debug_dir is not None else output_dir
has_debug_output = debug_output_dir is not None
# Build node mapping: each detection gets in/out nodes to allow merges/splits
# Node 0: source, Node 1: sink, Nodes 2+: detection_in/detection_out pairs
source_node = 0
sink_node = 1
next_node_id = 2
detection_to_nodes = {} # (frame, det_idx) -> (node_in, node_out)
node_in_to_detection = {}
node_out_to_detection = {}
for frame in frames:
df = detections_by_frame[frame].reset_index(drop=True)
for det_idx in range(len(df)):
node_in = next_node_id
node_out = next_node_id + 1
next_node_id += 2
detection_to_nodes[(frame, det_idx)] = (node_in, node_out)
node_in_to_detection[node_in] = (frame, det_idx)
node_out_to_detection[node_out] = (frame, det_idx)
total_nodes = next_node_id
if verbose:
print(f" Building flow network: {total_nodes} nodes ({len(detection_to_nodes)} detections)")
# Initialize min-cost flow solver
smcf = min_cost_flow.SimpleMinCostFlow()
# Scale costs to integers (OR-Tools requires integer costs)
COST_SCALE = 1000
# Add edges
edge_count = 0
link_count = 0
# 1. Source -> detection_in (birth cost), capacity 1 to prevent double births
for node_in in node_in_to_detection.keys():
smcf.add_arc_with_capacity_and_unit_cost(
source_node, node_in,
capacity=1,
unit_cost=int(self.mcf_birth_cost * COST_SCALE)
)
edge_count += 1
# 2. detection_in -> detection_out (allow two flows; second pays merge/split penalty)
# Apply quality bonus (negative cost) based on p_good to encourage including quality detections
quality_bonus_stats = {
'with_p_good': 0,
'applied': 0,
'total_bonus': 0.0,
'sum_p_good': 0.0,
'min_p_good': None,
'max_p_good': None,
}
for node_in, (frame, det_idx) in node_in_to_detection.items():
node_out = detection_to_nodes[(frame, det_idx)][1]
# Compute quality bonus based on p_good (0..1)
quality_bonus = 0.0
det = detections_by_frame[frame].iloc[det_idx]
p_good = None
if 'p_good' in det.index and pd.notna(det['p_good']):
p_good = float(det['p_good'])
else:
p_vals = []
if 'ch1_p_good' in det.index and pd.notna(det['ch1_p_good']):
p_vals.append(float(det['ch1_p_good']))
if 'ch2_p_good' in det.index and pd.notna(det['ch2_p_good']):
p_vals.append(float(det['ch2_p_good']))
if p_vals:
p_good = float(min(p_vals))
if p_good is not None:
p_good = float(np.clip(p_good, 0.0, 1.0))
quality_bonus_stats['with_p_good'] += 1
quality_bonus_stats['sum_p_good'] += p_good
if quality_bonus_stats['min_p_good'] is None:
quality_bonus_stats['min_p_good'] = p_good
quality_bonus_stats['max_p_good'] = p_good
else:
quality_bonus_stats['min_p_good'] = min(quality_bonus_stats['min_p_good'], p_good)
quality_bonus_stats['max_p_good'] = max(quality_bonus_stats['max_p_good'], p_good)
quality_bonus = -1.0 * self.mcf_quality_bonus_weight * p_good
if p_good > 0 and self.mcf_quality_bonus_weight > 0:
quality_bonus_stats['applied'] += 1
quality_bonus_stats['total_bonus'] += abs(quality_bonus)
# Single capacity transit arc with quality bonus + flat detection bonus (both negative costs = bonuses)
transit_cost = quality_bonus - self.mcf_detection_bonus
smcf.add_arc_with_capacity_and_unit_cost(
node_in, node_out,
capacity=1,
unit_cost=int(transit_cost * COST_SCALE)
)
edge_count += 1
# 3. detection_out -> sink (death cost), capacity 1 to prevent double deaths
for node_out in node_out_to_detection.keys():
smcf.add_arc_with_capacity_and_unit_cost(
node_out, sink_node,
capacity=1,
unit_cost=int(self.mcf_death_cost * COST_SCALE)
)
edge_count += 1
# 4. detection_out -> detection_in links (across frames)
for frame in frames:
df_curr = detections_by_frame[frame].reset_index(drop=True)
# Look ahead up to mcf_max_gap_frames + 1
for gap in range(1, self.mcf_max_gap_frames + 2):
next_frame = frame + gap
if next_frame not in detections_by_frame:
continue
df_next = detections_by_frame[next_frame].reset_index(drop=True)
# Build k-NN structure for next frame if limiting neighbors
max_k = mcf_max_neighbors if mcf_max_neighbors is not None else self.mcf_max_neighbors
use_knn = SCIPY_AVAILABLE and max_k is not None and max_k > 0
next_coords = None
next_tree = None
if use_knn:
if is_3d and self.pixel_size_z is not None and 'z' in df_next.columns:
next_coords = np.stack([
df_next['x'].to_numpy() * self.pixel_size_xy,
df_next['y'].to_numpy() * self.pixel_size_xy,
df_next['z'].to_numpy() * self.pixel_size_z
], axis=1)
else:
next_coords = np.stack([
df_next['x'].to_numpy() * self.pixel_size_xy,
df_next['y'].to_numpy() * self.pixel_size_xy
], axis=1)
next_tree = cKDTree(next_coords)
# Compute all possible links from current to next frame
for curr_idx, curr_det in df_curr.iterrows():
curr_node_out = detection_to_nodes[(frame, curr_idx)][1]
candidate_indices = range(len(df_next))
if use_knn:
if is_3d and self.pixel_size_z is not None and 'z' in df_curr.columns:
curr_coord = [
curr_det['x'] * self.pixel_size_xy,
curr_det['y'] * self.pixel_size_xy,
curr_det['z'] * self.pixel_size_z
]
else:
curr_coord = [
curr_det['x'] * self.pixel_size_xy,
curr_det['y'] * self.pixel_size_xy
]
dists, idxs = next_tree.query(curr_coord, k=min(max_k, len(df_next)))
if np.isscalar(idxs):
idxs = [int(idxs)]
else:
idxs = [int(i) for i in np.atleast_1d(idxs)]
candidate_indices = idxs
for next_idx in candidate_indices:
next_det = df_next.iloc[next_idx]
next_node_in = detection_to_nodes[(next_frame, next_idx)][0]
# Extract previous displacement from curr_det (if available)
prev_displacement = None
if 'prev_disp_x_um' in curr_det.index and pd.notna(curr_det.get('prev_disp_x_um')):
if is_3d and self.pixel_size_z is not None and 'prev_disp_z_um' in curr_det.index:
prev_displacement = np.array([
curr_det['prev_disp_x_um'],
curr_det['prev_disp_y_um'],
curr_det['prev_disp_z_um']
])
else:
prev_displacement = np.array([
curr_det['prev_disp_x_um'],
curr_det['prev_disp_y_um']
])
# Compute link cost using MCF-specific parameters (with motion smoothness if available)
cost = self._compute_link_cost(
curr_det, next_det, is_3d, gap,
use_mcf_costs=True,
prev_displacement=prev_displacement
)
if cost < np.inf:
smcf.add_arc_with_capacity_and_unit_cost(
curr_node_out, next_node_in,
capacity=1,
unit_cost=int(cost * COST_SCALE)
)
edge_count += 1
link_count += 1
# 5. Bypass arc: allow unused flow budget to go directly source->sink at cost 0
# This means tracks are only created when they beat the "do nothing" baseline
num_detections = len(detection_to_nodes)
bypass_cost = 0.0 # True bypass: no cost for unused flow
smcf.add_arc_with_capacity_and_unit_cost(
source_node, sink_node,
capacity=num_detections, # Upper bound on flow
unit_cost=int(bypass_cost * COST_SCALE)
)
edge_count += 1
# Always print link statistics
print(f" [MCF] Added {edge_count} total edges ({link_count} detection links, {edge_count - link_count} birth/death)")
print(f" [MCF] Birth/death cost: {self.mcf_birth_cost}/{self.mcf_death_cost}, max_displacement: {self.mcf_max_displacement_um} µm")
# Print quality bonus statistics
if quality_bonus_stats['with_p_good'] > 0:
mean_p_good = quality_bonus_stats['sum_p_good'] / quality_bonus_stats['with_p_good']
print(
f" [MCF] Quality bonus: weight={self.mcf_quality_bonus_weight:.2f}, "
f"applied to {quality_bonus_stats['applied']}/{quality_bonus_stats['with_p_good']} detections"
)
if self.verbosity_level >= 2:
print(
f" [MCF] p_good mean={mean_p_good:.3f}, "
f"min={quality_bonus_stats['min_p_good']:.3f}, max={quality_bonus_stats['max_p_good']:.3f}"
)
print(f" [MCF] Total bonus applied: {quality_bonus_stats['total_bonus']:.2f}")
elif self.mcf_quality_bonus_weight > 0:
print(f" [MCF] Quality bonus: weight={self.mcf_quality_bonus_weight:.2f}, but no p_good data found in detections")
# Print flat detection bonus
if self.mcf_detection_bonus > 0:
total_flat_bonus = self.mcf_detection_bonus * len(node_in_to_detection)
print(f" [MCF] Flat detection bonus: {self.mcf_detection_bonus:.2f} per detection, total={total_flat_bonus:.2f} across {len(node_in_to_detection)} detections")
# Flow budget: set to total number of detections (upper bound)
# The solver will only create tracks when they beat the bypass cost (0)
flow_budget = len(detection_to_nodes)
print(f" [MCF] Flow budget: {flow_budget} (detections: {flow_budget})")
print(f" [MCF] Bypass arc: cost={bypass_cost}, capacity={flow_budget}")
print(f" [MCF] Only tracks with negative total cost will be created (beat bypass)")
# Set supply/demand
# Source: supply = flow budget (upper bound on all flow)
# Sink: demand = -flow budget (must consume all supply)
# All detection nodes: supply/demand = 0 (flow through)
smcf.set_node_supply(source_node, flow_budget)
smcf.set_node_supply(sink_node, -flow_budget)
# Solve min-cost flow
if verbose:
print(f" Solving min-cost flow...")
try:
status = smcf.solve()
except Exception as e:
print(f" ERROR: Min-cost flow solver crashed: {e}")
import traceback
traceback.print_exc()
return pd.DataFrame()
# Status codes: OPTIMAL=0, NOT_SOLVED=1, FEASIBLE=2, INFEASIBLE=3, UNBALANCED=4, BAD_RESULT=5, BAD_COST_RANGE=6
status_names = {0: 'OPTIMAL', 1: 'NOT_SOLVED', 2: 'FEASIBLE', 3: 'INFEASIBLE', 4: 'UNBALANCED', 5: 'BAD_RESULT', 6: 'BAD_COST_RANGE'}
status_name = status_names.get(status, f'UNKNOWN({status})')
if status != smcf.OPTIMAL:
print(f" WARNING: Min-cost flow solver returned status: {status_name}")
print(f" This usually indicates:")
if status == 3: # INFEASIBLE
print(f" - Network is INFEASIBLE (no valid solution exists)")
print(f" - Check that supply ({flow_budget}) = demand (-{flow_budget})")
elif status == 4: # UNBALANCED
print(f" - Network is UNBALANCED (supply != demand)")
print(f" - Supply at source: {flow_budget}, Demand at sink: {flow_budget}")
elif status == 6: # BAD_COST_RANGE
print(f" - BAD_COST_RANGE (costs too large for integer solver)")
print(f" - Cost scale: {COST_SCALE}, birth/death costs: {self.birth_cost}/{self.death_cost}")
print(f" Returning empty result (will trigger greedy fallback in caller)")
return pd.DataFrame()
if verbose:
print(f" Optimal solution found! Total cost: {smcf.optimal_cost() / COST_SCALE:.2f}")
# Count bypass flow (source -> sink)
bypass_flow = 0
for i in range(smcf.num_arcs()):
if smcf.flow(i) > 0 and smcf.tail(i) == source_node and smcf.head(i) == sink_node:
bypass_flow = smcf.flow(i)
break
used_flow = flow_budget - bypass_flow
print(f" [MCF] Flow allocation: {used_flow}/{flow_budget} used for tracks, {bypass_flow} through bypass ({100*bypass_flow/flow_budget:.1f}%)")
# Extract tracks from flow solution
# Build adjacency from flow: detection -> [next_detection, ...]
# Use defaultdict to handle splits (one detection -> multiple next detections)
from collections import defaultdict
flow_graph = defaultdict(list) # (frame, idx) -> [(next_frame, next_idx), ...]
link_distances = []
# Death diagnostics: distance to nearest detection in next frame
death_to_nextframe_dists: List[float] = []
death_debug_records: List[dict] = []
# Cost component tracking for debug visualizations
cost_breakdowns = [] # List of {distance, distance_cost, total_cost, ambiguity_factor, smoothness_penalty}
ambiguity_factors = [] # Per-detection ambiguity factors
quality_bonuses = [] # Per-detection quality bonuses
smoothness_penalties = [] # Per-link smoothness penalties
for i in range(smcf.num_arcs()):
if smcf.flow(i) > 0:
tail = smcf.tail(i)
head = smcf.head(i)
if tail in node_out_to_detection and head in node_in_to_detection:
curr_det = node_out_to_detection[tail]
next_det = node_in_to_detection[head]
flow_graph[curr_det].append(next_det)
# Get detection data
f1, i1 = curr_det
f2, i2 = next_det
d1 = detections_by_frame[f1].iloc[i1]
d2 = detections_by_frame[f2].iloc[i2]
gap = f2 - f1
# Compute distance for diagnostics
if is_3d and self.pixel_size_z is not None and 'z' in d1 and 'z' in d2:
dist = np.sqrt(
((d1['x'] - d2['x']) * self.pixel_size_xy) ** 2 +
((d1['y'] - d2['y']) * self.pixel_size_xy) ** 2 +
((d1['z'] - d2['z']) * self.pixel_size_z) ** 2
)
else:
dist = np.sqrt(
((d1['x'] - d2['x']) * self.pixel_size_xy) ** 2 +
((d1['y'] - d2['y']) * self.pixel_size_xy) ** 2
)
link_distances.append(dist)
# Collect cost breakdown for debug visualizations
if debug_images:
# Extract previous displacement from d1 (if available)
prev_displacement = None
if 'prev_disp_x_um' in d1.index and pd.notna(d1.get('prev_disp_x_um')):
if is_3d and self.pixel_size_z is not None and 'prev_disp_z_um' in d1.index:
prev_displacement = np.array([
d1['prev_disp_x_um'],
d1['prev_disp_y_um'],
d1['prev_disp_z_um']
])
else:
prev_displacement = np.array([
d1['prev_disp_x_um'],
d1['prev_disp_y_um']
])
# Compute cost with breakdown
_, breakdown = self._compute_link_cost(
d1, d2, is_3d, gap,
use_mcf_costs=True,
prev_displacement=prev_displacement,
return_breakdown=True
)
cost_breakdowns.append(breakdown)
# Collect smoothness penalties
if breakdown['smoothness_penalty'] > 0:
smoothness_penalties.append(breakdown['smoothness_penalty'])
# Identify track starts: source -> detection_in with flow > 0
start_detections = []
for i in range(smcf.num_arcs()):
if smcf.flow(i) > 0 and smcf.tail(i) == source_node and smcf.head(i) in node_in_to_detection:
start_detections.append(node_in_to_detection[smcf.head(i)])
# Analyze split/merge events and double-used detections
# Count flow through each detection's in/out nodes
flow_in = {} # detection -> total flow entering detection_in
flow_out = {} # detection -> total flow leaving detection_out
transit_flow = {} # detection -> total flow through detection_in->detection_out
for i in range(smcf.num_arcs()):
if smcf.flow(i) > 0:
tail = smcf.tail(i)
head = smcf.head(i)
flow = smcf.flow(i)
# Count flow entering detection_in nodes (from source or from previous detections)
if head in node_in_to_detection:
det = node_in_to_detection[head]
flow_in[det] = flow_in.get(det, 0) + flow
# Count flow leaving detection_out nodes (to sink or to next detections)
if tail in node_out_to_detection:
det = node_out_to_detection[tail]
flow_out[det] = flow_out.get(det, 0) + flow
# Count flow through transit edges (detection_in -> detection_out)
if tail in node_in_to_detection and head in node_out_to_detection:
det_tail = node_in_to_detection[tail]
det_head = node_out_to_detection[head]
if det_tail == det_head: # Same detection
transit_flow[det_tail] = transit_flow.get(det_tail, 0) + flow
# Note: Split/merge detection is now done post-hoc on the output tracks
# (see detect_split_merge_events function)
# Identify deaths: detections that participated in flow but have no outgoing link
# Flatten adjacency lists to get all destination detections
all_destinations = set()
for next_dets in flow_graph.values():
all_destinations.update(next_dets)
used_detections = set(flow_graph.keys()) | all_destinations | set(start_detections)
dead_detections = used_detections - set(flow_graph.keys())
for (frame, idx) in dead_detections:
# Look only at immediate next frame
next_frame = frame + 1
if next_frame not in detections_by_frame:
continue
df_next = detections_by_frame[next_frame]
if len(df_next) == 0:
continue
det = detections_by_frame[frame].iloc[idx]
# Compute distances to all detections in next frame
if is_3d and self.pixel_size_z is not None and 'z' in det and 'z' in df_next.columns:
dx = (df_next['x'].to_numpy() - det['x']) * self.pixel_size_xy
dy = (df_next['y'].to_numpy() - det['y']) * self.pixel_size_xy
dz = (df_next['z'].to_numpy() - det['z']) * self.pixel_size_z
dists = np.sqrt(dx * dx + dy * dy + dz * dz)
else:
dx = (df_next['x'].to_numpy() - det['x']) * self.pixel_size_xy
dy = (df_next['y'].to_numpy() - det['y']) * self.pixel_size_xy
dists = np.sqrt(dx * dx + dy * dy)
if len(dists) > 0 and np.isfinite(dists).any():
min_dist = float(np.nanmin(dists))
if np.isfinite(min_dist):
death_to_nextframe_dists.append(min_dist)
# Re-run gate/pruning to count link candidates and min cost
max_k = mcf_max_neighbors if mcf_max_neighbors is not None else self.mcf_max_neighbors
gate = (dists <= (self.mcf_max_displacement_um if self.mcf_max_displacement_um is not None else self.max_displacement_um))
candidate_idxs = [i for i, keep in enumerate(gate) if keep]
if max_k is not None and max_k > 0 and candidate_idxs:
candidate_idxs = sorted(candidate_idxs, key=lambda j: dists[j])[:max_k]
min_link_cost = np.inf
for j in candidate_idxs:
next_det = df_next.iloc[j]
cost = self._compute_link_cost(det, next_det, is_3d, frame_gap=1, use_mcf_costs=True)
if np.isfinite(cost) and cost < min_link_cost:
min_link_cost = cost
death_debug_records.append({
'frame': frame,
'candidates': len(candidate_idxs),
'min_dist_um': float(np.nanmin(dists)) if len(dists) else np.nan,
'min_link_cost': float(min_link_cost) if np.isfinite(min_link_cost) else np.inf
})
# Build tracks by following flow from each starting detection
# NOTE: With splits/merges, this creates one track per "branch"
# A detection with 2 outgoing links will spawn 2 separate tracks
all_results = []
next_track_id = 0
visited = set()
track_costs = {} # track_id -> total cost (not computed with adjacency list approach)
# For now, we'll extract simple linear tracks by taking first branch at splits
# A proper split/merge-aware extraction would require more complex logic
for start_det in start_detections:
if start_det in visited:
continue
track_id = next_track_id
next_track_id += 1
curr_det = start_det
while curr_det is not None:
visited.add(curr_det)
frame, det_idx = curr_det
det = detections_by_frame[frame].iloc[det_idx]
result = det.to_dict()
result['track_id'] = track_id
result['frame'] = frame
all_results.append(result)
# Collect ambiguity and quality data for debug visualizations
if debug_images:
# Collect ambiguity factors (gap 1 is most common)
if 'ambiguity_gap1' in det.index and pd.notna(det.get('ambiguity_gap1')):
ambiguity_factors.append(float(det['ambiguity_gap1']))
# Collect quality bonuses (same logic as main MCF graph)
p_good = None
if 'p_good' in det.index and pd.notna(det['p_good']):
p_good = float(det['p_good'])
else:
p_vals = []
if 'ch1_p_good' in det.index and pd.notna(det['ch1_p_good']):
p_vals.append(float(det['ch1_p_good']))
if 'ch2_p_good' in det.index and pd.notna(det['ch2_p_good']):
p_vals.append(float(det['ch2_p_good']))
if p_vals:
p_good = float(min(p_vals))
if p_good is not None and p_good > 0:
p_good = float(np.clip(p_good, 0.0, 1.0))
quality_bonus = -1.0 * self.mcf_quality_bonus_weight * p_good
quality_bonuses.append(quality_bonus)
# Get next detection(s) from adjacency list
next_dets = flow_graph.get(curr_det, [])
if len(next_dets) == 0:
# End of track
curr_det = None
elif len(next_dets) == 1:
# Linear continuation
curr_det = next_dets[0]
else:
# Split: take first unvisited branch, others will be separate tracks
curr_det = None
for next_det in next_dets:
if next_det not in visited:
curr_det = next_det
break
if verbose:
print(f" Extracted {next_track_id} tracks from flow solution")
if link_distances:
link_dists = np.array(link_distances)
print(f" Link distance stats (µm): min={link_dists.min():.3f}, mean={link_dists.mean():.3f}, max={link_dists.max():.3f}")
# Simple histogram bins
bins = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.5, 2.0]
hist, edges = np.histogram(link_dists, bins=bins)
print(" Link distance histogram (µm):")
for count, start, end in zip(hist, edges[:-1], edges[1:]):
print(f" {start:.1f}-{end:.1f}: {count}")
# Save comprehensive debug visualizations if debug_images is enabled
if debug_images and has_debug_output:
try:
save_debug_visualizations(
detections_by_frame=detections_by_frame,
all_results=all_results,
link_distances=link_distances,
start_detections=start_detections,
death_to_nextframe_dists=death_to_nextframe_dists,
track_costs=track_costs,
is_3d=is_3d,
output_dir=debug_output_dir,
track_type=track_type,
type_id=type_id,
verbose=verbose,
pixel_size_xy=self.pixel_size_xy,
max_displacement_um=self.max_displacement_um,
mcf_max_displacement_um=self.mcf_max_displacement_um,
mcf_cost_distance=self.mcf_cost_distance,
mcf_crowding_multiplier=self.mcf_crowding_multiplier,
cost_breakdowns=cost_breakdowns,
ambiguity_factors=ambiguity_factors,
quality_bonuses=quality_bonuses,
smoothness_penalties=smoothness_penalties,
debug_dir=debug_dir
)
except Exception as e:
print(f" Warning: Could not save debug visualizations: {e}")
import traceback
traceback.print_exc()
if verbose and death_debug_records:
total_deaths = len(death_debug_records)
zero_candidates = sum(1 for r in death_debug_records if r['candidates'] == 0)
finite_costs = [r['min_link_cost'] for r in death_debug_records if np.isfinite(r['min_link_cost'])]
med_cost = np.median(finite_costs) if finite_costs else None
print(f" [MCF] Death debug: deaths={total_deaths}, zero_candidates={zero_candidates}, finite_min_costs={len(finite_costs)}, median_min_cost={med_cost}")
if not all_results:
return pd.DataFrame()
return pd.DataFrame(all_results)

mincostflow.py at commit 08be6e0, under MIT · at the source

Overview

Authors: Krassimira A Garbett¹, James P Allen¹, Jaybree M Lopez¹, Cassandra M Smith¹, Richard C Sando¹

ORCID iDs: Krassimira A Garbett, James P Allen, Jaybree M Lopez, Cassandra M Smith, Richard C Sando

Department of Pharmacology, Vanderbilt Brain Institute, Vanderbilt University Nashville United States

Institutions: Vanderbilt University (United States)

Journal: eLife, volume 14, article RP107635

Dates: published online 27 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.7554/elife.107635 · PMID 42507759 · PMCID PMC13405638 · OpenAlex W4413760595

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: mouse (organism), cellular / molecular (subfield)

Methods: Statistics, Evoked potentials, fMRI & imaging, Smoothing, state filtering, decompositions

Keywords: synapse, synapse formation, E/I balance, time-lapse imaging, live imaging, Mouse

MeSH: Synapses*, Animals, Mice, Neurodevelopment (* major topic)

Journal subjects: Cell Biology, Neuroscience

Topic: Neuroscience and Neuropharmacology Research (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex

Funding: National Institute of Mental Health (R00MH117235, DP2MH140134); Alfred P. Sloan Foundation (FG-2022-19037)

Citations: not cited yet (Europe PMC); 48 references in the paper

Research resources: Goat anti-mouse Alexa Fluor 546 RRID:AB_141370, HA-Tag (C29F4) rabbit mAb RRID:AB_1549585, RRID:AB_2138173, GFP polyclonal antibody rabbit RRID:AB_221569, Polyclonal anti-GABARα1 rabbit RRID:AB_2232180, RRID:AB_2301751, Goat anti-mouse Alexa Fluor 488 RRID:AB_2534069, Goat anti-rabbit Alexa Fluor 546 RRID:AB_2534077, Goat anti-mouse Alexa Fluor 647 RRID:AB_2535805, Goat anti-rabbit Alexa Fluor 647 RRID:AB_2535813, Goat anti-chicken Alexa Fluor 647 RRID:AB_2535866, Goat anti-guinea pig Alexa Fluor 647 RRID:AB_2535867, Goat anti-rabbit Alexa Fluor 488 RRID:AB_2576217, Monoclonal anti-Gephyrin mouse RRID:AB_2619837, Polyclonal anti-Shank 2 guinea pig RRID:AB_2619861, RRID:AB_2716712, Polyclonal anti-Homer 1 rabbit RRID:AB_887730, Polyclonal anti-Syn1/2 rabbit RRID:AB_887804, pAAV SEP-GluA2 TKIT (plasmid) RRID:Addgene_169442, HEK293T RRID:CVCL_1926, RRID:IMSR_JAX:000664, GraphPad Prism 8.0, 9.0 RRID:SCR_002798, ImageJ RRID:SCR_003070, Adobe Illustrator RRID:SCR_010279, Adobe Photoshop RRID:SCR_014199, SnapGene RRID:SCR_015052, NIS-Elements AR RRID:SCR_027181

Abstract

Diverse synaptic connections self-organize into neural circuits during brain development. A balance between excitatory and inhibitory synaptic function is required for information processing by these neural circuits. Despite the importance of this balance, the interplay between excitatory and inhibitory synaptic assembly during circuit establishment remains unclear due to a lack of means to monitor both processes simultaneously. Here, we develop imaging and analysis methods to visualize and track excitatory and inhibitory synapses. By applying these approaches, we find that despite continual dynamics, excitatory and inhibitory synaptic density remain at steady-state levels during synapse maturation. These results indicate balanced excitatory and inhibitory synapse assembly, despite continual synaptic turnover.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 11 matches between paragraphs and lines of code.

datajamesallen/Garbett_2026_eLife

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 08be6e0626df4dd4c9a9869b350dce35b416808a, 10 July 2026

Languages: Python (67)

Size: 76 files, 67 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, license file, environment (requirements.txt, setup.py, ConsensusDrift/setup.py)

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (45 files), pandas (23 files), SciPy (23 files), tifffile (16 files), Matplotlib (13 files), scikit-image (12 files), scikit-learn (4 files), ImageJ / Fiji (3 files), Pillow (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

69 files

ConsensusDrift/__init__.py, Python, 13 lines
ConsensusDrift/__main__.py, Python, 4 lines
ConsensusDrift/apply.py, Python, 463 lines
ConsensusDrift/cli.py, Python, 593 lines
ConsensusDrift/debug.py, Python, 646 lines
ConsensusDrift/ome_io.py, Python, 290 lines
ConsensusDrift/registration.py, Python, 2,698 lines, 1 match
ConsensusDrift/setup.py, Python, 42 lines
ConsensusDrift/smoothing.py, Python, 38 lines
ConsensusDrift/tiff_io.py, Python, 324 lines
ConsensusDrift/transforms.py, Python, 249 lines
ConsensusDrift/utils.py, Python, 92 lines
scripts/add_track_corner_markers.py, Python, 295 lines
scripts/analyze_tracks.py, Python, 2,981 lines, 1 match
scripts/categorize_tracks.py, Python, 1,864 lines
scripts/imagej_draw_lines.py, Python, 360 lines
scripts/imagej_overlay_merge.py, Python, 238 lines
scripts/imagej_puncta_preprocess.py, Python, 988 lines, 1 match
scripts/nd2_preprocess.py, Python, 1,120 lines
scripts/render_category_figures.py, Python, 193 lines
setup.py, Python, 76 lines
unisynapse/__init__.py, Python, 41 lines
unisynapse/__main__.py, Python, 23 lines
unisynapse/analysis/__init__.py, Python, 38 lines
unisynapse/analysis/codistribution.py, Python, 579 lines
unisynapse/analysis/skeleton_pruning.py, Python, 732 lines
unisynapse/cli.py, Python, 834 lines
unisynapse/config/__init__.py, Python, 41 lines
unisynapse/config/loader.py, Python, 398 lines
unisynapse/config/schema.py, Python, 620 lines, 1 match
unisynapse/detect/__init__.py, Python, 73 lines
unisynapse/detect/blob.py, Python, 475 lines
unisynapse/detect/filtering.py, Python, 376 lines
unisynapse/detect/gmm_filter.py, Python, 862 lines
unisynapse/detect/measurement.py, Python, 554 lines
unisynapse/detect/merging.py, Python, 243 lines
unisynapse/detect/refinement.py, Python, 260 lines
unisynapse/detect/segmentation.py, Python, 665 lines, 1 match
unisynapse/detect/tail_adaptive_filter.py, Python, 1,181 lines
unisynapse/detect/temporal.py, Python, 539 lines
unisynapse/image/__init__.py, Python, 48 lines
unisynapse/image/masks.py, Python, 1,471 lines, 1 match
unisynapse/image/morphology.py, Python, 360 lines
unisynapse/image/neurite_temporal.py, Python, 80 lines
unisynapse/image/transforms.py, Python, 249 lines
unisynapse/io/__init__.py, Python, 41 lines
unisynapse/io/ome.py, Python, 290 lines
unisynapse/io/output.py, Python, 400 lines
unisynapse/io/tiff.py, Python, 324 lines
unisynapse/pair/__init__.py, Python, 22 lines
unisynapse/pair/hungarian.py, Python, 850 lines, 1 match
unisynapse/pipelines/__init__.py, Python, 17 lines
unisynapse/pipelines/codistribution.py, Python, 1,227 lines
unisynapse/pipelines/detect.py, Python, 1,876 lines
unisynapse/pipelines/filter.py, Python, 1,307 lines
unisynapse/pipelines/neurite.py, Python, 1,048 lines
unisynapse/pipelines/pair.py, Python, 1,223 lines
unisynapse/pipelines/run.py, Python, 322 lines
unisynapse/pipelines/track.py, Python, 1,309 lines
unisynapse/track/__init__.py, Python, 26 lines
unisynapse/track/api.py, Python, 159 lines
unisynapse/track/mincostflow.py, Python, 2,591 lines, 4 matches
unisynapse/track/post_processing.py, Python, 314 lines
unisynapse/track/visualization.py, Python, 303 lines
unisynapse/viz/__init__.py, Python, 22 lines
unisynapse/viz/debug.py, Python, 224 lines
unisynapse/viz/overlays.py, Python, 2,055 lines
LICENSE, License, 21 lines
README.md, Text, 284 lines

At the source: github.com/datajamesallen/Garbett_2026_eLife

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
67 scripts, each with its path and the digest of its content;
11 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

doi:10.5061/dryad.5qfttdzp5, at Dryad; found in “Data availability”

Data availability

Code generated in this study is available at the following deposit: https://github.com/datajamesallen/Garbett_2026_eLife (copy archived at Allen, 2026). Live imaging data generated in this study is available at Dryad: https://doi.org/10.5061/dryad.5qfttdzp5. All materials generated in this study will be shared by the lead contact upon request.

The following dataset was generated:

GarbettKA AllenJP LopezJM SmithCM SandoRC Dryad Digital Repository10.5061/dryad.5qfttdzp52026Data from: Stable excitatory-inhibitory synapse balance despite dynamic turnoverPMC1340563842507759

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 5 authors, 6 keywords, 4 MeSH terms, 2 funders, 47 references, 27 RRIDs.

Cite

This paper

Garbett, K. A., Allen, J. P., Lopez, J. M., Smith, C. M., & Sando, R. C. (2026). Stable excitatory-inhibitory synapse balance despite dynamic turnover. eLife, 14, RP107635. https://doi.org/10.7554/elife.107635

BibTeX

@article{garbett2026stable, author = {Garbett, Krassimira A and Allen, James P and Lopez, Jaybree M and Smith, Cassandra M and Sando, Richard C}, title = {{Stable excitatory-inhibitory synapse balance despite dynamic turnover}}, journal = {eLife}, year = {2026}, month = jul, volume = {14}, pages = {RP107635}, publisher = {eLife Sciences Publications, Ltd}, issn = {2050-084X}, doi = {10.7554/elife.107635}, url = {https://doi.org/10.7554/elife.107635}, pmid = {42507759}, pmcid = {PMC13405638} }

RIS

TY - JOUR AU - Garbett, Krassimira A AU - Allen, James P AU - Lopez, Jaybree M AU - Smith, Cassandra M AU - Sando, Richard C TI - Stable excitatory-inhibitory synapse balance despite dynamic turnover T2 - eLife J2 - eLife PY - 2026 DA - 2026/07/27 VL - 14 SP - RP107635 SN - 2050-084X PB - eLife Sciences Publications, Ltd DO - 10.7554/elife.107635 UR - https://doi.org/10.7554/elife.107635 LA - en ER -

CSL-JSON

{ "id": "10.7554/elife.107635", "type": "article-journal", "title": "Stable excitatory-inhibitory synapse balance despite dynamic turnover", "container-title": "eLife", "author": [ { "family": "Garbett", "given": "Krassimira A" }, { "family": "Allen", "given": "James P" }, { "family": "Lopez", "given": "Jaybree M" }, { "family": "Smith", "given": "Cassandra M" }, { "family": "Sando", "given": "Richard C" } ], "container-title-short": "eLife", "volume": "14", "page": "RP107635", "DOI": "10.7554/elife.107635", "PMID": "42507759", "PMCID": "PMC13405638", "ISSN": "2050-084X", "publisher": "eLife Sciences Publications, Ltd", "URL": "https://doi.org/10.7554/elife.107635", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 27 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.