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_erosion, center_of_mass
   9from scipy.ndimage import label as ndimage_label
  10from skimage.measure import mesh_surface_area, regionprops
  11from skimage.morphology import 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
  17
  18
  19# ---------------------------------------------------------------------------
  20# Internal helpers
  21# ---------------------------------------------------------------------------
  22
  23def _to_sampling(voxel_size: Union[float, Dict[str, float]], ndim: int) -> np.ndarray:
  24    axes = ("z", "y", "x") if ndim == 3 else ("y", "x")
  25    if isinstance(voxel_size, dict):
  26        return np.array([voxel_size[ax] for ax in axes[:ndim]], dtype=float)
  27    return np.full(ndim, float(voxel_size))
  28
  29
  30def _voxel_radius(thickness_nm: float, voxel_size: Union[float, Dict[str, float]], ndim: int) -> int:
  31    return max(1, int(round(thickness_nm / float(np.mean(_to_sampling(voxel_size, ndim))))))
  32
  33
  34def _voxel_radius_xy(thickness_nm: float, voxel_size: Union[float, Dict[str, float]]) -> int:
  35    """Membrane radius in XY pixels — uses only the Y and X voxel sizes."""
  36    if isinstance(voxel_size, dict):
  37        xy_nm = (voxel_size["y"] + voxel_size["x"]) / 2.0
  38    else:
  39        xy_nm = float(voxel_size)
  40    return max(1, int(round(thickness_nm / xy_nm)))
  41
  42
  43def _gap_radius(
  44    voxel_size: Union[float, Dict[str, float]],
  45    membrane_thickness_nm: float,
  46    border_gap_nm: Optional[float],
  47    ndim: int,
  48) -> int:
  49    """Border-zone / mesh-trim radius in voxels; ``border_gap_nm`` defaults to ``membrane_thickness_nm``."""
  50    gap_nm = border_gap_nm if border_gap_nm is not None else membrane_thickness_nm
  51    return _voxel_radius(gap_nm, voxel_size, ndim)
  52
  53
  54def _border_zone(shape: tuple, radius: int) -> np.ndarray:
  55    """Boolean mask that is True within `radius` voxels of any face of the volume."""
  56    mask = np.zeros(shape, dtype=bool)
  57    for ax in range(len(shape)):
  58        idx_lo = [slice(None)] * len(shape)
  59        idx_hi = [slice(None)] * len(shape)
  60        idx_lo[ax] = slice(0, radius)
  61        idx_hi[ax] = slice(shape[ax] - radius, None)
  62        mask[tuple(idx_lo)] = True
  63        mask[tuple(idx_hi)] = True
  64    return mask
  65
  66
  67def _surface_mesh(
  68    mask: np.ndarray,
  69    sampling: np.ndarray,
  70    closed_faces: Optional[np.ndarray] = None,
  71) -> Optional[Tuple[np.ndarray, np.ndarray]]:
  72    """Triangle-mesh surface of a binary mask via marching cubes.
  73
  74    Each side of the array is padded by one background voxel before meshing so that objects touching
  75    the array edge yield a closed surface there. ``marching_cubes`` only emits triangles for cube
  76    cells that exist inside the array, so leaving a side *unpadded* omits that boundary face — an open
  77    mesh at that face — while interior surfaces still close. ``closed_faces`` chooses this per side:
  78    pad+close a face where there is genuine background beyond it, leave open a face where the object
  79    is clipped by the volume boundary (membrane presence unknown there).
  80
  81    Vertices are returned in the mask's own **unpadded** index frame (nm): a mask voxel at array index
  82    ``(z, y, x)`` maps to physical coordinates ``index * sampling`` regardless of the padding, so
  83    callers snapping voxel coordinates onto the mesh use ``index * sampling`` with no offset.
  84
  85    Args:
  86        mask: Binary segmentation.
  87        sampling: Voxel size per axis (nm), in array (z, y, x) order.
  88        closed_faces: Optional ``(ndim, 2)`` boolean array; ``closed_faces[a, s]`` True closes
  89            (pads) side ``s`` (0 = low, 1 = high) of axis ``a``, False leaves it open. Defaults to
  90            all True (fully closed watertight surface).
  91
  92    Returns:
  93        (vertices, faces) with vertices in nm (unpadded mask frame), or None if the mask is empty.
  94    """
  95    binary = mask.astype(bool)
  96    if not binary.any():
  97        return None
  98    if closed_faces is None:
  99        closed_faces = np.ones((binary.ndim, 2), dtype=bool)
 100    else:
 101        closed_faces = np.asarray(closed_faces, dtype=bool)
 102    pad_width = [(int(closed_faces[a, 0]), int(closed_faces[a, 1])) for a in range(binary.ndim)]
 103    padded = np.pad(binary.astype(np.float32), pad_width)
 104    verts, faces, _, _ = marching_cubes(padded, level=0.5, spacing=tuple(float(s) for s in sampling))
 105    pad_before = np.array([pw[0] for pw in pad_width], dtype=float)
 106    verts = verts - pad_before * np.asarray(sampling, dtype=float)
 107    return verts, faces
 108
 109
 110def _surface_area(
 111    mask: np.ndarray,
 112    sampling: np.ndarray,
 113    closed_faces: Optional[np.ndarray] = None,
 114) -> float:
 115    """Surface area (nm^2) of a binary mask via marching cubes.
 116
 117    Args:
 118        mask: Binary segmentation.
 119        sampling: Voxel size per axis (nm), in array (z, y, x) order.
 120        closed_faces: Optional per-side padding spec forwarded to :func:`_surface_mesh` — leave a
 121            clipped volume-boundary face open so its fabricated cap is not counted as surface area.
 122            Defaults to a fully closed (watertight) surface.
 123
 124    Returns:
 125        Surface area in nm^2, or NaN if the mask is empty.
 126    """
 127    mesh = _surface_mesh(mask, sampling, closed_faces=closed_faces)
 128    if mesh is None:
 129        return np.nan
 130    return float(mesh_surface_area(*mesh))
 131
 132
 133def _open_trimmed_mesh(
 134    mask: np.ndarray,
 135    sampling: np.ndarray,
 136    gap_radius: int,
 137    boundary: np.ndarray,
 138) -> Optional[Tuple[np.ndarray, np.ndarray]]:
 139    """Surface mesh of ``mask`` trimmed to the certain region and left OPEN at volume-boundary faces.
 140
 141    Near a clipped volume face the segmentation is cut off and membrane presence is unknown, so the
 142    surface must neither flare into that region nor be capped there — a cap is a fabricated flat disk
 143    that lets geodesics shortcut straight across instead of wrapping around the tube wall. For each
 144    ``(axis, side)`` flagged True in ``boundary`` (an ``(ndim, 2)`` bool of volume-boundary faces),
 145    ``gap_radius`` voxels are cropped off that face so the trim plane becomes an array boundary, and
 146    that face is then left open (marching cubes omits it). Interior faces stay closed.
 147
 148    Args:
 149        mask: Binary segmentation.
 150        sampling: Voxel size per axis (nm), in array (z, y, x) order.
 151        gap_radius: Border-zone width in voxels cropped off each flagged face (matches the membrane's
 152            border-gap trim).
 153        boundary: ``(ndim, 2)`` bool; True where the face is at the volume boundary (crop + open).
 154
 155    Returns:
 156        (vertices, faces) with vertices in nm in the original (uncropped) ``mask`` index frame, or None
 157        if the trimmed mask is empty.
 158    """
 159    ndim = mask.ndim
 160    boundary = np.asarray(boundary, dtype=bool)
 161    lo = [gap_radius if boundary[a, 0] else 0 for a in range(ndim)]
 162    hi = [mask.shape[a] - (gap_radius if boundary[a, 1] else 0) for a in range(ndim)]
 163    mesh = _surface_mesh(
 164        mask[tuple(slice(lo[a], hi[a]) for a in range(ndim))], sampling, closed_faces=~boundary
 165    )
 166    if mesh is None:
 167        return None
 168    verts, faces = mesh
 169    verts = verts + np.array(lo, dtype=float) * np.asarray(sampling, dtype=float)
 170    return verts, faces
 171
 172
 173def _medial_axis_thickness_nm(mask: np.ndarray, sampling: np.ndarray) -> float:
 174    """Local thickness (nm) of a mask via the distance transform, with no mesh generation.
 175
 176    The medial axis is approximated by the local maxima of the interior EDT; the thickness is
 177    ``2 × mean(EDT)`` there (the EDT at the medial axis is the half-thickness). This is the same
 178    estimator used by :func:`compute_crista_morphology`'s ``medial_axis`` branch, factored out so
 179    the distance-based (``method="fast"``) surface-area estimates can reuse it.
 180
 181    Args:
 182        mask: Binary segmentation.
 183        sampling: Voxel size per axis (nm), in array (z, y, x) order.
 184
 185    Returns:
 186        Mean local thickness in nm, or NaN if the mask is empty.
 187    """
 188    binary = mask.astype(bool)
 189    if not binary.any():
 190        return np.nan
 191    dist = distance_transform(binary, sampling=tuple(float(s) for s in sampling), number_of_threads=1)
 192    ridges = local_maxima(dist) & binary
 193    ridge_dists = dist[ridges]
 194    return float(2.0 * np.mean(ridge_dists)) if ridge_dists.size > 0 else np.nan
 195
 196
 197def _available_memory_bytes() -> int:
 198    """Best-effort available RAM in bytes (used to keep parallel working sets from OOMing)."""
 199    try:
 200        import psutil
 201        return int(psutil.virtual_memory().available)
 202    except Exception:
 203        try:
 204            return int(os.sysconf("SC_AVPHYS_PAGES") * os.sysconf("SC_PAGE_SIZE"))
 205        except Exception:
 206            return 4 * 1024 ** 3
 207
 208
 209def _bounded_workers(n_jobs: int, per_worker_bytes: int, fraction: float = 0.5) -> int:
 210    """Resolve n_jobs to a worker count whose combined working set fits in memory.
 211
 212    n_jobs: 1 = serial, -1 = all cores, else that many. The result is additionally capped so
 213    ``workers * per_worker_bytes <= fraction * available_RAM`` (at least 1).
 214    """
 215    workers = os.cpu_count() if n_jobs == -1 else max(1, int(n_jobs))
 216    if per_worker_bytes > 0:
 217        budget = int(_available_memory_bytes() * fraction)
 218        workers = min(workers, max(1, budget // int(per_worker_bytes)))
 219    return int(max(1, workers))
 220
 221
 222# ---------------------------------------------------------------------------
 223# Membrane approximation
 224# ---------------------------------------------------------------------------
 225
 226def approximate_membrane(
 227    mito_segmentation: np.ndarray,
 228    voxel_size: Union[float, Dict[str, float]],
 229    membrane_thickness_nm: float = 8.0,
 230    border_gap_nm: Optional[float] = None,
 231    n_jobs: int = 1,
 232    membrane_mode: str = "slice_2d",
 233    return_lumen: bool = False,
 234) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:
 235    """Approximate the mitochondrial membrane as the outer shell of the segmentation.
 236
 237    Two shell constructions are available via ``membrane_mode``:
 238
 239    - ``"slice_2d"`` (default): erode each Z-slice **independently in 2D** by an XY disk of radius
 240      ``round(thickness / xy_voxel)`` and keep ``slice & ~eroded``. A mitochondrion that changes shape
 241      rapidly along Z does not bleed into neighbouring slices, and the per-slice erosions are
 242      parallelised over Z (``n_jobs``). A separable Z-only erosion (radius ``round(thickness /
 243      z_voxel)``, no XY coupling) then adds **Z-caps** where a mito column truly ends in Z; ends
 244      clipped by a volume Z-face are left uncapped (``border_value=1`` + the ``border_gap`` trim). The
 245      XY shell can still fragment across slices. (2D inputs get a single 2D erosion.)
 246    - ``"shell_3d"``: a full 3D morphological erosion, ``mito & ~erode3d(mito, k)`` with
 247      ``k = round(thickness / mean_voxel)`` iterations of a 3×3×3 structuring element, per instance on
 248      its padded bounding box. A single **connected** shell including the Z-caps (no per-slice
 249      fragmentation), at a higher cost; thickness acts in all axes.
 250
 251    The eroded interior is the "lumen"; its surface is the single-wall mesh used by the geodesic
 252    backend and the display, so the mesh follows the chosen mode.
 253
 254    Membrane voxels within ``border_gap_nm`` of any volume face are removed so clipped mito edges are
 255    not treated as membrane. The lumen is NOT trimmed here — the mesh is trimmed to the certain region
 256    (and left open there) at mesh time by :func:`_open_trimmed_mesh`, which requires the untrimmed
 257    interior to produce an open cut rather than a fabricated cap.
 258
 259    Implementation notes: ``"slice_2d"`` erodes each Z-slice on the mito XY bbox with a
 260    ``membrane_radius`` margin, so the cropped ``border_value=1`` erosion matches eroding the full
 261    slice (empty slices are skipped), then adds Z-caps via a separable Z-only line erosion (which
 262    inspects only the same column, so no XY-shape bleed); ``border_value=1`` leaves ends clipped by a
 263    volume Z-face uncapped, and the ``border_gap`` removal clears anything near a face, so only true
 264    ends are capped. ``"shell_3d"`` erodes the *merged* binary in each instance's padded bbox so
 265    instances that share a boundary are handled together.
 266
 267    Args:
 268        mito_segmentation: Instance label array (background = 0).
 269        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 270        membrane_thickness_nm: Thickness of the membrane shell in nm.
 271        border_gap_nm: Distance from each volume face within which membrane voxels are
 272            suppressed. Defaults to membrane_thickness_nm when None.
 273        n_jobs: Workers for the per-Z-slice erosion (``"slice_2d"`` only): 1 = serial, -1 = all cores.
 274            Results are identical regardless of n_jobs.
 275        membrane_mode: ``"slice_2d"`` (default, per-slice 2D, z-parallel) or ``"shell_3d"``
 276            (connected 3D shell).
 277        return_lumen: If True, also return the eroded-mito interior ("lumen") mask — the single-wall
 278            surface source for the geodesic/display mesh. It is the plain eroded interior (not
 279            ``mito & ~membrane``, which would re-include the outer shell where the membrane is
 280            border-trimmed); trimming to the certain region happens at mesh time.
 281
 282    Returns:
 283        membrane_mask: Binary mask of the mitochondrial membrane (outer shell), with border-adjacent
 284            voxels zeroed out. If ``return_lumen`` is True, returns ``(membrane_mask, lumen_mask)``
 285            where ``lumen_mask`` is the (untrimmed) eroded interior described above.
 286    """
 287    if membrane_mode not in ("slice_2d", "shell_3d"):
 288        raise ValueError(f"membrane_mode must be 'slice_2d' or 'shell_3d', got {membrane_mode!r}")
 289    ndim = mito_segmentation.ndim
 290    mito_binary = mito_segmentation > 0
 291
 292    # NOTE (possible simplification, deferred): the "shell_3d" branch below builds the shell with an
 293    # iterated 3x3x3 erosion per instance. It could likely be a single anisotropic distance transform
 294    # instead — membrane = mito & (distance_transform(mito, sampling) <= thickness), lumen = the rest —
 295    # which is simpler and handles anisotropy directly. This would NOT replace "slice_2d": a distance
 296    # transform couples all axes, so it cannot reproduce slice_2d's per-Z-slice-independent erosion,
 297    # whose whole purpose is to stop the shell bleeding across slices in XY. Worth investigating.
 298    if membrane_mode == "shell_3d":
 299        sampling = _to_sampling(voxel_size, ndim)
 300        k = max(1, int(round(float(membrane_thickness_nm) / float(np.mean(sampling)))))
 301        struct = np.ones((3,) * ndim, dtype=bool)
 302        membrane_mask = np.zeros(mito_segmentation.shape, dtype=bool)
 303        lumen_mask = np.zeros(mito_segmentation.shape, dtype=bool)
 304        for prop in regionprops(mito_segmentation):
 305            bbox = prop.bbox
 306            sl = tuple(
 307                slice(max(0, bbox[i] - k), min(mito_segmentation.shape[i], bbox[i + ndim] + k))
 308                for i in range(ndim)
 309            )
 310            sub = mito_binary[sl]
 311            eroded = binary_erosion(sub, structure=struct, iterations=k, border_value=1)
 312            cur = mito_segmentation[sl] == prop.label
 313            membrane_mask[sl] |= cur & ~eroded
 314            lumen_mask[sl] |= cur & eroded
 315    elif ndim == 3:
 316        membrane_radius = _voxel_radius_xy(membrane_thickness_nm, voxel_size)
 317        struct = disk(membrane_radius)
 318        membrane_mask = np.zeros_like(mito_binary)
 319        lumen_mask = np.zeros_like(mito_binary)
 320        coords = np.argwhere(mito_binary)
 321        if coords.size:
 322            zmin, ymin, xmin = coords.min(axis=0)
 323            zmax, ymax, xmax = coords.max(axis=0) + 1
 324            m = membrane_radius
 325            y0, y1 = max(0, ymin - m), min(mito_binary.shape[1], ymax + m)
 326            x0, x1 = max(0, xmin - m), min(mito_binary.shape[2], xmax + m)
 327
 328            def _erode_slice(z):
 329                sl = mito_binary[z, y0:y1, x0:x1]
 330                if not sl.any():
 331                    return z, None
 332                eroded = binary_erosion(sl, structure=struct, border_value=1)
 333                return z, (sl & ~eroded, eroded)
 334
 335            z_range = range(int(zmin), int(zmax))
 336            if n_jobs == 1:
 337                results = [_erode_slice(z) for z in z_range]
 338            else:
 339                n_workers = mp.cpu_count() if n_jobs == -1 else n_jobs
 340                with futures.ThreadPoolExecutor(n_workers) as tp:
 341                    results = list(tp.map(_erode_slice, z_range))
 342            for z, res in results:
 343                if res is not None:
 344                    mem_sl, lum_sl = res
 345                    membrane_mask[z, y0:y1, x0:x1] = mem_sl
 346                    lumen_mask[z, y0:y1, x0:x1] = lum_sl
 347
 348            k_z = max(1, int(round(float(membrane_thickness_nm) / float(_to_sampling(voxel_size, ndim)[0]))))
 349            z0m, z1m = max(0, int(zmin) - k_z), min(mito_binary.shape[0], int(zmax) + k_z)
 350            sub = mito_binary[z0m:z1m, y0:y1, x0:x1]
 351            z_eroded = binary_erosion(sub, structure=np.ones((2 * k_z + 1, 1, 1), dtype=bool), border_value=1)
 352            membrane_mask[z0m:z1m, y0:y1, x0:x1] |= sub & ~z_eroded
 353            lumen_mask[z0m:z1m, y0:y1, x0:x1] &= z_eroded
 354    else:
 355        membrane_radius = _voxel_radius(membrane_thickness_nm, voxel_size, ndim)
 356        eroded = binary_erosion(mito_binary, structure=disk(membrane_radius), border_value=1)
 357        membrane_mask = mito_binary & ~eroded
 358        lumen_mask = mito_binary & eroded
 359
 360    gap_radius = _gap_radius(voxel_size, membrane_thickness_nm, border_gap_nm, ndim)
 361    membrane_mask &= ~_border_zone(mito_segmentation.shape, gap_radius)
 362    if return_lumen:
 363        return membrane_mask.astype(bool), lumen_mask.astype(bool)
 364    return membrane_mask.astype(bool)
 365
 366
 367# ---------------------------------------------------------------------------
 368# Orientation
 369# ---------------------------------------------------------------------------
 370
 371def compute_crista_orientation(
 372    crista_mask: np.ndarray,
 373    voxel_size: Union[float, Dict[str, float]],
 374    neighborhood_size_nm: float = 30.0,
 375) -> np.ndarray:
 376    """Compute the per-voxel crista orientation anisotropy via the structure tensor.
 377
 378    Uses ``bioimage_cpp.filters.structure_tensor_eigenvalues`` (a fast C++ routine). Only the
 379    anisotropy is produced (the principal directions / eigenvectors are not computed).
 380
 381    The structure tensor's outer/integration sigma is ``neighborhood_size_nm`` per axis (in voxels);
 382    the inner (derivative) sigma must be > 0, so a minimal 1-voxel scale is used. Eigenvalues are
 383    non-negative in theory, but the solver emits tiny negatives for near-rank-deficient tensors
 384    (degenerate sheets/tubes), so they are clamped to 0 and the ratio is taken as
 385    ``max/min`` over the trailing axis — order-agnostic and sign-safe, so a tiny negative minor
 386    eigenvalue cannot flip the denominator and blow the ratio up.
 387
 388    Args:
 389        crista_mask: Binary crista segmentation.
 390        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 391        neighborhood_size_nm: Gaussian integration radius in nm for tensor averaging (the
 392            structure tensor's outer/integration scale).
 393
 394    Returns:
 395        anisotropy: (...) — λ_max / (λ_min + ε) per voxel. High values indicate a strongly
 396            directional crista (e.g. parallel lamellae); low values indicate isotropic or
 397            tubular/disordered morphology. Magnitude only (rotation-invariant).
 398    """
 399    ndim = crista_mask.ndim
 400    sampling = _to_sampling(voxel_size, ndim)
 401
 402    outer_sigma = [float(s) for s in (neighborhood_size_nm / sampling)]
 403    inner_sigma = 1.0
 404    evals = structure_tensor_eigenvalues(crista_mask.astype(np.float32), inner_sigma, outer_sigma)
 405    evals = np.clip(evals, 0.0, None)
 406    anisotropy = evals.max(axis=-1) / (evals.min(axis=-1) + 1e-10)
 407    return anisotropy.astype(np.float32)
 408
 409
 410def _scale_voxel_size(
 411    voxel_size: Union[float, Dict[str, float]], factor: float
 412) -> Union[float, Dict[str, float]]:
 413    """Multiply a voxel size (scalar or z/y/x dict) by ``factor``, preserving its type."""
 414    if isinstance(voxel_size, dict):
 415        return {ax: voxel_size[ax] * factor for ax in voxel_size}
 416    return float(voxel_size) * factor
 417
 418
 419def _downsampled_orientation_anisotropy(
 420    crista_mask: np.ndarray,
 421    voxel_size: Union[float, Dict[str, float]],
 422    factor: int = 2,
 423) -> float:
 424    """Mean crista orientation anisotropy computed on a downsampled crop (fast, approximate).
 425
 426    The crista mask is **nearest-neighbour** downsampled by ``factor`` per axis (strided subsampling,
 427    i.e. every ``factor``-th voxel), and :func:`compute_crista_orientation` is run on that coarser grid
 428    at the correspondingly scaled voxel size (so the physical structure-tensor neighbourhood is
 429    unchanged). This is ~``factor**ndim`` times cheaper than the full-resolution structure tensor — the
 430    dominant cost of the analysis. Nearest-neighbour keeps the segmentation binary; block-mean
 431    (``downscale_local_mean``) would blur it into meaningless partial-occupancy values.
 432
 433    Because thin cristae (only a few voxels across) lose structure when downsampled, the returned
 434    anisotropy is a *relative* indicator only: it preserves the ordering between mitochondria but is
 435    not comparable in magnitude to the full-resolution ``method="exact"`` value.
 436
 437    Args:
 438        crista_mask: Binary crista segmentation.
 439        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 440        factor: Integer downsampling factor per axis.
 441
 442    Returns:
 443        Mean anisotropy over the downsampled crista region, or NaN if it vanishes when downsampled.
 444    """
 445    ndim = crista_mask.ndim
 446    sub = crista_mask[(slice(None, None, factor),) * ndim]
 447    region = sub.astype(bool)
 448    if not region.any():
 449        return np.nan
 450    anisotropy = compute_crista_orientation(sub, _scale_voxel_size(voxel_size, factor))
 451    return float(np.mean(anisotropy[region]))
 452
 453
 454# ---------------------------------------------------------------------------
 455# Proximity
 456# ---------------------------------------------------------------------------
 457
 458def compute_crista_proximity(
 459    crista_mask: np.ndarray,
 460    membrane_mask: np.ndarray,
 461    voxel_size: Union[float, Dict[str, float]],
 462    membrane_distance: Optional[np.ndarray] = None,
 463) -> Tuple[np.ndarray, Dict[str, float]]:
 464    """Distance from each crista voxel to the nearest membrane voxel (nm).
 465
 466    Args:
 467        crista_mask: Binary crista segmentation.
 468        membrane_mask: Binary membrane mask (OM or IMM).
 469        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 470        membrane_distance: Optional precomputed per-voxel distance to the nearest membrane
 471            voxel (nm), i.e. ``distance_transform(~membrane_mask, sampling=...)``. When
 472            given, the distance transform is not recomputed (used to avoid redundant work in
 473            :func:`compute_mito_crista_statistics`).
 474
 475    Returns:
 476        distance_map: Per-voxel distance to membrane (nm); zero outside crista.
 477        summary_stats: min_nm, median_nm, max_nm.
 478    """
 479    sampling = _to_sampling(voxel_size, crista_mask.ndim)
 480    if membrane_distance is None:
 481        dist = distance_transform(~membrane_mask.astype(bool), sampling=sampling.tolist(), number_of_threads=1)
 482    else:
 483        dist = membrane_distance
 484    crista_dists = dist[crista_mask.astype(bool)]
 485
 486    if crista_dists.size == 0:
 487        summary: Dict[str, float] = {"min_nm": np.nan, "median_nm": np.nan, "max_nm": np.nan}
 488    else:
 489        summary = {
 490            "min_nm": float(crista_dists.min()),
 491            "median_nm": float(np.median(crista_dists)),
 492            "max_nm": float(crista_dists.max()),
 493        }
 494
 495    distance_map = np.zeros(crista_mask.shape, dtype=np.float32)
 496    distance_map[crista_mask.astype(bool)] = crista_dists
 497    return distance_map, summary
 498
 499
 500# ---------------------------------------------------------------------------
 501# Contact sites
 502# ---------------------------------------------------------------------------
 503
 504def detect_contact_sites(
 505    crista_mask: np.ndarray,
 506    membrane_mask: np.ndarray,
 507    voxel_size: Union[float, Dict[str, float]],
 508) -> Tuple[np.ndarray, Dict[str, float]]:
 509    """Detect crista-membrane contact sites as the direct overlap of the two masks.
 510
 511    Contact = crista voxels that are also membrane voxels (the pure intersection of the crista
 512    mask and the mitochondrial membrane band). No dilation/erosion is applied here, so the
 513    detected junctions correspond exactly to the visible overlap of the two layers; connected
 514    overlaps are grouped into junctions with 26-connectivity in 3D.
 515
 516    Args:
 517        crista_mask: Binary crista segmentation.
 518        membrane_mask: Binary mitochondrial membrane mask (as displayed).
 519        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 520
 521    Returns:
 522        contact_labels: Integer array (same shape as the input) where each connected
 523            junction has a unique ID (0 = background, 1..n = junctions). Contact voxel
 524            coordinates are recoverable via ``np.argwhere(contact_labels > 0)``.
 525        summary: contact_voxel_count, crista_junction_count, contact_volume_nm3.
 526    """
 527    ndim = crista_mask.ndim
 528    sampling = _to_sampling(voxel_size, ndim)
 529    voxel_vol = float(np.prod(sampling))
 530
 531    contact_mask = crista_mask.astype(bool) & membrane_mask.astype(bool)
 532
 533    connectivity_struct = np.ones(ndim * (3,), dtype=bool)
 534    contact_labels, n_regions = ndimage_label(contact_mask, structure=connectivity_struct)
 535    contact_voxel_count = int(np.count_nonzero(contact_labels))
 536
 537    return contact_labels, {
 538        "contact_voxel_count": contact_voxel_count,
 539        "crista_junction_count": int(n_regions),
 540        "contact_volume_nm3": float(contact_voxel_count) * voxel_vol,
 541    }
 542
 543
 544_JUNCTION_DISTANCE_NAN = {
 545    "junction_count": 0,
 546    "mean_nn_junction_distance_nm": np.nan,
 547    "median_nn_junction_distance_nm": np.nan,
 548    "junction_clustering_index": np.nan,
 549}
 550
 551
 552def _junction_matrix_mesh(
 553    centroids: np.ndarray,
 554    sampling: np.ndarray,
 555    vertices: np.ndarray,
 556    faces: np.ndarray,
 557    n_jobs: int = 1,
 558) -> Optional[np.ndarray]:
 559    """Pairwise junction geodesic distances (nm) along a triangle-mesh surface.
 560
 561    Snaps each junction centroid to its nearest mesh vertex and calls
 562    ``bioimage_cpp.distance.geodesic_distances_mesh`` (passing all junction vertices as sources, so it
 563    returns the full pairwise matrix directly). The mesh (see :func:`_surface_mesh`) returns vertices
 564    in the unpadded mask index frame, so voxel centroids map to it as ``index * sampling`` with no
 565    offset. Returns None (junction distances become NaN) when the mesh is empty.
 566
 567    Args:
 568        centroids: (n, ndim) junction centroids in voxel coordinates.
 569        sampling: Voxel size per axis (nm), array (z, y, x) order.
 570        vertices: Mesh vertices (n_vertices, 3) in nm (unpadded mask frame).
 571        faces: Mesh triangle indices (n_faces, 3).
 572        n_jobs: Forwarded to the C++ solver's ``number_of_threads``: -1/0 map to 0 (the solver's
 573            "use hardware_concurrency"), otherwise that many threads.
 574
 575    Returns:
 576        (n, n) geodesic distance matrix in nm (0 diagonal, NaN for disconnected pairs), or None.
 577        Disconnected pairs come back from the solver as ``+inf`` and are converted to NaN.
 578    """
 579    verts = np.ascontiguousarray(vertices, dtype=np.float64)
 580    tris = np.ascontiguousarray(faces, dtype=np.int64)
 581    if verts.shape[0] == 0 or tris.shape[0] == 0:
 582        return None
 583    from scipy.spatial import cKDTree
 584
 585    points = np.asarray(centroids, dtype=float) * sampling
 586    _, vertex_ids = cKDTree(verts).query(points)
 587    vertex_ids = np.atleast_1d(np.asarray(vertex_ids, dtype=np.int64))
 588    n_threads = 0 if n_jobs in (-1, 0) else max(1, int(n_jobs))
 589    dm = np.asarray(
 590        geodesic_distances_mesh(verts, tris, vertex_ids, number_of_threads=n_threads), dtype=float
 591    )
 592    dm[~np.isfinite(dm)] = np.nan
 593    np.fill_diagonal(dm, 0.0)
 594    return dm
 595
 596
 597def compute_junction_distances(
 598    contact_labels: np.ndarray,
 599    membrane_mask: np.ndarray,
 600    voxel_size: Union[float, Dict[str, float]],
 601    surface_area_nm2: Optional[float] = None,
 602    n_jobs: int = 1,
 603    mesh_vertices: Optional[np.ndarray] = None,
 604    mesh_faces: Optional[np.ndarray] = None,
 605) -> Tuple[np.ndarray, Dict[str, float]]:
 606    """Geodesic distances between crista-membrane junctions along the eroded-mito surface mesh.
 607
 608    Each junction (a connected component in ``contact_labels``) is reduced to its centroid, snapped to
 609    the nearest vertex of a triangle mesh, and pairwise surface geodesics are computed with
 610    ``bioimage_cpp.distance.geodesic_distances_mesh``. The mesh is the **eroded-mito (lumen) surface**
 611    passed in as ``mesh_vertices``/``mesh_faces`` by :func:`_single_mito_row` (a clean, single-wall
 612    surface at the membrane's inner edge). If no mesh is supplied — or no usable surface mesh exists
 613    (empty membrane / degenerate mesh) — the junction distances are NaN. (There is no membrane-band
 614    fallback mesh: the metric is defined on the lumen surface, and meshing the thick membrane band
 615    would give a different, capped double-wall surface.)
 616
 617    A Clark-Evans nearest-neighbour index summarises whether the junctions are clustered.
 618
 619    Args:
 620        contact_labels: Integer junction label array (0 = background, 1..n = junctions),
 621            e.g. the first return value of :func:`detect_contact_sites`.
 622        membrane_mask: Binary mitochondrial membrane mask the junctions sit on. Only used for the
 623            empty-membrane early-out (no membrane → NaN); it is not meshed.
 624        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 625        surface_area_nm2: Membrane/mito surface area used as the reference area for the
 626            Clark-Evans expectation. If None or non-positive, the clustering index is NaN.
 627        n_jobs: 1 = serial, -1 = all cores (forwarded to the mesh solver's thread count).
 628        mesh_vertices: Optional (n_vertices, 3) mesh vertices in nm (unpadded mask frame; see
 629            :func:`_surface_mesh`) — the eroded-mito (lumen) surface. If omitted, the junction
 630            distances are NaN.
 631        mesh_faces: Optional (n_faces, 3) triangle indices matching ``mesh_vertices``.
 632
 633    Returns:
 634        distance_matrix: (n, n) geodesic distances in nm between junctions; the diagonal is
 635            0 and unreachable pairs (disconnected fragments) are NaN. Empty for fewer than two
 636            junctions.
 637        summary: junction_count, mean_nn_junction_distance_nm, median_nn_junction_distance_nm,
 638            junction_clustering_index (Clark-Evans R: < 1 clustered, ~ 1 random, > 1 dispersed).
 639
 640    Notes:
 641        The Clark-Evans expected nearest-neighbour distance uses the standard 2-D planar
 642        approximation ``0.5 * sqrt(A / n)`` with ``A = surface_area_nm2``.
 643
 644        Each junction's nearest-neighbour distance is the smallest distance to a *reachable* other
 645        junction: self (diagonal) and unreachable (NaN) pairs are set to +inf before the per-row
 646        minimum, and rows with no reachable neighbour (min stays +inf) are dropped.
 647    """
 648    ndim = contact_labels.ndim
 649    sampling = _to_sampling(voxel_size, ndim)
 650    membrane = membrane_mask.astype(bool)
 651
 652    labels = [lbl for lbl in np.unique(contact_labels) if lbl != 0]
 653    n = len(labels)
 654    if n < 2 or not membrane.any():
 655        summary = dict(_JUNCTION_DISTANCE_NAN)
 656        summary["junction_count"] = n
 657        return np.zeros((n, n), dtype=float), summary
 658
 659    centroids = np.atleast_2d(
 660        np.asarray(center_of_mass(contact_labels > 0, labels=contact_labels, index=labels), dtype=float)
 661    )
 662
 663    if mesh_vertices is not None and mesh_faces is not None and len(mesh_faces) > 0:
 664        distance_matrix = _junction_matrix_mesh(centroids, sampling, mesh_vertices, mesh_faces, n_jobs)
 665    else:
 666        distance_matrix = None
 667
 668    if distance_matrix is None:
 669        summary = dict(_JUNCTION_DISTANCE_NAN)
 670        summary["junction_count"] = n
 671        return np.full((n, n), np.nan, dtype=float), summary
 672
 673    dm = distance_matrix.copy()
 674    np.fill_diagonal(dm, np.inf)
 675    dm[~np.isfinite(dm)] = np.inf
 676    row_min = dm.min(axis=1)
 677    nn_distances = row_min[np.isfinite(row_min)]
 678
 679    mean_nn = float(np.mean(nn_distances)) if nn_distances.size else np.nan
 680    median_nn = float(np.median(nn_distances)) if nn_distances.size else np.nan
 681
 682    clustering_index = np.nan
 683    if surface_area_nm2 is not None and surface_area_nm2 > 0 and np.isfinite(mean_nn):
 684        expected_nn = 0.5 * np.sqrt(float(surface_area_nm2) / n)
 685        if expected_nn > 0:
 686            clustering_index = mean_nn / expected_nn
 687
 688    return distance_matrix, {
 689        "junction_count": n,
 690        "mean_nn_junction_distance_nm": mean_nn,
 691        "median_nn_junction_distance_nm": median_nn,
 692        "junction_clustering_index": clustering_index,
 693    }
 694
 695
 696# ---------------------------------------------------------------------------
 697# Morphology
 698# ---------------------------------------------------------------------------
 699
 700def compute_crista_morphology(
 701    crista_mask: np.ndarray,
 702    voxel_size: Union[float, Dict[str, float]],
 703    method: str = "both",
 704) -> Dict[str, float]:
 705    """Compute crista shape metrics from binary mask.
 706
 707    Args:
 708        crista_mask: Binary crista segmentation.
 709        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 710        method: "area" | "medial_axis" | "both".
 711
 712    Returns:
 713        Dict with cristae_surface_area_nm2 (area/both) and avg_thickness_nm (medial_axis/both).
 714        avg_thickness_nm is 2 × mean distance-transform value at skeleton voxels.
 715    """
 716    if method not in ("area", "medial_axis", "both"):
 717        raise ValueError(f"method must be 'area', 'medial_axis', or 'both', got {method!r}")
 718
 719    sampling = _to_sampling(voxel_size, crista_mask.ndim)
 720    result: Dict[str, float] = {}
 721
 722    if method in ("area", "both"):
 723        result["cristae_surface_area_nm2"] = _surface_area(crista_mask, sampling)
 724
 725    if method in ("medial_axis", "both"):
 726        result["avg_thickness_nm"] = _medial_axis_thickness_nm(crista_mask, sampling)
 727
 728    return result
 729
 730
 731# ---------------------------------------------------------------------------
 732# Per-mitochondrion statistics
 733# ---------------------------------------------------------------------------
 734
 735def _single_mito_row(
 736    label: int,
 737    bbox: tuple,
 738    mito_crop: np.ndarray,
 739    crista_crop: np.ndarray,
 740    membrane_crop: np.ndarray,
 741    voxel_size: Union[float, Dict[str, float]],
 742    sampling: np.ndarray,
 743    voxel_vol: float,
 744    vol_shape: tuple,
 745    border_radius: int,
 746    method: str = "skip",
 747    inner_n_jobs: int = 1,
 748    lumen_crop: Optional[np.ndarray] = None,
 749) -> Dict[str, float]:
 750    """Compute the statistics row for a single mitochondrion instance.
 751
 752    Factored out of :func:`compute_mito_crista_statistics` so the per-mito work (which is
 753    independent between instances) can be parallelised. Takes the bbox-cropped label arrays
 754    (``mito_crop`` = label array cropped to ``bbox``; ``crista_crop``/``membrane_crop`` the
 755    matching binary crops), so a worker only touches its own mito's bbox region rather than the
 756    whole volume. ``inner_n_jobs`` is forwarded to the (parallelisable) junction-distance stage.
 757
 758    ``method`` controls only the crista orientation anisotropy — every other metric (marching-cubes
 759    surface areas, geodesic junction distances, EDT proximity/thickness) is computed identically for
 760    all modes. ``"skip"`` (the default) leaves the anisotropy NaN and is the fastest; ``"fast"``
 761    computes it on a 2× downsampled crop (~8× cheaper, a *relative* indicator only, not comparable to
 762    exact); ``"exact"`` uses the full-resolution structure tensor (the dominant cost).
 763
 764    ``lumen_crop`` is the (optional) bbox-cropped eroded lumen from :func:`approximate_membrane`
 765    (``return_lumen=True``) used for the junction geodesic mesh; without it the mesh falls back to
 766    ``mito_local & ~membrane_local`` (used only when a caller supplies their own membrane, and
 767    contaminated near clipped faces). Both the lumen geodesic mesh and the mito outer-surface-area mesh
 768    are trimmed/opened at faces where the mito is clipped by the volume boundary: the lumen via
 769    :func:`_open_trimmed_mesh` (trimmed to the certain region and left open, so geodesics do not
 770    shortcut across a cap), the mito surface via ``closed_faces`` (open, so a fabricated cap is not
 771    counted as membrane area). The membrane distance transform is freed before the (memory-heavy)
 772    orientation stage to cap peak memory.
 773    """
 774    ndim = mito_crop.ndim
 775    touches_border = any(
 776        bbox[i] < border_radius or bbox[i + ndim] > vol_shape[i] - border_radius
 777        for i in range(ndim)
 778    )
 779
 780    mito_local = mito_crop == label
 781    crista_local = crista_crop & mito_local
 782    membrane_local = membrane_crop & mito_local
 783
 784    mito_vol = float(mito_local.sum()) * voxel_vol
 785    crista_vol = float(crista_local.sum()) * voxel_vol
 786
 787    boundary = np.array(
 788        [[bbox[a] == 0, bbox[a + ndim] == vol_shape[a]] for a in range(ndim)], dtype=bool
 789    )
 790
 791    has_crista = crista_local.any()
 792    has_membrane = membrane_local.any()
 793    mito_surface = _surface_area(mito_local, sampling, closed_faces=~boundary)
 794
 795    if has_crista and has_membrane:
 796        contact_labels_local, contact_summary = detect_contact_sites(crista_local, membrane_local, voxel_size)
 797        membrane_distance = distance_transform(~membrane_local, sampling=sampling.tolist(), number_of_threads=1)
 798        _, proximity = compute_crista_proximity(
 799            crista_local, membrane_local, voxel_size, membrane_distance=membrane_distance
 800        )
 801        lumen_local = (lumen_crop & mito_local) if lumen_crop is not None else (mito_local & ~membrane_local)
 802        lumen_mesh = _open_trimmed_mesh(lumen_local, sampling, border_radius, boundary)
 803        mesh_verts, mesh_faces = lumen_mesh if lumen_mesh is not None else (None, None)
 804        _, junction_dist = compute_junction_distances(
 805            contact_labels_local, membrane_local, voxel_size,
 806            surface_area_nm2=mito_surface, n_jobs=inner_n_jobs,
 807            mesh_vertices=mesh_verts, mesh_faces=mesh_faces,
 808        )
 809        del membrane_distance, contact_labels_local
 810    else:
 811        contact_summary = {"contact_voxel_count": 0, "crista_junction_count": 0, "contact_volume_nm3": 0.0}
 812        proximity = {"median_nm": np.nan}
 813        junction_dist = dict(_JUNCTION_DISTANCE_NAN)
 814
 815    if has_crista:
 816        morph = compute_crista_morphology(crista_local, voxel_size)
 817        crista_surface = morph.get("cristae_surface_area_nm2", np.nan)
 818        avg_thickness_nm = morph.get("avg_thickness_nm", np.nan)
 819        if method == "skip":
 820            crista_orientation_anisotropy = np.nan
 821        elif method == "fast":
 822            crista_orientation_anisotropy = _downsampled_orientation_anisotropy(
 823                crista_local, voxel_size, factor=2
 824            )
 825        else:
 826            anisotropy = compute_crista_orientation(crista_local, voxel_size)
 827            crista_orientation_anisotropy = float(np.mean(anisotropy[crista_local]))
 828    else:
 829        crista_orientation_anisotropy = np.nan
 830        crista_surface = np.nan
 831        avg_thickness_nm = np.nan
 832
 833    if mito_surface and mito_surface > 0 and np.isfinite(crista_surface):
 834        crista_to_mito_surface_ratio = crista_surface / mito_surface
 835    else:
 836        crista_to_mito_surface_ratio = np.nan
 837
 838    return {
 839        "mito_label_id": int(label),
 840        "mito_touches_border": touches_border,
 841        "mito_volume_nm3": mito_vol,
 842        "crista_volume_nm3": crista_vol,
 843        "crista_fraction": crista_vol / mito_vol if mito_vol > 0 else np.nan,
 844        "contact_voxel_count": contact_summary["contact_voxel_count"],
 845        "crista_junction_count": contact_summary["crista_junction_count"],
 846        "contact_volume_nm3": contact_summary["contact_volume_nm3"],
 847        "avg_crista_to_membrane_nm": proximity["median_nm"],
 848        "mean_nn_junction_distance_nm": junction_dist["mean_nn_junction_distance_nm"],
 849        "median_nn_junction_distance_nm": junction_dist["median_nn_junction_distance_nm"],
 850        "junction_clustering_index": junction_dist["junction_clustering_index"],
 851        "crista_orientation_anisotropy": crista_orientation_anisotropy,
 852        "cristae_surface_area_nm2": crista_surface,
 853        "mito_surface_area_nm2": mito_surface,
 854        "crista_to_mito_surface_ratio": crista_to_mito_surface_ratio,
 855        "avg_thickness_nm": avg_thickness_nm,
 856    }
 857
 858
 859def compute_mito_crista_statistics(
 860    crista_mask: np.ndarray,
 861    mito_segmentation: np.ndarray,
 862    voxel_size: Union[float, Dict[str, float]],
 863    membrane_mask: Optional[np.ndarray] = None,
 864    membrane_thickness_nm: float = 8.0,
 865    border_gap_nm: Optional[float] = None,
 866    method: str = "skip",
 867    n_jobs: int = 1,
 868    verbose: bool = False,
 869    progress_callback: Optional[Callable[[int, int], None]] = None,
 870    membrane_mode: str = "slice_2d",
 871    lumen_mask: Optional[np.ndarray] = None,
 872) -> pd.DataFrame:
 873    """Compute all crista metrics organised by mitochondrial instance.
 874
 875    Args:
 876        crista_mask: Binary crista segmentation (global volume).
 877        mito_segmentation: Instance label array (background = 0).
 878        voxel_size: Voxel size in nm — scalar or dict with "z"/"y"/"x" keys.
 879        membrane_mask: Precomputed membrane mask; recomputed if None.
 880        membrane_thickness_nm: Membrane shell thickness used if membrane_mask is None.
 881        border_gap_nm: Border suppression distance passed to approximate_membrane;
 882            defaults to membrane_thickness_nm when None.
 883        method: How the crista orientation anisotropy is computed — this is the ONLY metric that
 884            differs between modes; surface areas (marching cubes), junction distances (geodesic
 885            along the membrane) and thickness/proximity (EDT) are computed identically for all of
 886            them. ``"skip"`` (default) does not compute orientation at all
 887            (``crista_orientation_anisotropy`` is NaN) and is the fastest — use it when only the other
 888            metrics are needed. ``"fast"`` computes the anisotropy on a 2× downsampled crista crop
 889            (~8× cheaper — the structure tensor is by far the dominant cost); the resulting value is a
 890            *relative* indicator that preserves the ordering between mitochondria but is systematically
 891            different in magnitude from the full-resolution value and is NOT comparable to
 892            ``method="exact"``. ``"exact"`` computes the anisotropy from the full-resolution structure
 893            tensor (use it when the magnitude must be precise).
 894        n_jobs: Number of workers for processing mitochondria in parallel (they are
 895            independent). 1 (default) runs serially; other values use a ``concurrent.futures``
 896            thread pool (-1 = all cores). Results are identical regardless of n_jobs.
 897        verbose: If True, show a terminal tqdm progress bar over mitochondria.
 898        progress_callback: Optional callable invoked once per completed mitochondrion with
 899            (completed_count, total_count) — e.g. to drive a napari progress bar. It is
 900            always called from the calling thread (the futures are consumed here as they
 901            complete), so GUI updates from it need no cross-thread marshaling.
 902        membrane_mode: How the membrane shell is built when ``membrane_mask`` is None —
 903            ``"slice_2d"`` (default, per-Z-slice 2D erosion, z-parallel) or ``"shell_3d"`` (connected
 904            3D shell). See :func:`approximate_membrane`.
 905        lumen_mask: Optional eroded-mito interior matching ``membrane_mask``, i.e. the second return
 906            value of ``approximate_membrane(..., return_lumen=True)``. It is the clean single-wall
 907            surface the junction geodesics run along. Only used when ``membrane_mask`` is also
 908            supplied (when the membrane is built here, the matching lumen is derived automatically);
 909            when neither is available the geodesic mesh falls back to ``mito & ~membrane``, which is
 910            contaminated by the membrane's border-gap suppression near clipped volume faces.
 911
 912    The junction nearest-neighbour distances are geodesics along the eroded-mito surface mesh
 913    (``bioimage_cpp.distance.geodesic_distances_mesh``); for a mito with no usable mesh (empty
 914    membrane / degenerate mesh) those columns are NaN.
 915
 916    Implementation notes: each mito is pre-cropped to its bounding box by basic slicing (views, so
 917    cropping is memory-free). Parallelism is adaptive and single-level (never oversubscribed): with
 918    many mitochondria the work is parallelised *across* them on a ``concurrent.futures``
 919    ``ThreadPoolExecutor`` — the heavy per-mito stages (structure tensor, EDT, geodesics) are
 920    GIL-releasing C++, so threads scale them — with each worker's inner stages kept single-threaded
 921    (the EDT/geodesic solvers are called with ``number_of_threads=1``); with few mitochondria they run
 922    serially and each mito's junction-distance stage gets all cores. The concurrent worker count is
 923    additionally capped so the combined per-mito working set (tensor components + label crops,
 924    ~40 bytes/voxel of the largest mito) fits in RAM. Results stream in as they complete and are
 925    finally sorted by label for an n_jobs-independent ordering.
 926
 927    Returns:
 928        DataFrame with one row per mito instance:
 929        label | mito_volume_nm3 | crista_volume_nm3 | crista_fraction |
 930        contact_voxel_count | crista_junction_count | contact_volume_nm3 |
 931        avg_crista_to_membrane_nm | mean_nn_junction_distance_nm | median_nn_junction_distance_nm |
 932        junction_clustering_index | crista_orientation_anisotropy | cristae_surface_area_nm2 |
 933        mito_surface_area_nm2 | crista_to_mito_surface_ratio | avg_thickness_nm.
 934        cristae_surface_area_nm2 is the crista surface area; crista_to_mito_surface_ratio is
 935        crista surface / mitochondrial outer-membrane surface (can exceed 1 for folded cristae).
 936        The *_nn_junction_distance_nm columns are geodesic nearest-neighbour distances between
 937        crista-membrane junctions along the membrane; junction_clustering_index is a Clark-Evans
 938        index (< 1 clustered, ~ 1 random, > 1 dispersed). ``crista_orientation_anisotropy`` is
 939        computed at full resolution for ``method="exact"``, on a downsampled crop (relative-only,
 940        not comparable) for ``method="fast"``, and left NaN for ``method="skip"``.
 941    """
 942    if method not in ("fast", "exact", "skip"):
 943        raise ValueError(f"method must be 'fast', 'exact', or 'skip', got {method!r}")
 944    if membrane_mask is None:
 945        membrane_mask, lumen_mask = approximate_membrane(
 946            mito_segmentation, voxel_size, membrane_thickness_nm, border_gap_nm,
 947            n_jobs=n_jobs, membrane_mode=membrane_mode, return_lumen=True,
 948        )
 949
 950    ndim = mito_segmentation.ndim
 951    sampling = _to_sampling(voxel_size, ndim)
 952    voxel_vol = float(np.prod(sampling))
 953    crista_binary = crista_mask.astype(bool)
 954    vol_shape = mito_segmentation.shape
 955    border_radius = _gap_radius(voxel_size, membrane_thickness_nm, border_gap_nm, ndim)
 956
 957    tasks = []
 958    for prop in regionprops(mito_segmentation):
 959        bbox = prop.bbox
 960        slices = tuple(slice(bbox[i], bbox[i + ndim]) for i in range(ndim))
 961        lumen_crop = None if lumen_mask is None else lumen_mask[slices]
 962        tasks.append((
 963            int(prop.label), bbox,
 964            mito_segmentation[slices], crista_binary[slices], membrane_mask[slices],
 965            lumen_crop,
 966        ))
 967    total = len(tasks)
 968
 969    def _run(task, inner_n_jobs):
 970        label, bbox, mito_crop, crista_crop, membrane_crop, lumen_crop = task
 971        return _single_mito_row(
 972            label, bbox, mito_crop, crista_crop, membrane_crop,
 973            voxel_size, sampling, voxel_vol, vol_shape, border_radius,
 974            method=method, inner_n_jobs=inner_n_jobs, lumen_crop=lumen_crop,
 975        )
 976
 977    n_workers = os.cpu_count() if n_jobs == -1 else max(1, n_jobs)
 978    across = n_workers > 1 and total >= n_workers
 979
 980    rows = []
 981
 982    def _consume(results):
 983        for i, row in enumerate(
 984            tqdm(results, total=total, desc="Cristae analysis", disable=not verbose), start=1
 985        ):
 986            rows.append(row)
 987            if progress_callback is not None:
 988                progress_callback(i, total)
 989
 990    if across:
 991        max_voxels = max(int(task[2].size) for task in tasks)
 992        across_workers = _bounded_workers(n_jobs, per_worker_bytes=max_voxels * 40)
 993        # Parallelise across mitochondria with a thread pool (the heavy per-mito stages — structure
 994        # tensor, EDT, geodesics — are GIL-releasing C++). Each worker's inner stages run
 995        # single-threaded (``_run(task, 1)`` passes ``number_of_threads=1`` down to the EDT/geodesic
 996        # solvers) so the across-mito threads do not oversubscribe the cores. Results stream in as they
 997        # complete (``as_completed``) to drive the progress bar; rows are label-sorted below.
 998        with futures.ThreadPoolExecutor(across_workers) as tp:
 999            submitted = [tp.submit(_run, task, 1) for task in tasks]
1000            _consume(future.result() for future in futures.as_completed(submitted))
1001    else:
1002        _consume(_run(task, n_workers) for task in tasks)
1003
1004    rows.sort(key=lambda row: row["mito_label_id"])
1005    return pd.DataFrame(rows)