synapse_net.cristae_analysis
1import multiprocessing as mp 2import os 3from concurrent import futures 4from typing import Callable, Dict, Optional, Tuple, Union 5 6import numpy as np 7import pandas as pd 8from scipy.ndimage import binary_dilation, binary_erosion, center_of_mass 9from scipy.ndimage import label as ndimage_label 10from skimage.measure import mesh_surface_area, regionprops 11from skimage.morphology import ball, disk, local_maxima 12from tqdm import tqdm 13 14from bioimage_cpp.distance import distance_transform, geodesic_distances_mesh 15from bioimage_cpp.filters import structure_tensor_eigenvalues 16from bioimage_cpp.mesh import marching_cubes 17from bioimage_cpp.skeleton import teasar 18 19 20# --------------------------------------------------------------------------- 21# Internal helpers 22# --------------------------------------------------------------------------- 23 24def _to_sampling(voxel_size: Union[float, Dict[str, float]], ndim: int) -> np.ndarray: 25 axes = ("z", "y", "x") if ndim == 3 else ("y", "x") 26 if isinstance(voxel_size, dict): 27 return np.array([voxel_size[ax] for ax in axes[:ndim]], dtype=float) 28 return np.full(ndim, float(voxel_size)) 29 30 31def _voxel_radius(thickness_nm: float, voxel_size: Union[float, Dict[str, float]], ndim: int) -> int: 32 return max(1, int(round(thickness_nm / float(np.mean(_to_sampling(voxel_size, ndim)))))) 33 34 35def _voxel_radius_xy(thickness_nm: float, voxel_size: Union[float, Dict[str, float]]) -> int: 36 """Membrane radius in XY pixels — uses only the Y and X voxel sizes.""" 37 if isinstance(voxel_size, dict): 38 xy_nm = (voxel_size["y"] + voxel_size["x"]) / 2.0 39 else: 40 xy_nm = float(voxel_size) 41 return max(1, int(round(thickness_nm / xy_nm))) 42 43 44def _gap_radius( 45 voxel_size: Union[float, Dict[str, float]], 46 membrane_thickness_nm: float, 47 border_gap_nm: Optional[float], 48 ndim: int, 49) -> int: 50 """Border-zone / mesh-trim radius in voxels; ``border_gap_nm`` defaults to ``membrane_thickness_nm``.""" 51 gap_nm = border_gap_nm if border_gap_nm is not None else membrane_thickness_nm 52 return _voxel_radius(gap_nm, voxel_size, ndim) 53 54 55def _border_zone( 56 shape: tuple, radius: Union[int, np.ndarray], boundary: Optional[np.ndarray] = None 57) -> np.ndarray: 58 """Boolean mask that is True within ``radius`` voxels of a volume face. 59 60 This is the region where the segmentation is cut off by the field of view, so membrane presence is 61 *unknown* rather than absent — see the trim in :func:`approximate_membrane`. 62 63 **The zone is a global fact about the volume**, which is why a caller working on a bounding-box 64 crop cannot describe it with a per-face boolean. A crop whose low face sits one voxel inside the 65 volume still overlaps the zone across its first ``radius - 1`` voxels, yet an exact-face test 66 ("is this bbox face at the volume face?") reports no overlap at all and suppresses nothing. Pass 67 **per-face radii** instead. Measured over 4000 random ``(vol_shape, bbox, radius)`` combinations 68 against the ground truth ``_border_zone(vol_shape, radius)[crop_slices]``, per-face radii disagree 69 in 0 cases and the boolean form in 1250 (31%). 70 71 Args: 72 shape: Shape of the array to build the mask for. 73 radius: Width of the zone in voxels — either a scalar, applied to every face, or an 74 ``(ndim, 2)`` array of per-face widths where ``radius[a, s]`` is side ``s`` (0 = low, 75 1 = high) of axis ``a``. **Pass per-face radii whenever ``shape`` is a bounding-box 76 crop**: for a crop at ``bbox`` in a volume of shape ``vol_shape`` they are 77 ``max(0, radius - bbox[a])`` and ``max(0, radius - (vol_shape[a] - bbox_hi[a]))``, which 78 give back the full radius at a genuine volume face and 0 at a face clear of the zone — 79 so the boolean behaviour below falls out of them for free. 80 boundary: Optional ``(ndim, 2)`` bool zeroing the radius on faces that are not genuine volume 81 faces; defaults to all True. This answers the *mesh* question ("is the object clipped by 82 the volume here?", see :func:`_open_trimmed_mesh`), not the border-zone one, and it is 83 redundant once per-face radii are given — an interior face already gets radius 0. **Do 84 not pass both**: the radii are multiplied by it, so a per-face radius on a face this 85 marks False is silently zeroed. 86 87 Returns: 88 Boolean mask of the border zone. 89 """ 90 ndim = len(shape) 91 radii = np.broadcast_to(np.asarray(radius, dtype=int), (ndim, 2)) 92 if boundary is not None: 93 radii = radii * np.asarray(boundary, dtype=bool) 94 95 mask = np.zeros(shape, dtype=bool) 96 for ax in range(ndim): 97 lo, hi = int(radii[ax, 0]), int(radii[ax, 1]) 98 if lo: 99 idx_lo = [slice(None)] * ndim 100 idx_lo[ax] = slice(0, lo) 101 mask[tuple(idx_lo)] = True 102 if hi: 103 idx_hi = [slice(None)] * ndim 104 # max(0, ...): a radius wider than the axis would give a negative start, which counts 105 # from the end and leaves the low voxels outside the zone. 106 idx_hi[ax] = slice(max(0, shape[ax] - hi), None) 107 mask[tuple(idx_hi)] = True 108 return mask 109 110 111def _surface_mesh( 112 mask: np.ndarray, 113 sampling: np.ndarray, 114 closed_faces: Optional[np.ndarray] = None, 115) -> Optional[Tuple[np.ndarray, np.ndarray]]: 116 """Triangle-mesh surface of a binary mask via marching cubes. 117 118 Each side of the array is padded by one background voxel before meshing so that objects touching 119 the array edge yield a closed surface there. ``marching_cubes`` only emits triangles for cube 120 cells that exist inside the array, so leaving a side *unpadded* omits that boundary face — an open 121 mesh at that face — while interior surfaces still close. ``closed_faces`` chooses this per side: 122 pad+close a face where there is genuine background beyond it, leave open a face where the object 123 is clipped by the volume boundary (membrane presence unknown there). 124 125 Vertices are returned in the mask's own **unpadded** index frame (nm): a mask voxel at array index 126 ``(z, y, x)`` maps to physical coordinates ``index * sampling`` regardless of the padding, so 127 callers snapping voxel coordinates onto the mesh use ``index * sampling`` with no offset. 128 129 Args: 130 mask: Binary segmentation. 131 sampling: Voxel size per axis (nm), in array (z, y, x) order. 132 closed_faces: Optional ``(ndim, 2)`` boolean array; ``closed_faces[a, s]`` True closes 133 (pads) side ``s`` (0 = low, 1 = high) of axis ``a``, False leaves it open. Defaults to 134 all True (fully closed watertight surface). 135 136 Returns: 137 (vertices, faces) with vertices in nm (unpadded mask frame), or None if the mask is empty. 138 """ 139 binary = mask.astype(bool) 140 if not binary.any(): 141 return None 142 if closed_faces is None: 143 closed_faces = np.ones((binary.ndim, 2), dtype=bool) 144 else: 145 closed_faces = np.asarray(closed_faces, dtype=bool) 146 pad_width = [(int(closed_faces[a, 0]), int(closed_faces[a, 1])) for a in range(binary.ndim)] 147 padded = np.pad(binary.astype(np.float32), pad_width) 148 verts, faces, _, _ = marching_cubes(padded, level=0.5, spacing=tuple(float(s) for s in sampling)) 149 pad_before = np.array([pw[0] for pw in pad_width], dtype=float) 150 verts = verts - pad_before * np.asarray(sampling, dtype=float) 151 return verts, faces 152 153 154def _surface_area( 155 mask: np.ndarray, 156 sampling: np.ndarray, 157 closed_faces: Optional[np.ndarray] = None, 158) -> float: 159 """Surface area (nm^2) of a binary mask via marching cubes. 160 161 Args: 162 mask: Binary segmentation. 163 sampling: Voxel size per axis (nm), in array (z, y, x) order. 164 closed_faces: Optional per-side padding spec forwarded to :func:`_surface_mesh` — leave a 165 clipped volume-boundary face open so its fabricated cap is not counted as surface area. 166 Defaults to a fully closed (watertight) surface. 167 168 Returns: 169 Surface area in nm^2, or NaN if the mask is empty. 170 """ 171 mesh = _surface_mesh(mask, sampling, closed_faces=closed_faces) 172 if mesh is None: 173 return np.nan 174 return float(mesh_surface_area(*mesh)) 175 176 177def _open_trimmed_mesh( 178 mask: np.ndarray, 179 sampling: np.ndarray, 180 gap_radius: int, 181 boundary: np.ndarray, 182) -> Optional[Tuple[np.ndarray, np.ndarray]]: 183 """Surface mesh of ``mask`` trimmed to the certain region and left OPEN at volume-boundary faces. 184 185 Near a clipped volume face the segmentation is cut off and membrane presence is unknown, so the 186 surface must neither flare into that region nor be capped there — a cap is a fabricated flat disk 187 that lets geodesics shortcut straight across instead of wrapping around the tube wall. For each 188 ``(axis, side)`` flagged True in ``boundary`` (an ``(ndim, 2)`` bool of volume-boundary faces), 189 ``gap_radius`` voxels are cropped off that face so the trim plane becomes an array boundary, and 190 that face is then left open (marching cubes omits it). Interior faces stay closed. 191 192 Args: 193 mask: Binary segmentation. 194 sampling: Voxel size per axis (nm), in array (z, y, x) order. 195 gap_radius: Border-zone width in voxels cropped off each flagged face (matches the membrane's 196 border-gap trim). 197 boundary: ``(ndim, 2)`` bool; True where the face is at the volume boundary (crop + open). 198 199 Returns: 200 (vertices, faces) with vertices in nm in the original (uncropped) ``mask`` index frame, or None 201 if the trimmed mask is empty. 202 """ 203 ndim = mask.ndim 204 boundary = np.asarray(boundary, dtype=bool) 205 lo = [gap_radius if boundary[a, 0] else 0 for a in range(ndim)] 206 hi = [mask.shape[a] - (gap_radius if boundary[a, 1] else 0) for a in range(ndim)] 207 mesh = _surface_mesh( 208 mask[tuple(slice(lo[a], hi[a]) for a in range(ndim))], sampling, closed_faces=~boundary 209 ) 210 if mesh is None: 211 return None 212 verts, faces = mesh 213 verts = verts + np.array(lo, dtype=float) * np.asarray(sampling, dtype=float) 214 return verts, faces 215 216 217def _medial_axis_thickness_nm(mask: np.ndarray, sampling: np.ndarray) -> float: 218 """Local thickness (nm) of a mask via the distance transform, with no mesh generation. 219 220 The medial axis is approximated by the local maxima of the interior EDT; the thickness is 221 ``2 × mean(EDT)`` there (the EDT at the medial axis is the half-thickness). This is the same 222 estimator used by :func:`compute_crista_morphology`'s ``medial_axis`` branch, factored out so 223 the distance-based (``method="fast"``) surface-area estimates can reuse it. 224 225 Args: 226 mask: Binary segmentation. 227 sampling: Voxel size per axis (nm), in array (z, y, x) order. 228 229 Returns: 230 Mean local thickness in nm, or NaN if the mask is empty. 231 """ 232 binary = mask.astype(bool) 233 if not binary.any(): 234 return np.nan 235 dist = distance_transform(binary, sampling=tuple(float(s) for s in sampling), number_of_threads=1) 236 ridges = local_maxima(dist) & binary 237 ridge_dists = dist[ridges] 238 return float(2.0 * np.mean(ridge_dists)) if ridge_dists.size > 0 else np.nan 239 240 241def _available_memory_bytes() -> int: 242 """Best-effort available RAM in bytes (used to keep parallel working sets from OOMing).""" 243 try: 244 import psutil 245 return int(psutil.virtual_memory().available) 246 except Exception: 247 try: 248 return int(os.sysconf("SC_AVPHYS_PAGES") * os.sysconf("SC_PAGE_SIZE")) 249 except Exception: 250 return 4 * 1024 ** 3 251 252 253def _bounded_workers(n_jobs: int, per_worker_bytes: int, fraction: float = 0.5) -> int: 254 """Resolve n_jobs to a worker count whose combined working set fits in memory. 255 256 n_jobs: 1 = serial, -1 = all cores, else that many. The result is additionally capped so 257 ``workers * per_worker_bytes <= fraction * available_RAM`` (at least 1). 258 """ 259 workers = os.cpu_count() if n_jobs == -1 else max(1, int(n_jobs)) 260 if per_worker_bytes > 0: 261 budget = int(_available_memory_bytes() * fraction) 262 workers = min(workers, max(1, budget // int(per_worker_bytes))) 263 return int(max(1, workers)) 264 265 266# --------------------------------------------------------------------------- 267# Membrane approximation 268# --------------------------------------------------------------------------- 269 270def approximate_membrane( 271 mito_segmentation: np.ndarray, 272 voxel_size: Union[float, Dict[str, float]], 273 membrane_thickness_nm: float = 8.0, 274 border_gap_nm: Optional[float] = None, 275 n_jobs: int = 1, 276 membrane_mode: str = "slice_2d", 277 return_lumen: bool = False, 278) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]: 279 """Approximate the mitochondrial membrane as the outer shell of the segmentation. 280 281 Two shell constructions are available via ``membrane_mode``: 282 283 - ``"slice_2d"`` (default): erode each Z-slice **independently in 2D** by an XY disk of radius 284 ``round(thickness / xy_voxel)`` and keep ``slice & ~eroded``. A mitochondrion that changes shape 285 rapidly along Z does not bleed into neighbouring slices, and the per-slice erosions are 286 parallelised over Z (``n_jobs``). A separable Z-only erosion (radius ``round(thickness / 287 z_voxel)``, no XY coupling) then adds **Z-caps** where a mito column truly ends in Z; ends 288 clipped by a volume Z-face are left uncapped (``border_value=1`` + the ``border_gap`` trim). The 289 XY shell can still fragment across slices. (2D inputs get a single 2D erosion.) 290 - ``"shell_3d"``: a full 3D morphological erosion, ``mito & ~erode3d(mito, k)`` with 291 ``k = round(thickness / mean_voxel)`` iterations of a 3×3×3 structuring element, per instance on 292 its padded bounding box. A single **connected** shell including the Z-caps (no per-slice 293 fragmentation), at a higher cost; thickness acts in all axes. 294 295 The eroded interior is the "lumen"; its surface is the single-wall mesh used by the geodesic 296 backend and the display, so the mesh follows the chosen mode. 297 298 Membrane voxels within ``border_gap_nm`` of any volume face are removed so clipped mito edges are 299 not treated as membrane. The lumen is NOT trimmed here — the mesh is trimmed to the certain region 300 (and left open there) at mesh time by :func:`_open_trimmed_mesh`, which requires the untrimmed 301 interior to produce an open cut rather than a fabricated cap. 302 303 Implementation notes: ``"slice_2d"`` erodes each Z-slice on the mito XY bbox with a 304 ``membrane_radius`` margin, so the cropped ``border_value=1`` erosion matches eroding the full 305 slice (empty slices are skipped), then adds Z-caps via a separable Z-only line erosion (which 306 inspects only the same column, so no XY-shape bleed); ``border_value=1`` leaves ends clipped by a 307 volume Z-face uncapped, and the ``border_gap`` removal clears anything near a face, so only true 308 ends are capped. ``"shell_3d"`` erodes the *merged* binary in each instance's padded bbox so 309 instances that share a boundary are handled together. 310 311 Args: 312 mito_segmentation: Instance label array (background = 0). 313 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 314 membrane_thickness_nm: Thickness of the membrane shell in nm. 315 border_gap_nm: Distance from each volume face within which membrane voxels are 316 suppressed. Defaults to membrane_thickness_nm when None. 317 n_jobs: Workers for the per-Z-slice erosion (``"slice_2d"`` only): 1 = serial, -1 = all cores. 318 Results are identical regardless of n_jobs. 319 membrane_mode: ``"slice_2d"`` (default, per-slice 2D, z-parallel) or ``"shell_3d"`` 320 (connected 3D shell). 321 return_lumen: If True, also return the eroded-mito interior ("lumen") mask — the single-wall 322 surface source for the geodesic/display mesh. It is the plain eroded interior (not 323 ``mito & ~membrane``, which would re-include the outer shell where the membrane is 324 border-trimmed); trimming to the certain region happens at mesh time. 325 326 Returns: 327 membrane_mask: Binary mask of the mitochondrial membrane (outer shell), with border-adjacent 328 voxels zeroed out. If ``return_lumen`` is True, returns ``(membrane_mask, lumen_mask)`` 329 where ``lumen_mask`` is the (untrimmed) eroded interior described above. 330 """ 331 if membrane_mode not in ("slice_2d", "shell_3d"): 332 raise ValueError(f"membrane_mode must be 'slice_2d' or 'shell_3d', got {membrane_mode!r}") 333 ndim = mito_segmentation.ndim 334 mito_binary = mito_segmentation > 0 335 336 # NOTE (possible simplification, deferred): the "shell_3d" branch below builds the shell with an 337 # iterated 3x3x3 erosion per instance. It could likely be a single anisotropic distance transform 338 # instead — membrane = mito & (distance_transform(mito, sampling) <= thickness), lumen = the rest — 339 # which is simpler and handles anisotropy directly. This would NOT replace "slice_2d": a distance 340 # transform couples all axes, so it cannot reproduce slice_2d's per-Z-slice-independent erosion, 341 # whose whole purpose is to stop the shell bleeding across slices in XY. Worth investigating. 342 if membrane_mode == "shell_3d": 343 sampling = _to_sampling(voxel_size, ndim) 344 k = max(1, int(round(float(membrane_thickness_nm) / float(np.mean(sampling))))) 345 struct = np.ones((3,) * ndim, dtype=bool) 346 membrane_mask = np.zeros(mito_segmentation.shape, dtype=bool) 347 lumen_mask = np.zeros(mito_segmentation.shape, dtype=bool) 348 for prop in regionprops(mito_segmentation): 349 bbox = prop.bbox 350 sl = tuple( 351 slice(max(0, bbox[i] - k), min(mito_segmentation.shape[i], bbox[i + ndim] + k)) 352 for i in range(ndim) 353 ) 354 sub = mito_binary[sl] 355 eroded = binary_erosion(sub, structure=struct, iterations=k, border_value=1) 356 cur = mito_segmentation[sl] == prop.label 357 membrane_mask[sl] |= cur & ~eroded 358 lumen_mask[sl] |= cur & eroded 359 elif ndim == 3: 360 membrane_radius = _voxel_radius_xy(membrane_thickness_nm, voxel_size) 361 struct = disk(membrane_radius) 362 membrane_mask = np.zeros_like(mito_binary) 363 lumen_mask = np.zeros_like(mito_binary) 364 coords = np.argwhere(mito_binary) 365 if coords.size: 366 zmin, ymin, xmin = coords.min(axis=0) 367 zmax, ymax, xmax = coords.max(axis=0) + 1 368 m = membrane_radius 369 y0, y1 = max(0, ymin - m), min(mito_binary.shape[1], ymax + m) 370 x0, x1 = max(0, xmin - m), min(mito_binary.shape[2], xmax + m) 371 372 def _erode_slice(z): 373 sl = mito_binary[z, y0:y1, x0:x1] 374 if not sl.any(): 375 return z, None 376 eroded = binary_erosion(sl, structure=struct, border_value=1) 377 return z, (sl & ~eroded, eroded) 378 379 z_range = range(int(zmin), int(zmax)) 380 if n_jobs == 1: 381 results = [_erode_slice(z) for z in z_range] 382 else: 383 n_workers = mp.cpu_count() if n_jobs == -1 else n_jobs 384 with futures.ThreadPoolExecutor(n_workers) as tp: 385 results = list(tp.map(_erode_slice, z_range)) 386 for z, res in results: 387 if res is not None: 388 mem_sl, lum_sl = res 389 membrane_mask[z, y0:y1, x0:x1] = mem_sl 390 lumen_mask[z, y0:y1, x0:x1] = lum_sl 391 392 k_z = max(1, int(round(float(membrane_thickness_nm) / float(_to_sampling(voxel_size, ndim)[0])))) 393 z0m, z1m = max(0, int(zmin) - k_z), min(mito_binary.shape[0], int(zmax) + k_z) 394 sub = mito_binary[z0m:z1m, y0:y1, x0:x1] 395 z_eroded = binary_erosion(sub, structure=np.ones((2 * k_z + 1, 1, 1), dtype=bool), border_value=1) 396 membrane_mask[z0m:z1m, y0:y1, x0:x1] |= sub & ~z_eroded 397 lumen_mask[z0m:z1m, y0:y1, x0:x1] &= z_eroded 398 else: 399 membrane_radius = _voxel_radius(membrane_thickness_nm, voxel_size, ndim) 400 eroded = binary_erosion(mito_binary, structure=disk(membrane_radius), border_value=1) 401 membrane_mask = mito_binary & ~eroded 402 lumen_mask = mito_binary & eroded 403 404 gap_radius = _gap_radius(voxel_size, membrane_thickness_nm, border_gap_nm, ndim) 405 membrane_mask &= ~_border_zone(mito_segmentation.shape, gap_radius) 406 if return_lumen: 407 return membrane_mask.astype(bool), lumen_mask.astype(bool) 408 return membrane_mask.astype(bool) 409 410 411# --------------------------------------------------------------------------- 412# Orientation 413# --------------------------------------------------------------------------- 414 415def compute_crista_orientation( 416 crista_mask: np.ndarray, 417 voxel_size: Union[float, Dict[str, float]], 418 neighborhood_size_nm: float = 30.0, 419) -> np.ndarray: 420 """Compute the per-voxel crista orientation anisotropy via the structure tensor. 421 422 Uses ``bioimage_cpp.filters.structure_tensor_eigenvalues`` (a fast C++ routine). Only the 423 anisotropy is produced (the principal directions / eigenvectors are not computed). 424 425 The structure tensor's outer/integration sigma is ``neighborhood_size_nm`` per axis (in voxels); 426 the inner (derivative) sigma must be > 0, so a minimal 1-voxel scale is used. Eigenvalues are 427 non-negative in theory, but the solver emits tiny negatives for near-rank-deficient tensors 428 (degenerate sheets/tubes), so they are clamped to 0 and the ratio is taken as 429 ``max/min`` over the trailing axis — order-agnostic and sign-safe, so a tiny negative minor 430 eigenvalue cannot flip the denominator and blow the ratio up. 431 432 Args: 433 crista_mask: Binary crista segmentation. 434 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 435 neighborhood_size_nm: Gaussian integration radius in nm for tensor averaging (the 436 structure tensor's outer/integration scale). 437 438 Returns: 439 anisotropy: (...) — λ_max / (λ_min + ε) per voxel. High values indicate a strongly 440 directional crista (e.g. parallel lamellae); low values indicate isotropic or 441 tubular/disordered morphology. Magnitude only (rotation-invariant). 442 """ 443 ndim = crista_mask.ndim 444 sampling = _to_sampling(voxel_size, ndim) 445 446 outer_sigma = [float(s) for s in (neighborhood_size_nm / sampling)] 447 inner_sigma = 1.0 448 evals = structure_tensor_eigenvalues(crista_mask.astype(np.float32), inner_sigma, outer_sigma) 449 evals = np.clip(evals, 0.0, None) 450 anisotropy = evals.max(axis=-1) / (evals.min(axis=-1) + 1e-10) 451 return anisotropy.astype(np.float32) 452 453 454def _scale_voxel_size( 455 voxel_size: Union[float, Dict[str, float]], factor: float 456) -> Union[float, Dict[str, float]]: 457 """Multiply a voxel size (scalar or z/y/x dict) by ``factor``, preserving its type.""" 458 if isinstance(voxel_size, dict): 459 return {ax: voxel_size[ax] * factor for ax in voxel_size} 460 return float(voxel_size) * factor 461 462 463def _downsampled_orientation_anisotropy( 464 crista_mask: np.ndarray, 465 voxel_size: Union[float, Dict[str, float]], 466 factor: int = 2, 467) -> float: 468 """Mean crista orientation anisotropy computed on a downsampled crop (fast, approximate). 469 470 The crista mask is **nearest-neighbour** downsampled by ``factor`` per axis (strided subsampling, 471 i.e. every ``factor``-th voxel), and :func:`compute_crista_orientation` is run on that coarser grid 472 at the correspondingly scaled voxel size (so the physical structure-tensor neighbourhood is 473 unchanged). This is ~``factor**ndim`` times cheaper than the full-resolution structure tensor — the 474 dominant cost of the analysis. Nearest-neighbour keeps the segmentation binary; block-mean 475 (``downscale_local_mean``) would blur it into meaningless partial-occupancy values. 476 477 Because thin cristae (only a few voxels across) lose structure when downsampled, the returned 478 anisotropy is a *relative* indicator only: it preserves the ordering between mitochondria but is 479 not comparable in magnitude to the full-resolution ``method="exact"`` value. 480 481 Args: 482 crista_mask: Binary crista segmentation. 483 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 484 factor: Integer downsampling factor per axis. 485 486 Returns: 487 Mean anisotropy over the downsampled crista region, or NaN if it vanishes when downsampled. 488 """ 489 ndim = crista_mask.ndim 490 sub = crista_mask[(slice(None, None, factor),) * ndim] 491 region = sub.astype(bool) 492 if not region.any(): 493 return np.nan 494 anisotropy = compute_crista_orientation(sub, _scale_voxel_size(voxel_size, factor)) 495 return float(np.mean(anisotropy[region])) 496 497 498# --------------------------------------------------------------------------- 499# Proximity 500# --------------------------------------------------------------------------- 501 502def compute_crista_proximity( 503 crista_mask: np.ndarray, 504 membrane_mask: np.ndarray, 505 voxel_size: Union[float, Dict[str, float]], 506 membrane_distance: Optional[np.ndarray] = None, 507) -> Tuple[np.ndarray, Dict[str, float]]: 508 """Distance from each crista voxel to the nearest membrane voxel (nm). 509 510 Args: 511 crista_mask: Binary crista segmentation. 512 membrane_mask: Binary membrane mask (OM or IMM). 513 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 514 membrane_distance: Optional precomputed per-voxel distance to the nearest membrane 515 voxel (nm), i.e. ``distance_transform(~membrane_mask, sampling=...)``. When 516 given, the distance transform is not recomputed (used to avoid redundant work in 517 :func:`compute_mito_crista_statistics`). 518 519 Returns: 520 distance_map: Per-voxel distance to membrane (nm); zero outside crista. 521 summary_stats: min_nm, median_nm, max_nm. 522 """ 523 sampling = _to_sampling(voxel_size, crista_mask.ndim) 524 if membrane_distance is None: 525 dist = distance_transform(~membrane_mask.astype(bool), sampling=sampling.tolist(), number_of_threads=1) 526 else: 527 dist = membrane_distance 528 crista_dists = dist[crista_mask.astype(bool)] 529 530 if crista_dists.size == 0: 531 summary: Dict[str, float] = {"min_nm": np.nan, "median_nm": np.nan, "max_nm": np.nan} 532 else: 533 summary = { 534 "min_nm": float(crista_dists.min()), 535 "median_nm": float(np.median(crista_dists)), 536 "max_nm": float(crista_dists.max()), 537 } 538 539 distance_map = np.zeros(crista_mask.shape, dtype=np.float32) 540 distance_map[crista_mask.astype(bool)] = crista_dists 541 return distance_map, summary 542 543 544# --------------------------------------------------------------------------- 545# Contact sites 546# --------------------------------------------------------------------------- 547 548def detect_contact_sites( 549 crista_mask: np.ndarray, 550 membrane_mask: np.ndarray, 551 voxel_size: Union[float, Dict[str, float]], 552) -> Tuple[np.ndarray, Dict[str, float]]: 553 """Detect crista-membrane contact sites as the direct overlap of the two masks. 554 555 Contact = crista voxels that are also membrane voxels (the pure intersection of the crista 556 mask and the mitochondrial membrane band). No dilation/erosion is applied here, so the 557 detected junctions correspond exactly to the visible overlap of the two layers; connected 558 overlaps are grouped into junctions with 26-connectivity in 3D. 559 560 Args: 561 crista_mask: Binary crista segmentation. 562 membrane_mask: Binary mitochondrial membrane mask (as displayed). 563 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 564 565 Returns: 566 contact_labels: Integer array (same shape as the input) where each connected 567 junction has a unique ID (0 = background, 1..n = junctions). Contact voxel 568 coordinates are recoverable via ``np.argwhere(contact_labels > 0)``. 569 summary: contact_voxel_count, crista_junction_count, contact_volume_nm3. 570 """ 571 ndim = crista_mask.ndim 572 sampling = _to_sampling(voxel_size, ndim) 573 voxel_vol = float(np.prod(sampling)) 574 575 contact_mask = crista_mask.astype(bool) & membrane_mask.astype(bool) 576 577 connectivity_struct = np.ones(ndim * (3,), dtype=bool) 578 contact_labels, n_regions = ndimage_label(contact_mask, structure=connectivity_struct) 579 contact_voxel_count = int(np.count_nonzero(contact_labels)) 580 581 return contact_labels, { 582 "contact_voxel_count": contact_voxel_count, 583 "crista_junction_count": int(n_regions), 584 "contact_volume_nm3": float(contact_voxel_count) * voxel_vol, 585 } 586 587 588def _solver_threads(n_jobs: int) -> int: 589 """Map the ``n_jobs`` convention (-1/0 = all cores) onto the C++ solvers' ``number_of_threads``.""" 590 return 0 if n_jobs in (-1, 0) else max(1, int(n_jobs)) 591 592 593def _given(**options): 594 """Drop the options left as None, so the callee's own default applies instead. 595 596 Lets an intermediate function forward an optional tuning parameter without having to know, or 597 restate, the default the consuming function declares. ``None`` therefore means "unset" all the way 598 down rather than being resolved into a number at every layer. 599 """ 600 return {name: value for name, value in options.items() if value is not None} 601 602 603def _skeleton_graph(vertices: np.ndarray, edges: np.ndarray, min_skeleton_nm: float): 604 """Build the TEASAR skeleton as a weighted graph and drop speck components. 605 606 Every segmentation speck contributes its own miniature skeleton, and each of those contributes 607 termini, so a noisy crista mask produces far more free ends than it has cristae. Filtering whole 608 components by total arclength removes them at the root. Keyed on nm rather than vertex count so the 609 result does not depend on the voxel size. 610 611 Args: 612 vertices: (n, 3) skeleton vertex coordinates in nm. 613 edges: (m, 2) integer vertex-index pairs. 614 min_skeleton_nm: Components with a total edge length below this are dropped, as are isolated 615 single vertices. 0 keeps the graph as TEASAR produced it. 616 617 Returns: 618 A ``networkx.Graph`` over the surviving vertex indices, each edge weighted by its length in nm. 619 Node ids index into ``vertices`` unchanged, so the caller can look coordinates up directly. 620 """ 621 import networkx as nx 622 623 graph = nx.Graph() 624 graph.add_nodes_from(range(len(vertices))) 625 for node_a, node_b in edges: 626 node_a, node_b = int(node_a), int(node_b) 627 if node_a == node_b: # a self-loop carries no length and would only distort the degrees 628 continue 629 graph.add_edge(node_a, node_b, weight=float(np.linalg.norm(vertices[node_a] - vertices[node_b]))) 630 if min_skeleton_nm <= 0: 631 return graph 632 633 for component in list(nx.connected_components(graph)): 634 subgraph = graph.subgraph(component) 635 total = sum(data["weight"] for _, _, data in subgraph.edges(data=True)) 636 if len(component) == 1 or total < min_skeleton_nm: 637 graph.remove_nodes_from(list(component)) 638 return graph 639 640 641def _merge_termini(graph, vertices: np.ndarray, terminus_nodes, merge_nm: float): 642 """Collapse each cluster of nearby termini on one skeleton component to a single representative. 643 644 This is what actually tames the terminus count. TEASAR spans a lamella with a caterpillar of short 645 side branches, so one crista end appears as a fan of degree-1 nodes a few nm apart, all describing 646 the same end. Clustering is single-linkage over pairs within ``merge_nm``, **restricted to pairs on 647 the same connected component of** ``graph``. That restriction is essential rather than cosmetic: 648 densely packed cristae sit a few nm apart, so unrestricted single-linkage chains termini across 649 separate cristae and merges the whole field. Measured on a 40-lamella mask at 6 nm spacing, an 8 nm 650 radius gives 5 clusters for the entire volume without the restriction and 200 with it. 651 652 The kept representative is the cluster member nearest the cluster centroid — a real skeleton 653 vertex, not the centroid itself, so a terminus always lies on the skeleton and inside the crista. 654 655 Args: 656 graph: The skeleton graph from :func:`_skeleton_graph`. 657 vertices: (n, 3) vertex coordinates in nm. 658 terminus_nodes: Iterable of candidate terminus node ids (degree-1 nodes of ``graph``). 659 merge_nm: Cluster radius in nm. 0 disables merging. 660 661 Returns: 662 A sorted list of the surviving terminus node ids. 663 """ 664 import networkx as nx 665 666 terminus_nodes = list(terminus_nodes) 667 if merge_nm <= 0 or len(terminus_nodes) < 2: 668 return sorted(terminus_nodes) 669 670 from scipy.spatial import cKDTree 671 672 points = vertices[terminus_nodes] 673 component_of = {} 674 for index, component in enumerate(nx.connected_components(graph)): 675 for node in component: 676 component_of[node] = index 677 component_ids = np.array([component_of[node] for node in terminus_nodes]) 678 679 clusters = nx.Graph() 680 clusters.add_nodes_from(range(len(terminus_nodes))) 681 for i, j in cKDTree(points).query_pairs(merge_nm): 682 if component_ids[i] == component_ids[j]: 683 clusters.add_edge(i, j) 684 685 kept = [] 686 for cluster in nx.connected_components(clusters): 687 members = list(cluster) 688 centroid = points[members].mean(axis=0) 689 kept.append(terminus_nodes[members[int(np.argmin(np.linalg.norm(points[members] - centroid, axis=1)))]]) 690 return sorted(kept) 691 692 693def compute_crista_skeleton( 694 crista_mask: np.ndarray, 695 voxel_size: Union[float, Dict[str, float]], 696 n_jobs: int = 1, 697 min_skeleton_nm: float = 10.0, 698 terminus_merge_nm: float = 4.0, 699 return_edges: bool = False, 700): 701 """The crista centerline skeleton and its termini, for inspection and display. 702 703 Wraps ``bioimage_cpp.skeleton.teasar`` and then **cleans the resulting graph up**, because the raw 704 TEASAR output is not usable as a set of crista ends: it spans a lamella with a caterpillar of short 705 side branches whose tips line the sheet rim, every segmentation speck adds its own miniature 706 skeleton, and counting ``degree <= 1`` also counts isolated vertices. On a tomogram-scale 707 40-lamella mask that produces 10080 termini where roughly 80 are real. 708 709 The cleanup is :func:`_skeleton_graph` (drop speck components) followed by :func:`_merge_termini` 710 (collapse each fan of nearby termini on one component), which brings the same mask to 480, and 711 ``degree == 1`` instead of ``<= 1``. 712 713 **Why no leaf-branch pruning.** Cutting short leaf branches off a surviving component is the obvious 714 first idea, and it was measured to be actively harmful. TEASAR's medial axis of a lamella is a 715 caterpillar — a main path with many short side leaves — and at the sheet's real end the main path 716 itself arrives as a short leaf off a nearby branch node. A length threshold therefore deletes the 717 real crista end: on the test lamella spanning y 13-46 it left termini only at y 45-46, losing the 718 y=13 junction entirely. It also fragmented components, so fewer termini could be merged afterwards 719 (10080 -> 960 termini with pruning, against 480 without it on the same tomogram-scale mask). 720 Component filtering plus the per-component merge is both simpler and strictly better. 721 722 :func:`detect_junctions_skeleton` keys its terminus filter on exactly these points, so plotting 723 them is the way to see why a junction was accepted or rejected — that is what the napari widget's 724 **Show Crista Skeleton** option displays. 725 726 Args: 727 crista_mask: Binary crista segmentation (3D). 728 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 729 n_jobs: Forwarded to TEASAR's ``number_of_threads`` (-1/0 = all cores). 730 min_skeleton_nm: Skeleton components whose total length is below this (nm) are dropped as 731 specks, along with their termini — so a crista smaller than this cannot contribute a 732 terminus at all. 0 keeps every component TEASAR produced. 733 terminus_merge_nm: Termini within this distance (nm) of each other **on the same skeleton 734 component** collapse to one; 0 disables merging. This is what actually tames the terminus 735 count (4640 -> 960 on the mask above), since component filtering alone leaves the rim fan 736 intact. The value is the scale of the ragged fan at one crista end, a few voxels — **not** 737 the crista size. It is deliberately not larger even though larger radii look tidier: the 738 junction detector accepts a region only if it lies within ``terminus_nm`` of *some* 739 terminus, so collapsing a whole sheet rim to one point costs real junctions at the other 740 end of that rim. A 16 nm radius reaches the cosmetic ideal of two termini per lamella and 741 reduces a small sheet to a single terminus, which is exactly that failure. 742 return_edges: If True, also return the surviving edges, so a caller can draw the skeleton as 743 connected segments rather than a point cloud (the widget's Vectors layer). 744 745 Returns: 746 (vertices, is_terminus), or (vertices, is_terminus, edges) when ``return_edges`` is True. 747 ``vertices`` is (n, 3) in **nm**, array (z, y, x) order (divide by the voxel size for array 748 indices) and covers only the vertices of surviving components; ``is_terminus`` is the matching 749 boolean mask; ``edges`` is (m, 2) of indices into ``vertices``. All empty if the skeleton is 750 empty. 751 752 Raises: 753 ValueError: If the input is not 3D — TEASAR has no 2D implementation. 754 """ 755 if crista_mask.ndim != 3: 756 raise ValueError( 757 "compute_crista_skeleton requires a 3D volume (bioimage_cpp.skeleton.teasar has no 2D " 758 f"implementation), got ndim={crista_mask.ndim}." 759 ) 760 sampling = _to_sampling(voxel_size, 3) 761 vertices, edges, _ = teasar( 762 crista_mask.astype(bool), spacing=tuple(float(s) for s in sampling), 763 number_of_threads=_solver_threads(n_jobs), 764 ) 765 empty = (np.zeros((0, 3), dtype=float), np.zeros(0, dtype=bool)) 766 if len(vertices) == 0: 767 return (*empty, np.zeros((0, 2), dtype=int)) if return_edges else empty 768 769 vertices = np.asarray(vertices, dtype=float) 770 graph = _skeleton_graph(vertices, edges, min_skeleton_nm) 771 if graph.number_of_nodes() == 0: 772 return (*empty, np.zeros((0, 2), dtype=int)) if return_edges else empty 773 774 termini = set(_merge_termini(graph, vertices, [n for n in graph.nodes if graph.degree(n) == 1], 775 terminus_merge_nm)) 776 777 kept = sorted(graph.nodes) 778 remap = {node: index for index, node in enumerate(kept)} 779 out_vertices = vertices[kept] 780 is_terminus = np.array([node in termini for node in kept], dtype=bool) 781 if not return_edges: 782 return out_vertices, is_terminus 783 out_edges = np.array([[remap[a], remap[b]] for a, b in graph.edges], dtype=int).reshape(-1, 2) 784 return out_vertices, is_terminus, out_edges 785 786 787def _inner_surface_distance( 788 lumen_mask: np.ndarray, sampling: np.ndarray, n_jobs: int = 1 789) -> np.ndarray: 790 """Distance in nm to the **inner boundary membrane surface**, for localising junctions. 791 792 The membrane produced by :func:`approximate_membrane` is a *band* several nm thick, and 793 ``distance_transform(~band)`` is identically 0 on every voxel of it. That makes the band useless as 794 a reference for *where* inside itself a contact sits: a crista penetrating an 8 nm band has band 795 distance 0.0 across all of its in-band voxels, while its distance to the inner surface spreads over 796 1.5-7.5 nm. Measuring from the surface instead is what gives a junction a position and a size. 797 798 The surface used is the outermost layer of the lumen, i.e. the same single-wall geometry the 799 junction geodesics already run on (:func:`_open_trimmed_mesh`), so junction positions and junction 800 spacing finally refer to one surface. 801 802 Args: 803 lumen_mask: The eroded-mito interior from ``approximate_membrane(..., return_lumen=True)``. 804 sampling: Voxel size per axis in nm, as produced by :func:`_to_sampling`. 805 n_jobs: Thread budget for the distance transform (-1/0 = all cores). 806 807 Returns: 808 Float array of the same shape, the distance in nm to the nearest inner-surface voxel. 809 """ 810 inner = lumen_mask & ~binary_erosion(lumen_mask, border_value=1) 811 return distance_transform( 812 ~inner, sampling=sampling.tolist(), number_of_threads=_solver_threads(n_jobs) 813 ) 814 815 816def detect_junctions_skeleton( 817 crista_mask: np.ndarray, 818 membrane_mask: np.ndarray, 819 voxel_size: Union[float, Dict[str, float]], 820 max_extension_nm: float = 8.0, 821 min_extension_nm: float = 0.0, 822 terminus_nm: float = 20.0, 823 min_junction_volume_nm3: float = 50.0, 824 border_radius: int = 0, 825 boundary: Optional[np.ndarray] = None, 826 n_jobs: int = 1, 827 membrane_distance: Optional[np.ndarray] = None, 828 lumen_mask: Optional[np.ndarray] = None, 829 footprint_nm: Optional[float] = None, 830 min_skeleton_nm: Optional[float] = None, 831 terminus_merge_nm: Optional[float] = None, 832) -> Tuple[np.ndarray, Dict[str, float]]: 833 """Detect crista-membrane junctions as crista regions that reach close to the membrane. 834 835 The gap-tolerant alternative to :func:`detect_contact_sites`. Where the overlap detector requires 836 the crista mask to physically intersect the membrane band, this one accepts a crista region that 837 comes within ``max_extension_nm`` of it, so a crista segmented a few nm short of the inner boundary 838 membrane still registers. A TEASAR skeleton (``bioimage_cpp.skeleton.teasar``) then gates each 839 region on being near a crista *terminus*, which is what separates a crista ending at the membrane 840 from one running alongside it. 841 842 A junction is one 26-connected component of ``crista & (distance_to_membrane <= 843 max_extension_nm)``. Taking identity from the contact geometry this way is what makes the result 844 stable, and it is worth recording why, because the two obvious alternatives both fail: 845 846 - **Deriving identity from per-endpoint hits plus a merge radius does not work.** Measured on a real 847 mitochondrion with three junctions, no radius returns three: the count steps 10, 8, 4, 2 as the 848 radius grows, because a small radius fragments one junction while a large one fuses two that are 849 only 13.9 nm apart. 850 - **Extending each skeleton end along its tangent does not find the junctions.** On the same data 851 the direction from a skeleton end to its nearest real contact was 142-146 degrees away from that 852 end's tangent — pointing backwards — for all three. A crista-membrane contact is a *rim* feature 853 while a skeleton end is a *centerline* feature, and for a sheet meeting the membrane obliquely 854 their directions are unrelated. Widening the ray to a +/-60 degree cone changed nothing; only an 855 omnidirectional search found all three, which is a proximity test in disguise. That is this 856 function. 857 858 Because a component of the near-membrane mask is a subset of the crista mask, two disconnected 859 cristae can never be merged into one junction, and no tuning parameter governs that. The dual also 860 holds: the terminus filter is keyed per crista, so one crista can never be validated by another 861 one's end (see ``terminus_nm``). 862 863 **Known limitation — this mode over-detects on densely packed cristae.** Proximity is not the same 864 as junction: on a real mitochondrion with many cristae (TS_PS_01 mito 1, 0.8681 nm voxels, 8 nm 865 membrane) it reports 21 junctions of which only 2 involve any literal crista-membrane contact; the 866 other 19 are cristae merely passing within 8 nm of the inner boundary membrane. The false positives 867 are full-sized (up to ~1750 nm3), so ``min_junction_volume_nm3`` does not remove them. **Five** 868 discriminators have now been measured against real data and none separates the two populations: 869 region elongation (real junctions are 1.5-2.4 elongated too), region axis versus the membrane normal 870 (75-90 degrees for every region, because a contact patch spreads along the membrane by nature), 871 the rate at which membrane distance drops toward the terminus (0.42-0.73 for every region), a 872 minimum size, and the crista **sheet normal** versus the membrane normal (below). Treat the count as 873 an upper bound on a dense mitochondrion and inspect the result -- the widget's 874 **Show Crista Skeleton** layers exist for exactly that. Use ``"overlap"`` when only junctions with 875 actual contact should count. 876 877 **The sheet-normal discriminator was implemented, measured and removed.** The idea: a crista is a 878 lamella, so compare its sheet normal (the gradient of the mask smoothed at the sheet thickness) to 879 the membrane normal (the gradient of the reference distance field). A crista running parallel to the 880 membrane should give ``|cos| -> 1`` and one meeting it end-on ``|cos| -> 0``. On synthetic geometry 881 it behaved exactly so, 0.00 against 0.70. On the ``cutout_mito2`` cristae with three hand-verified 882 contacts it **inverted**: sampled at each region's closest-approach voxel the two contacting regions 883 scored 0.786 and 0.881 while the two non-contacting ones scored 0.000 and 0.766, so a threshold of 884 0.7 removed all three real junctions and kept the false positive. Region means do not separate 885 either (0.582/0.606 against 0.558/0.842). Part of the reason is that the measure is ill-posed 886 precisely where one wants to sample it: at a closest-approach voxel the distance field can be 887 locally flat, and a vanishing gradient normalises to a meaningless direction -- the exact 0.000 888 above is that artefact. Reproduce with ``scripts/cooper/measure_terminus_alignment.py``. Do not 889 re-propose it without new evidence from that script. 890 891 Args: 892 crista_mask: Binary crista segmentation. Must be 3D — TEASAR has no 2D implementation. 893 membrane_mask: Binary mitochondrial membrane mask, e.g. from :func:`approximate_membrane`. 894 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 895 max_extension_nm: How far (nm) a crista may fall short of the membrane and still count. Set it 896 to about one membrane thickness; an over-large value starts flagging cristae that merely 897 pass near the boundary, and it also fuses junctions whose near-membrane regions touch. 898 min_extension_nm: Smallest gap (nm) that counts, measured as a region's *closest* approach. 899 0 (the default) accepts regions already overlapping the membrane; raise it to isolate 900 non-overlapping junctions only. 901 terminus_nm: A region is kept only if it lies within this distance (nm) of a skeleton end **of 902 its own 26-connected crista**, so a crista running alongside the membrane is rejected. 903 Keying on the region's own crista is what stops a neighbouring crista's end from validating 904 it; with a single pooled terminus tree a region whose own skeleton was dropped by 905 ``min_skeleton_nm`` could still pass. On real data this changes nothing — measured 4 -> 4 906 on ``cutout_mito2``, whose crista mask merges into just 5 connected components — so do not 907 "simplify" it back out on the grounds that it makes no difference there. The leak needs a 908 genuinely isolated crista, which sparse and synthetic masks readily produce; see 909 ``test_a_crista_cannot_borrow_another_cristas_terminus``. The 20 nm default is generous by design: 910 the measured terminus distances of real junctions were 0.0, 0.0 and 8.6 nm, so it rejects 911 only cristae running well clear of any end, such as one sliding along the membrane. Pass 912 ``inf`` to disable the filter and skip TEASAR entirely. 913 min_junction_volume_nm3: Regions smaller than this are dropped. This removes specks only — real 914 data produces 1-, 3- and 12-voxel regions that are not junctions on any reading, while real 915 junctions measure hundreds to thousands of nm^3. It is **not** a remedy for the 916 over-detection described above, whose false positives are full-sized. 917 border_radius: Width in voxels of the volume-border zone to exclude, matching the membrane's 918 own border-gap trim (``_gap_radius``). 0 (the default) disables the exclusion and 919 reproduces the misaligned behaviour, so **every production caller passes it**: a crista 920 within this distance of a clipped face would otherwise be matched against membrane that 921 :func:`approximate_membrane` deliberately removed as unknown — the distance transform would 922 measure straight across the deleted region to the nearest surviving membrane voxel and 923 assert a junction against a membrane it was told nothing about. Regions are trimmed rather 924 than discarded, so a crista entering the unknown zone still counts wherever else it 925 genuinely reaches the membrane. May also be an ``(ndim, 2)`` array of **per-face** radii, 926 which is what a caller passing a bounding-box crop must do — the zone is a global fact and 927 a crop starting one voxel inside the volume still overlaps it. See :func:`_border_zone`. 928 boundary: Optional ``(ndim, 2)`` bool forwarded to :func:`_border_zone`, zeroing the radius on 929 faces that are not genuine volume faces; None means all faces. It is the *mesh* question, 930 and per-face ``border_radius`` already encodes it, so **pass one or the other, never 931 both** — the radii are multiplied by this mask. 932 n_jobs: Forwarded to TEASAR's ``number_of_threads`` (-1/0 = all cores). 933 membrane_distance: Optional precomputed ``distance_transform(~membrane_mask)`` in nm, to avoid 934 recomputing it when the caller already has one (see :func:`_single_mito_row`). Must match 935 ``membrane_mask`` and use the same sampling. Used only as the fallback reference when no 936 ``lumen_mask`` is given. 937 lumen_mask: The eroded-mito interior from ``approximate_membrane(..., return_lumen=True)``. 938 When given, distances are measured to the **inner boundary membrane surface** derived from 939 it rather than to the membrane band, which is what lets a junction be localised at all — 940 see :func:`_inner_surface_distance`. Pass only a genuine lumen: ``mito & ~membrane`` is not 941 one near a clipped face, where it re-includes the mito's outer shell. 942 footprint_nm: How far (nm) beyond a region's closest approach the painted label extends. This 943 sets the junction label's thickness and therefore its centroid, but never the count. None 944 (the default) means **one voxel diagonal**, which cannot be written as a literal here since 945 it depends on ``voxel_size``: it is the tightest tolerance that still captures a contact 946 patch lying oblique to the grid. Measured at 1.5 nm isotropic voxels it keeps the footprint 947 at 6% of the candidate region, where 4 nm would already take 56% of it. 948 min_skeleton_nm: Forwarded to :func:`compute_crista_skeleton` — skeleton components shorter 949 than this are dropped, so a crista smaller than this contributes no terminus and cannot 950 score a junction (which holds only because the terminus trees are per-crista; see 951 ``terminus_nm``). None leaves that function's own default in force. 952 terminus_merge_nm: Forwarded to :func:`compute_crista_skeleton` — the radius within which 953 termini on one component collapse to a single representative. None leaves that function's 954 own default in force. 955 956 Returns: 957 contact_labels: Integer array (same shape as the input) where each junction has a unique ID 958 (0 = background, 1..n = junctions) — interchangeable with the first return value of 959 :func:`detect_contact_sites`, so it feeds :func:`compute_junction_distances` and the napari 960 junction layer unchanged. The labelled region is the junction's **closest-approach 961 footprint**: the voxels of the candidate region within ``footprint_nm`` of its nearest 962 approach to the reference surface, not the whole region. It is therefore a thin patch at 963 the membrane rather than a slab of crista (measured: 80 voxels against 1360 for the region), 964 it lies inside the crista and hence inside the mitochondrion, and its centroid is a usable 965 junction position for the geodesic stage. 966 summary: contact_voxel_count, crista_junction_count, contact_volume_nm3, 967 mean_junction_extension_nm (the mean over junctions of each region's closest approach to 968 the reference surface; 0 where the crista reaches it). Only ``crista_junction_count`` and 969 ``mean_junction_extension_nm`` are mode-specific; ``contact_voxel_count`` and 970 ``contact_volume_nm3`` still describe the genuine crista-membrane overlap, so those two 971 columns mean the same thing in both modes. 972 973 Raises: 974 ValueError: If the input is not 3D, the voxel size is not positive on every axis, or the 975 extension range is negative / inverted. 976 """ 977 if crista_mask.ndim != 3: 978 raise ValueError( 979 "detect_junctions_skeleton requires a 3D volume (bioimage_cpp.skeleton.teasar has no 2D " 980 f"implementation), got ndim={crista_mask.ndim}. Use junction_mode='overlap' for 2D data." 981 ) 982 if min_extension_nm < 0 or max_extension_nm < min_extension_nm: 983 raise ValueError( 984 "need 0 <= min_extension_nm <= max_extension_nm, got " 985 f"min_extension_nm={min_extension_nm}, max_extension_nm={max_extension_nm}" 986 ) 987 988 sampling = _to_sampling(voxel_size, 3) 989 if not np.all(sampling > 0): 990 raise ValueError(f"voxel_size must be positive on every axis, got {voxel_size!r}") 991 voxel_vol = float(np.prod(sampling)) 992 crista = crista_mask.astype(bool) 993 membrane = membrane_mask.astype(bool) 994 995 overlap_voxels = int(np.count_nonzero(crista & membrane)) 996 summary = { 997 "contact_voxel_count": overlap_voxels, 998 "crista_junction_count": 0, 999 "contact_volume_nm3": float(overlap_voxels) * voxel_vol, 1000 "mean_junction_extension_nm": np.nan, 1001 } 1002 labels = np.zeros(crista.shape, dtype=np.int32) 1003 if not crista.any() or not membrane.any(): 1004 return labels, summary 1005 1006 if membrane_distance is None: 1007 membrane_distance = distance_transform( 1008 ~membrane, sampling=sampling.tolist(), number_of_threads=_solver_threads(n_jobs) 1009 ) 1010 reference_distance = membrane_distance 1011 if lumen_mask is not None and lumen_mask.any(): 1012 reference_distance = _inner_surface_distance(lumen_mask.astype(bool), sampling, n_jobs) 1013 if footprint_nm is None: 1014 footprint_nm = float(np.linalg.norm(sampling)) 1015 1016 near = crista & (reference_distance <= max_extension_nm) 1017 if np.any(np.asarray(border_radius) >= 1): 1018 near &= ~_border_zone(crista.shape, border_radius, boundary) 1019 if not near.any(): 1020 return labels, summary 1021 1022 region_labels, n_regions = ndimage_label(near, structure=np.ones(3 * (3,), dtype=bool)) 1023 1024 # One KD-tree PER CRISTA, not one for all of them: a region is validated only by a terminus of 1025 # its own 26-connected crista. A single pooled tree let an unrelated crista ending nearby accept a 1026 # region whose own skeleton was dropped by ``min_skeleton_nm`` — measured on the test lamella, a 1027 # sheet alone gave 1 junction, an 18-voxel blob alone 0, and the two together 2. 1028 terminus_trees = None 1029 if np.isfinite(terminus_nm): 1030 vertices, is_terminus = compute_crista_skeleton( 1031 crista, voxel_size, n_jobs=n_jobs, 1032 **_given(min_skeleton_nm=min_skeleton_nm, terminus_merge_nm=terminus_merge_nm), 1033 ) 1034 endpoints = vertices[is_terminus] 1035 if len(endpoints) == 0: 1036 return labels, summary 1037 from scipy.spatial import cKDTree 1038 1039 crista_components, _ = ndimage_label(crista, structure=np.ones(3 * (3,), dtype=bool)) 1040 # TEASAR returns physical coordinates, i.e. index * spacing exactly (measured residual 0.0 at 1041 # isotropic and 3.6e-15 at anisotropic spacing), so rounding recovers the index and always 1042 # lands on a foreground voxel — a terminus can never be attributed to component 0. 1043 endpoint_component = crista_components[tuple(np.rint(endpoints / sampling).astype(int).T)] 1044 terminus_trees = { 1045 int(component): cKDTree(endpoints[endpoint_component == component]) 1046 for component in np.unique(endpoint_component) 1047 } 1048 # ``near`` is a subset of ``crista`` under the same connectivity, so each region lies in 1049 # exactly one component. Collapse the lookup to one entry per region and free the int32 1050 # volume (120 MB on a 401x257x290 tomogram) before the loop. 1051 region_component = np.zeros(n_regions + 1, dtype=np.int32) 1052 region_component[region_labels[near]] = crista_components[near] 1053 del crista_components 1054 1055 assigned = 0 1056 gaps = [] 1057 for region in range(1, n_regions + 1): 1058 region_mask = region_labels == region 1059 if float(np.count_nonzero(region_mask)) * voxel_vol < min_junction_volume_nm3: 1060 continue 1061 region_distance = reference_distance[region_mask] 1062 gap = float(region_distance.min()) 1063 if not (min_extension_nm <= gap <= max_extension_nm): 1064 continue 1065 if terminus_trees is not None: 1066 tree = terminus_trees.get(int(region_component[region])) 1067 if tree is None: # this crista contributed no terminus of its own 1068 continue 1069 # Capped query: past terminus_nm the KD-tree returns inf rather than a real distance. 1070 points = np.argwhere(region_mask) * sampling 1071 hit = tree.query(points, distance_upper_bound=terminus_nm)[0] 1072 if not np.isfinite(hit).any(): 1073 continue 1074 assigned += 1 1075 gaps.append(gap) 1076 labels[region_mask & (reference_distance <= gap + footprint_nm)] = assigned 1077 1078 summary["crista_junction_count"] = assigned 1079 if gaps: 1080 summary["mean_junction_extension_nm"] = float(np.mean(gaps)) 1081 return labels, summary 1082 1083 1084def detect_junctions( 1085 crista_mask: np.ndarray, 1086 membrane_mask: np.ndarray, 1087 voxel_size: Union[float, Dict[str, float]], 1088 junction_mode: str = "overlap", 1089 max_extension_nm: float = 8.0, 1090 n_jobs: int = 1, 1091 **kwargs, 1092) -> Tuple[np.ndarray, Dict[str, float]]: 1093 """Detect crista-membrane junctions with the chosen algorithm. 1094 1095 Args: 1096 crista_mask: Binary crista segmentation. 1097 membrane_mask: Binary mitochondrial membrane mask. 1098 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1099 junction_mode: ``"overlap"`` (default) counts connected components of the direct 1100 crista-membrane intersection via :func:`detect_contact_sites`; ``"skeleton"`` counts 1101 crista regions that come within ``max_extension_nm`` of the membrane near a crista 1102 terminus, via :func:`detect_junctions_skeleton`. 1103 max_extension_nm: How far (nm) a crista may fall short of the membrane; ``"skeleton"`` only. 1104 n_jobs: Thread budget; ``"skeleton"`` only. 1105 **kwargs: Further :func:`detect_junctions_skeleton` options (``min_extension_nm``, 1106 ``terminus_nm``, ``min_junction_volume_nm3``, ``border_radius``, ``boundary``, 1107 ``membrane_distance``, ``lumen_mask``, ``footprint_nm``, 1108 ``min_skeleton_nm``, ``terminus_merge_nm``); ignored in ``"overlap"`` mode, which needs no 1109 border parameter because it is border-safe by construction. Any of them passed as ``None`` 1110 is dropped here (:func:`_given`), so a caller that treats ``None`` as "unset" — the widget's 1111 zeroed spin-boxes, the CLI's unset flags, the per-mito path — gets the detector's own 1112 default rather than having to restate it. 1113 1114 Returns: 1115 (contact_labels, summary) as documented on the two backends. ``summary`` always carries 1116 ``mean_junction_extension_nm`` (NaN in ``"overlap"`` mode, which measures no gap) so callers 1117 can read it without branching on the mode. 1118 1119 Raises: 1120 ValueError: If ``junction_mode`` is not one of the two supported values. 1121 """ 1122 if junction_mode == "overlap": 1123 contact_labels, summary = detect_contact_sites(crista_mask, membrane_mask, voxel_size) 1124 summary["mean_junction_extension_nm"] = np.nan 1125 return contact_labels, summary 1126 if junction_mode == "skeleton": 1127 return detect_junctions_skeleton( 1128 crista_mask, membrane_mask, voxel_size, 1129 max_extension_nm=max_extension_nm, n_jobs=n_jobs, **_given(**kwargs), 1130 ) 1131 raise ValueError(f"junction_mode must be 'overlap' or 'skeleton', got {junction_mode!r}") 1132 1133 1134 1135def detect_junctions_per_mito( 1136 crista_mask: np.ndarray, 1137 mito_segmentation: np.ndarray, 1138 membrane_mask: np.ndarray, 1139 voxel_size: Union[float, Dict[str, float]], 1140 lumen_mask: Optional[np.ndarray] = None, 1141 border_radius: int = 0, 1142 n_jobs: int = 1, 1143 **kwargs, 1144) -> Tuple[np.ndarray, Dict[str, float]]: 1145 """Detect junctions the way the statistics table does: **per mitochondrial instance**. 1146 1147 :func:`detect_junctions` run once over a whole volume and the same function run per instance are 1148 not the same measurement, and in ``"skeleton"`` mode they routinely disagree. Restricting a crista 1149 to one mitochondrion clips it, and clipping moves its skeleton endpoints: a crista tube running 1150 through a mitochondrion and out the other side has its global endpoints far away from the membrane 1151 (no junction), while the clipped tube ends *at* the mitochondrial boundary (two junctions). Any 1152 caller that displays junctions next to numbers from 1153 :func:`compute_mito_crista_statistics` must therefore detect them this way, or the picture and the 1154 table disagree at the same settings. 1155 1156 It is also **cheaper** than the global pass, not more expensive, which is unintuitive enough to be 1157 worth recording: the cost is dominated by two whole-volume distance transforms, and mitochondrial 1158 bounding boxes sum to a fraction of a tomogram (measured 0.06 s against 0.33 s, a 0.19x ratio, on 1159 4 box mitos whose bboxes covered 14% of an 80x160x160 volume — same 30 junctions either way). 1160 TEASAR is per connected component either way, so splitting the mask by instance multiplies no 1161 skeletonisation work. 1162 1163 The per-instance masking and the crop-aware border radii here mirror :func:`_single_mito_row` 1164 exactly; ``test_per_mito_detector_matches_the_table`` pins the two together. 1165 1166 Args: 1167 crista_mask: Binary crista segmentation (global volume). 1168 mito_segmentation: Instance label array (background = 0). 1169 membrane_mask: Binary mitochondrial membrane mask, e.g. from :func:`approximate_membrane`. 1170 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1171 lumen_mask: Optional eroded-mito interior matching ``membrane_mask``; restricted to each 1172 instance, as in :func:`_single_mito_row`. 1173 border_radius: Volume-border zone width in voxels (``_gap_radius``). Converted to crop-aware 1174 per-face radii for each bounding box — see :func:`_border_zone`. 1175 n_jobs: Thread budget forwarded to the detector. 1176 **kwargs: Further :func:`detect_junctions` options (``junction_mode``, ``max_extension_nm``, 1177 ``terminus_nm``, ...). 1178 1179 Returns: 1180 (junction_labels, summary) with the same keys as :func:`detect_junctions`. Ids are 1181 per-instance and therefore repeat between mitochondria; the painted voxels do not, since each 1182 lies inside its own mito. ``mean_junction_extension_nm`` is averaged over the instances that 1183 reported one. 1184 """ 1185 ndim = mito_segmentation.ndim 1186 crista_binary = crista_mask.astype(bool) 1187 membrane_binary = membrane_mask.astype(bool) 1188 vol_shape = mito_segmentation.shape 1189 1190 labels = np.zeros(vol_shape, dtype=np.int32) 1191 summary = { 1192 "contact_voxel_count": 0, "crista_junction_count": 0, "contact_volume_nm3": 0.0, 1193 "mean_junction_extension_nm": np.nan, 1194 } 1195 extensions = [] 1196 1197 for prop in regionprops(mito_segmentation): 1198 bbox = prop.bbox 1199 slices = tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim)) 1200 mito_local = mito_segmentation[slices] == prop.label 1201 crista_local = crista_binary[slices] & mito_local 1202 membrane_local = membrane_binary[slices] & mito_local 1203 if not crista_local.any() or not membrane_local.any(): 1204 continue 1205 border_radii = np.array( 1206 [[max(0, border_radius - bbox[a]), 1207 max(0, border_radius - (vol_shape[a] - bbox[a + ndim]))] for a in range(ndim)], dtype=int 1208 ) 1209 local_labels, local_summary = detect_junctions( 1210 crista_local, membrane_local, voxel_size, n_jobs=n_jobs, 1211 border_radius=border_radii, 1212 lumen_mask=None if lumen_mask is None else lumen_mask[slices] & mito_local, 1213 **kwargs, 1214 ) 1215 written = local_labels > 0 1216 labels[slices][written] = local_labels[written] # masked: bounding boxes overlap 1217 summary["contact_voxel_count"] += local_summary["contact_voxel_count"] 1218 summary["contact_volume_nm3"] += local_summary["contact_volume_nm3"] 1219 summary["crista_junction_count"] += local_summary["crista_junction_count"] 1220 if np.isfinite(local_summary["mean_junction_extension_nm"]): 1221 extensions.append(local_summary["mean_junction_extension_nm"]) 1222 1223 if extensions: 1224 summary["mean_junction_extension_nm"] = float(np.mean(extensions)) 1225 return labels, summary 1226 1227 1228_JUNCTION_DISTANCE_NAN = { 1229 "junction_count": 0, 1230 "mean_nn_junction_distance_nm": np.nan, 1231 "median_nn_junction_distance_nm": np.nan, 1232 "junction_clustering_index": np.nan, 1233} 1234 1235 1236def _junction_matrix_mesh( 1237 centroids: np.ndarray, 1238 sampling: np.ndarray, 1239 vertices: np.ndarray, 1240 faces: np.ndarray, 1241 n_jobs: int = 1, 1242) -> Optional[np.ndarray]: 1243 """Pairwise junction geodesic distances (nm) along a triangle-mesh surface. 1244 1245 Snaps each junction centroid to its nearest mesh vertex and calls 1246 ``bioimage_cpp.distance.geodesic_distances_mesh`` (passing all junction vertices as sources, so it 1247 returns the full pairwise matrix directly). The mesh (see :func:`_surface_mesh`) returns vertices 1248 in the unpadded mask index frame, so voxel centroids map to it as ``index * sampling`` with no 1249 offset. Returns None (junction distances become NaN) when the mesh is empty. 1250 1251 Args: 1252 centroids: (n, ndim) junction centroids in voxel coordinates. 1253 sampling: Voxel size per axis (nm), array (z, y, x) order. 1254 vertices: Mesh vertices (n_vertices, 3) in nm (unpadded mask frame). 1255 faces: Mesh triangle indices (n_faces, 3). 1256 n_jobs: Forwarded to the C++ solver's ``number_of_threads``: -1/0 map to 0 (the solver's 1257 "use hardware_concurrency"), otherwise that many threads. 1258 1259 Returns: 1260 (n, n) geodesic distance matrix in nm (0 diagonal, NaN for disconnected pairs), or None. 1261 Disconnected pairs come back from the solver as ``+inf`` and are converted to NaN. 1262 """ 1263 verts = np.ascontiguousarray(vertices, dtype=np.float64) 1264 tris = np.ascontiguousarray(faces, dtype=np.int64) 1265 if verts.shape[0] == 0 or tris.shape[0] == 0: 1266 return None 1267 from scipy.spatial import cKDTree 1268 1269 points = np.asarray(centroids, dtype=float) * sampling 1270 _, vertex_ids = cKDTree(verts).query(points) 1271 vertex_ids = np.atleast_1d(np.asarray(vertex_ids, dtype=np.int64)) 1272 dm = np.asarray( 1273 geodesic_distances_mesh( 1274 verts, tris, vertex_ids, number_of_threads=_solver_threads(n_jobs) 1275 ), 1276 dtype=float, 1277 ) 1278 dm[~np.isfinite(dm)] = np.nan 1279 np.fill_diagonal(dm, 0.0) 1280 return dm 1281 1282 1283def compute_junction_distances( 1284 contact_labels: np.ndarray, 1285 membrane_mask: np.ndarray, 1286 voxel_size: Union[float, Dict[str, float]], 1287 surface_area_nm2: Optional[float] = None, 1288 n_jobs: int = 1, 1289 mesh_vertices: Optional[np.ndarray] = None, 1290 mesh_faces: Optional[np.ndarray] = None, 1291) -> Tuple[np.ndarray, Dict[str, float]]: 1292 """Geodesic distances between crista-membrane junctions along the eroded-mito surface mesh. 1293 1294 Each junction (a connected component in ``contact_labels``) is reduced to its centroid, snapped to 1295 the nearest vertex of a triangle mesh, and pairwise surface geodesics are computed with 1296 ``bioimage_cpp.distance.geodesic_distances_mesh``. The mesh is the **eroded-mito (lumen) surface** 1297 passed in as ``mesh_vertices``/``mesh_faces`` by :func:`_single_mito_row` (a clean, single-wall 1298 surface at the membrane's inner edge). If no mesh is supplied — or no usable surface mesh exists 1299 (empty membrane / degenerate mesh) — the junction distances are NaN. (There is no membrane-band 1300 fallback mesh: the metric is defined on the lumen surface, and meshing the thick membrane band 1301 would give a different, capped double-wall surface.) 1302 1303 A Clark-Evans nearest-neighbour index summarises whether the junctions are clustered. 1304 1305 Args: 1306 contact_labels: Integer junction label array (0 = background, 1..n = junctions), 1307 e.g. the first return value of :func:`detect_contact_sites`. 1308 membrane_mask: Binary mitochondrial membrane mask the junctions sit on. Only used for the 1309 empty-membrane early-out (no membrane → NaN); it is not meshed. 1310 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1311 surface_area_nm2: Membrane/mito surface area used as the reference area for the 1312 Clark-Evans expectation. If None or non-positive, the clustering index is NaN. 1313 n_jobs: 1 = serial, -1 = all cores (forwarded to the mesh solver's thread count). 1314 mesh_vertices: Optional (n_vertices, 3) mesh vertices in nm (unpadded mask frame; see 1315 :func:`_surface_mesh`) — the eroded-mito (lumen) surface. If omitted, the junction 1316 distances are NaN. 1317 mesh_faces: Optional (n_faces, 3) triangle indices matching ``mesh_vertices``. 1318 1319 Returns: 1320 distance_matrix: (n, n) geodesic distances in nm between junctions; the diagonal is 1321 0 and unreachable pairs (disconnected fragments) are NaN. Empty for fewer than two 1322 junctions. 1323 summary: junction_count, mean_nn_junction_distance_nm, median_nn_junction_distance_nm, 1324 junction_clustering_index (Clark-Evans R: < 1 clustered, ~ 1 random, > 1 dispersed). 1325 1326 Notes: 1327 The Clark-Evans expected nearest-neighbour distance uses the standard 2-D planar 1328 approximation ``0.5 * sqrt(A / n)`` with ``A = surface_area_nm2``. 1329 1330 Each junction's nearest-neighbour distance is the smallest distance to a *reachable* other 1331 junction: self (diagonal) and unreachable (NaN) pairs are set to +inf before the per-row 1332 minimum, and rows with no reachable neighbour (min stays +inf) are dropped. 1333 """ 1334 ndim = contact_labels.ndim 1335 sampling = _to_sampling(voxel_size, ndim) 1336 membrane = membrane_mask.astype(bool) 1337 1338 labels = [lbl for lbl in np.unique(contact_labels) if lbl != 0] 1339 n = len(labels) 1340 if n < 2 or not membrane.any(): 1341 summary = dict(_JUNCTION_DISTANCE_NAN) 1342 summary["junction_count"] = n 1343 return np.zeros((n, n), dtype=float), summary 1344 1345 centroids = np.atleast_2d( 1346 np.asarray(center_of_mass(contact_labels > 0, labels=contact_labels, index=labels), dtype=float) 1347 ) 1348 1349 if mesh_vertices is not None and mesh_faces is not None and len(mesh_faces) > 0: 1350 distance_matrix = _junction_matrix_mesh(centroids, sampling, mesh_vertices, mesh_faces, n_jobs) 1351 else: 1352 distance_matrix = None 1353 1354 if distance_matrix is None: 1355 summary = dict(_JUNCTION_DISTANCE_NAN) 1356 summary["junction_count"] = n 1357 return np.full((n, n), np.nan, dtype=float), summary 1358 1359 dm = distance_matrix.copy() 1360 np.fill_diagonal(dm, np.inf) 1361 dm[~np.isfinite(dm)] = np.inf 1362 row_min = dm.min(axis=1) 1363 nn_distances = row_min[np.isfinite(row_min)] 1364 1365 mean_nn = float(np.mean(nn_distances)) if nn_distances.size else np.nan 1366 median_nn = float(np.median(nn_distances)) if nn_distances.size else np.nan 1367 1368 clustering_index = np.nan 1369 if surface_area_nm2 is not None and surface_area_nm2 > 0 and np.isfinite(mean_nn): 1370 expected_nn = 0.5 * np.sqrt(float(surface_area_nm2) / n) 1371 if expected_nn > 0: 1372 clustering_index = mean_nn / expected_nn 1373 1374 return distance_matrix, { 1375 "junction_count": n, 1376 "mean_nn_junction_distance_nm": mean_nn, 1377 "median_nn_junction_distance_nm": median_nn, 1378 "junction_clustering_index": clustering_index, 1379 } 1380 1381 1382# --------------------------------------------------------------------------- 1383# Morphology 1384# --------------------------------------------------------------------------- 1385 1386def compute_crista_morphology( 1387 crista_mask: np.ndarray, 1388 voxel_size: Union[float, Dict[str, float]], 1389 method: str = "both", 1390) -> Dict[str, float]: 1391 """Compute crista shape metrics from binary mask. 1392 1393 Args: 1394 crista_mask: Binary crista segmentation. 1395 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1396 method: "area" | "medial_axis" | "both". 1397 1398 Returns: 1399 Dict with cristae_surface_area_nm2 (area/both) and avg_thickness_nm (medial_axis/both). 1400 avg_thickness_nm is 2 × mean distance-transform value at skeleton voxels. 1401 """ 1402 if method not in ("area", "medial_axis", "both"): 1403 raise ValueError(f"method must be 'area', 'medial_axis', or 'both', got {method!r}") 1404 1405 sampling = _to_sampling(voxel_size, crista_mask.ndim) 1406 result: Dict[str, float] = {} 1407 1408 if method in ("area", "both"): 1409 result["cristae_surface_area_nm2"] = _surface_area(crista_mask, sampling) 1410 1411 if method in ("medial_axis", "both"): 1412 result["avg_thickness_nm"] = _medial_axis_thickness_nm(crista_mask, sampling) 1413 1414 return result 1415 1416 1417# --------------------------------------------------------------------------- 1418# Per-mitochondrion statistics 1419# --------------------------------------------------------------------------- 1420 1421def _single_mito_row( 1422 label: int, 1423 bbox: tuple, 1424 mito_crop: np.ndarray, 1425 crista_crop: np.ndarray, 1426 membrane_crop: np.ndarray, 1427 voxel_size: Union[float, Dict[str, float]], 1428 sampling: np.ndarray, 1429 voxel_vol: float, 1430 vol_shape: tuple, 1431 border_radius: int, 1432 method: str = "skip", 1433 inner_n_jobs: int = 1, 1434 lumen_crop: Optional[np.ndarray] = None, 1435 junction_mode: str = "overlap", 1436 max_extension_nm: float = 8.0, 1437 terminus_nm: Optional[float] = None, 1438 min_junction_volume_nm3: Optional[float] = None, 1439 footprint_nm: Optional[float] = None, 1440 min_skeleton_nm: Optional[float] = None, 1441 terminus_merge_nm: Optional[float] = None, 1442 junction_out: Optional[np.ndarray] = None, 1443) -> Dict[str, float]: 1444 """Compute the statistics row for a single mitochondrion instance. 1445 1446 Factored out of :func:`compute_mito_crista_statistics` so the per-mito work (which is 1447 independent between instances) can be parallelised. Takes the bbox-cropped label arrays 1448 (``mito_crop`` = label array cropped to ``bbox``; ``crista_crop``/``membrane_crop`` the 1449 matching binary crops), so a worker only touches its own mito's bbox region rather than the 1450 whole volume. ``inner_n_jobs`` is forwarded to the (parallelisable) junction-distance stage. 1451 1452 ``method`` controls only the crista orientation anisotropy — every other metric (marching-cubes 1453 surface areas, geodesic junction distances, EDT proximity/thickness) is computed identically for 1454 all modes. ``"skip"`` (the default) leaves the anisotropy NaN and is the fastest; ``"fast"`` 1455 computes it on a 2× downsampled crop (~8× cheaper, a *relative* indicator only, not comparable to 1456 exact); ``"exact"`` uses the full-resolution structure tensor (the dominant cost). 1457 1458 ``junction_mode`` selects the junction detector (see :func:`detect_junctions`), with 1459 ``max_extension_nm`` its gap tolerance and ``terminus_nm`` its crista-terminus filter; these only 1460 affect ``crista_junction_count`` and ``mean_junction_extension_nm``. The junction *distances* are 1461 computed the same way in either mode, since both detectors return the same kind of labeled 1462 junction array. ``membrane_distance`` is computed once here and shared between the junction 1463 detector and the crista-proximity stage. 1464 1465 ``lumen_crop`` is the (optional) bbox-cropped eroded lumen from :func:`approximate_membrane` 1466 (``return_lumen=True``) used for the junction geodesic mesh and for the inner boundary membrane 1467 area in ``imm_surface_area_nm2``; without it the mesh falls back to 1468 ``mito_local & ~membrane_local`` (used only when a caller supplies their own membrane, and 1469 contaminated near clipped faces). That fallback is deliberately **not** handed to the junction 1470 detector as ``lumen_mask``: near a clipped face :func:`approximate_membrane` has deleted the 1471 membrane, so ``mito & ~membrane`` re-includes the mito's *outer* shell there, which would put the 1472 inner boundary membrane on the outside of the mitochondrion. Without a real lumen the detector is 1473 given ``None`` and falls back to the membrane band, which is border-safe. For the same reason the 1474 fallback is only taken when a membrane band exists at all: with an empty band it would degenerate 1475 to ``mito_local``, reporting the *outer* surface as the inner boundary membrane, so the lumen is 1476 left ``None`` and ``imm_surface_area_nm2`` is NaN instead. 1477 Both the lumen geodesic mesh and the mito outer-surface-area mesh 1478 are trimmed/opened at faces where the mito is clipped by the volume boundary: the lumen via 1479 :func:`_open_trimmed_mesh` (trimmed to the certain region and left open, so geodesics do not 1480 shortcut across a cap), the mito surface via ``closed_faces`` (open, so a fabricated cap is not 1481 counted as membrane area). The membrane distance transform is freed before the (memory-heavy) 1482 orientation stage to cap peak memory. 1483 """ 1484 ndim = mito_crop.ndim 1485 touches_border = any( 1486 bbox[i] < border_radius or bbox[i + ndim] > vol_shape[i] - border_radius 1487 for i in range(ndim) 1488 ) 1489 1490 mito_local = mito_crop == label 1491 crista_local = crista_crop & mito_local 1492 membrane_local = membrane_crop & mito_local 1493 1494 mito_vol = float(mito_local.sum()) * voxel_vol 1495 crista_vol = float(crista_local.sum()) * voxel_vol 1496 1497 # Three DIFFERENT border questions are asked below; they need three different answers, and 1498 # conflating the last two is what made the junction exclusion silently inert on an offset crop. 1499 # 1. Does this mito enter the unknown zone at all? -> `touches_border` above, global-aware. 1500 # 2. Is this bbox face a genuine volume face? -> `boundary`, an exact-face test. Right for 1501 # the meshes: it asks whether the object is CLIPPED here, and a mito whose bbox starts one 1502 # voxel in genuinely ends there, so its surface is real and must be closed, not trimmed. 1503 # 3. Where INSIDE this crop is the membrane unknown? -> `border_radii`. The zone is a global 1504 # fact (`approximate_membrane` trims it on the whole volume), so a crop starting one voxel 1505 # inside the volume still overlaps it while `boundary` reports no volume face at all. 1506 boundary = np.array( 1507 [[bbox[a] == 0, bbox[a + ndim] == vol_shape[a]] for a in range(ndim)], dtype=bool 1508 ) 1509 border_radii = np.array( 1510 [[max(0, border_radius - bbox[a]), 1511 max(0, border_radius - (vol_shape[a] - bbox[a + ndim]))] for a in range(ndim)], dtype=int 1512 ) 1513 1514 has_crista = crista_local.any() 1515 has_membrane = membrane_local.any() 1516 mito_surface = _surface_area(mito_local, sampling, closed_faces=~boundary) 1517 1518 if lumen_crop is not None: 1519 lumen_local = lumen_crop & mito_local 1520 elif has_membrane: 1521 lumen_local = mito_local & ~membrane_local 1522 else: 1523 lumen_local = None # the fallback would be mito_local itself, i.e. the OUTER surface. 1524 # ponytail: the lumen spans flat across each crista junction mouth (cristae are never subtracted 1525 # from it), so the inner boundary membrane area is over-counted by ~one small disk per junction. 1526 lumen_mesh = ( 1527 _open_trimmed_mesh(lumen_local, sampling, border_radius, boundary) 1528 if lumen_local is not None else None 1529 ) 1530 mesh_verts, mesh_faces = lumen_mesh if lumen_mesh is not None else (None, None) 1531 ibm_surface = float(mesh_surface_area(*lumen_mesh)) if lumen_mesh is not None else np.nan 1532 1533 if has_crista and has_membrane: 1534 membrane_distance = distance_transform(~membrane_local, sampling=sampling.tolist(), number_of_threads=1) 1535 contact_labels_local, contact_summary = detect_junctions( 1536 crista_local, membrane_local, voxel_size, 1537 junction_mode=junction_mode, max_extension_nm=max_extension_nm, 1538 terminus_nm=terminus_nm, min_junction_volume_nm3=min_junction_volume_nm3, 1539 # per-face radii already encode `boundary`; passing both would zero them where they matter 1540 border_radius=border_radii, 1541 n_jobs=inner_n_jobs, membrane_distance=membrane_distance, lumen_mask=lumen_local, 1542 footprint_nm=footprint_nm, 1543 min_skeleton_nm=min_skeleton_nm, terminus_merge_nm=terminus_merge_nm, 1544 ) 1545 _, proximity = compute_crista_proximity( 1546 crista_local, membrane_local, voxel_size, membrane_distance=membrane_distance 1547 ) 1548 _, junction_dist = compute_junction_distances( 1549 contact_labels_local, membrane_local, voxel_size, 1550 surface_area_nm2=mito_surface, n_jobs=inner_n_jobs, 1551 mesh_vertices=mesh_verts, mesh_faces=mesh_faces, 1552 ) 1553 if junction_out is not None: 1554 # A basic-slicing VIEW of the caller's volume, so this writes through with no copy. The 1555 # write is masked rather than a plain assignment because bounding boxes overlap: a plain 1556 # one would zero a neighbour's junctions. Ids are per-instance (1..n) and so repeat across 1557 # mitochondria, but the painted voxels cannot: they lie inside crista_local, hence inside 1558 # mito_local, which is disjoint between instances — which also makes this thread-safe. 1559 written = contact_labels_local > 0 1560 junction_out[written] = contact_labels_local[written] 1561 del membrane_distance, contact_labels_local 1562 else: 1563 contact_summary = { 1564 "contact_voxel_count": 0, "crista_junction_count": 0, "contact_volume_nm3": 0.0, 1565 "mean_junction_extension_nm": np.nan, 1566 } 1567 proximity = {"median_nm": np.nan} 1568 junction_dist = dict(_JUNCTION_DISTANCE_NAN) 1569 1570 if has_crista: 1571 morph = compute_crista_morphology(crista_local, voxel_size) 1572 crista_surface = morph.get("cristae_surface_area_nm2", np.nan) 1573 avg_thickness_nm = morph.get("avg_thickness_nm", np.nan) 1574 if method == "skip": 1575 crista_orientation_anisotropy = np.nan 1576 elif method == "fast": 1577 crista_orientation_anisotropy = _downsampled_orientation_anisotropy( 1578 crista_local, voxel_size, factor=2 1579 ) 1580 else: 1581 anisotropy = compute_crista_orientation(crista_local, voxel_size) 1582 crista_orientation_anisotropy = float(np.mean(anisotropy[crista_local])) 1583 else: 1584 crista_orientation_anisotropy = np.nan 1585 crista_surface = np.nan 1586 avg_thickness_nm = np.nan 1587 1588 if mito_surface and mito_surface > 0 and np.isfinite(crista_surface): 1589 crista_to_mito_surface_ratio = crista_surface / mito_surface 1590 else: 1591 crista_to_mito_surface_ratio = np.nan 1592 1593 # Inner mitochondrial membrane = inner boundary membrane + crista membrane. A crista mask is a 1594 # slab, so its closed mesh already wraps both leaflets — that is the crista membrane area. 1595 # ponytail: cristae contribute 0.0 when absent, not NaN — a crista-less mito still has an IBM. 1596 # ponytail: crista area is meshed fully closed while the lumen is border-trimmed and left open; 1597 # they only disagree when the mito is clipped AND cristae reach the clipped face. 1598 imm_surface = ibm_surface + (crista_surface if np.isfinite(crista_surface) else 0.0) 1599 1600 return { 1601 "mito_label_id": int(label), 1602 "mito_touches_border": touches_border, 1603 "mito_volume_nm3": mito_vol, 1604 "crista_volume_nm3": crista_vol, 1605 "crista_fraction": crista_vol / mito_vol if mito_vol > 0 else np.nan, 1606 "contact_voxel_count": contact_summary["contact_voxel_count"], 1607 "crista_junction_count": contact_summary["crista_junction_count"], 1608 "contact_volume_nm3": contact_summary["contact_volume_nm3"], 1609 "mean_junction_extension_nm": contact_summary["mean_junction_extension_nm"], 1610 "avg_crista_to_membrane_nm": proximity["median_nm"], 1611 "mean_nn_junction_distance_nm": junction_dist["mean_nn_junction_distance_nm"], 1612 "median_nn_junction_distance_nm": junction_dist["median_nn_junction_distance_nm"], 1613 "junction_clustering_index": junction_dist["junction_clustering_index"], 1614 "crista_orientation_anisotropy": crista_orientation_anisotropy, 1615 "cristae_surface_area_nm2": crista_surface, 1616 "mito_surface_area_nm2": mito_surface, 1617 "crista_to_mito_surface_ratio": crista_to_mito_surface_ratio, 1618 "imm_surface_area_nm2": imm_surface, 1619 "imm_surface_per_mito_volume": imm_surface / mito_vol if mito_vol > 0 else np.nan, 1620 "imm_surface_per_crista_volume": imm_surface / crista_vol if crista_vol > 0 else np.nan, 1621 "avg_thickness_nm": avg_thickness_nm, 1622 } 1623 1624 1625def compute_mito_crista_statistics( 1626 crista_mask: np.ndarray, 1627 mito_segmentation: np.ndarray, 1628 voxel_size: Union[float, Dict[str, float]], 1629 membrane_mask: Optional[np.ndarray] = None, 1630 membrane_thickness_nm: float = 8.0, 1631 border_gap_nm: Optional[float] = None, 1632 method: str = "skip", 1633 n_jobs: int = 1, 1634 verbose: bool = False, 1635 progress_callback: Optional[Callable[[int, int], None]] = None, 1636 membrane_mode: str = "slice_2d", 1637 lumen_mask: Optional[np.ndarray] = None, 1638 junction_mode: str = "overlap", 1639 max_extension_nm: Optional[float] = None, 1640 terminus_nm: Optional[float] = None, 1641 min_junction_volume_nm3: Optional[float] = None, 1642 footprint_nm: Optional[float] = None, 1643 min_skeleton_nm: Optional[float] = None, 1644 terminus_merge_nm: Optional[float] = None, 1645 return_junction_labels: bool = False, 1646) -> Union[pd.DataFrame, Tuple[pd.DataFrame, np.ndarray]]: 1647 """Compute all crista metrics organised by mitochondrial instance. 1648 1649 Args: 1650 crista_mask: Binary crista segmentation (global volume). 1651 mito_segmentation: Instance label array (background = 0). 1652 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1653 membrane_mask: Precomputed membrane mask; recomputed if None. 1654 membrane_thickness_nm: Membrane shell thickness used if membrane_mask is None. 1655 border_gap_nm: Border suppression distance passed to approximate_membrane; 1656 defaults to membrane_thickness_nm when None. 1657 method: How the crista orientation anisotropy is computed — this is the ONLY metric that 1658 differs between modes; surface areas (marching cubes), junction distances (geodesic 1659 along the membrane) and thickness/proximity (EDT) are computed identically for all of 1660 them. ``"skip"`` (default) does not compute orientation at all 1661 (``crista_orientation_anisotropy`` is NaN) and is the fastest — use it when only the other 1662 metrics are needed. ``"fast"`` computes the anisotropy on a 2× downsampled crista crop 1663 (~8× cheaper — the structure tensor is by far the dominant cost); the resulting value is a 1664 *relative* indicator that preserves the ordering between mitochondria but is systematically 1665 different in magnitude from the full-resolution value and is NOT comparable to 1666 ``method="exact"``. ``"exact"`` computes the anisotropy from the full-resolution structure 1667 tensor (use it when the magnitude must be precise). 1668 n_jobs: Number of workers for processing mitochondria in parallel (they are 1669 independent). 1 (default) runs serially; other values use a ``concurrent.futures`` 1670 thread pool (-1 = all cores). Results are identical regardless of n_jobs. 1671 verbose: If True, show a terminal tqdm progress bar over mitochondria. 1672 progress_callback: Optional callable invoked once per completed mitochondrion with 1673 (completed_count, total_count) — e.g. to drive a napari progress bar. It is 1674 always called from the calling thread (the futures are consumed here as they 1675 complete), so GUI updates from it need no cross-thread marshaling. 1676 membrane_mode: How the membrane shell is built when ``membrane_mask`` is None — 1677 ``"slice_2d"`` (default, per-Z-slice 2D erosion, z-parallel) or ``"shell_3d"`` (connected 1678 3D shell). See :func:`approximate_membrane`. 1679 lumen_mask: Optional eroded-mito interior matching ``membrane_mask``, i.e. the second return 1680 value of ``approximate_membrane(..., return_lumen=True)``. It is the clean single-wall 1681 surface the junction geodesics run along. Only used when ``membrane_mask`` is also 1682 supplied (when the membrane is built here, the matching lumen is derived automatically); 1683 when neither is available the geodesic mesh falls back to ``mito & ~membrane``, which is 1684 contaminated by the membrane's border-gap suppression near clipped volume faces. 1685 In ``junction_mode="skeleton"`` the lumen is **also** the junction reference surface (see 1686 :func:`_inner_surface_distance`), so supplying it changes ``crista_junction_count`` and 1687 ``mean_junction_extension_nm`` — it is no longer a display-only input. The 1688 ``mito & ~membrane`` fallback is deliberately *not* used for that purpose. 1689 junction_mode: Which junction detector fills ``crista_junction_count`` — ``"overlap"`` 1690 (default) counts connected components of the direct crista-membrane intersection, while 1691 ``"skeleton"`` counts crista regions reaching within ``max_extension_nm`` of the inner 1692 boundary membrane near a crista terminus. See :func:`detect_junctions`. ``"skeleton"`` 1693 additionally fills ``mean_junction_extension_nm`` and requires 3D input. Every other 1694 column is identical between the two modes. 1695 max_extension_nm: How far (nm) a crista may fall short of the inner boundary membrane surface 1696 and still count (``junction_mode="skeleton"`` only). Defaults to 1697 ``membrane_thickness_nm`` when None, following ``border_gap_nm``. 1698 terminus_nm: A near-membrane crista region counts as a junction only if it lies within this 1699 distance (nm) of a crista terminus (``junction_mode="skeleton"`` only), which rejects a 1700 crista running alongside the membrane. Pass ``inf`` to disable the filter. See 1701 :func:`detect_junctions_skeleton`. 1702 min_junction_volume_nm3: Smallest junction volume (nm^3) that counts 1703 (``junction_mode="skeleton"`` only). Removes specks only — see the known limitation in 1704 :func:`detect_junctions_skeleton`. Measured on the candidate region, not on the painted 1705 footprint. 1706 footprint_nm: How thick (nm) the painted junction label is around each region's closest 1707 approach to the membrane (``junction_mode="skeleton"`` only). Affects the junction label 1708 array and hence the junction centroids, never the count. None means one voxel diagonal. 1709 min_skeleton_nm: Skeleton components shorter than this (nm) are dropped as specks before the 1710 terminus filter runs (``junction_mode="skeleton"`` only). A crista whose entire skeleton is 1711 shorter than this has no termini and so cannot score a junction. 1712 terminus_merge_nm: Termini within this distance (nm) on the same skeleton component collapse to 1713 one (``junction_mode="skeleton"`` only). See :func:`compute_crista_skeleton`. 1714 return_junction_labels: If True, also return the junction label array these counts were 1715 computed from, so a caller can *display* exactly what the table reports instead of 1716 re-detecting junctions itself and getting a different answer. Each worker writes its own 1717 instance's labels into a view of the output volume, so this costs one int32 volume and no 1718 extra computation. Ids are per-instance (``1..crista_junction_count``) and therefore repeat 1719 between mitochondria; the painted voxels never do, since each lies inside its own mito. 1720 1721 Every ``junction_mode="skeleton"`` tuning parameter above is ``None`` by default, meaning "leave the 1722 detector's own default in force" — the concrete values live in :func:`detect_junctions_skeleton` and 1723 :func:`compute_crista_skeleton` rather than being restated here. 1724 1725 The junction nearest-neighbour distances are geodesics along the eroded-mito surface mesh 1726 (``bioimage_cpp.distance.geodesic_distances_mesh``); for a mito with no usable mesh (empty 1727 membrane / degenerate mesh) those columns are NaN. 1728 1729 Implementation notes: each mito is pre-cropped to its bounding box by basic slicing (views, so 1730 cropping is memory-free). Parallelism is adaptive and single-level (never oversubscribed): with 1731 many mitochondria the work is parallelised *across* them on a ``concurrent.futures`` 1732 ``ThreadPoolExecutor`` — the heavy per-mito stages (structure tensor, EDT, geodesics) are 1733 GIL-releasing C++, so threads scale them — with each worker's inner stages kept single-threaded 1734 (the EDT/geodesic solvers are called with ``number_of_threads=1``); with few mitochondria they run 1735 serially and each mito's junction-distance stage gets all cores. The concurrent worker count is 1736 additionally capped so the combined per-mito working set (tensor components + label crops, 1737 ~40 bytes/voxel of the largest mito) fits in RAM. Results stream in as they complete and are 1738 finally sorted by label for an n_jobs-independent ordering. 1739 1740 Returns: 1741 DataFrame with one row per mito instance: 1742 label | mito_volume_nm3 | crista_volume_nm3 | crista_fraction | 1743 contact_voxel_count | crista_junction_count | contact_volume_nm3 | 1744 mean_junction_extension_nm | 1745 avg_crista_to_membrane_nm | mean_nn_junction_distance_nm | median_nn_junction_distance_nm | 1746 junction_clustering_index | crista_orientation_anisotropy | cristae_surface_area_nm2 | 1747 mito_surface_area_nm2 | crista_to_mito_surface_ratio | imm_surface_area_nm2 | 1748 imm_surface_per_mito_volume | imm_surface_per_crista_volume | avg_thickness_nm. 1749 cristae_surface_area_nm2 is the crista surface area; crista_to_mito_surface_ratio is 1750 crista surface / mitochondrial outer-membrane surface (can exceed 1 for folded cristae). 1751 imm_surface_area_nm2 is the inner mitochondrial membrane area — the inner boundary membrane 1752 (the lumen surface) plus cristae_surface_area_nm2 — and the imm_surface_per_*_volume columns 1753 divide it by mito_volume_nm3 and crista_volume_nm3 respectively (nm^-1, "cristae surface 1754 density"). The inner-boundary-membrane area alone is 1755 imm_surface_area_nm2 - cristae_surface_area_nm2. 1756 When ``return_junction_labels`` is True the return value is instead 1757 ``(DataFrame, junction_labels)``, the second an int32 array of the input shape. 1758 1759 The *_nn_junction_distance_nm columns are geodesic nearest-neighbour distances between 1760 crista-membrane junctions along the membrane; junction_clustering_index is a Clark-Evans 1761 index (< 1 clustered, ~ 1 random, > 1 dispersed). ``mean_junction_extension_nm`` is the mean 1762 gap each skeleton end had to bridge to reach the membrane and is NaN unless 1763 ``junction_mode="skeleton"``. ``crista_orientation_anisotropy`` is 1764 computed at full resolution for ``method="exact"``, on a downsampled crop (relative-only, 1765 not comparable) for ``method="fast"``, and left NaN for ``method="skip"``. 1766 """ 1767 if method not in ("fast", "exact", "skip"): 1768 raise ValueError(f"method must be 'fast', 'exact', or 'skip', got {method!r}") 1769 if junction_mode not in ("overlap", "skeleton"): 1770 raise ValueError(f"junction_mode must be 'overlap' or 'skeleton', got {junction_mode!r}") 1771 if max_extension_nm is None: 1772 max_extension_nm = membrane_thickness_nm 1773 if membrane_mask is None: 1774 membrane_mask, lumen_mask = approximate_membrane( 1775 mito_segmentation, voxel_size, membrane_thickness_nm, border_gap_nm, 1776 n_jobs=n_jobs, membrane_mode=membrane_mode, return_lumen=True, 1777 ) 1778 1779 ndim = mito_segmentation.ndim 1780 sampling = _to_sampling(voxel_size, ndim) 1781 voxel_vol = float(np.prod(sampling)) 1782 crista_binary = crista_mask.astype(bool) 1783 vol_shape = mito_segmentation.shape 1784 border_radius = _gap_radius(voxel_size, membrane_thickness_nm, border_gap_nm, ndim) 1785 1786 junction_labels = np.zeros(vol_shape, dtype=np.int32) if return_junction_labels else None 1787 1788 tasks = [] 1789 for prop in regionprops(mito_segmentation): 1790 bbox = prop.bbox 1791 slices = tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim)) 1792 lumen_crop = None if lumen_mask is None else lumen_mask[slices] 1793 tasks.append(( 1794 int(prop.label), bbox, 1795 mito_segmentation[slices], crista_binary[slices], membrane_mask[slices], 1796 lumen_crop, 1797 # a basic-slicing view, so each worker writes its junctions straight into the output 1798 None if junction_labels is None else junction_labels[slices], 1799 )) 1800 total = len(tasks) 1801 1802 def _run(task, inner_n_jobs): 1803 label, bbox, mito_crop, crista_crop, membrane_crop, lumen_crop, junction_out = task 1804 return _single_mito_row( 1805 label, bbox, mito_crop, crista_crop, membrane_crop, 1806 voxel_size, sampling, voxel_vol, vol_shape, border_radius, 1807 method=method, inner_n_jobs=inner_n_jobs, lumen_crop=lumen_crop, 1808 junction_mode=junction_mode, max_extension_nm=max_extension_nm, 1809 terminus_nm=terminus_nm, min_junction_volume_nm3=min_junction_volume_nm3, 1810 footprint_nm=footprint_nm, 1811 min_skeleton_nm=min_skeleton_nm, terminus_merge_nm=terminus_merge_nm, 1812 junction_out=junction_out, 1813 ) 1814 1815 n_workers = os.cpu_count() if n_jobs == -1 else max(1, n_jobs) 1816 across = n_workers > 1 and total >= n_workers 1817 1818 rows = [] 1819 1820 def _consume(results): 1821 for i, row in enumerate( 1822 tqdm(results, total=total, desc="Cristae analysis", disable=not verbose), start=1 1823 ): 1824 rows.append(row) 1825 if progress_callback is not None: 1826 progress_callback(i, total) 1827 1828 if across: 1829 max_voxels = max(int(task[2].size) for task in tasks) 1830 across_workers = _bounded_workers(n_jobs, per_worker_bytes=max_voxels * 40) 1831 # Parallelise across mitochondria with a thread pool (the heavy per-mito stages — structure 1832 # tensor, EDT, geodesics — are GIL-releasing C++). Each worker's inner stages run 1833 # single-threaded (``_run(task, 1)`` passes ``number_of_threads=1`` down to the EDT/geodesic 1834 # solvers) so the across-mito threads do not oversubscribe the cores. Results stream in as they 1835 # complete (``as_completed``) to drive the progress bar; rows are label-sorted below. 1836 with futures.ThreadPoolExecutor(across_workers) as tp: 1837 submitted = [tp.submit(_run, task, 1) for task in tasks] 1838 _consume(future.result() for future in futures.as_completed(submitted)) 1839 else: 1840 _consume(_run(task, n_workers) for task in tasks) 1841 1842 rows.sort(key=lambda row: row["mito_label_id"]) 1843 table = pd.DataFrame(rows) 1844 return (table, junction_labels) if return_junction_labels else table
def
approximate_membrane( mito_segmentation: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], membrane_thickness_nm: float = 8.0, border_gap_nm: Optional[float] = None, n_jobs: int = 1, membrane_mode: str = 'slice_2d', return_lumen: bool = False) -> Union[numpy.ndarray, Tuple[numpy.ndarray, numpy.ndarray]]:
271def approximate_membrane( 272 mito_segmentation: np.ndarray, 273 voxel_size: Union[float, Dict[str, float]], 274 membrane_thickness_nm: float = 8.0, 275 border_gap_nm: Optional[float] = None, 276 n_jobs: int = 1, 277 membrane_mode: str = "slice_2d", 278 return_lumen: bool = False, 279) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]: 280 """Approximate the mitochondrial membrane as the outer shell of the segmentation. 281 282 Two shell constructions are available via ``membrane_mode``: 283 284 - ``"slice_2d"`` (default): erode each Z-slice **independently in 2D** by an XY disk of radius 285 ``round(thickness / xy_voxel)`` and keep ``slice & ~eroded``. A mitochondrion that changes shape 286 rapidly along Z does not bleed into neighbouring slices, and the per-slice erosions are 287 parallelised over Z (``n_jobs``). A separable Z-only erosion (radius ``round(thickness / 288 z_voxel)``, no XY coupling) then adds **Z-caps** where a mito column truly ends in Z; ends 289 clipped by a volume Z-face are left uncapped (``border_value=1`` + the ``border_gap`` trim). The 290 XY shell can still fragment across slices. (2D inputs get a single 2D erosion.) 291 - ``"shell_3d"``: a full 3D morphological erosion, ``mito & ~erode3d(mito, k)`` with 292 ``k = round(thickness / mean_voxel)`` iterations of a 3×3×3 structuring element, per instance on 293 its padded bounding box. A single **connected** shell including the Z-caps (no per-slice 294 fragmentation), at a higher cost; thickness acts in all axes. 295 296 The eroded interior is the "lumen"; its surface is the single-wall mesh used by the geodesic 297 backend and the display, so the mesh follows the chosen mode. 298 299 Membrane voxels within ``border_gap_nm`` of any volume face are removed so clipped mito edges are 300 not treated as membrane. The lumen is NOT trimmed here — the mesh is trimmed to the certain region 301 (and left open there) at mesh time by :func:`_open_trimmed_mesh`, which requires the untrimmed 302 interior to produce an open cut rather than a fabricated cap. 303 304 Implementation notes: ``"slice_2d"`` erodes each Z-slice on the mito XY bbox with a 305 ``membrane_radius`` margin, so the cropped ``border_value=1`` erosion matches eroding the full 306 slice (empty slices are skipped), then adds Z-caps via a separable Z-only line erosion (which 307 inspects only the same column, so no XY-shape bleed); ``border_value=1`` leaves ends clipped by a 308 volume Z-face uncapped, and the ``border_gap`` removal clears anything near a face, so only true 309 ends are capped. ``"shell_3d"`` erodes the *merged* binary in each instance's padded bbox so 310 instances that share a boundary are handled together. 311 312 Args: 313 mito_segmentation: Instance label array (background = 0). 314 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 315 membrane_thickness_nm: Thickness of the membrane shell in nm. 316 border_gap_nm: Distance from each volume face within which membrane voxels are 317 suppressed. Defaults to membrane_thickness_nm when None. 318 n_jobs: Workers for the per-Z-slice erosion (``"slice_2d"`` only): 1 = serial, -1 = all cores. 319 Results are identical regardless of n_jobs. 320 membrane_mode: ``"slice_2d"`` (default, per-slice 2D, z-parallel) or ``"shell_3d"`` 321 (connected 3D shell). 322 return_lumen: If True, also return the eroded-mito interior ("lumen") mask — the single-wall 323 surface source for the geodesic/display mesh. It is the plain eroded interior (not 324 ``mito & ~membrane``, which would re-include the outer shell where the membrane is 325 border-trimmed); trimming to the certain region happens at mesh time. 326 327 Returns: 328 membrane_mask: Binary mask of the mitochondrial membrane (outer shell), with border-adjacent 329 voxels zeroed out. If ``return_lumen`` is True, returns ``(membrane_mask, lumen_mask)`` 330 where ``lumen_mask`` is the (untrimmed) eroded interior described above. 331 """ 332 if membrane_mode not in ("slice_2d", "shell_3d"): 333 raise ValueError(f"membrane_mode must be 'slice_2d' or 'shell_3d', got {membrane_mode!r}") 334 ndim = mito_segmentation.ndim 335 mito_binary = mito_segmentation > 0 336 337 # NOTE (possible simplification, deferred): the "shell_3d" branch below builds the shell with an 338 # iterated 3x3x3 erosion per instance. It could likely be a single anisotropic distance transform 339 # instead — membrane = mito & (distance_transform(mito, sampling) <= thickness), lumen = the rest — 340 # which is simpler and handles anisotropy directly. This would NOT replace "slice_2d": a distance 341 # transform couples all axes, so it cannot reproduce slice_2d's per-Z-slice-independent erosion, 342 # whose whole purpose is to stop the shell bleeding across slices in XY. Worth investigating. 343 if membrane_mode == "shell_3d": 344 sampling = _to_sampling(voxel_size, ndim) 345 k = max(1, int(round(float(membrane_thickness_nm) / float(np.mean(sampling))))) 346 struct = np.ones((3,) * ndim, dtype=bool) 347 membrane_mask = np.zeros(mito_segmentation.shape, dtype=bool) 348 lumen_mask = np.zeros(mito_segmentation.shape, dtype=bool) 349 for prop in regionprops(mito_segmentation): 350 bbox = prop.bbox 351 sl = tuple( 352 slice(max(0, bbox[i] - k), min(mito_segmentation.shape[i], bbox[i + ndim] + k)) 353 for i in range(ndim) 354 ) 355 sub = mito_binary[sl] 356 eroded = binary_erosion(sub, structure=struct, iterations=k, border_value=1) 357 cur = mito_segmentation[sl] == prop.label 358 membrane_mask[sl] |= cur & ~eroded 359 lumen_mask[sl] |= cur & eroded 360 elif ndim == 3: 361 membrane_radius = _voxel_radius_xy(membrane_thickness_nm, voxel_size) 362 struct = disk(membrane_radius) 363 membrane_mask = np.zeros_like(mito_binary) 364 lumen_mask = np.zeros_like(mito_binary) 365 coords = np.argwhere(mito_binary) 366 if coords.size: 367 zmin, ymin, xmin = coords.min(axis=0) 368 zmax, ymax, xmax = coords.max(axis=0) + 1 369 m = membrane_radius 370 y0, y1 = max(0, ymin - m), min(mito_binary.shape[1], ymax + m) 371 x0, x1 = max(0, xmin - m), min(mito_binary.shape[2], xmax + m) 372 373 def _erode_slice(z): 374 sl = mito_binary[z, y0:y1, x0:x1] 375 if not sl.any(): 376 return z, None 377 eroded = binary_erosion(sl, structure=struct, border_value=1) 378 return z, (sl & ~eroded, eroded) 379 380 z_range = range(int(zmin), int(zmax)) 381 if n_jobs == 1: 382 results = [_erode_slice(z) for z in z_range] 383 else: 384 n_workers = mp.cpu_count() if n_jobs == -1 else n_jobs 385 with futures.ThreadPoolExecutor(n_workers) as tp: 386 results = list(tp.map(_erode_slice, z_range)) 387 for z, res in results: 388 if res is not None: 389 mem_sl, lum_sl = res 390 membrane_mask[z, y0:y1, x0:x1] = mem_sl 391 lumen_mask[z, y0:y1, x0:x1] = lum_sl 392 393 k_z = max(1, int(round(float(membrane_thickness_nm) / float(_to_sampling(voxel_size, ndim)[0])))) 394 z0m, z1m = max(0, int(zmin) - k_z), min(mito_binary.shape[0], int(zmax) + k_z) 395 sub = mito_binary[z0m:z1m, y0:y1, x0:x1] 396 z_eroded = binary_erosion(sub, structure=np.ones((2 * k_z + 1, 1, 1), dtype=bool), border_value=1) 397 membrane_mask[z0m:z1m, y0:y1, x0:x1] |= sub & ~z_eroded 398 lumen_mask[z0m:z1m, y0:y1, x0:x1] &= z_eroded 399 else: 400 membrane_radius = _voxel_radius(membrane_thickness_nm, voxel_size, ndim) 401 eroded = binary_erosion(mito_binary, structure=disk(membrane_radius), border_value=1) 402 membrane_mask = mito_binary & ~eroded 403 lumen_mask = mito_binary & eroded 404 405 gap_radius = _gap_radius(voxel_size, membrane_thickness_nm, border_gap_nm, ndim) 406 membrane_mask &= ~_border_zone(mito_segmentation.shape, gap_radius) 407 if return_lumen: 408 return membrane_mask.astype(bool), lumen_mask.astype(bool) 409 return membrane_mask.astype(bool)
Approximate the mitochondrial membrane as the outer shell of the segmentation.
Two shell constructions are available via membrane_mode:
"slice_2d"(default): erode each Z-slice independently in 2D by an XY disk of radiusround(thickness / xy_voxel)and keepslice & ~eroded. A mitochondrion that changes shape rapidly along Z does not bleed into neighbouring slices, and the per-slice erosions are parallelised over Z (n_jobs). A separable Z-only erosion (radiusround(thickness / z_voxel), no XY coupling) then adds Z-caps where a mito column truly ends in Z; ends clipped by a volume Z-face are left uncapped (border_value=1+ theborder_gaptrim). The XY shell can still fragment across slices. (2D inputs get a single 2D erosion.)"shell_3d": a full 3D morphological erosion,mito & ~erode3d(mito, k)withk = round(thickness / mean_voxel)iterations of a 3×3×3 structuring element, per instance on its padded bounding box. A single connected shell including the Z-caps (no per-slice fragmentation), at a higher cost; thickness acts in all axes.
The eroded interior is the "lumen"; its surface is the single-wall mesh used by the geodesic backend and the display, so the mesh follows the chosen mode.
Membrane voxels within border_gap_nm of any volume face are removed so clipped mito edges are
not treated as membrane. The lumen is NOT trimmed here — the mesh is trimmed to the certain region
(and left open there) at mesh time by _open_trimmed_mesh(), which requires the untrimmed
interior to produce an open cut rather than a fabricated cap.
Implementation notes: "slice_2d" erodes each Z-slice on the mito XY bbox with a
membrane_radius margin, so the cropped border_value=1 erosion matches eroding the full
slice (empty slices are skipped), then adds Z-caps via a separable Z-only line erosion (which
inspects only the same column, so no XY-shape bleed); border_value=1 leaves ends clipped by a
volume Z-face uncapped, and the border_gap removal clears anything near a face, so only true
ends are capped. "shell_3d" erodes the merged binary in each instance's padded bbox so
instances that share a boundary are handled together.
Arguments:
- mito_segmentation: Instance label array (background = 0).
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- membrane_thickness_nm: Thickness of the membrane shell in nm.
- border_gap_nm: Distance from each volume face within which membrane voxels are suppressed. Defaults to membrane_thickness_nm when None.
- n_jobs: Workers for the per-Z-slice erosion (
"slice_2d"only): 1 = serial, -1 = all cores. Results are identical regardless of n_jobs. - membrane_mode:
"slice_2d"(default, per-slice 2D, z-parallel) or"shell_3d"(connected 3D shell). - return_lumen: If True, also return the eroded-mito interior ("lumen") mask — the single-wall
surface source for the geodesic/display mesh. It is the plain eroded interior (not
mito & ~membrane, which would re-include the outer shell where the membrane is border-trimmed); trimming to the certain region happens at mesh time.
Returns:
membrane_mask: Binary mask of the mitochondrial membrane (outer shell), with border-adjacent voxels zeroed out. If
return_lumenis True, returns(membrane_mask, lumen_mask)wherelumen_maskis the (untrimmed) eroded interior described above.
def
compute_crista_orientation( crista_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], neighborhood_size_nm: float = 30.0) -> numpy.ndarray:
416def compute_crista_orientation( 417 crista_mask: np.ndarray, 418 voxel_size: Union[float, Dict[str, float]], 419 neighborhood_size_nm: float = 30.0, 420) -> np.ndarray: 421 """Compute the per-voxel crista orientation anisotropy via the structure tensor. 422 423 Uses ``bioimage_cpp.filters.structure_tensor_eigenvalues`` (a fast C++ routine). Only the 424 anisotropy is produced (the principal directions / eigenvectors are not computed). 425 426 The structure tensor's outer/integration sigma is ``neighborhood_size_nm`` per axis (in voxels); 427 the inner (derivative) sigma must be > 0, so a minimal 1-voxel scale is used. Eigenvalues are 428 non-negative in theory, but the solver emits tiny negatives for near-rank-deficient tensors 429 (degenerate sheets/tubes), so they are clamped to 0 and the ratio is taken as 430 ``max/min`` over the trailing axis — order-agnostic and sign-safe, so a tiny negative minor 431 eigenvalue cannot flip the denominator and blow the ratio up. 432 433 Args: 434 crista_mask: Binary crista segmentation. 435 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 436 neighborhood_size_nm: Gaussian integration radius in nm for tensor averaging (the 437 structure tensor's outer/integration scale). 438 439 Returns: 440 anisotropy: (...) — λ_max / (λ_min + ε) per voxel. High values indicate a strongly 441 directional crista (e.g. parallel lamellae); low values indicate isotropic or 442 tubular/disordered morphology. Magnitude only (rotation-invariant). 443 """ 444 ndim = crista_mask.ndim 445 sampling = _to_sampling(voxel_size, ndim) 446 447 outer_sigma = [float(s) for s in (neighborhood_size_nm / sampling)] 448 inner_sigma = 1.0 449 evals = structure_tensor_eigenvalues(crista_mask.astype(np.float32), inner_sigma, outer_sigma) 450 evals = np.clip(evals, 0.0, None) 451 anisotropy = evals.max(axis=-1) / (evals.min(axis=-1) + 1e-10) 452 return anisotropy.astype(np.float32)
Compute the per-voxel crista orientation anisotropy via the structure tensor.
Uses bioimage_cpp.filters.structure_tensor_eigenvalues (a fast C++ routine). Only the
anisotropy is produced (the principal directions / eigenvectors are not computed).
The structure tensor's outer/integration sigma is neighborhood_size_nm per axis (in voxels);
the inner (derivative) sigma must be > 0, so a minimal 1-voxel scale is used. Eigenvalues are
non-negative in theory, but the solver emits tiny negatives for near-rank-deficient tensors
(degenerate sheets/tubes), so they are clamped to 0 and the ratio is taken as
max/min over the trailing axis — order-agnostic and sign-safe, so a tiny negative minor
eigenvalue cannot flip the denominator and blow the ratio up.
Arguments:
- crista_mask: Binary crista segmentation.
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- neighborhood_size_nm: Gaussian integration radius in nm for tensor averaging (the structure tensor's outer/integration scale).
Returns:
anisotropy: (...) — λ_max / (λ_min + ε) per voxel. High values indicate a strongly directional crista (e.g. parallel lamellae); low values indicate isotropic or tubular/disordered morphology. Magnitude only (rotation-invariant).
def
compute_crista_proximity( crista_mask: numpy.ndarray, membrane_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], membrane_distance: Optional[numpy.ndarray] = None) -> Tuple[numpy.ndarray, Dict[str, float]]:
503def compute_crista_proximity( 504 crista_mask: np.ndarray, 505 membrane_mask: np.ndarray, 506 voxel_size: Union[float, Dict[str, float]], 507 membrane_distance: Optional[np.ndarray] = None, 508) -> Tuple[np.ndarray, Dict[str, float]]: 509 """Distance from each crista voxel to the nearest membrane voxel (nm). 510 511 Args: 512 crista_mask: Binary crista segmentation. 513 membrane_mask: Binary membrane mask (OM or IMM). 514 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 515 membrane_distance: Optional precomputed per-voxel distance to the nearest membrane 516 voxel (nm), i.e. ``distance_transform(~membrane_mask, sampling=...)``. When 517 given, the distance transform is not recomputed (used to avoid redundant work in 518 :func:`compute_mito_crista_statistics`). 519 520 Returns: 521 distance_map: Per-voxel distance to membrane (nm); zero outside crista. 522 summary_stats: min_nm, median_nm, max_nm. 523 """ 524 sampling = _to_sampling(voxel_size, crista_mask.ndim) 525 if membrane_distance is None: 526 dist = distance_transform(~membrane_mask.astype(bool), sampling=sampling.tolist(), number_of_threads=1) 527 else: 528 dist = membrane_distance 529 crista_dists = dist[crista_mask.astype(bool)] 530 531 if crista_dists.size == 0: 532 summary: Dict[str, float] = {"min_nm": np.nan, "median_nm": np.nan, "max_nm": np.nan} 533 else: 534 summary = { 535 "min_nm": float(crista_dists.min()), 536 "median_nm": float(np.median(crista_dists)), 537 "max_nm": float(crista_dists.max()), 538 } 539 540 distance_map = np.zeros(crista_mask.shape, dtype=np.float32) 541 distance_map[crista_mask.astype(bool)] = crista_dists 542 return distance_map, summary
Distance from each crista voxel to the nearest membrane voxel (nm).
Arguments:
- crista_mask: Binary crista segmentation.
- membrane_mask: Binary membrane mask (OM or IMM).
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- membrane_distance: Optional precomputed per-voxel distance to the nearest membrane
voxel (nm), i.e.
distance_transform(~membrane_mask, sampling=...). When given, the distance transform is not recomputed (used to avoid redundant work incompute_mito_crista_statistics()).
Returns:
distance_map: Per-voxel distance to membrane (nm); zero outside crista. summary_stats: min_nm, median_nm, max_nm.
def
detect_contact_sites( crista_mask: numpy.ndarray, membrane_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]]) -> Tuple[numpy.ndarray, Dict[str, float]]:
549def detect_contact_sites( 550 crista_mask: np.ndarray, 551 membrane_mask: np.ndarray, 552 voxel_size: Union[float, Dict[str, float]], 553) -> Tuple[np.ndarray, Dict[str, float]]: 554 """Detect crista-membrane contact sites as the direct overlap of the two masks. 555 556 Contact = crista voxels that are also membrane voxels (the pure intersection of the crista 557 mask and the mitochondrial membrane band). No dilation/erosion is applied here, so the 558 detected junctions correspond exactly to the visible overlap of the two layers; connected 559 overlaps are grouped into junctions with 26-connectivity in 3D. 560 561 Args: 562 crista_mask: Binary crista segmentation. 563 membrane_mask: Binary mitochondrial membrane mask (as displayed). 564 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 565 566 Returns: 567 contact_labels: Integer array (same shape as the input) where each connected 568 junction has a unique ID (0 = background, 1..n = junctions). Contact voxel 569 coordinates are recoverable via ``np.argwhere(contact_labels > 0)``. 570 summary: contact_voxel_count, crista_junction_count, contact_volume_nm3. 571 """ 572 ndim = crista_mask.ndim 573 sampling = _to_sampling(voxel_size, ndim) 574 voxel_vol = float(np.prod(sampling)) 575 576 contact_mask = crista_mask.astype(bool) & membrane_mask.astype(bool) 577 578 connectivity_struct = np.ones(ndim * (3,), dtype=bool) 579 contact_labels, n_regions = ndimage_label(contact_mask, structure=connectivity_struct) 580 contact_voxel_count = int(np.count_nonzero(contact_labels)) 581 582 return contact_labels, { 583 "contact_voxel_count": contact_voxel_count, 584 "crista_junction_count": int(n_regions), 585 "contact_volume_nm3": float(contact_voxel_count) * voxel_vol, 586 }
Detect crista-membrane contact sites as the direct overlap of the two masks.
Contact = crista voxels that are also membrane voxels (the pure intersection of the crista mask and the mitochondrial membrane band). No dilation/erosion is applied here, so the detected junctions correspond exactly to the visible overlap of the two layers; connected overlaps are grouped into junctions with 26-connectivity in 3D.
Arguments:
- crista_mask: Binary crista segmentation.
- membrane_mask: Binary mitochondrial membrane mask (as displayed).
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
Returns:
contact_labels: Integer array (same shape as the input) where each connected junction has a unique ID (0 = background, 1..n = junctions). Contact voxel coordinates are recoverable via
np.argwhere(contact_labels > 0). summary: contact_voxel_count, crista_junction_count, contact_volume_nm3.
def
compute_crista_skeleton( crista_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], n_jobs: int = 1, min_skeleton_nm: float = 10.0, terminus_merge_nm: float = 4.0, return_edges: bool = False):
694def compute_crista_skeleton( 695 crista_mask: np.ndarray, 696 voxel_size: Union[float, Dict[str, float]], 697 n_jobs: int = 1, 698 min_skeleton_nm: float = 10.0, 699 terminus_merge_nm: float = 4.0, 700 return_edges: bool = False, 701): 702 """The crista centerline skeleton and its termini, for inspection and display. 703 704 Wraps ``bioimage_cpp.skeleton.teasar`` and then **cleans the resulting graph up**, because the raw 705 TEASAR output is not usable as a set of crista ends: it spans a lamella with a caterpillar of short 706 side branches whose tips line the sheet rim, every segmentation speck adds its own miniature 707 skeleton, and counting ``degree <= 1`` also counts isolated vertices. On a tomogram-scale 708 40-lamella mask that produces 10080 termini where roughly 80 are real. 709 710 The cleanup is :func:`_skeleton_graph` (drop speck components) followed by :func:`_merge_termini` 711 (collapse each fan of nearby termini on one component), which brings the same mask to 480, and 712 ``degree == 1`` instead of ``<= 1``. 713 714 **Why no leaf-branch pruning.** Cutting short leaf branches off a surviving component is the obvious 715 first idea, and it was measured to be actively harmful. TEASAR's medial axis of a lamella is a 716 caterpillar — a main path with many short side leaves — and at the sheet's real end the main path 717 itself arrives as a short leaf off a nearby branch node. A length threshold therefore deletes the 718 real crista end: on the test lamella spanning y 13-46 it left termini only at y 45-46, losing the 719 y=13 junction entirely. It also fragmented components, so fewer termini could be merged afterwards 720 (10080 -> 960 termini with pruning, against 480 without it on the same tomogram-scale mask). 721 Component filtering plus the per-component merge is both simpler and strictly better. 722 723 :func:`detect_junctions_skeleton` keys its terminus filter on exactly these points, so plotting 724 them is the way to see why a junction was accepted or rejected — that is what the napari widget's 725 **Show Crista Skeleton** option displays. 726 727 Args: 728 crista_mask: Binary crista segmentation (3D). 729 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 730 n_jobs: Forwarded to TEASAR's ``number_of_threads`` (-1/0 = all cores). 731 min_skeleton_nm: Skeleton components whose total length is below this (nm) are dropped as 732 specks, along with their termini — so a crista smaller than this cannot contribute a 733 terminus at all. 0 keeps every component TEASAR produced. 734 terminus_merge_nm: Termini within this distance (nm) of each other **on the same skeleton 735 component** collapse to one; 0 disables merging. This is what actually tames the terminus 736 count (4640 -> 960 on the mask above), since component filtering alone leaves the rim fan 737 intact. The value is the scale of the ragged fan at one crista end, a few voxels — **not** 738 the crista size. It is deliberately not larger even though larger radii look tidier: the 739 junction detector accepts a region only if it lies within ``terminus_nm`` of *some* 740 terminus, so collapsing a whole sheet rim to one point costs real junctions at the other 741 end of that rim. A 16 nm radius reaches the cosmetic ideal of two termini per lamella and 742 reduces a small sheet to a single terminus, which is exactly that failure. 743 return_edges: If True, also return the surviving edges, so a caller can draw the skeleton as 744 connected segments rather than a point cloud (the widget's Vectors layer). 745 746 Returns: 747 (vertices, is_terminus), or (vertices, is_terminus, edges) when ``return_edges`` is True. 748 ``vertices`` is (n, 3) in **nm**, array (z, y, x) order (divide by the voxel size for array 749 indices) and covers only the vertices of surviving components; ``is_terminus`` is the matching 750 boolean mask; ``edges`` is (m, 2) of indices into ``vertices``. All empty if the skeleton is 751 empty. 752 753 Raises: 754 ValueError: If the input is not 3D — TEASAR has no 2D implementation. 755 """ 756 if crista_mask.ndim != 3: 757 raise ValueError( 758 "compute_crista_skeleton requires a 3D volume (bioimage_cpp.skeleton.teasar has no 2D " 759 f"implementation), got ndim={crista_mask.ndim}." 760 ) 761 sampling = _to_sampling(voxel_size, 3) 762 vertices, edges, _ = teasar( 763 crista_mask.astype(bool), spacing=tuple(float(s) for s in sampling), 764 number_of_threads=_solver_threads(n_jobs), 765 ) 766 empty = (np.zeros((0, 3), dtype=float), np.zeros(0, dtype=bool)) 767 if len(vertices) == 0: 768 return (*empty, np.zeros((0, 2), dtype=int)) if return_edges else empty 769 770 vertices = np.asarray(vertices, dtype=float) 771 graph = _skeleton_graph(vertices, edges, min_skeleton_nm) 772 if graph.number_of_nodes() == 0: 773 return (*empty, np.zeros((0, 2), dtype=int)) if return_edges else empty 774 775 termini = set(_merge_termini(graph, vertices, [n for n in graph.nodes if graph.degree(n) == 1], 776 terminus_merge_nm)) 777 778 kept = sorted(graph.nodes) 779 remap = {node: index for index, node in enumerate(kept)} 780 out_vertices = vertices[kept] 781 is_terminus = np.array([node in termini for node in kept], dtype=bool) 782 if not return_edges: 783 return out_vertices, is_terminus 784 out_edges = np.array([[remap[a], remap[b]] for a, b in graph.edges], dtype=int).reshape(-1, 2) 785 return out_vertices, is_terminus, out_edges
The crista centerline skeleton and its termini, for inspection and display.
Wraps bioimage_cpp.skeleton.teasar and then cleans the resulting graph up, because the raw
TEASAR output is not usable as a set of crista ends: it spans a lamella with a caterpillar of short
side branches whose tips line the sheet rim, every segmentation speck adds its own miniature
skeleton, and counting degree <= 1 also counts isolated vertices. On a tomogram-scale
40-lamella mask that produces 10080 termini where roughly 80 are real.
The cleanup is _skeleton_graph() (drop speck components) followed by _merge_termini()
(collapse each fan of nearby termini on one component), which brings the same mask to 480, and
degree == 1 instead of <= 1.
Why no leaf-branch pruning. Cutting short leaf branches off a surviving component is the obvious first idea, and it was measured to be actively harmful. TEASAR's medial axis of a lamella is a caterpillar — a main path with many short side leaves — and at the sheet's real end the main path itself arrives as a short leaf off a nearby branch node. A length threshold therefore deletes the real crista end: on the test lamella spanning y 13-46 it left termini only at y 45-46, losing the y=13 junction entirely. It also fragmented components, so fewer termini could be merged afterwards (10080 -> 960 termini with pruning, against 480 without it on the same tomogram-scale mask). Component filtering plus the per-component merge is both simpler and strictly better.
detect_junctions_skeleton() keys its terminus filter on exactly these points, so plotting
them is the way to see why a junction was accepted or rejected — that is what the napari widget's
Show Crista Skeleton option displays.
Arguments:
- crista_mask: Binary crista segmentation (3D).
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- n_jobs: Forwarded to TEASAR's
number_of_threads(-1/0 = all cores). - min_skeleton_nm: Skeleton components whose total length is below this (nm) are dropped as specks, along with their termini — so a crista smaller than this cannot contribute a terminus at all. 0 keeps every component TEASAR produced.
- terminus_merge_nm: Termini within this distance (nm) of each other on the same skeleton
component collapse to one; 0 disables merging. This is what actually tames the terminus
count (4640 -> 960 on the mask above), since component filtering alone leaves the rim fan
intact. The value is the scale of the ragged fan at one crista end, a few voxels — not
the crista size. It is deliberately not larger even though larger radii look tidier: the
junction detector accepts a region only if it lies within
terminus_nmof some terminus, so collapsing a whole sheet rim to one point costs real junctions at the other end of that rim. A 16 nm radius reaches the cosmetic ideal of two termini per lamella and reduces a small sheet to a single terminus, which is exactly that failure. - return_edges: If True, also return the surviving edges, so a caller can draw the skeleton as connected segments rather than a point cloud (the widget's Vectors layer).
Returns:
(vertices, is_terminus), or (vertices, is_terminus, edges) when
return_edgesis True.verticesis (n, 3) in nm, array (z, y, x) order (divide by the voxel size for array indices) and covers only the vertices of surviving components;is_terminusis the matching boolean mask;edgesis (m, 2) of indices intovertices. All empty if the skeleton is empty.
Raises:
- ValueError: If the input is not 3D — TEASAR has no 2D implementation.
def
detect_junctions_skeleton( crista_mask: numpy.ndarray, membrane_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], max_extension_nm: float = 8.0, min_extension_nm: float = 0.0, terminus_nm: float = 20.0, min_junction_volume_nm3: float = 50.0, border_radius: int = 0, boundary: Optional[numpy.ndarray] = None, n_jobs: int = 1, membrane_distance: Optional[numpy.ndarray] = None, lumen_mask: Optional[numpy.ndarray] = None, footprint_nm: Optional[float] = None, min_skeleton_nm: Optional[float] = None, terminus_merge_nm: Optional[float] = None) -> Tuple[numpy.ndarray, Dict[str, float]]:
817def detect_junctions_skeleton( 818 crista_mask: np.ndarray, 819 membrane_mask: np.ndarray, 820 voxel_size: Union[float, Dict[str, float]], 821 max_extension_nm: float = 8.0, 822 min_extension_nm: float = 0.0, 823 terminus_nm: float = 20.0, 824 min_junction_volume_nm3: float = 50.0, 825 border_radius: int = 0, 826 boundary: Optional[np.ndarray] = None, 827 n_jobs: int = 1, 828 membrane_distance: Optional[np.ndarray] = None, 829 lumen_mask: Optional[np.ndarray] = None, 830 footprint_nm: Optional[float] = None, 831 min_skeleton_nm: Optional[float] = None, 832 terminus_merge_nm: Optional[float] = None, 833) -> Tuple[np.ndarray, Dict[str, float]]: 834 """Detect crista-membrane junctions as crista regions that reach close to the membrane. 835 836 The gap-tolerant alternative to :func:`detect_contact_sites`. Where the overlap detector requires 837 the crista mask to physically intersect the membrane band, this one accepts a crista region that 838 comes within ``max_extension_nm`` of it, so a crista segmented a few nm short of the inner boundary 839 membrane still registers. A TEASAR skeleton (``bioimage_cpp.skeleton.teasar``) then gates each 840 region on being near a crista *terminus*, which is what separates a crista ending at the membrane 841 from one running alongside it. 842 843 A junction is one 26-connected component of ``crista & (distance_to_membrane <= 844 max_extension_nm)``. Taking identity from the contact geometry this way is what makes the result 845 stable, and it is worth recording why, because the two obvious alternatives both fail: 846 847 - **Deriving identity from per-endpoint hits plus a merge radius does not work.** Measured on a real 848 mitochondrion with three junctions, no radius returns three: the count steps 10, 8, 4, 2 as the 849 radius grows, because a small radius fragments one junction while a large one fuses two that are 850 only 13.9 nm apart. 851 - **Extending each skeleton end along its tangent does not find the junctions.** On the same data 852 the direction from a skeleton end to its nearest real contact was 142-146 degrees away from that 853 end's tangent — pointing backwards — for all three. A crista-membrane contact is a *rim* feature 854 while a skeleton end is a *centerline* feature, and for a sheet meeting the membrane obliquely 855 their directions are unrelated. Widening the ray to a +/-60 degree cone changed nothing; only an 856 omnidirectional search found all three, which is a proximity test in disguise. That is this 857 function. 858 859 Because a component of the near-membrane mask is a subset of the crista mask, two disconnected 860 cristae can never be merged into one junction, and no tuning parameter governs that. The dual also 861 holds: the terminus filter is keyed per crista, so one crista can never be validated by another 862 one's end (see ``terminus_nm``). 863 864 **Known limitation — this mode over-detects on densely packed cristae.** Proximity is not the same 865 as junction: on a real mitochondrion with many cristae (TS_PS_01 mito 1, 0.8681 nm voxels, 8 nm 866 membrane) it reports 21 junctions of which only 2 involve any literal crista-membrane contact; the 867 other 19 are cristae merely passing within 8 nm of the inner boundary membrane. The false positives 868 are full-sized (up to ~1750 nm3), so ``min_junction_volume_nm3`` does not remove them. **Five** 869 discriminators have now been measured against real data and none separates the two populations: 870 region elongation (real junctions are 1.5-2.4 elongated too), region axis versus the membrane normal 871 (75-90 degrees for every region, because a contact patch spreads along the membrane by nature), 872 the rate at which membrane distance drops toward the terminus (0.42-0.73 for every region), a 873 minimum size, and the crista **sheet normal** versus the membrane normal (below). Treat the count as 874 an upper bound on a dense mitochondrion and inspect the result -- the widget's 875 **Show Crista Skeleton** layers exist for exactly that. Use ``"overlap"`` when only junctions with 876 actual contact should count. 877 878 **The sheet-normal discriminator was implemented, measured and removed.** The idea: a crista is a 879 lamella, so compare its sheet normal (the gradient of the mask smoothed at the sheet thickness) to 880 the membrane normal (the gradient of the reference distance field). A crista running parallel to the 881 membrane should give ``|cos| -> 1`` and one meeting it end-on ``|cos| -> 0``. On synthetic geometry 882 it behaved exactly so, 0.00 against 0.70. On the ``cutout_mito2`` cristae with three hand-verified 883 contacts it **inverted**: sampled at each region's closest-approach voxel the two contacting regions 884 scored 0.786 and 0.881 while the two non-contacting ones scored 0.000 and 0.766, so a threshold of 885 0.7 removed all three real junctions and kept the false positive. Region means do not separate 886 either (0.582/0.606 against 0.558/0.842). Part of the reason is that the measure is ill-posed 887 precisely where one wants to sample it: at a closest-approach voxel the distance field can be 888 locally flat, and a vanishing gradient normalises to a meaningless direction -- the exact 0.000 889 above is that artefact. Reproduce with ``scripts/cooper/measure_terminus_alignment.py``. Do not 890 re-propose it without new evidence from that script. 891 892 Args: 893 crista_mask: Binary crista segmentation. Must be 3D — TEASAR has no 2D implementation. 894 membrane_mask: Binary mitochondrial membrane mask, e.g. from :func:`approximate_membrane`. 895 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 896 max_extension_nm: How far (nm) a crista may fall short of the membrane and still count. Set it 897 to about one membrane thickness; an over-large value starts flagging cristae that merely 898 pass near the boundary, and it also fuses junctions whose near-membrane regions touch. 899 min_extension_nm: Smallest gap (nm) that counts, measured as a region's *closest* approach. 900 0 (the default) accepts regions already overlapping the membrane; raise it to isolate 901 non-overlapping junctions only. 902 terminus_nm: A region is kept only if it lies within this distance (nm) of a skeleton end **of 903 its own 26-connected crista**, so a crista running alongside the membrane is rejected. 904 Keying on the region's own crista is what stops a neighbouring crista's end from validating 905 it; with a single pooled terminus tree a region whose own skeleton was dropped by 906 ``min_skeleton_nm`` could still pass. On real data this changes nothing — measured 4 -> 4 907 on ``cutout_mito2``, whose crista mask merges into just 5 connected components — so do not 908 "simplify" it back out on the grounds that it makes no difference there. The leak needs a 909 genuinely isolated crista, which sparse and synthetic masks readily produce; see 910 ``test_a_crista_cannot_borrow_another_cristas_terminus``. The 20 nm default is generous by design: 911 the measured terminus distances of real junctions were 0.0, 0.0 and 8.6 nm, so it rejects 912 only cristae running well clear of any end, such as one sliding along the membrane. Pass 913 ``inf`` to disable the filter and skip TEASAR entirely. 914 min_junction_volume_nm3: Regions smaller than this are dropped. This removes specks only — real 915 data produces 1-, 3- and 12-voxel regions that are not junctions on any reading, while real 916 junctions measure hundreds to thousands of nm^3. It is **not** a remedy for the 917 over-detection described above, whose false positives are full-sized. 918 border_radius: Width in voxels of the volume-border zone to exclude, matching the membrane's 919 own border-gap trim (``_gap_radius``). 0 (the default) disables the exclusion and 920 reproduces the misaligned behaviour, so **every production caller passes it**: a crista 921 within this distance of a clipped face would otherwise be matched against membrane that 922 :func:`approximate_membrane` deliberately removed as unknown — the distance transform would 923 measure straight across the deleted region to the nearest surviving membrane voxel and 924 assert a junction against a membrane it was told nothing about. Regions are trimmed rather 925 than discarded, so a crista entering the unknown zone still counts wherever else it 926 genuinely reaches the membrane. May also be an ``(ndim, 2)`` array of **per-face** radii, 927 which is what a caller passing a bounding-box crop must do — the zone is a global fact and 928 a crop starting one voxel inside the volume still overlaps it. See :func:`_border_zone`. 929 boundary: Optional ``(ndim, 2)`` bool forwarded to :func:`_border_zone`, zeroing the radius on 930 faces that are not genuine volume faces; None means all faces. It is the *mesh* question, 931 and per-face ``border_radius`` already encodes it, so **pass one or the other, never 932 both** — the radii are multiplied by this mask. 933 n_jobs: Forwarded to TEASAR's ``number_of_threads`` (-1/0 = all cores). 934 membrane_distance: Optional precomputed ``distance_transform(~membrane_mask)`` in nm, to avoid 935 recomputing it when the caller already has one (see :func:`_single_mito_row`). Must match 936 ``membrane_mask`` and use the same sampling. Used only as the fallback reference when no 937 ``lumen_mask`` is given. 938 lumen_mask: The eroded-mito interior from ``approximate_membrane(..., return_lumen=True)``. 939 When given, distances are measured to the **inner boundary membrane surface** derived from 940 it rather than to the membrane band, which is what lets a junction be localised at all — 941 see :func:`_inner_surface_distance`. Pass only a genuine lumen: ``mito & ~membrane`` is not 942 one near a clipped face, where it re-includes the mito's outer shell. 943 footprint_nm: How far (nm) beyond a region's closest approach the painted label extends. This 944 sets the junction label's thickness and therefore its centroid, but never the count. None 945 (the default) means **one voxel diagonal**, which cannot be written as a literal here since 946 it depends on ``voxel_size``: it is the tightest tolerance that still captures a contact 947 patch lying oblique to the grid. Measured at 1.5 nm isotropic voxels it keeps the footprint 948 at 6% of the candidate region, where 4 nm would already take 56% of it. 949 min_skeleton_nm: Forwarded to :func:`compute_crista_skeleton` — skeleton components shorter 950 than this are dropped, so a crista smaller than this contributes no terminus and cannot 951 score a junction (which holds only because the terminus trees are per-crista; see 952 ``terminus_nm``). None leaves that function's own default in force. 953 terminus_merge_nm: Forwarded to :func:`compute_crista_skeleton` — the radius within which 954 termini on one component collapse to a single representative. None leaves that function's 955 own default in force. 956 957 Returns: 958 contact_labels: Integer array (same shape as the input) where each junction has a unique ID 959 (0 = background, 1..n = junctions) — interchangeable with the first return value of 960 :func:`detect_contact_sites`, so it feeds :func:`compute_junction_distances` and the napari 961 junction layer unchanged. The labelled region is the junction's **closest-approach 962 footprint**: the voxels of the candidate region within ``footprint_nm`` of its nearest 963 approach to the reference surface, not the whole region. It is therefore a thin patch at 964 the membrane rather than a slab of crista (measured: 80 voxels against 1360 for the region), 965 it lies inside the crista and hence inside the mitochondrion, and its centroid is a usable 966 junction position for the geodesic stage. 967 summary: contact_voxel_count, crista_junction_count, contact_volume_nm3, 968 mean_junction_extension_nm (the mean over junctions of each region's closest approach to 969 the reference surface; 0 where the crista reaches it). Only ``crista_junction_count`` and 970 ``mean_junction_extension_nm`` are mode-specific; ``contact_voxel_count`` and 971 ``contact_volume_nm3`` still describe the genuine crista-membrane overlap, so those two 972 columns mean the same thing in both modes. 973 974 Raises: 975 ValueError: If the input is not 3D, the voxel size is not positive on every axis, or the 976 extension range is negative / inverted. 977 """ 978 if crista_mask.ndim != 3: 979 raise ValueError( 980 "detect_junctions_skeleton requires a 3D volume (bioimage_cpp.skeleton.teasar has no 2D " 981 f"implementation), got ndim={crista_mask.ndim}. Use junction_mode='overlap' for 2D data." 982 ) 983 if min_extension_nm < 0 or max_extension_nm < min_extension_nm: 984 raise ValueError( 985 "need 0 <= min_extension_nm <= max_extension_nm, got " 986 f"min_extension_nm={min_extension_nm}, max_extension_nm={max_extension_nm}" 987 ) 988 989 sampling = _to_sampling(voxel_size, 3) 990 if not np.all(sampling > 0): 991 raise ValueError(f"voxel_size must be positive on every axis, got {voxel_size!r}") 992 voxel_vol = float(np.prod(sampling)) 993 crista = crista_mask.astype(bool) 994 membrane = membrane_mask.astype(bool) 995 996 overlap_voxels = int(np.count_nonzero(crista & membrane)) 997 summary = { 998 "contact_voxel_count": overlap_voxels, 999 "crista_junction_count": 0, 1000 "contact_volume_nm3": float(overlap_voxels) * voxel_vol, 1001 "mean_junction_extension_nm": np.nan, 1002 } 1003 labels = np.zeros(crista.shape, dtype=np.int32) 1004 if not crista.any() or not membrane.any(): 1005 return labels, summary 1006 1007 if membrane_distance is None: 1008 membrane_distance = distance_transform( 1009 ~membrane, sampling=sampling.tolist(), number_of_threads=_solver_threads(n_jobs) 1010 ) 1011 reference_distance = membrane_distance 1012 if lumen_mask is not None and lumen_mask.any(): 1013 reference_distance = _inner_surface_distance(lumen_mask.astype(bool), sampling, n_jobs) 1014 if footprint_nm is None: 1015 footprint_nm = float(np.linalg.norm(sampling)) 1016 1017 near = crista & (reference_distance <= max_extension_nm) 1018 if np.any(np.asarray(border_radius) >= 1): 1019 near &= ~_border_zone(crista.shape, border_radius, boundary) 1020 if not near.any(): 1021 return labels, summary 1022 1023 region_labels, n_regions = ndimage_label(near, structure=np.ones(3 * (3,), dtype=bool)) 1024 1025 # One KD-tree PER CRISTA, not one for all of them: a region is validated only by a terminus of 1026 # its own 26-connected crista. A single pooled tree let an unrelated crista ending nearby accept a 1027 # region whose own skeleton was dropped by ``min_skeleton_nm`` — measured on the test lamella, a 1028 # sheet alone gave 1 junction, an 18-voxel blob alone 0, and the two together 2. 1029 terminus_trees = None 1030 if np.isfinite(terminus_nm): 1031 vertices, is_terminus = compute_crista_skeleton( 1032 crista, voxel_size, n_jobs=n_jobs, 1033 **_given(min_skeleton_nm=min_skeleton_nm, terminus_merge_nm=terminus_merge_nm), 1034 ) 1035 endpoints = vertices[is_terminus] 1036 if len(endpoints) == 0: 1037 return labels, summary 1038 from scipy.spatial import cKDTree 1039 1040 crista_components, _ = ndimage_label(crista, structure=np.ones(3 * (3,), dtype=bool)) 1041 # TEASAR returns physical coordinates, i.e. index * spacing exactly (measured residual 0.0 at 1042 # isotropic and 3.6e-15 at anisotropic spacing), so rounding recovers the index and always 1043 # lands on a foreground voxel — a terminus can never be attributed to component 0. 1044 endpoint_component = crista_components[tuple(np.rint(endpoints / sampling).astype(int).T)] 1045 terminus_trees = { 1046 int(component): cKDTree(endpoints[endpoint_component == component]) 1047 for component in np.unique(endpoint_component) 1048 } 1049 # ``near`` is a subset of ``crista`` under the same connectivity, so each region lies in 1050 # exactly one component. Collapse the lookup to one entry per region and free the int32 1051 # volume (120 MB on a 401x257x290 tomogram) before the loop. 1052 region_component = np.zeros(n_regions + 1, dtype=np.int32) 1053 region_component[region_labels[near]] = crista_components[near] 1054 del crista_components 1055 1056 assigned = 0 1057 gaps = [] 1058 for region in range(1, n_regions + 1): 1059 region_mask = region_labels == region 1060 if float(np.count_nonzero(region_mask)) * voxel_vol < min_junction_volume_nm3: 1061 continue 1062 region_distance = reference_distance[region_mask] 1063 gap = float(region_distance.min()) 1064 if not (min_extension_nm <= gap <= max_extension_nm): 1065 continue 1066 if terminus_trees is not None: 1067 tree = terminus_trees.get(int(region_component[region])) 1068 if tree is None: # this crista contributed no terminus of its own 1069 continue 1070 # Capped query: past terminus_nm the KD-tree returns inf rather than a real distance. 1071 points = np.argwhere(region_mask) * sampling 1072 hit = tree.query(points, distance_upper_bound=terminus_nm)[0] 1073 if not np.isfinite(hit).any(): 1074 continue 1075 assigned += 1 1076 gaps.append(gap) 1077 labels[region_mask & (reference_distance <= gap + footprint_nm)] = assigned 1078 1079 summary["crista_junction_count"] = assigned 1080 if gaps: 1081 summary["mean_junction_extension_nm"] = float(np.mean(gaps)) 1082 return labels, summary
Detect crista-membrane junctions as crista regions that reach close to the membrane.
The gap-tolerant alternative to detect_contact_sites(). Where the overlap detector requires
the crista mask to physically intersect the membrane band, this one accepts a crista region that
comes within max_extension_nm of it, so a crista segmented a few nm short of the inner boundary
membrane still registers. A TEASAR skeleton (bioimage_cpp.skeleton.teasar) then gates each
region on being near a crista terminus, which is what separates a crista ending at the membrane
from one running alongside it.
A junction is one 26-connected component of crista & (distance_to_membrane <=
max_extension_nm). Taking identity from the contact geometry this way is what makes the result
stable, and it is worth recording why, because the two obvious alternatives both fail:
- Deriving identity from per-endpoint hits plus a merge radius does not work. Measured on a real mitochondrion with three junctions, no radius returns three: the count steps 10, 8, 4, 2 as the radius grows, because a small radius fragments one junction while a large one fuses two that are only 13.9 nm apart.
- Extending each skeleton end along its tangent does not find the junctions. On the same data the direction from a skeleton end to its nearest real contact was 142-146 degrees away from that end's tangent — pointing backwards — for all three. A crista-membrane contact is a rim feature while a skeleton end is a centerline feature, and for a sheet meeting the membrane obliquely their directions are unrelated. Widening the ray to a +/-60 degree cone changed nothing; only an omnidirectional search found all three, which is a proximity test in disguise. That is this function.
Because a component of the near-membrane mask is a subset of the crista mask, two disconnected
cristae can never be merged into one junction, and no tuning parameter governs that. The dual also
holds: the terminus filter is keyed per crista, so one crista can never be validated by another
one's end (see terminus_nm).
Known limitation — this mode over-detects on densely packed cristae. Proximity is not the same
as junction: on a real mitochondrion with many cristae (TS_PS_01 mito 1, 0.8681 nm voxels, 8 nm
membrane) it reports 21 junctions of which only 2 involve any literal crista-membrane contact; the
other 19 are cristae merely passing within 8 nm of the inner boundary membrane. The false positives
are full-sized (up to ~1750 nm3), so min_junction_volume_nm3 does not remove them. Five
discriminators have now been measured against real data and none separates the two populations:
region elongation (real junctions are 1.5-2.4 elongated too), region axis versus the membrane normal
(75-90 degrees for every region, because a contact patch spreads along the membrane by nature),
the rate at which membrane distance drops toward the terminus (0.42-0.73 for every region), a
minimum size, and the crista sheet normal versus the membrane normal (below). Treat the count as
an upper bound on a dense mitochondrion and inspect the result -- the widget's
Show Crista Skeleton layers exist for exactly that. Use "overlap" when only junctions with
actual contact should count.
The sheet-normal discriminator was implemented, measured and removed. The idea: a crista is a
lamella, so compare its sheet normal (the gradient of the mask smoothed at the sheet thickness) to
the membrane normal (the gradient of the reference distance field). A crista running parallel to the
membrane should give |cos| -> 1 and one meeting it end-on |cos| -> 0. On synthetic geometry
it behaved exactly so, 0.00 against 0.70. On the cutout_mito2 cristae with three hand-verified
contacts it inverted: sampled at each region's closest-approach voxel the two contacting regions
scored 0.786 and 0.881 while the two non-contacting ones scored 0.000 and 0.766, so a threshold of
0.7 removed all three real junctions and kept the false positive. Region means do not separate
either (0.582/0.606 against 0.558/0.842). Part of the reason is that the measure is ill-posed
precisely where one wants to sample it: at a closest-approach voxel the distance field can be
locally flat, and a vanishing gradient normalises to a meaningless direction -- the exact 0.000
above is that artefact. Reproduce with scripts/cooper/measure_terminus_alignment.py. Do not
re-propose it without new evidence from that script.
Arguments:
- crista_mask: Binary crista segmentation. Must be 3D — TEASAR has no 2D implementation.
- membrane_mask: Binary mitochondrial membrane mask, e.g. from
approximate_membrane(). - voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- max_extension_nm: How far (nm) a crista may fall short of the membrane and still count. Set it to about one membrane thickness; an over-large value starts flagging cristae that merely pass near the boundary, and it also fuses junctions whose near-membrane regions touch.
- min_extension_nm: Smallest gap (nm) that counts, measured as a region's closest approach. 0 (the default) accepts regions already overlapping the membrane; raise it to isolate non-overlapping junctions only.
- terminus_nm: A region is kept only if it lies within this distance (nm) of a skeleton end of
its own 26-connected crista, so a crista running alongside the membrane is rejected.
Keying on the region's own crista is what stops a neighbouring crista's end from validating
it; with a single pooled terminus tree a region whose own skeleton was dropped by
min_skeleton_nmcould still pass. On real data this changes nothing — measured 4 -> 4 oncutout_mito2, whose crista mask merges into just 5 connected components — so do not "simplify" it back out on the grounds that it makes no difference there. The leak needs a genuinely isolated crista, which sparse and synthetic masks readily produce; seetest_a_crista_cannot_borrow_another_cristas_terminus. The 20 nm default is generous by design: the measured terminus distances of real junctions were 0.0, 0.0 and 8.6 nm, so it rejects only cristae running well clear of any end, such as one sliding along the membrane. Passinfto disable the filter and skip TEASAR entirely. - min_junction_volume_nm3: Regions smaller than this are dropped. This removes specks only — real data produces 1-, 3- and 12-voxel regions that are not junctions on any reading, while real junctions measure hundreds to thousands of nm^3. It is not a remedy for the over-detection described above, whose false positives are full-sized.
- border_radius: Width in voxels of the volume-border zone to exclude, matching the membrane's
own border-gap trim (
_gap_radius). 0 (the default) disables the exclusion and reproduces the misaligned behaviour, so every production caller passes it: a crista within this distance of a clipped face would otherwise be matched against membrane thatapproximate_membrane()deliberately removed as unknown — the distance transform would measure straight across the deleted region to the nearest surviving membrane voxel and assert a junction against a membrane it was told nothing about. Regions are trimmed rather than discarded, so a crista entering the unknown zone still counts wherever else it genuinely reaches the membrane. May also be an(ndim, 2)array of per-face radii, which is what a caller passing a bounding-box crop must do — the zone is a global fact and a crop starting one voxel inside the volume still overlaps it. See_border_zone(). - boundary: Optional
(ndim, 2)bool forwarded to_border_zone(), zeroing the radius on faces that are not genuine volume faces; None means all faces. It is the mesh question, and per-faceborder_radiusalready encodes it, so pass one or the other, never both — the radii are multiplied by this mask. - n_jobs: Forwarded to TEASAR's
number_of_threads(-1/0 = all cores). - membrane_distance: Optional precomputed
distance_transform(~membrane_mask)in nm, to avoid recomputing it when the caller already has one (see_single_mito_row()). Must matchmembrane_maskand use the same sampling. Used only as the fallback reference when nolumen_maskis given. - lumen_mask: The eroded-mito interior from
approximate_membrane(..., return_lumen=True). When given, distances are measured to the inner boundary membrane surface derived from it rather than to the membrane band, which is what lets a junction be localised at all — see_inner_surface_distance(). Pass only a genuine lumen:mito & ~membraneis not one near a clipped face, where it re-includes the mito's outer shell. - footprint_nm: How far (nm) beyond a region's closest approach the painted label extends. This
sets the junction label's thickness and therefore its centroid, but never the count. None
(the default) means one voxel diagonal, which cannot be written as a literal here since
it depends on
voxel_size: it is the tightest tolerance that still captures a contact patch lying oblique to the grid. Measured at 1.5 nm isotropic voxels it keeps the footprint at 6% of the candidate region, where 4 nm would already take 56% of it. - min_skeleton_nm: Forwarded to
compute_crista_skeleton()— skeleton components shorter than this are dropped, so a crista smaller than this contributes no terminus and cannot score a junction (which holds only because the terminus trees are per-crista; seeterminus_nm). None leaves that function's own default in force. - terminus_merge_nm: Forwarded to
compute_crista_skeleton()— the radius within which termini on one component collapse to a single representative. None leaves that function's own default in force.
Returns:
contact_labels: Integer array (same shape as the input) where each junction has a unique ID (0 = background, 1..n = junctions) — interchangeable with the first return value of
detect_contact_sites(), so it feedscompute_junction_distances()and the napari junction layer unchanged. The labelled region is the junction's closest-approach footprint: the voxels of the candidate region withinfootprint_nmof its nearest approach to the reference surface, not the whole region. It is therefore a thin patch at the membrane rather than a slab of crista (measured: 80 voxels against 1360 for the region), it lies inside the crista and hence inside the mitochondrion, and its centroid is a usable junction position for the geodesic stage. summary: contact_voxel_count, crista_junction_count, contact_volume_nm3, mean_junction_extension_nm (the mean over junctions of each region's closest approach to the reference surface; 0 where the crista reaches it). Onlycrista_junction_countandmean_junction_extension_nmare mode-specific;contact_voxel_countandcontact_volume_nm3still describe the genuine crista-membrane overlap, so those two columns mean the same thing in both modes.
Raises:
- ValueError: If the input is not 3D, the voxel size is not positive on every axis, or the extension range is negative / inverted.
def
detect_junctions( crista_mask: numpy.ndarray, membrane_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], junction_mode: str = 'overlap', max_extension_nm: float = 8.0, n_jobs: int = 1, **kwargs) -> Tuple[numpy.ndarray, Dict[str, float]]:
1085def detect_junctions( 1086 crista_mask: np.ndarray, 1087 membrane_mask: np.ndarray, 1088 voxel_size: Union[float, Dict[str, float]], 1089 junction_mode: str = "overlap", 1090 max_extension_nm: float = 8.0, 1091 n_jobs: int = 1, 1092 **kwargs, 1093) -> Tuple[np.ndarray, Dict[str, float]]: 1094 """Detect crista-membrane junctions with the chosen algorithm. 1095 1096 Args: 1097 crista_mask: Binary crista segmentation. 1098 membrane_mask: Binary mitochondrial membrane mask. 1099 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1100 junction_mode: ``"overlap"`` (default) counts connected components of the direct 1101 crista-membrane intersection via :func:`detect_contact_sites`; ``"skeleton"`` counts 1102 crista regions that come within ``max_extension_nm`` of the membrane near a crista 1103 terminus, via :func:`detect_junctions_skeleton`. 1104 max_extension_nm: How far (nm) a crista may fall short of the membrane; ``"skeleton"`` only. 1105 n_jobs: Thread budget; ``"skeleton"`` only. 1106 **kwargs: Further :func:`detect_junctions_skeleton` options (``min_extension_nm``, 1107 ``terminus_nm``, ``min_junction_volume_nm3``, ``border_radius``, ``boundary``, 1108 ``membrane_distance``, ``lumen_mask``, ``footprint_nm``, 1109 ``min_skeleton_nm``, ``terminus_merge_nm``); ignored in ``"overlap"`` mode, which needs no 1110 border parameter because it is border-safe by construction. Any of them passed as ``None`` 1111 is dropped here (:func:`_given`), so a caller that treats ``None`` as "unset" — the widget's 1112 zeroed spin-boxes, the CLI's unset flags, the per-mito path — gets the detector's own 1113 default rather than having to restate it. 1114 1115 Returns: 1116 (contact_labels, summary) as documented on the two backends. ``summary`` always carries 1117 ``mean_junction_extension_nm`` (NaN in ``"overlap"`` mode, which measures no gap) so callers 1118 can read it without branching on the mode. 1119 1120 Raises: 1121 ValueError: If ``junction_mode`` is not one of the two supported values. 1122 """ 1123 if junction_mode == "overlap": 1124 contact_labels, summary = detect_contact_sites(crista_mask, membrane_mask, voxel_size) 1125 summary["mean_junction_extension_nm"] = np.nan 1126 return contact_labels, summary 1127 if junction_mode == "skeleton": 1128 return detect_junctions_skeleton( 1129 crista_mask, membrane_mask, voxel_size, 1130 max_extension_nm=max_extension_nm, n_jobs=n_jobs, **_given(**kwargs), 1131 ) 1132 raise ValueError(f"junction_mode must be 'overlap' or 'skeleton', got {junction_mode!r}")
Detect crista-membrane junctions with the chosen algorithm.
Arguments:
- crista_mask: Binary crista segmentation.
- membrane_mask: Binary mitochondrial membrane mask.
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- junction_mode:
"overlap"(default) counts connected components of the direct crista-membrane intersection viadetect_contact_sites();"skeleton"counts crista regions that come withinmax_extension_nmof the membrane near a crista terminus, viadetect_junctions_skeleton(). - max_extension_nm: How far (nm) a crista may fall short of the membrane;
"skeleton"only. - n_jobs: Thread budget;
"skeleton"only. - **kwargs: Further
detect_junctions_skeleton()options (min_extension_nm,terminus_nm,min_junction_volume_nm3,border_radius,boundary,membrane_distance,lumen_mask,footprint_nm,min_skeleton_nm,terminus_merge_nm); ignored in"overlap"mode, which needs no border parameter because it is border-safe by construction. Any of them passed asNoneis dropped here (_given()), so a caller that treatsNoneas "unset" — the widget's zeroed spin-boxes, the CLI's unset flags, the per-mito path — gets the detector's own default rather than having to restate it.
Returns:
(contact_labels, summary) as documented on the two backends.
summaryalways carriesmean_junction_extension_nm(NaN in"overlap"mode, which measures no gap) so callers can read it without branching on the mode.
Raises:
- ValueError: If
junction_modeis not one of the two supported values.
def
detect_junctions_per_mito( crista_mask: numpy.ndarray, mito_segmentation: numpy.ndarray, membrane_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], lumen_mask: Optional[numpy.ndarray] = None, border_radius: int = 0, n_jobs: int = 1, **kwargs) -> Tuple[numpy.ndarray, Dict[str, float]]:
1136def detect_junctions_per_mito( 1137 crista_mask: np.ndarray, 1138 mito_segmentation: np.ndarray, 1139 membrane_mask: np.ndarray, 1140 voxel_size: Union[float, Dict[str, float]], 1141 lumen_mask: Optional[np.ndarray] = None, 1142 border_radius: int = 0, 1143 n_jobs: int = 1, 1144 **kwargs, 1145) -> Tuple[np.ndarray, Dict[str, float]]: 1146 """Detect junctions the way the statistics table does: **per mitochondrial instance**. 1147 1148 :func:`detect_junctions` run once over a whole volume and the same function run per instance are 1149 not the same measurement, and in ``"skeleton"`` mode they routinely disagree. Restricting a crista 1150 to one mitochondrion clips it, and clipping moves its skeleton endpoints: a crista tube running 1151 through a mitochondrion and out the other side has its global endpoints far away from the membrane 1152 (no junction), while the clipped tube ends *at* the mitochondrial boundary (two junctions). Any 1153 caller that displays junctions next to numbers from 1154 :func:`compute_mito_crista_statistics` must therefore detect them this way, or the picture and the 1155 table disagree at the same settings. 1156 1157 It is also **cheaper** than the global pass, not more expensive, which is unintuitive enough to be 1158 worth recording: the cost is dominated by two whole-volume distance transforms, and mitochondrial 1159 bounding boxes sum to a fraction of a tomogram (measured 0.06 s against 0.33 s, a 0.19x ratio, on 1160 4 box mitos whose bboxes covered 14% of an 80x160x160 volume — same 30 junctions either way). 1161 TEASAR is per connected component either way, so splitting the mask by instance multiplies no 1162 skeletonisation work. 1163 1164 The per-instance masking and the crop-aware border radii here mirror :func:`_single_mito_row` 1165 exactly; ``test_per_mito_detector_matches_the_table`` pins the two together. 1166 1167 Args: 1168 crista_mask: Binary crista segmentation (global volume). 1169 mito_segmentation: Instance label array (background = 0). 1170 membrane_mask: Binary mitochondrial membrane mask, e.g. from :func:`approximate_membrane`. 1171 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1172 lumen_mask: Optional eroded-mito interior matching ``membrane_mask``; restricted to each 1173 instance, as in :func:`_single_mito_row`. 1174 border_radius: Volume-border zone width in voxels (``_gap_radius``). Converted to crop-aware 1175 per-face radii for each bounding box — see :func:`_border_zone`. 1176 n_jobs: Thread budget forwarded to the detector. 1177 **kwargs: Further :func:`detect_junctions` options (``junction_mode``, ``max_extension_nm``, 1178 ``terminus_nm``, ...). 1179 1180 Returns: 1181 (junction_labels, summary) with the same keys as :func:`detect_junctions`. Ids are 1182 per-instance and therefore repeat between mitochondria; the painted voxels do not, since each 1183 lies inside its own mito. ``mean_junction_extension_nm`` is averaged over the instances that 1184 reported one. 1185 """ 1186 ndim = mito_segmentation.ndim 1187 crista_binary = crista_mask.astype(bool) 1188 membrane_binary = membrane_mask.astype(bool) 1189 vol_shape = mito_segmentation.shape 1190 1191 labels = np.zeros(vol_shape, dtype=np.int32) 1192 summary = { 1193 "contact_voxel_count": 0, "crista_junction_count": 0, "contact_volume_nm3": 0.0, 1194 "mean_junction_extension_nm": np.nan, 1195 } 1196 extensions = [] 1197 1198 for prop in regionprops(mito_segmentation): 1199 bbox = prop.bbox 1200 slices = tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim)) 1201 mito_local = mito_segmentation[slices] == prop.label 1202 crista_local = crista_binary[slices] & mito_local 1203 membrane_local = membrane_binary[slices] & mito_local 1204 if not crista_local.any() or not membrane_local.any(): 1205 continue 1206 border_radii = np.array( 1207 [[max(0, border_radius - bbox[a]), 1208 max(0, border_radius - (vol_shape[a] - bbox[a + ndim]))] for a in range(ndim)], dtype=int 1209 ) 1210 local_labels, local_summary = detect_junctions( 1211 crista_local, membrane_local, voxel_size, n_jobs=n_jobs, 1212 border_radius=border_radii, 1213 lumen_mask=None if lumen_mask is None else lumen_mask[slices] & mito_local, 1214 **kwargs, 1215 ) 1216 written = local_labels > 0 1217 labels[slices][written] = local_labels[written] # masked: bounding boxes overlap 1218 summary["contact_voxel_count"] += local_summary["contact_voxel_count"] 1219 summary["contact_volume_nm3"] += local_summary["contact_volume_nm3"] 1220 summary["crista_junction_count"] += local_summary["crista_junction_count"] 1221 if np.isfinite(local_summary["mean_junction_extension_nm"]): 1222 extensions.append(local_summary["mean_junction_extension_nm"]) 1223 1224 if extensions: 1225 summary["mean_junction_extension_nm"] = float(np.mean(extensions)) 1226 return labels, summary
Detect junctions the way the statistics table does: per mitochondrial instance.
detect_junctions() run once over a whole volume and the same function run per instance are
not the same measurement, and in "skeleton" mode they routinely disagree. Restricting a crista
to one mitochondrion clips it, and clipping moves its skeleton endpoints: a crista tube running
through a mitochondrion and out the other side has its global endpoints far away from the membrane
(no junction), while the clipped tube ends at the mitochondrial boundary (two junctions). Any
caller that displays junctions next to numbers from
compute_mito_crista_statistics() must therefore detect them this way, or the picture and the
table disagree at the same settings.
It is also cheaper than the global pass, not more expensive, which is unintuitive enough to be worth recording: the cost is dominated by two whole-volume distance transforms, and mitochondrial bounding boxes sum to a fraction of a tomogram (measured 0.06 s against 0.33 s, a 0.19x ratio, on 4 box mitos whose bboxes covered 14% of an 80x160x160 volume — same 30 junctions either way). TEASAR is per connected component either way, so splitting the mask by instance multiplies no skeletonisation work.
The per-instance masking and the crop-aware border radii here mirror _single_mito_row()
exactly; test_per_mito_detector_matches_the_table pins the two together.
Arguments:
- crista_mask: Binary crista segmentation (global volume).
- mito_segmentation: Instance label array (background = 0).
- membrane_mask: Binary mitochondrial membrane mask, e.g. from
approximate_membrane(). - voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- lumen_mask: Optional eroded-mito interior matching
membrane_mask; restricted to each instance, as in_single_mito_row(). - border_radius: Volume-border zone width in voxels (
_gap_radius). Converted to crop-aware per-face radii for each bounding box — see_border_zone(). - n_jobs: Thread budget forwarded to the detector.
- **kwargs: Further
detect_junctions()options (junction_mode,max_extension_nm,terminus_nm, ...).
Returns:
(junction_labels, summary) with the same keys as
detect_junctions(). Ids are per-instance and therefore repeat between mitochondria; the painted voxels do not, since each lies inside its own mito.mean_junction_extension_nmis averaged over the instances that reported one.
def
compute_junction_distances( contact_labels: numpy.ndarray, membrane_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], surface_area_nm2: Optional[float] = None, n_jobs: int = 1, mesh_vertices: Optional[numpy.ndarray] = None, mesh_faces: Optional[numpy.ndarray] = None) -> Tuple[numpy.ndarray, Dict[str, float]]:
1284def compute_junction_distances( 1285 contact_labels: np.ndarray, 1286 membrane_mask: np.ndarray, 1287 voxel_size: Union[float, Dict[str, float]], 1288 surface_area_nm2: Optional[float] = None, 1289 n_jobs: int = 1, 1290 mesh_vertices: Optional[np.ndarray] = None, 1291 mesh_faces: Optional[np.ndarray] = None, 1292) -> Tuple[np.ndarray, Dict[str, float]]: 1293 """Geodesic distances between crista-membrane junctions along the eroded-mito surface mesh. 1294 1295 Each junction (a connected component in ``contact_labels``) is reduced to its centroid, snapped to 1296 the nearest vertex of a triangle mesh, and pairwise surface geodesics are computed with 1297 ``bioimage_cpp.distance.geodesic_distances_mesh``. The mesh is the **eroded-mito (lumen) surface** 1298 passed in as ``mesh_vertices``/``mesh_faces`` by :func:`_single_mito_row` (a clean, single-wall 1299 surface at the membrane's inner edge). If no mesh is supplied — or no usable surface mesh exists 1300 (empty membrane / degenerate mesh) — the junction distances are NaN. (There is no membrane-band 1301 fallback mesh: the metric is defined on the lumen surface, and meshing the thick membrane band 1302 would give a different, capped double-wall surface.) 1303 1304 A Clark-Evans nearest-neighbour index summarises whether the junctions are clustered. 1305 1306 Args: 1307 contact_labels: Integer junction label array (0 = background, 1..n = junctions), 1308 e.g. the first return value of :func:`detect_contact_sites`. 1309 membrane_mask: Binary mitochondrial membrane mask the junctions sit on. Only used for the 1310 empty-membrane early-out (no membrane → NaN); it is not meshed. 1311 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1312 surface_area_nm2: Membrane/mito surface area used as the reference area for the 1313 Clark-Evans expectation. If None or non-positive, the clustering index is NaN. 1314 n_jobs: 1 = serial, -1 = all cores (forwarded to the mesh solver's thread count). 1315 mesh_vertices: Optional (n_vertices, 3) mesh vertices in nm (unpadded mask frame; see 1316 :func:`_surface_mesh`) — the eroded-mito (lumen) surface. If omitted, the junction 1317 distances are NaN. 1318 mesh_faces: Optional (n_faces, 3) triangle indices matching ``mesh_vertices``. 1319 1320 Returns: 1321 distance_matrix: (n, n) geodesic distances in nm between junctions; the diagonal is 1322 0 and unreachable pairs (disconnected fragments) are NaN. Empty for fewer than two 1323 junctions. 1324 summary: junction_count, mean_nn_junction_distance_nm, median_nn_junction_distance_nm, 1325 junction_clustering_index (Clark-Evans R: < 1 clustered, ~ 1 random, > 1 dispersed). 1326 1327 Notes: 1328 The Clark-Evans expected nearest-neighbour distance uses the standard 2-D planar 1329 approximation ``0.5 * sqrt(A / n)`` with ``A = surface_area_nm2``. 1330 1331 Each junction's nearest-neighbour distance is the smallest distance to a *reachable* other 1332 junction: self (diagonal) and unreachable (NaN) pairs are set to +inf before the per-row 1333 minimum, and rows with no reachable neighbour (min stays +inf) are dropped. 1334 """ 1335 ndim = contact_labels.ndim 1336 sampling = _to_sampling(voxel_size, ndim) 1337 membrane = membrane_mask.astype(bool) 1338 1339 labels = [lbl for lbl in np.unique(contact_labels) if lbl != 0] 1340 n = len(labels) 1341 if n < 2 or not membrane.any(): 1342 summary = dict(_JUNCTION_DISTANCE_NAN) 1343 summary["junction_count"] = n 1344 return np.zeros((n, n), dtype=float), summary 1345 1346 centroids = np.atleast_2d( 1347 np.asarray(center_of_mass(contact_labels > 0, labels=contact_labels, index=labels), dtype=float) 1348 ) 1349 1350 if mesh_vertices is not None and mesh_faces is not None and len(mesh_faces) > 0: 1351 distance_matrix = _junction_matrix_mesh(centroids, sampling, mesh_vertices, mesh_faces, n_jobs) 1352 else: 1353 distance_matrix = None 1354 1355 if distance_matrix is None: 1356 summary = dict(_JUNCTION_DISTANCE_NAN) 1357 summary["junction_count"] = n 1358 return np.full((n, n), np.nan, dtype=float), summary 1359 1360 dm = distance_matrix.copy() 1361 np.fill_diagonal(dm, np.inf) 1362 dm[~np.isfinite(dm)] = np.inf 1363 row_min = dm.min(axis=1) 1364 nn_distances = row_min[np.isfinite(row_min)] 1365 1366 mean_nn = float(np.mean(nn_distances)) if nn_distances.size else np.nan 1367 median_nn = float(np.median(nn_distances)) if nn_distances.size else np.nan 1368 1369 clustering_index = np.nan 1370 if surface_area_nm2 is not None and surface_area_nm2 > 0 and np.isfinite(mean_nn): 1371 expected_nn = 0.5 * np.sqrt(float(surface_area_nm2) / n) 1372 if expected_nn > 0: 1373 clustering_index = mean_nn / expected_nn 1374 1375 return distance_matrix, { 1376 "junction_count": n, 1377 "mean_nn_junction_distance_nm": mean_nn, 1378 "median_nn_junction_distance_nm": median_nn, 1379 "junction_clustering_index": clustering_index, 1380 }
Geodesic distances between crista-membrane junctions along the eroded-mito surface mesh.
Each junction (a connected component in contact_labels) is reduced to its centroid, snapped to
the nearest vertex of a triangle mesh, and pairwise surface geodesics are computed with
bioimage_cpp.distance.geodesic_distances_mesh. The mesh is the eroded-mito (lumen) surface
passed in as mesh_vertices/mesh_faces by _single_mito_row() (a clean, single-wall
surface at the membrane's inner edge). If no mesh is supplied — or no usable surface mesh exists
(empty membrane / degenerate mesh) — the junction distances are NaN. (There is no membrane-band
fallback mesh: the metric is defined on the lumen surface, and meshing the thick membrane band
would give a different, capped double-wall surface.)
A Clark-Evans nearest-neighbour index summarises whether the junctions are clustered.
Arguments:
- contact_labels: Integer junction label array (0 = background, 1..n = junctions),
e.g. the first return value of
detect_contact_sites(). - membrane_mask: Binary mitochondrial membrane mask the junctions sit on. Only used for the empty-membrane early-out (no membrane → NaN); it is not meshed.
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- surface_area_nm2: Membrane/mito surface area used as the reference area for the Clark-Evans expectation. If None or non-positive, the clustering index is NaN.
- n_jobs: 1 = serial, -1 = all cores (forwarded to the mesh solver's thread count).
- mesh_vertices: Optional (n_vertices, 3) mesh vertices in nm (unpadded mask frame; see
_surface_mesh()) — the eroded-mito (lumen) surface. If omitted, the junction distances are NaN. - mesh_faces: Optional (n_faces, 3) triangle indices matching
mesh_vertices.
Returns:
distance_matrix: (n, n) geodesic distances in nm between junctions; the diagonal is 0 and unreachable pairs (disconnected fragments) are NaN. Empty for fewer than two junctions. summary: junction_count, mean_nn_junction_distance_nm, median_nn_junction_distance_nm, junction_clustering_index (Clark-Evans R: < 1 clustered, ~ 1 random, > 1 dispersed).
Notes:
The Clark-Evans expected nearest-neighbour distance uses the standard 2-D planar approximation
0.5 * sqrt(A / n)withA = surface_area_nm2.Each junction's nearest-neighbour distance is the smallest distance to a reachable other junction: self (diagonal) and unreachable (NaN) pairs are set to +inf before the per-row minimum, and rows with no reachable neighbour (min stays +inf) are dropped.
def
compute_crista_morphology( crista_mask: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], method: str = 'both') -> Dict[str, float]:
1387def compute_crista_morphology( 1388 crista_mask: np.ndarray, 1389 voxel_size: Union[float, Dict[str, float]], 1390 method: str = "both", 1391) -> Dict[str, float]: 1392 """Compute crista shape metrics from binary mask. 1393 1394 Args: 1395 crista_mask: Binary crista segmentation. 1396 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1397 method: "area" | "medial_axis" | "both". 1398 1399 Returns: 1400 Dict with cristae_surface_area_nm2 (area/both) and avg_thickness_nm (medial_axis/both). 1401 avg_thickness_nm is 2 × mean distance-transform value at skeleton voxels. 1402 """ 1403 if method not in ("area", "medial_axis", "both"): 1404 raise ValueError(f"method must be 'area', 'medial_axis', or 'both', got {method!r}") 1405 1406 sampling = _to_sampling(voxel_size, crista_mask.ndim) 1407 result: Dict[str, float] = {} 1408 1409 if method in ("area", "both"): 1410 result["cristae_surface_area_nm2"] = _surface_area(crista_mask, sampling) 1411 1412 if method in ("medial_axis", "both"): 1413 result["avg_thickness_nm"] = _medial_axis_thickness_nm(crista_mask, sampling) 1414 1415 return result
Compute crista shape metrics from binary mask.
Arguments:
- crista_mask: Binary crista segmentation.
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- method: "area" | "medial_axis" | "both".
Returns:
Dict with cristae_surface_area_nm2 (area/both) and avg_thickness_nm (medial_axis/both). avg_thickness_nm is 2 × mean distance-transform value at skeleton voxels.
def
compute_mito_crista_statistics( crista_mask: numpy.ndarray, mito_segmentation: numpy.ndarray, voxel_size: Union[float, Dict[str, float]], membrane_mask: Optional[numpy.ndarray] = None, membrane_thickness_nm: float = 8.0, border_gap_nm: Optional[float] = None, method: str = 'skip', n_jobs: int = 1, verbose: bool = False, progress_callback: Optional[Callable[[int, int], NoneType]] = None, membrane_mode: str = 'slice_2d', lumen_mask: Optional[numpy.ndarray] = None, junction_mode: str = 'overlap', max_extension_nm: Optional[float] = None, terminus_nm: Optional[float] = None, min_junction_volume_nm3: Optional[float] = None, footprint_nm: Optional[float] = None, min_skeleton_nm: Optional[float] = None, terminus_merge_nm: Optional[float] = None, return_junction_labels: bool = False) -> Union[pandas.DataFrame, Tuple[pandas.DataFrame, numpy.ndarray]]:
1626def compute_mito_crista_statistics( 1627 crista_mask: np.ndarray, 1628 mito_segmentation: np.ndarray, 1629 voxel_size: Union[float, Dict[str, float]], 1630 membrane_mask: Optional[np.ndarray] = None, 1631 membrane_thickness_nm: float = 8.0, 1632 border_gap_nm: Optional[float] = None, 1633 method: str = "skip", 1634 n_jobs: int = 1, 1635 verbose: bool = False, 1636 progress_callback: Optional[Callable[[int, int], None]] = None, 1637 membrane_mode: str = "slice_2d", 1638 lumen_mask: Optional[np.ndarray] = None, 1639 junction_mode: str = "overlap", 1640 max_extension_nm: Optional[float] = None, 1641 terminus_nm: Optional[float] = None, 1642 min_junction_volume_nm3: Optional[float] = None, 1643 footprint_nm: Optional[float] = None, 1644 min_skeleton_nm: Optional[float] = None, 1645 terminus_merge_nm: Optional[float] = None, 1646 return_junction_labels: bool = False, 1647) -> Union[pd.DataFrame, Tuple[pd.DataFrame, np.ndarray]]: 1648 """Compute all crista metrics organised by mitochondrial instance. 1649 1650 Args: 1651 crista_mask: Binary crista segmentation (global volume). 1652 mito_segmentation: Instance label array (background = 0). 1653 voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys. 1654 membrane_mask: Precomputed membrane mask; recomputed if None. 1655 membrane_thickness_nm: Membrane shell thickness used if membrane_mask is None. 1656 border_gap_nm: Border suppression distance passed to approximate_membrane; 1657 defaults to membrane_thickness_nm when None. 1658 method: How the crista orientation anisotropy is computed — this is the ONLY metric that 1659 differs between modes; surface areas (marching cubes), junction distances (geodesic 1660 along the membrane) and thickness/proximity (EDT) are computed identically for all of 1661 them. ``"skip"`` (default) does not compute orientation at all 1662 (``crista_orientation_anisotropy`` is NaN) and is the fastest — use it when only the other 1663 metrics are needed. ``"fast"`` computes the anisotropy on a 2× downsampled crista crop 1664 (~8× cheaper — the structure tensor is by far the dominant cost); the resulting value is a 1665 *relative* indicator that preserves the ordering between mitochondria but is systematically 1666 different in magnitude from the full-resolution value and is NOT comparable to 1667 ``method="exact"``. ``"exact"`` computes the anisotropy from the full-resolution structure 1668 tensor (use it when the magnitude must be precise). 1669 n_jobs: Number of workers for processing mitochondria in parallel (they are 1670 independent). 1 (default) runs serially; other values use a ``concurrent.futures`` 1671 thread pool (-1 = all cores). Results are identical regardless of n_jobs. 1672 verbose: If True, show a terminal tqdm progress bar over mitochondria. 1673 progress_callback: Optional callable invoked once per completed mitochondrion with 1674 (completed_count, total_count) — e.g. to drive a napari progress bar. It is 1675 always called from the calling thread (the futures are consumed here as they 1676 complete), so GUI updates from it need no cross-thread marshaling. 1677 membrane_mode: How the membrane shell is built when ``membrane_mask`` is None — 1678 ``"slice_2d"`` (default, per-Z-slice 2D erosion, z-parallel) or ``"shell_3d"`` (connected 1679 3D shell). See :func:`approximate_membrane`. 1680 lumen_mask: Optional eroded-mito interior matching ``membrane_mask``, i.e. the second return 1681 value of ``approximate_membrane(..., return_lumen=True)``. It is the clean single-wall 1682 surface the junction geodesics run along. Only used when ``membrane_mask`` is also 1683 supplied (when the membrane is built here, the matching lumen is derived automatically); 1684 when neither is available the geodesic mesh falls back to ``mito & ~membrane``, which is 1685 contaminated by the membrane's border-gap suppression near clipped volume faces. 1686 In ``junction_mode="skeleton"`` the lumen is **also** the junction reference surface (see 1687 :func:`_inner_surface_distance`), so supplying it changes ``crista_junction_count`` and 1688 ``mean_junction_extension_nm`` — it is no longer a display-only input. The 1689 ``mito & ~membrane`` fallback is deliberately *not* used for that purpose. 1690 junction_mode: Which junction detector fills ``crista_junction_count`` — ``"overlap"`` 1691 (default) counts connected components of the direct crista-membrane intersection, while 1692 ``"skeleton"`` counts crista regions reaching within ``max_extension_nm`` of the inner 1693 boundary membrane near a crista terminus. See :func:`detect_junctions`. ``"skeleton"`` 1694 additionally fills ``mean_junction_extension_nm`` and requires 3D input. Every other 1695 column is identical between the two modes. 1696 max_extension_nm: How far (nm) a crista may fall short of the inner boundary membrane surface 1697 and still count (``junction_mode="skeleton"`` only). Defaults to 1698 ``membrane_thickness_nm`` when None, following ``border_gap_nm``. 1699 terminus_nm: A near-membrane crista region counts as a junction only if it lies within this 1700 distance (nm) of a crista terminus (``junction_mode="skeleton"`` only), which rejects a 1701 crista running alongside the membrane. Pass ``inf`` to disable the filter. See 1702 :func:`detect_junctions_skeleton`. 1703 min_junction_volume_nm3: Smallest junction volume (nm^3) that counts 1704 (``junction_mode="skeleton"`` only). Removes specks only — see the known limitation in 1705 :func:`detect_junctions_skeleton`. Measured on the candidate region, not on the painted 1706 footprint. 1707 footprint_nm: How thick (nm) the painted junction label is around each region's closest 1708 approach to the membrane (``junction_mode="skeleton"`` only). Affects the junction label 1709 array and hence the junction centroids, never the count. None means one voxel diagonal. 1710 min_skeleton_nm: Skeleton components shorter than this (nm) are dropped as specks before the 1711 terminus filter runs (``junction_mode="skeleton"`` only). A crista whose entire skeleton is 1712 shorter than this has no termini and so cannot score a junction. 1713 terminus_merge_nm: Termini within this distance (nm) on the same skeleton component collapse to 1714 one (``junction_mode="skeleton"`` only). See :func:`compute_crista_skeleton`. 1715 return_junction_labels: If True, also return the junction label array these counts were 1716 computed from, so a caller can *display* exactly what the table reports instead of 1717 re-detecting junctions itself and getting a different answer. Each worker writes its own 1718 instance's labels into a view of the output volume, so this costs one int32 volume and no 1719 extra computation. Ids are per-instance (``1..crista_junction_count``) and therefore repeat 1720 between mitochondria; the painted voxels never do, since each lies inside its own mito. 1721 1722 Every ``junction_mode="skeleton"`` tuning parameter above is ``None`` by default, meaning "leave the 1723 detector's own default in force" — the concrete values live in :func:`detect_junctions_skeleton` and 1724 :func:`compute_crista_skeleton` rather than being restated here. 1725 1726 The junction nearest-neighbour distances are geodesics along the eroded-mito surface mesh 1727 (``bioimage_cpp.distance.geodesic_distances_mesh``); for a mito with no usable mesh (empty 1728 membrane / degenerate mesh) those columns are NaN. 1729 1730 Implementation notes: each mito is pre-cropped to its bounding box by basic slicing (views, so 1731 cropping is memory-free). Parallelism is adaptive and single-level (never oversubscribed): with 1732 many mitochondria the work is parallelised *across* them on a ``concurrent.futures`` 1733 ``ThreadPoolExecutor`` — the heavy per-mito stages (structure tensor, EDT, geodesics) are 1734 GIL-releasing C++, so threads scale them — with each worker's inner stages kept single-threaded 1735 (the EDT/geodesic solvers are called with ``number_of_threads=1``); with few mitochondria they run 1736 serially and each mito's junction-distance stage gets all cores. The concurrent worker count is 1737 additionally capped so the combined per-mito working set (tensor components + label crops, 1738 ~40 bytes/voxel of the largest mito) fits in RAM. Results stream in as they complete and are 1739 finally sorted by label for an n_jobs-independent ordering. 1740 1741 Returns: 1742 DataFrame with one row per mito instance: 1743 label | mito_volume_nm3 | crista_volume_nm3 | crista_fraction | 1744 contact_voxel_count | crista_junction_count | contact_volume_nm3 | 1745 mean_junction_extension_nm | 1746 avg_crista_to_membrane_nm | mean_nn_junction_distance_nm | median_nn_junction_distance_nm | 1747 junction_clustering_index | crista_orientation_anisotropy | cristae_surface_area_nm2 | 1748 mito_surface_area_nm2 | crista_to_mito_surface_ratio | imm_surface_area_nm2 | 1749 imm_surface_per_mito_volume | imm_surface_per_crista_volume | avg_thickness_nm. 1750 cristae_surface_area_nm2 is the crista surface area; crista_to_mito_surface_ratio is 1751 crista surface / mitochondrial outer-membrane surface (can exceed 1 for folded cristae). 1752 imm_surface_area_nm2 is the inner mitochondrial membrane area — the inner boundary membrane 1753 (the lumen surface) plus cristae_surface_area_nm2 — and the imm_surface_per_*_volume columns 1754 divide it by mito_volume_nm3 and crista_volume_nm3 respectively (nm^-1, "cristae surface 1755 density"). The inner-boundary-membrane area alone is 1756 imm_surface_area_nm2 - cristae_surface_area_nm2. 1757 When ``return_junction_labels`` is True the return value is instead 1758 ``(DataFrame, junction_labels)``, the second an int32 array of the input shape. 1759 1760 The *_nn_junction_distance_nm columns are geodesic nearest-neighbour distances between 1761 crista-membrane junctions along the membrane; junction_clustering_index is a Clark-Evans 1762 index (< 1 clustered, ~ 1 random, > 1 dispersed). ``mean_junction_extension_nm`` is the mean 1763 gap each skeleton end had to bridge to reach the membrane and is NaN unless 1764 ``junction_mode="skeleton"``. ``crista_orientation_anisotropy`` is 1765 computed at full resolution for ``method="exact"``, on a downsampled crop (relative-only, 1766 not comparable) for ``method="fast"``, and left NaN for ``method="skip"``. 1767 """ 1768 if method not in ("fast", "exact", "skip"): 1769 raise ValueError(f"method must be 'fast', 'exact', or 'skip', got {method!r}") 1770 if junction_mode not in ("overlap", "skeleton"): 1771 raise ValueError(f"junction_mode must be 'overlap' or 'skeleton', got {junction_mode!r}") 1772 if max_extension_nm is None: 1773 max_extension_nm = membrane_thickness_nm 1774 if membrane_mask is None: 1775 membrane_mask, lumen_mask = approximate_membrane( 1776 mito_segmentation, voxel_size, membrane_thickness_nm, border_gap_nm, 1777 n_jobs=n_jobs, membrane_mode=membrane_mode, return_lumen=True, 1778 ) 1779 1780 ndim = mito_segmentation.ndim 1781 sampling = _to_sampling(voxel_size, ndim) 1782 voxel_vol = float(np.prod(sampling)) 1783 crista_binary = crista_mask.astype(bool) 1784 vol_shape = mito_segmentation.shape 1785 border_radius = _gap_radius(voxel_size, membrane_thickness_nm, border_gap_nm, ndim) 1786 1787 junction_labels = np.zeros(vol_shape, dtype=np.int32) if return_junction_labels else None 1788 1789 tasks = [] 1790 for prop in regionprops(mito_segmentation): 1791 bbox = prop.bbox 1792 slices = tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim)) 1793 lumen_crop = None if lumen_mask is None else lumen_mask[slices] 1794 tasks.append(( 1795 int(prop.label), bbox, 1796 mito_segmentation[slices], crista_binary[slices], membrane_mask[slices], 1797 lumen_crop, 1798 # a basic-slicing view, so each worker writes its junctions straight into the output 1799 None if junction_labels is None else junction_labels[slices], 1800 )) 1801 total = len(tasks) 1802 1803 def _run(task, inner_n_jobs): 1804 label, bbox, mito_crop, crista_crop, membrane_crop, lumen_crop, junction_out = task 1805 return _single_mito_row( 1806 label, bbox, mito_crop, crista_crop, membrane_crop, 1807 voxel_size, sampling, voxel_vol, vol_shape, border_radius, 1808 method=method, inner_n_jobs=inner_n_jobs, lumen_crop=lumen_crop, 1809 junction_mode=junction_mode, max_extension_nm=max_extension_nm, 1810 terminus_nm=terminus_nm, min_junction_volume_nm3=min_junction_volume_nm3, 1811 footprint_nm=footprint_nm, 1812 min_skeleton_nm=min_skeleton_nm, terminus_merge_nm=terminus_merge_nm, 1813 junction_out=junction_out, 1814 ) 1815 1816 n_workers = os.cpu_count() if n_jobs == -1 else max(1, n_jobs) 1817 across = n_workers > 1 and total >= n_workers 1818 1819 rows = [] 1820 1821 def _consume(results): 1822 for i, row in enumerate( 1823 tqdm(results, total=total, desc="Cristae analysis", disable=not verbose), start=1 1824 ): 1825 rows.append(row) 1826 if progress_callback is not None: 1827 progress_callback(i, total) 1828 1829 if across: 1830 max_voxels = max(int(task[2].size) for task in tasks) 1831 across_workers = _bounded_workers(n_jobs, per_worker_bytes=max_voxels * 40) 1832 # Parallelise across mitochondria with a thread pool (the heavy per-mito stages — structure 1833 # tensor, EDT, geodesics — are GIL-releasing C++). Each worker's inner stages run 1834 # single-threaded (``_run(task, 1)`` passes ``number_of_threads=1`` down to the EDT/geodesic 1835 # solvers) so the across-mito threads do not oversubscribe the cores. Results stream in as they 1836 # complete (``as_completed``) to drive the progress bar; rows are label-sorted below. 1837 with futures.ThreadPoolExecutor(across_workers) as tp: 1838 submitted = [tp.submit(_run, task, 1) for task in tasks] 1839 _consume(future.result() for future in futures.as_completed(submitted)) 1840 else: 1841 _consume(_run(task, n_workers) for task in tasks) 1842 1843 rows.sort(key=lambda row: row["mito_label_id"]) 1844 table = pd.DataFrame(rows) 1845 return (table, junction_labels) if return_junction_labels else table
Compute all crista metrics organised by mitochondrial instance.
Arguments:
- crista_mask: Binary crista segmentation (global volume).
- mito_segmentation: Instance label array (background = 0).
- voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
- membrane_mask: Precomputed membrane mask; recomputed if None.
- membrane_thickness_nm: Membrane shell thickness used if membrane_mask is None.
- border_gap_nm: Border suppression distance passed to approximate_membrane; defaults to membrane_thickness_nm when None.
- method: How the crista orientation anisotropy is computed — this is the ONLY metric that
differs between modes; surface areas (marching cubes), junction distances (geodesic
along the membrane) and thickness/proximity (EDT) are computed identically for all of
them.
"skip"(default) does not compute orientation at all (crista_orientation_anisotropyis NaN) and is the fastest — use it when only the other metrics are needed."fast"computes the anisotropy on a 2× downsampled crista crop (~8× cheaper — the structure tensor is by far the dominant cost); the resulting value is a relative indicator that preserves the ordering between mitochondria but is systematically different in magnitude from the full-resolution value and is NOT comparable tomethod="exact"."exact"computes the anisotropy from the full-resolution structure tensor (use it when the magnitude must be precise). - n_jobs: Number of workers for processing mitochondria in parallel (they are
independent). 1 (default) runs serially; other values use a
concurrent.futuresthread pool (-1 = all cores). Results are identical regardless of n_jobs. - verbose: If True, show a terminal tqdm progress bar over mitochondria.
- progress_callback: Optional callable invoked once per completed mitochondrion with (completed_count, total_count) — e.g. to drive a napari progress bar. It is always called from the calling thread (the futures are consumed here as they complete), so GUI updates from it need no cross-thread marshaling.
- membrane_mode: How the membrane shell is built when
membrane_maskis None —"slice_2d"(default, per-Z-slice 2D erosion, z-parallel) or"shell_3d"(connected 3D shell). Seeapproximate_membrane(). - lumen_mask: Optional eroded-mito interior matching
membrane_mask, i.e. the second return value ofapproximate_membrane(..., return_lumen=True). It is the clean single-wall surface the junction geodesics run along. Only used whenmembrane_maskis also supplied (when the membrane is built here, the matching lumen is derived automatically); when neither is available the geodesic mesh falls back tomito & ~membrane, which is contaminated by the membrane's border-gap suppression near clipped volume faces. Injunction_mode="skeleton"the lumen is also the junction reference surface (see_inner_surface_distance()), so supplying it changescrista_junction_countandmean_junction_extension_nm— it is no longer a display-only input. Themito & ~membranefallback is deliberately not used for that purpose. - junction_mode: Which junction detector fills
crista_junction_count—"overlap"(default) counts connected components of the direct crista-membrane intersection, while"skeleton"counts crista regions reaching withinmax_extension_nmof the inner boundary membrane near a crista terminus. Seedetect_junctions()."skeleton"additionally fillsmean_junction_extension_nmand requires 3D input. Every other column is identical between the two modes. - max_extension_nm: How far (nm) a crista may fall short of the inner boundary membrane surface
and still count (
junction_mode="skeleton"only). Defaults tomembrane_thickness_nmwhen None, followingborder_gap_nm. - terminus_nm: A near-membrane crista region counts as a junction only if it lies within this
distance (nm) of a crista terminus (
junction_mode="skeleton"only), which rejects a crista running alongside the membrane. Passinfto disable the filter. Seedetect_junctions_skeleton(). - min_junction_volume_nm3: Smallest junction volume (nm^3) that counts
(
junction_mode="skeleton"only). Removes specks only — see the known limitation indetect_junctions_skeleton(). Measured on the candidate region, not on the painted footprint. - footprint_nm: How thick (nm) the painted junction label is around each region's closest
approach to the membrane (
junction_mode="skeleton"only). Affects the junction label array and hence the junction centroids, never the count. None means one voxel diagonal. - min_skeleton_nm: Skeleton components shorter than this (nm) are dropped as specks before the
terminus filter runs (
junction_mode="skeleton"only). A crista whose entire skeleton is shorter than this has no termini and so cannot score a junction. - terminus_merge_nm: Termini within this distance (nm) on the same skeleton component collapse to
one (
junction_mode="skeleton"only). Seecompute_crista_skeleton(). - return_junction_labels: If True, also return the junction label array these counts were
computed from, so a caller can display exactly what the table reports instead of
re-detecting junctions itself and getting a different answer. Each worker writes its own
instance's labels into a view of the output volume, so this costs one int32 volume and no
extra computation. Ids are per-instance (
1..crista_junction_count) and therefore repeat between mitochondria; the painted voxels never do, since each lies inside its own mito.
Every junction_mode="skeleton" tuning parameter above is None by default, meaning "leave the
detector's own default in force" — the concrete values live in detect_junctions_skeleton() and
compute_crista_skeleton() rather than being restated here.
The junction nearest-neighbour distances are geodesics along the eroded-mito surface mesh
(bioimage_cpp.distance.geodesic_distances_mesh); for a mito with no usable mesh (empty
membrane / degenerate mesh) those columns are NaN.
Implementation notes: each mito is pre-cropped to its bounding box by basic slicing (views, so
cropping is memory-free). Parallelism is adaptive and single-level (never oversubscribed): with
many mitochondria the work is parallelised across them on a concurrent.futures
ThreadPoolExecutor — the heavy per-mito stages (structure tensor, EDT, geodesics) are
GIL-releasing C++, so threads scale them — with each worker's inner stages kept single-threaded
(the EDT/geodesic solvers are called with number_of_threads=1); with few mitochondria they run
serially and each mito's junction-distance stage gets all cores. The concurrent worker count is
additionally capped so the combined per-mito working set (tensor components + label crops,
~40 bytes/voxel of the largest mito) fits in RAM. Results stream in as they complete and are
finally sorted by label for an n_jobs-independent ordering.
Returns:
DataFrame with one row per mito instance: label | mito_volume_nm3 | crista_volume_nm3 | crista_fraction | contact_voxel_count | crista_junction_count | contact_volume_nm3 | mean_junction_extension_nm | avg_crista_to_membrane_nm | mean_nn_junction_distance_nm | median_nn_junction_distance_nm | junction_clustering_index | crista_orientation_anisotropy | cristae_surface_area_nm2 | mito_surface_area_nm2 | crista_to_mito_surface_ratio | imm_surface_area_nm2 | imm_surface_per_mito_volume | imm_surface_per_crista_volume | avg_thickness_nm. cristae_surface_area_nm2 is the crista surface area; crista_to_mito_surface_ratio is crista surface / mitochondrial outer-membrane surface (can exceed 1 for folded cristae). imm_surface_area_nm2 is the inner mitochondrial membrane area — the inner boundary membrane (the lumen surface) plus cristae_surface_area_nm2 — and the imm_surface_per_*_volume columns divide it by mito_volume_nm3 and crista_volume_nm3 respectively (nm^-1, "cristae surface density"). The inner-boundary-membrane area alone is imm_surface_area_nm2 - cristae_surface_area_nm2. When
return_junction_labelsis True the return value is instead(DataFrame, junction_labels), the second an int32 array of the input shape.The *_nn_junction_distance_nm columns are geodesic nearest-neighbour distances between crista-membrane junctions along the membrane; junction_clustering_index is a Clark-Evans index (< 1 clustered, ~ 1 random, > 1 dispersed).
mean_junction_extension_nmis the mean gap each skeleton end had to bridge to reach the membrane and is NaN unlessjunction_mode="skeleton".crista_orientation_anisotropyis computed at full resolution formethod="exact", on a downsampled crop (relative-only, not comparable) formethod="fast", and left NaN formethod="skip".