OSCR

Quantitative Assessment of Mesoscale Cellular Order and Organization in the Mouse Hippocampus.

Code ↔ Paper

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

The 3 matches
  1. [1] § Results › Intensity Processing and Segmentation of Nuclei in High‐Density Brain Areas ↔ code/FunctionsClean.py, lines 68 to 107 · score 0.60 · background illumination, Gaussian filter, dilation, fitting · Report a mismatch
  2. [2] § Material and Methods › Image Acquisition ↔ code/cbhs3d_onefile.py, lines 95 to 96 · score 0.54 · Zeiss LSM, 0.1135 um · Report a mismatch
  3. [3] § Results › Intensity Processing and Segmentation of Nuclei in High‐Density Brain Areas ↔ code/cbhs3d_onefile.py, lines 357 to 395 · score 0.50 · Gaussian filter, noise, signals, dilation, erosion, components · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,250 lines · 55 KB · no license · 2 matches

  1. #!/usr/bin/env python
  2. r"""CBHS 3D pipeline, single file: raw TIFF z-stack -> 3D segmented nuclear masks + shape metrics.
  3. Implements the README "Current pipeline" (the chain that produced ``M_synth3_yuki_*`` and
  4. ``M_synth3_yuki_anchored``) and the README "Shape metrics" section, with no imports from
  5. the ``cbhs3d`` package -- every function is inlined here.
  6. raw z-slices --(1) cellpose 2D, independent-----> M_indep.tiff
  7. M_indep + raw --(2) synthetic DAPI v3-----------> synth_dapi_v3.tiff
  8. synth --(3) StarDist 3D (own env)-------> M_stardist_v3.tiff
  9. M_stardist --(4) Yukimatsu refine (reduced)--> M_synth3_yuki_reduced.tiff
  10. + M_indep --(5) anchor 3D->2D + finish------> M_synth3_yuki_reduced_anchored.tiff
  11. final masks --(6) shape metrics---------------> per_nucleus.csv, per_slice.csv
  12. Stage 3 needs TensorFlow, which the analysis environment does not have. Pass
  13. ``--stardist-python`` and this file re-invokes *itself* with that interpreter for stage 3
  14. only (``--stardist-stage``); omit it and stage 3 runs in-process instead. Stages are
  15. cached: an existing output file skips its stage unless ``--force``.
  16. PROVENANCE (original script -> section below)
  17. Sec. 2 Cellpose front-ends ............ Cellpose3DSegment.py, Compare3DvsMean2D.cellpose_segment
  18. Sec. 5a round_caps, split_touching_shape_aware, relabel_compact_1n
  19. Yukimatsu.py
  20. Sec. 5b cleanup_labels ................. YukimatsuOnSynth2Stardist.cleanup_morphology
  21. VARIANTS ....................... YukimatsuOnSynth3.py
  22. Sec. 5c link_2d_to_3d .................. HybridSegment.build_hybrid
  23. --- included with explicit approval, not from the four scripts above ---
  24. Sec. 3 synthetic DAPI ................. BuildSyntheticDapi.py / _v2.py / _v3.py, PipelineV2.py st.1
  25. Sec. 4 run_stardist3d ................. Stardist3DSegment.py, Stardist3DOnSynth*.py, PipelineV2 st.2
  26. Sec. 5d anchor_to_2d, anchor_and_finish AnchorTo2D.py
  27. Sec. 6 metrics_3d, metrics_2d_* ....... Compare3DvsMean2D.py
  28. per_slice_table ................ SliceSampling3Dvs2D.build_per_slice
  29. Sec. 1 z-stack loading, MIPs, save .... library glue (no original-script equivalent)
  30. NOT included: audit_fragmentation (AuditHybridYuki.py) -- not needed by this chain.
  31. Usage
  32. python cbhs3d_onefile.py --raw-dir "<...IHCFinal_6wETHO_4_1 c=0.tiff_files>" \
  33. --out-dir out \
  34. --cellpose-weights <cpsam_20250704_112416> \
  35. --stardist-model-dir <Ferme-et-al_2025/stardist> \
  36. --stardist-python C:\Users\Admin\.conda\envs\stardist-gpu\python.exe \
  37. --synth v3 --variant reduced --anchor --metrics
  38. Environments
  39. analysis: numpy, scipy, scikit-image, pandas, tifffile, tqdm, cellpose, torch (GPU rounding)
  40. stardist-gpu: Python 3.10, TF 2.10, CUDA 11.2, stardist, csbdeep (stage 3 only)
  41. """
  42. import argparse
  43. import glob
  44. import json
  45. import math
  46. import os
  47. import re
  48. import subprocess
  49. import sys
  50. os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
  51. import numpy as np
  52. import tifffile
  53. from scipy.ndimage import (binary_dilation, binary_erosion, binary_fill_holes,
  54. distance_transform_edt, find_objects, gaussian_filter,
  55. label as cc_label, zoom)
  56. from scipy.spatial import ConvexHull
  57. from skimage.feature import peak_local_max
  58. from skimage.measure import marching_cubes, mesh_surface_area, regionprops
  59. from skimage.segmentation import watershed
  60. from tqdm import tqdm
  61. try: # only Sec. 6 needs pandas; the stardist env may lack it
  62. import pandas as pd
  63. except ImportError:
  64. pd = None
  65. try: # re-exported at scipy.spatial only from SciPy 1.8
  66. from scipy.spatial import QhullError
  67. except ImportError:
  68. from scipy.spatial.qhull import QhullError
  69. try: # torch missing or broken -> CPU morphology path
  70. import torch
  71. import torch.nn.functional as F
  72. _HAS_CUDA = torch.cuda.is_available()
  73. except Exception:
  74. torch = None
  75. F = None
  76. _HAS_CUDA = False
  77. # ============================================================================
  78. # 0. CONFIG -- calibration and parameter presets (config.py)
  79. # Every preset is a plain dict; copy and edit rather than mutating in place.
  80. # ============================================================================
  81. # (Z, Y, X) voxel size in micrometres -- Zeiss LSM, 40x, 0.33 um z-step
  82. VOXEL_UM = (0.33, 0.1135, 0.1135)
  83. # Cellpose 2D (custom cpsam model)
  84. CELLPOSE = dict(diameter=45, flow_threshold=0.6, min_size=30, stitch_threshold=0.4)
  85. # StarDist 3D (Ferme et al. 2025, Zenodo 10.5281/zenodo.15120155, model_A)
  86. STARDIST = dict(model_name="model_A", n_tiles=(2, 8, 8),
  87. prob_thresh=None, nms_thresh=None, normalize_pct=(1.0, 99.8))
  88. # Synthetic DAPI (2D masks -> 3D-segmentable volume)
  89. SYNTH_V1 = dict(kind="binary", intensity_interior=255, intensity_rim=160, rim_width_vox=2)
  90. SYNTH_V2 = dict(kind="masked", xy_erosion_vox=2, z_dilation=True, xy_smooth_sigma=0.7,
  91. pct_lo=2.0, pct_hi=98.0, noise_sigma=5.0, seed=42)
  92. SYNTH_V3 = dict(SYNTH_V2, xy_erosion_vox=1) # reduced-rim variant (current)
  93. # Yukimatsu cap rounding
  94. ROUNDING_DEFAULT = dict(
  95. z_boost=2.6, max_z_scale=4, gaussian_sigma=1.1,
  96. close_radius_xy_um=0.35, close_radius_z_um=0.90, open_radius_um=0.30,
  97. virtual_z_margin_um=1.2, imagej_interp=True, z_push_um=0.30,
  98. )
  99. # reduced rounding: no z-push, smaller Gaussian and z-closing radii, no opening (synth v3 variant)
  100. ROUNDING_REDUCED = dict(ROUNDING_DEFAULT, z_push_um=0.0, gaussian_sigma=0.6,
  101. close_radius_xy_um=0.35, close_radius_z_um=0.45,
  102. open_radius_um=0.0, virtual_z_margin_um=0.4)
  103. # rounding applied after anchoring to the 2D masks; closes residual slice steps
  104. ROUNDING_POST_ANCHOR = dict(ROUNDING_DEFAULT, z_push_um=0.25, gaussian_sigma=0.6,
  105. close_radius_xy_um=0.20, close_radius_z_um=0.70,
  106. open_radius_um=0.0, virtual_z_margin_um=1.0)
  107. # rounding with larger radii, applied to Cellpose-stitched input (Compare3DvsMean2D M_final)
  108. ROUNDING_STITCHED = dict(ROUNDING_DEFAULT, gaussian_sigma=2.0,
  109. close_radius_xy_um=0.80, close_radius_z_um=1.50,
  110. open_radius_um=0.70, virtual_z_margin_um=1.5, z_push_um=0.45)
  111. # Yukimatsu shape-aware split (intensity-free, DT markers + neck/shape guards)
  112. SPLIT_DEFAULT = dict(min_peak_sep_um=1.0, max_peaks=5, min_rel_dt_at_peak=0.22,
  113. neck_rel=0.35, min_child_voxels=2600, min_child_vol_ratio=0.25,
  114. elong_merge_if_ratio_below=1.15)
  115. SPLIT_STITCHED = dict(min_peak_sep_um=0.8, max_peaks=8, min_rel_dt_at_peak=0.18,
  116. neck_rel=0.25, min_child_voxels=1500, min_child_vol_ratio=0.18,
  117. elong_merge_if_ratio_below=1.10)
  118. # Per-label cleanup (largest CC + hole fill + ellipsoid opening + size filter)
  119. CLEANUP_DEFAULT = dict(open_radius_um=1.2, min_voxels=200)
  120. CLEANUP_REDUCED = dict(open_radius_um=0.5, min_voxels=200)
  121. CLEANUP_ANCHOR = dict(open_radius_um=0.6, min_voxels=200) # removes shelves 1-3 slices thick
  122. CLEANUP_NO_OPENING = dict(open_radius_um=0.0, min_voxels=200)
  123. # Named end-to-end variants (name -> rounding, cleanup); None = no rounding.
  124. # Former names: default = published, reduced = gentle, no_rounding = minimal.
  125. VARIANTS = {
  126. "default": (ROUNDING_DEFAULT, CLEANUP_DEFAULT),
  127. "reduced": (ROUNDING_REDUCED, CLEANUP_REDUCED),
  128. "no_rounding": (None, CLEANUP_NO_OPENING),
  129. }
  130. # 2D <-> 3D linking
  131. HYBRID_OVERLAP_THR = 0.30 # Cellpose blob must have >=30 % of its pixels in one 3D label
  132. ANCHOR_PAD_PX = 20 # Cellpose blob may exceed the 3D bbox by this many pixels
  133. # ============================================================================
  134. # 1. I/O (io.py -- library glue)
  135. # ============================================================================
  136. Z_INDEX_PATTERN = re.compile(r"b0z(\d+)c", re.IGNORECASE)
  137. def load_zstack_from_folder(folder, pattern=Z_INDEX_PATTERN):
  138. """Stack per-z TIFFs (Zeiss export ``..._b0z<N>c0...tiff``) into a (Z, Y, X) volume.
  139. RGB/RGBA exports with R==G==B are reduced to channel 0."""
  140. files = [f for f in glob.glob(os.path.join(folder, "*.tif*"))
  141. if "_metadata" not in os.path.basename(f)]
  142. indexed = []
  143. for f in files:
  144. m = pattern.search(os.path.basename(f))
  145. if m:
  146. indexed.append((int(m.group(1)), f))
  147. if not indexed:
  148. raise FileNotFoundError(f"no per-z tiffs matching {pattern.pattern!r} in {folder}")
  149. indexed.sort()
  150. slices = []
  151. for _, f in indexed:
  152. a = tifffile.imread(f)
  153. if a.ndim == 3 and a.shape[-1] in (3, 4):
  154. a = a[..., 0]
  155. slices.append(a)
  156. return np.stack(slices, axis=0)
  157. def imread_3d(path):
  158. arr = tifffile.imread(path)
  159. if arr.ndim != 3:
  160. raise ValueError(f"{os.path.basename(path)} is not 3D (ndim={arr.ndim})")
  161. return arr
  162. def load_raw(path):
  163. """A folder of per-z TIFFs or a single 3D TIFF."""
  164. return load_zstack_from_folder(path) if os.path.isdir(path) else imread_3d(path)
  165. def n_labels(masks):
  166. return int(np.unique(masks).size - (1 if (masks == 0).any() else 0))
  167. def relabel_compact_1n(masks, dtype=np.uint32):
  168. """Relabel to 1..N preserving label order; keeps all labels. [Yukimatsu.relabel_compact_1N]"""
  169. labels = np.unique(masks)
  170. labels = labels[labels > 0]
  171. if labels.size == 0:
  172. return masks.astype(dtype, copy=False)
  173. lut = np.zeros(int(masks.max()) + 1, dtype=dtype)
  174. lut[labels] = np.arange(1, labels.size + 1, dtype=dtype)
  175. return lut[masks]
  176. def ensure_unique(path):
  177. if not os.path.exists(path):
  178. return path
  179. stem, ext = os.path.splitext(path)
  180. k = 2
  181. while os.path.exists(f"{stem}__v{k}{ext}"):
  182. k += 1
  183. return f"{stem}__v{k}{ext}"
  184. def save_mips(stem, vol):
  185. """Write XY / XZ / YZ maximum-intensity projections next to ``stem``."""
  186. for name, axis in (("XY", 0), ("XZ", 1), ("YZ", 2)):
  187. tifffile.imwrite(ensure_unique(f"{stem}_MIP_{name}.tiff"),
  188. np.max(vol, axis=axis).astype(np.uint16, copy=False))
  189. def save_labels(path, masks):
  190. tifffile.imwrite(path, np.asarray(masks).astype(np.uint32, copy=False))
  191. # ============================================================================
  192. # 2. CELLPOSE FRONT-ENDS [Cellpose3DSegment.py, Compare3DvsMean2D.cellpose_segment]
  193. # Require the ``cellpose`` package (CUDA optional).
  194. # ============================================================================
  195. def _cellpose_model(weights_path, gpu=True):
  196. from cellpose import models
  197. return models.CellposeModel(gpu=gpu, pretrained_model=weights_path)
  198. def cellpose_segment_independent(stack, weights_path, diameter=CELLPOSE["diameter"],
  199. flow_threshold=CELLPOSE["flow_threshold"],
  200. min_size=CELLPOSE["min_size"], gpu=True):
  201. """2D per slice, no stitching. Labels are offset per slice so they are globally unique.
  202. This is stage 1 of the current pipeline."""
  203. model = _cellpose_model(weights_path, gpu)
  204. out = model.eval([stack[z] for z in range(stack.shape[0])], diameter=diameter,
  205. do_3D=False, flow_threshold=flow_threshold, min_size=min_size)
  206. per_slice = out[0] if isinstance(out, tuple) else out
  207. full = np.zeros(stack.shape, dtype=np.uint32)
  208. offset = 0
  209. for z, m in enumerate(per_slice):
  210. m = np.asarray(m, dtype=np.uint32)
  211. nz = m > 0
  212. if nz.any():
  213. m[nz] += offset
  214. offset = int(m.max())
  215. full[z] = m
  216. return full
  217. def cellpose_segment_stitched(stack, weights_path, stitch_threshold=CELLPOSE["stitch_threshold"],
  218. diameter=CELLPOSE["diameter"], flow_threshold=CELLPOSE["flow_threshold"],
  219. min_size=CELLPOSE["min_size"], gpu=True):
  220. """2D per slice + z-stitching -> labels consistent across z. Returns uint32 (Z, Y, X).
  221. This is the legacy (paper Fig. 3 F/G) stage 1."""
  222. model = _cellpose_model(weights_path, gpu)
  223. out = model.eval(stack, z_axis=0, diameter=diameter, do_3D=False,
  224. stitch_threshold=stitch_threshold, flow_threshold=flow_threshold,
  225. min_size=min_size)
  226. masks = out[0] if isinstance(out, tuple) else out
  227. return np.asarray(masks).astype(np.uint32)
  228. def cellpose_segment_3d(stack, weights_path, voxel_um=VOXEL_UM, diameter=CELLPOSE["diameter"],
  229. flow_threshold=CELLPOSE["flow_threshold"], min_size=CELLPOSE["min_size"],
  230. gpu=True):
  231. """Native Cellpose 3D (do_3D=True): XY/XZ/YZ flows combined. Slow (30-60 min on GPU)."""
  232. model = _cellpose_model(weights_path, gpu)
  233. anisotropy = voxel_um[0] / voxel_um[1]
  234. out = model.eval(stack, z_axis=0, diameter=diameter, do_3D=True,
  235. flow_threshold=flow_threshold, min_size=min_size, anisotropy=anisotropy)
  236. masks = out[0] if isinstance(out, tuple) else out
  237. return np.asarray(masks).astype(np.uint32)
  238. CELLPOSE_MODES = {
  239. "independent": cellpose_segment_independent, # current pipeline, stage 1
  240. "stitched": cellpose_segment_stitched, # legacy chain
  241. "3d": cellpose_segment_3d,
  242. }
  243. # ============================================================================
  244. # 3. SYNTHETIC DAPI [BuildSyntheticDapi.py / _v2.py / _v3.py, PipelineV2.py stage 1]
  245. # Volumes built from the independent 2D masks, so that a 3D segmenter (StarDist 3D)
  246. # can link slices while the in-plane shape stays that of the 2D masks.
  247. #
  248. # v1 binary interior + darker rim, hard edges -> build_synth_dapi_binary
  249. # v2 raw DAPI masked by 2-vox-eroded blobs + refinements -> build_synth_dapi_masked(2)
  250. # v3 as v2 but 1-vox erosion (reduced rim, current) -> build_synth_dapi_masked(1)
  251. # ============================================================================
  252. def z_dilate_labels(M):
  253. """Propagate every 2D blob one slice up and down, into empty voxels only.
  254. Recovers faint cap slices without bridging stacked nuclei."""
  255. Z = M.shape[0]
  256. out = M.copy()
  257. for z in range(Z - 1):
  258. fill = (out[z + 1] == 0) & (M[z] > 0)
  259. out[z + 1] = np.where(fill, M[z], out[z + 1])
  260. for z in range(Z - 1, 0, -1):
  261. fill = (out[z - 1] == 0) & (M[z] > 0)
  262. out[z - 1] = np.where(fill, M[z], out[z - 1])
  263. return out
  264. def erode_blobs_xy(M, iterations):
  265. """Per-blob 2D erosion; returns a boolean foreground with gaps between neighbours."""
  266. Z = M.shape[0]
  267. fg = np.zeros(M.shape, dtype=bool)
  268. for z in tqdm(range(Z), desc="erode xy", unit="slice", leave=False):
  269. sl2 = M[z]
  270. for idx, sl in enumerate(find_objects(sl2)):
  271. if sl is None:
  272. continue
  273. m = sl2[sl] == idx + 1
  274. if m.any():
  275. fg[z][sl] |= binary_erosion(m, iterations=iterations) if iterations > 0 else m
  276. return fg
  277. def build_synth_dapi_binary(M_indep, intensity_interior=SYNTH_V1["intensity_interior"],
  278. intensity_rim=SYNTH_V1["intensity_rim"],
  279. rim_width_vox=SYNTH_V1["rim_width_vox"]):
  280. """v1: interior=255, rim=160; the rim erosion doubles as the gap to neighbours."""
  281. Z = M_indep.shape[0]
  282. synth = np.zeros(M_indep.shape, dtype=np.uint8)
  283. for z in tqdm(range(Z), desc="synth v1", unit="slice", leave=False):
  284. sl2 = M_indep[z]
  285. for idx, sl in enumerate(find_objects(sl2)):
  286. if sl is None:
  287. continue
  288. m = sl2[sl] == idx + 1
  289. if not m.any():
  290. continue
  291. interior = binary_erosion(m, iterations=rim_width_vox)
  292. sub = synth[z][sl]
  293. sub[m & ~interior] = intensity_rim
  294. sub[interior] = intensity_interior
  295. return synth
  296. def build_synth_dapi_masked(raw, M_indep, xy_erosion_vox=SYNTH_V3["xy_erosion_vox"],
  297. z_dilation=SYNTH_V2["z_dilation"],
  298. xy_smooth_sigma=SYNTH_V2["xy_smooth_sigma"],
  299. pct_lo=SYNTH_V2["pct_lo"], pct_hi=SYNTH_V2["pct_hi"],
  300. noise_sigma=SYNTH_V2["noise_sigma"], seed=SYNTH_V2["seed"]):
  301. """v2/v3: raw DAPI inside eroded 2D blobs, per-3D-connected-component percentile
  302. normalisation (preserves the natural z-cap fade), light XY smoothing and additive
  303. Gaussian noise on signal voxels. Returns uint8."""
  304. if raw.shape != M_indep.shape:
  305. raise ValueError(f"shape mismatch raw {raw.shape} vs masks {M_indep.shape}")
  306. M = z_dilate_labels(M_indep) if z_dilation else M_indep
  307. fg = erode_blobs_xy(M, xy_erosion_vox)
  308. synth = np.zeros(raw.shape, dtype=np.float32)
  309. synth[fg] = raw[fg]
  310. cc, _ = cc_label(fg, structure=np.ones((3, 3, 3)))
  311. out = np.zeros_like(synth)
  312. for idx, sl in enumerate(tqdm(find_objects(cc), desc="normalise cc", unit="cc", leave=False)):
  313. if sl is None:
  314. continue
  315. m = cc[sl] == idx + 1
  316. if not m.any():
  317. continue
  318. vals = synth[sl][m]
  319. if vals.size < 5:
  320. out[sl][m] = vals
  321. continue
  322. p_lo, p_hi = np.percentile(vals, [pct_lo, pct_hi])
  323. if p_hi <= p_lo:
  324. out[sl][m] = vals
  325. continue
  326. out[sl][m] = np.clip((vals - p_lo) / (p_hi - p_lo) * 240.0 + 15.0, 0, 255)
  327. out = gaussian_filter(out, sigma=(0.0, xy_smooth_sigma, xy_smooth_sigma))
  328. if noise_sigma > 0:
  329. rng = np.random.default_rng(seed)
  330. out = out + rng.normal(0.0, noise_sigma, out.shape).astype(np.float32) * (out > 0)
  331. return np.clip(out, 0, 255).astype(np.uint8)
  332. def build_synth_dapi(raw, M_indep, version="v3"):
  333. """Dispatch by version name: 'v1' | 'v2' | 'v3'."""
  334. if version == "v1":
  335. return build_synth_dapi_binary(M_indep, **{k: v for k, v in SYNTH_V1.items() if k != "kind"})
  336. cfg = SYNTH_V2 if version == "v2" else SYNTH_V3
  337. return build_synth_dapi_masked(raw, M_indep, **{k: v for k, v in cfg.items() if k != "kind"})
  338. # ============================================================================
  339. # 4. STARDIST 3D [Stardist3DSegment.py, Stardist3DOnSynth*.py, PipelineV2.py stage 2]
  340. # Needs the ``stardist-gpu`` environment (Python 3.10, TF 2.10, CUDA 11.2,
  341. # stardist, csbdeep). Nothing else in this file imports TensorFlow, and the
  342. # imports below are function-local so the analysis env can load this module.
  343. # ============================================================================
  344. def run_stardist3d(stack, model_basedir, model_name=STARDIST["model_name"],
  345. n_tiles=STARDIST["n_tiles"], prob_thresh=STARDIST["prob_thresh"],
  346. nms_thresh=STARDIST["nms_thresh"], normalize_pct=STARDIST["normalize_pct"],
  347. verbose=True):
  348. """Percentile-normalise, pad Z to a multiple of 4 (UNet alignment), predict, crop back.
  349. Returns uint32 labels of the input shape."""
  350. os.environ.setdefault("TF_CPP_MIN_LOG_LEVEL", "2")
  351. os.environ.setdefault("TF_ENABLE_ONEDNN_OPTS", "0")
  352. from csbdeep.utils import normalize
  353. from stardist.models import StarDist3D
  354. x = normalize(stack.astype(np.float32, copy=False), *normalize_pct, axis=None)
  355. Z = x.shape[0]
  356. Z_pad = ((Z + 3) // 4) * 4
  357. if Z_pad != Z:
  358. x = np.pad(x, ((0, Z_pad - Z), (0, 0), (0, 0)), mode="edge")
  359. model = StarDist3D(None, name=model_name, basedir=model_basedir)
  360. kw = dict(n_tiles=n_tiles, verbose=verbose)
  361. if prob_thresh is not None:
  362. kw["prob_thresh"] = prob_thresh
  363. if nms_thresh is not None:
  364. kw["nms_thresh"] = nms_thresh
  365. labels, _ = model.predict_instances(x, **kw)
  366. return np.asarray(labels[:Z]).astype(np.uint32)
  367. # ============================================================================
  368. # 5. YUKIMATSU -- label-volume refinement for anisotropic nuclear masks
  369. #
  370. # (a) round_caps z-upsample, z-push, Gaussian smooth, ellipsoidal
  371. # close/open, max-pool back [Yukimatsu.py]
  372. # split_touching_shape_aware DT watershed with neck / size / elongation guards
  373. # [Yukimatsu.py]
  374. # (b) cleanup_labels largest CC + hole fill + ellipsoid opening + size
  375. # filter [YukimatsuOnSynth2Stardist.py]
  376. # (c) link_2d_to_3d each 2D blob takes the label of the 3D object it
  377. # overlaps most [HybridSegment.py]
  378. # (d) anchor_to_2d back-match: rebuild every 3D label slice by slice
  379. # from the 2D masks [AnchorTo2D.py]
  380. # (e) refine chain: compact -> round -> [split] -> cleanup ->
  381. # [anchor]
  382. #
  383. # All radii are in micrometres and converted with the (Z, Y, X) voxel size.
  384. # ============================================================================
  385. def _r(v, minimum=0):
  386. return int(max(minimum, round(float(v))))
  387. def ellipsoid_se(rz, ry, rx):
  388. rz_s, ry_s, rx_s = max(rz, 1), max(ry, 1), max(rx, 1)
  389. Z, Y, X = np.ogrid[-rz:rz + 1, -ry:ry + 1, -rx:rx + 1]
  390. return (Z * Z) / (rz_s * rz_s) + (Y * Y) / (ry_s * ry_s) + (X * X) / (rx_s * rx_s) <= 1.0
  391. def largest_cc(mask, connectivity=1):
  392. structure = np.ones((3, 3, 3)) if connectivity == 3 else None
  393. cc, n = cc_label(mask, structure=structure)
  394. if n <= 1:
  395. return mask
  396. counts = np.bincount(cc.ravel())
  397. counts[0] = 0
  398. return cc == int(np.argmax(counts))
  399. def _radii_vox(um, voxel_um, minimum=0):
  400. return tuple(_r(um / v, minimum) for v in voxel_um)
  401. # ---- GPU morphology (separable box approximation of the ellipsoid, as in the original) ----
  402. def _gauss1d_gpu(s, dev):
  403. if s <= 0:
  404. return torch.tensor([1.0], device=dev)
  405. r = max(1, int(3 * s + 0.5))
  406. x = torch.arange(-r, r + 1, device=dev, dtype=torch.float32)
  407. k = torch.exp(-(x * x) / (2 * s * s))
  408. return k / k.sum()
  409. def _blur3d_gpu(vol, s, dev):
  410. k = _gauss1d_gpu(s, dev)
  411. p = k.numel() // 2
  412. vol = F.conv3d(F.pad(vol, (0, 0, 0, 0, p, p), mode="replicate"), k.view(1, 1, -1, 1, 1))
  413. vol = F.conv3d(F.pad(vol, (0, 0, p, p, 0, 0), mode="replicate"), k.view(1, 1, 1, -1, 1))
  414. vol = F.conv3d(F.pad(vol, (p, p, 0, 0, 0, 0), mode="replicate"), k.view(1, 1, 1, 1, -1))
  415. return vol
  416. def _morph_gpu(vol, rz, ry, rx, op, dev):
  417. def pass1(v, r, axis):
  418. if r <= 0:
  419. return v
  420. shape = [1, 1, 1, 1, 1]
  421. shape[2 + axis] = 2 * r + 1
  422. k = torch.ones(shape, device=dev)
  423. pad = [0, 0, 0, 0, 0, 0]
  424. pad[2 * (2 - axis)] = pad[2 * (2 - axis) + 1] = r
  425. c = F.conv3d(F.pad(v, pad, mode="constant", value=0.0), k)
  426. return (c > 0).float() if op == "dilate" else (c >= (2 * r + 1) - 1e-6).float()
  427. v = pass1(vol, rz, 0)
  428. v = pass1(v, ry, 1)
  429. return pass1(v, rx, 2)
  430. def _round_one_gpu(sub, z_scale, tri, rz_push, sigma, close_r, open_r, dev):
  431. t = torch.from_numpy(sub.astype(np.float32))[None, None].to(dev, non_blocking=True)
  432. if z_scale > 1:
  433. t = F.interpolate(t, scale_factor=(z_scale, 1, 1),
  434. mode="trilinear" if tri else "nearest",
  435. align_corners=False if tri else None)
  436. if rz_push > 0:
  437. k = torch.ones((1, 1, 2 * rz_push + 1, 1, 1), device=dev)
  438. t = (F.conv3d(F.pad(t, (0, 0, 0, 0, rz_push, rz_push)), k) > 0).float()
  439. t = (_blur3d_gpu(t, sigma, dev) >= 0.5).float()
  440. t = _morph_gpu(t, *close_r, "dilate", dev)
  441. t = _morph_gpu(t, *close_r, "erode", dev)
  442. if max(open_r) > 0:
  443. t = _morph_gpu(t, *open_r, "erode", dev)
  444. t = _morph_gpu(t, *open_r, "dilate", dev)
  445. if z_scale > 1:
  446. t = F.max_pool3d(t, kernel_size=(z_scale, 1, 1), stride=(z_scale, 1, 1))
  447. return t.squeeze().detach().cpu().numpy() > 0.5
  448. def _round_one_cpu(sub, z_scale, tri, rz_push, sigma, close_r, open_r):
  449. if z_scale > 1:
  450. hi = zoom(sub.astype(np.float32), zoom=(z_scale, 1, 1), order=1 if tri else 0) >= 0.5
  451. else:
  452. hi = sub.astype(bool, copy=False)
  453. if rz_push > 0:
  454. hi = binary_dilation(hi, structure=np.ones((2 * rz_push + 1, 1, 1), dtype=bool))
  455. hi = gaussian_filter(hi.astype(np.float32), sigma=sigma, mode="nearest") >= 0.5
  456. se = ellipsoid_se(*close_r)
  457. hi = binary_erosion(binary_dilation(hi, structure=se), structure=se)
  458. if max(open_r) > 0:
  459. se2 = ellipsoid_se(*open_r)
  460. hi = binary_dilation(binary_erosion(hi, structure=se2), structure=se2)
  461. if z_scale <= 1:
  462. return hi
  463. Zhi = hi.shape[0]
  464. out = np.zeros(sub.shape, dtype=bool)
  465. for zi in range(sub.shape[0]):
  466. z0 = max(0, min(Zhi, int(round(zi * z_scale))))
  467. z1 = max(z0 + 1, min(Zhi, int(round((zi + 1) * z_scale))))
  468. out[zi] = np.any(hi[z0:z1], axis=0)
  469. return out
  470. # ---- 5a. cap rounding [Yukimatsu.round_caps] ----
  471. def round_caps(masks, voxel_um=VOXEL_UM, cfg=ROUNDING_DEFAULT, use_gpu=True,
  472. xy_pad=4, progress=True):
  473. """Per-label z-cap rounding. Each label is upsampled along Z (``z_boost`` x anisotropy,
  474. capped at ``max_z_scale``), optionally pushed along Z, Gaussian-smoothed, closed and
  475. opened with ellipsoidal elements, then max-pooled back to the original z grid so the
  476. rounded cap may spill into ``virtual_z_margin_um`` of neighbouring original slices.
  477. Larger labels are written first; overlaps are first-come."""
  478. masks = np.asarray(masks)
  479. if masks.max() == 0:
  480. return masks.copy()
  481. z_um, y_um, x_um = map(float, voxel_um)
  482. anis = max(1.0, z_um / min(y_um, x_um))
  483. z_scale = int(min(cfg.get("max_z_scale", 4), max(1, round(anis * cfg.get("z_boost", 2.0)))))
  484. z_hi = z_um / z_scale
  485. sigma = float(cfg.get("gaussian_sigma", 1.8))
  486. tri = bool(cfg.get("imagej_interp", True))
  487. close_r = (_r(cfg.get("close_radius_z_um", 0.90) / z_hi, 1),
  488. _r(cfg.get("close_radius_xy_um", 0.35) / y_um, 1),
  489. _r(cfg.get("close_radius_xy_um", 0.35) / x_um, 1))
  490. open_um = float(cfg.get("open_radius_um", 0.30))
  491. open_r = (_r(open_um / z_hi), _r(open_um / y_um), _r(open_um / x_um))
  492. rz_push = _r(cfg.get("z_push_um", 0.60) / z_hi)
  493. spill_um = max(float(cfg.get("virtual_z_margin_um", 1.2)), z_um + 1e-6)
  494. z_pad = max(1, int(np.ceil(spill_um / z_um)))
  495. gpu = bool(use_gpu and _HAS_CUDA)
  496. dev = torch.device("cuda") if gpu else None
  497. out = np.zeros_like(masks)
  498. sizes = np.bincount(masks.ravel())
  499. bboxes = find_objects(masks)
  500. labs = [i for i in range(1, sizes.size) if sizes[i] > 0 and bboxes[i - 1] is not None]
  501. labs.sort(key=lambda L: sizes[L], reverse=True)
  502. Z, Y, X = masks.shape
  503. it = tqdm(labs, desc=f"round caps ({'GPU' if gpu else 'CPU'})", unit="obj", leave=False) \
  504. if progress else labs
  505. for lab in it:
  506. sl = bboxes[lab - 1]
  507. z0, z1 = max(0, sl[0].start - 1 - z_pad), min(Z, sl[0].stop + 1 + z_pad)
  508. y0, y1 = max(0, sl[1].start - xy_pad), min(Y, sl[1].stop + xy_pad)
  509. x0, x1 = max(0, sl[2].start - xy_pad), min(X, sl[2].stop + xy_pad)
  510. slc = (slice(z0, z1), slice(y0, y1), slice(x0, x1))
  511. sub = masks[slc] == lab
  512. if not sub.any():
  513. continue
  514. if gpu:
  515. res = _round_one_gpu(sub, z_scale, tri, rz_push, sigma, close_r, open_r, dev)
  516. else:
  517. res = _round_one_cpu(sub, z_scale, tri, rz_push, sigma, close_r, open_r)
  518. write = res & (out[slc] == 0)
  519. out[slc][write] = lab
  520. return out
  521. # ---- 5a. shape-aware split [Yukimatsu.split_touching_shape_aware] ----
  522. def _peak_markers(sub, dt, voxel_um, cfg):
  523. fp = ellipsoid_se(*_radii_vox(cfg["min_peak_sep_um"], voxel_um, 1))
  524. peaks = peak_local_max(dt, labels=sub.astype(np.uint8), footprint=fp, exclude_border=False)
  525. if peaks.shape[0] > cfg["max_peaks"]:
  526. peaks = peaks[np.argsort(dt[tuple(peaks.T)])[-cfg["max_peaks"]:]]
  527. if peaks.size == 0:
  528. return None, None
  529. dt_max = float(dt[sub].max())
  530. if dt_max <= 0:
  531. return None, None
  532. peaks = peaks[dt[tuple(peaks.T)] >= cfg["min_rel_dt_at_peak"] * dt_max]
  533. if peaks.shape[0] < 2:
  534. return None, None
  535. markers = np.zeros(dt.shape, dtype=np.int32)
  536. for i, p in enumerate(peaks, start=1):
  537. markers[tuple(int(v) for v in p)] = i
  538. return markers, dt_max
  539. def _axes_ratio(mask):
  540. props = regionprops(mask.astype(np.uint8))
  541. if not props:
  542. return 1.0
  543. ev = np.sort(np.linalg.eigvalsh(props[0].inertia_tensor))
  544. return float(np.sqrt(ev[-1] / ev[0])) if ev[0] > 0 else 1.0
  545. def _merge_by_rules(parent, ws, dt, dt_max, cfg):
  546. out = ws.copy()
  547. total = int(parent.sum())
  548. changed = True
  549. while changed:
  550. changed = False
  551. ids = [i for i in np.unique(out) if i > 0]
  552. if len(ids) < 2:
  553. break
  554. counts = {i: int((out == i).sum()) for i in ids}
  555. dil = {i: binary_dilation(out == i) for i in ids}
  556. to_merge = set()
  557. for a in range(len(ids)):
  558. for b in range(a + 1, len(ids)):
  559. i, j = ids[a], ids[b]
  560. touch = dil[i] & (out == j)
  561. if not touch.any():
  562. continue
  563. if float(dt[touch].min()) >= cfg["neck_rel"] * dt_max:
  564. to_merge.add((i, j)) # thick boundary -> one nucleus
  565. continue
  566. if any(counts[c] < cfg["min_child_voxels"] or
  567. counts[c] < cfg["min_child_vol_ratio"] * total for c in (i, j)):
  568. to_merge.add((i, j)) # tiny fragment
  569. continue
  570. if (_axes_ratio(out == i) < cfg["elong_merge_if_ratio_below"] and
  571. _axes_ratio(out == j) < cfg["elong_merge_if_ratio_below"]):
  572. to_merge.add((i, j)) # both children round -> false split
  573. for i, j in to_merge:
  574. if (out == i).any() and (out == j).any():
  575. out[out == j] = i
  576. changed = True
  577. final = np.zeros_like(out)
  578. for k, old in enumerate([i for i in np.unique(out) if i > 0], start=1):
  579. final[out == old] = k
  580. return final
  581. def split_touching_shape_aware(masks, voxel_um=VOXEL_UM, cfg=None, progress=True):
  582. """Split labels whose distance transform has >=2 deep maxima separated by a thin neck.
  583. The largest child keeps the parent label; others get new ids appended at the end."""
  584. cfg = cfg or SPLIT_DEFAULT
  585. out = np.asarray(masks).copy()
  586. if out.max() == 0:
  587. return out
  588. bboxes = find_objects(out)
  589. next_id = int(out.max())
  590. it = range(1, len(bboxes) + 1)
  591. it = tqdm(it, desc="shape-aware split", unit="obj", leave=False) if progress else it
  592. for lab in it:
  593. sl = bboxes[lab - 1]
  594. if sl is None:
  595. continue
  596. sub = out[sl] == lab
  597. if not sub.any():
  598. continue
  599. dt = distance_transform_edt(sub, sampling=voxel_um)
  600. markers, dt_max = _peak_markers(sub, dt, voxel_um, cfg)
  601. if markers is None:
  602. continue
  603. ws = watershed(-dt, markers=markers, mask=sub)
  604. if ws.max() <= 1:
  605. continue
  606. ws = _merge_by_rules(sub, ws, dt, dt_max, cfg)
  607. if ws.max() <= 1:
  608. continue
  609. counts = np.bincount(ws.ravel())
  610. order = np.argsort(counts[1:])[::-1] + 1
  611. out[sl][sub] = 0
  612. out[sl][ws == order[0]] = lab
  613. for p in order[1:]:
  614. next_id += 1
  615. out[sl][ws == p] = next_id
  616. return out
  617. # ---- 5b. cleanup [YukimatsuOnSynth2Stardist.cleanup_morphology] ----
  618. def cleanup_labels(masks, voxel_um=VOXEL_UM, open_radius_um=CLEANUP_DEFAULT["open_radius_um"],
  619. min_voxels=CLEANUP_DEFAULT["min_voxels"], progress=True):
  620. """Per label: largest connected component -> hole fill -> ellipsoid opening
  621. (skipped when ``open_radius_um`` is 0/None) -> largest CC -> drop if < ``min_voxels``.
  622. Returns (uint32 volume, stats dict)."""
  623. masks = np.asarray(masks)
  624. se = None
  625. rz = ry = rx = 0
  626. if open_radius_um:
  627. rz, ry, rx = _radii_vox(open_radius_um, voxel_um)
  628. se = ellipsoid_se(rz, ry, rx)
  629. pad = max(rz, ry, rx) + 1
  630. Z, Y, X = masks.shape
  631. out = np.zeros(masks.shape, dtype=np.uint32)
  632. kept = dropped = 0
  633. bboxes = find_objects(masks)
  634. it = enumerate(bboxes)
  635. it = tqdm(it, total=len(bboxes), desc="cleanup", unit="obj", leave=False) if progress else it
  636. for idx, sl in it:
  637. if sl is None:
  638. continue
  639. lab = idx + 1
  640. slc = (slice(max(0, sl[0].start - pad), min(Z, sl[0].stop + pad)),
  641. slice(max(0, sl[1].start - pad), min(Y, sl[1].stop + pad)),
  642. slice(max(0, sl[2].start - pad), min(X, sl[2].stop + pad)))
  643. m = masks[slc] == lab
  644. if not m.any():
  645. continue
  646. m = binary_fill_holes(largest_cc(m))
  647. if se is not None:
  648. m = largest_cc(binary_dilation(binary_erosion(m, structure=se), structure=se))
  649. if m.sum() < min_voxels:
  650. dropped += 1
  651. continue
  652. write = m & (out[slc] == 0)
  653. out[slc][write] = lab
  654. kept += 1
  655. return out, dict(kept=kept, dropped=dropped, open_radius_um=open_radius_um,
  656. min_voxels=min_voxels)
  657. # ---- 5c. 2D -> 3D linking (hybrid) [HybridSegment.build_hybrid] ----
  658. def link_2d_to_3d(masks2d, masks3d, overlap_thr=HYBRID_OVERLAP_THR, progress=True):
  659. """Hybrid volume: each independent 2D blob takes the label of the 3D object covering
  660. the largest share of its pixels, if that share is >= ``overlap_thr``. Pixel boundaries
  661. come from the 2D masks, z-linking from the 3D masks. Returns (uint32 volume, stats).
  662. Not part of the current pipeline -- an alternative to 5d, claiming in the opposite
  663. direction. Reachable via ``--hybrid``."""
  664. if masks2d.shape != masks3d.shape:
  665. raise ValueError("shape mismatch")
  666. out = np.zeros(masks3d.shape, dtype=np.uint32)
  667. st = dict(blobs=0, assigned=0, no_overlap=0, low_overlap=0)
  668. Z = masks2d.shape[0]
  669. it = tqdm(range(Z), desc="link 2d->3d", unit="slice", leave=False) if progress else range(Z)
  670. for z in it:
  671. s2, s3 = masks2d[z], masks3d[z]
  672. for idx, sl in enumerate(find_objects(s2)):
  673. if sl is None:
  674. continue
  675. m = s2[sl] == idx + 1
  676. if not m.any():
  677. continue
  678. st["blobs"] += 1
  679. under = s3[sl][m]
  680. under = under[under > 0]
  681. if under.size == 0:
  682. st["no_overlap"] += 1
  683. continue
  684. counts = np.bincount(under)
  685. best = int(counts.argmax())
  686. if counts[best] / int(m.sum()) < overlap_thr:
  687. st["low_overlap"] += 1
  688. continue
  689. out[z][sl][m] = best
  690. st["assigned"] += 1
  691. st["labels_out"] = int(np.unique(out).size - 1)
  692. return out, st
  693. # ---- 5d. 3D -> 2D anchoring (back-matching) [AnchorTo2D.py] ----
  694. def anchor_to_2d(masks3d, masks2d, min_voxels=CLEANUP_ANCHOR["min_voxels"],
  695. pad=ANCHOR_PAD_PX, progress=True):
  696. """Back-match a 3D label volume to the 2D masks, slice by slice.
  697. For every 3D label and z-slice: if the label claims a 2D blob there (max-overlap claim,
  698. conflicts between 3D labels resolved by overlap count), the slice becomes exactly that
  699. 2D blob; if the 2D segmentation saw nothing (caps, gap slices) the 3D slice is kept,
  700. minus any 2D blobs claimed by other labels. Label ids are preserved.
  701. Returns (uint32 volume, stats)."""
  702. if masks2d.shape != masks3d.shape:
  703. raise ValueError("shape mismatch")
  704. objs = find_objects(masks3d)
  705. labels = [i + 1 for i, sl in enumerate(objs) if sl is not None]
  706. Z, Y, X = masks3d.shape
  707. claims = {} # (z, label2d) -> (count, label3d)
  708. it = tqdm(labels, desc="anchor: claims", unit="obj", leave=False) if progress else labels
  709. for lab in it:
  710. sl = objs[lab - 1]
  711. for z in range(sl[0].start, sl[0].stop):
  712. m3 = masks3d[z, sl[1], sl[2]] == lab
  713. if not m3.any():
  714. continue
  715. under = masks2d[z, sl[1], sl[2]][m3]
  716. under = under[under != 0]
  717. if under.size == 0:
  718. continue
  719. l2, c = np.unique(under, return_counts=True)
  720. for a, b in zip(l2, c):
  721. key = (z, int(a))
  722. if key not in claims or b > claims[key][0]:
  723. claims[key] = (int(b), lab)
  724. out = np.zeros(masks3d.shape, dtype=np.uint32)
  725. st = dict(kept=0, dropped=0, slices_anchored=0, slices_kept3d=0, claims=len(claims))
  726. it = tqdm(labels, desc="anchor: rebuild", unit="obj", leave=False) if progress else labels
  727. for lab in it:
  728. sl = objs[lab - 1]
  729. y0, y1 = max(0, sl[1].start - pad), min(Y, sl[1].stop + pad)
  730. x0, x1 = max(0, sl[2].start - pad), min(X, sl[2].stop + pad)
  731. result = {}
  732. vox = 0
  733. for z in range(sl[0].start, sl[0].stop):
  734. m3 = masks3d[z, y0:y1, x0:x1] == lab
  735. if not m3.any():
  736. continue
  737. cp = masks2d[z, y0:y1, x0:x1]
  738. under = np.unique(cp[m3])
  739. under = under[under != 0]
  740. mine = [int(l) for l in under if claims.get((z, int(l)), (0, -1))[1] == lab]
  741. foreign = [int(l) for l in under if claims.get((z, int(l)), (0, -1))[1] != lab]
  742. if mine:
  743. res = np.isin(cp, mine)
  744. st["slices_anchored"] += 1
  745. else:
  746. res = m3.copy()
  747. if foreign:
  748. res &= ~np.isin(cp, foreign)
  749. st["slices_kept3d"] += 1
  750. if res.any():
  751. result[z] = res
  752. vox += int(res.sum())
  753. if vox < min_voxels:
  754. st["dropped"] += 1
  755. continue
  756. for z, res in result.items():
  757. tgt = out[z, y0:y1, x0:x1]
  758. tgt[res & (tgt == 0)] = lab
  759. st["kept"] += 1
  760. return out, st
  761. def anchor_and_finish(masks3d, masks2d, voxel_um=VOXEL_UM, cleanup=CLEANUP_ANCHOR,
  762. rounding=ROUNDING_POST_ANCHOR, use_gpu=True, progress=True):
  763. """Anchoring followed by two finishing steps: an ellipsoid opening that removes shelves
  764. one to three slices thick, and a rounding pass (``ROUNDING_POST_ANCHOR``) that closes the
  765. residual slice steps."""
  766. out, st = anchor_to_2d(masks3d, masks2d, min_voxels=cleanup["min_voxels"], progress=progress)
  767. out, st_clean = cleanup_labels(out, voxel_um, progress=progress, **cleanup)
  768. if rounding is not None:
  769. out = round_caps(out, voxel_um, rounding, use_gpu=use_gpu, progress=progress)
  770. st.update({f"cleanup_{k}": v for k, v in st_clean.items()})
  771. return out.astype(np.uint32, copy=False), st
  772. # ---- 5e. convenience chain ----
  773. def refine(masks, voxel_um=VOXEL_UM, rounding=ROUNDING_DEFAULT, split=None,
  774. cleanup=CLEANUP_DEFAULT, anchor_to=None, anchor_cleanup=CLEANUP_ANCHOR,
  775. anchor_rounding=ROUNDING_POST_ANCHOR, use_gpu=True, progress=True):
  776. """compact labels -> round caps (``rounding``; None = skip) -> shape-aware split
  777. (``split`` cfg; None = skip) -> cleanup -> optional anchoring to 2D masks
  778. (``anchor_to`` = independent 2D label volume). Returns (uint32 volume, stats)."""
  779. vol = relabel_compact_1n(masks)
  780. stats = dict(labels_in=int(vol.max()))
  781. if rounding is not None:
  782. vol = round_caps(vol, voxel_um, rounding, use_gpu=use_gpu, progress=progress)
  783. if split is not None:
  784. vol = split_touching_shape_aware(vol, voxel_um, split, progress=progress)
  785. stats["labels_after_split"] = int(np.unique(vol).size - 1)
  786. vol, st = cleanup_labels(vol, voxel_um, progress=progress, **cleanup)
  787. stats.update({f"cleanup_{k}": v for k, v in st.items()})
  788. if anchor_to is not None:
  789. vol, st = anchor_and_finish(vol, anchor_to, voxel_um, anchor_cleanup, anchor_rounding,
  790. use_gpu=use_gpu, progress=progress)
  791. stats.update({f"anchor_{k}": v for k, v in st.items()})
  792. stats["labels_out"] = int(np.unique(vol).size - 1)
  793. return vol.astype(np.uint32, copy=False), stats
  794. # ============================================================================
  795. # 6. SHAPE METRICS -- per-nucleus 3D and per-slice 2D
  796. # [Compare3DvsMean2D.py metrics_3D / metrics_2D_sameMask / metrics_2D_indep;
  797. # SliceSampling3Dvs2D.build_per_slice]
  798. #
  799. # Naming follows the committed per_nucleus.csv: ``R3D`` = Wadell sphericity psi
  800. # (pi^(1/3) (6V)^(2/3) / A), ``S3D`` = min/max bounding-box extent;
  801. # ``roundness_2d`` = 4 pi A / P^2 (circularity), ``sphericity_2d`` = min/max 2D bbox.
  802. # ============================================================================
  803. def _require_pandas():
  804. if pd is None:
  805. raise RuntimeError("the shape-metric tables need pandas; pip install pandas")
  806. def metrics_3d(masks, voxel_um=VOXEL_UM, progress=True):
  807. """dict label -> {S3D, R3D, V_um3, A_um2, n_slices, z0, z1, y0, y1, x0, x1}."""
  808. zsp, ysp, xsp = voxel_um
  809. bboxes = find_objects(masks)
  810. out = {}
  811. it = enumerate(bboxes)
  812. it = tqdm(it, total=len(bboxes), desc="metrics_3d", unit="obj", leave=False) if progress else it
  813. for idx, sl in it:
  814. if sl is None:
  815. continue
  816. lab = idx + 1
  817. sub = masks[sl] == lab
  818. if not sub.any():
  819. continue
  820. ext = ((sl[0].stop - sl[0].start) * zsp, (sl[1].stop - sl[1].start) * ysp,
  821. (sl[2].stop - sl[2].start) * xsp)
  822. V = float(sub.sum() * zsp * ysp * xsp)
  823. try:
  824. verts, faces, _, _ = marching_cubes(np.pad(sub.astype(np.uint8), 1), level=0.5,
  825. spacing=voxel_um)
  826. A = float(mesh_surface_area(verts, faces))
  827. R3D = float(math.pi ** (1 / 3) * (6 * V) ** (2 / 3) / A) if A > 0 else float("nan")
  828. except (RuntimeError, ValueError):
  829. A, R3D = float("nan"), float("nan")
  830. out[lab] = dict(S3D=float(min(ext) / max(ext)), R3D=R3D, V_um3=V, A_um2=A,
  831. n_slices=int((sub.sum(axis=(1, 2)) > 0).sum()),
  832. z0=sl[0].start, z1=sl[0].stop, y0=sl[1].start, y1=sl[1].stop,
  833. x0=sl[2].start, x1=sl[2].stop)
  834. return out
  835. def blob_metrics_2d(blob):
  836. """(sphericity_2d, roundness_2d) of one boolean 2D blob; pixel units cancel."""
  837. if not blob.any():
  838. return float("nan"), float("nan")
  839. yy, xx = np.where(blob)
  840. by, bx = yy.max() - yy.min() + 1, xx.max() - xx.min() + 1
  841. sph = float(min(by, bx) / max(by, bx))
  842. props = regionprops(blob.astype(np.uint8))
  843. if not props or props[0].perimeter == 0:
  844. return sph, float("nan")
  845. return sph, float(4.0 * math.pi * props[0].area / props[0].perimeter ** 2)
  846. def metrics_2d_same_mask(masks, meta3d, progress=True):
  847. """Mean over the z cross-sections of the SAME 3D label (mode A)."""
  848. out = {}
  849. it = meta3d.items()
  850. it = tqdm(it, desc="metrics_2d_same_mask", unit="obj", leave=False) if progress else it
  851. for lab, m in it:
  852. sphs, rnds = [], []
  853. for z in range(m["z0"], m["z1"]):
  854. blob = masks[z, m["y0"]:m["y1"], m["x0"]:m["x1"]] == lab
  855. if not blob.any():
  856. continue
  857. s, r = blob_metrics_2d(blob)
  858. if not math.isnan(s):
  859. sphs.append(s)
  860. if not math.isnan(r):
  861. rnds.append(r)
  862. out[lab] = dict(S_2Dmean_sameMask=float(np.mean(sphs)) if sphs else float("nan"),
  863. R_2Dmean_sameMask=float(np.mean(rnds)) if rnds else float("nan"),
  864. n_used_sameMask=len(sphs))
  865. return out
  866. def per_slice_bboxes(masks2d):
  867. return [find_objects(masks2d[z]) for z in range(masks2d.shape[0])]
  868. def metrics_2d_independent(masks3d, masks2d, meta3d, bboxes2d=None, progress=True):
  869. """Mean over the best-overlap INDEPENDENT 2D blob in each slice of the label (mode B)."""
  870. bboxes2d = bboxes2d or per_slice_bboxes(masks2d)
  871. out = {}
  872. it = meta3d.items()
  873. it = tqdm(it, desc="metrics_2d_independent", unit="obj", leave=False) if progress else it
  874. for lab, m in it:
  875. sphs, rnds, n = [], [], 0
  876. for z in range(m["z0"], m["z1"]):
  877. this = masks3d[z, m["y0"]:m["y1"], m["x0"]:m["x1"]] == lab
  878. if not this.any():
  879. continue
  880. ids, counts = np.unique(masks2d[z, m["y0"]:m["y1"], m["x0"]:m["x1"]][this],
  881. return_counts=True)
  882. keep = ids != 0
  883. ids, counts = ids[keep], counts[keep]
  884. if ids.size == 0:
  885. continue
  886. best = int(ids[np.argmax(counts)])
  887. if best - 1 >= len(bboxes2d[z]) or bboxes2d[z][best - 1] is None:
  888. continue
  889. blob = masks2d[z][bboxes2d[z][best - 1]] == best
  890. if not blob.any():
  891. continue
  892. s, r = blob_metrics_2d(blob)
  893. if not math.isnan(s):
  894. sphs.append(s)
  895. n += 1
  896. if not math.isnan(r):
  897. rnds.append(r)
  898. out[lab] = dict(S_2Dmean_indep=float(np.mean(sphs)) if sphs else float("nan"),
  899. R_2Dmean_indep=float(np.mean(rnds)) if rnds else float("nan"),
  900. n_indep_matches=n)
  901. return out
  902. def per_nucleus_table(masks3d, masks2d=None, voxel_um=VOXEL_UM, progress=True):
  903. """DataFrame with one row per 3D label: 3D metrics + 2D-mean same-mask metrics and,
  904. when ``masks2d`` is given, 2D-mean independent metrics."""
  905. _require_pandas()
  906. meta = metrics_3d(masks3d, voxel_um, progress)
  907. same = metrics_2d_same_mask(masks3d, meta, progress)
  908. indep = metrics_2d_independent(masks3d, masks2d, meta, progress=progress) \
  909. if masks2d is not None else {}
  910. rows = []
  911. for lab in sorted(meta):
  912. row = dict(label=lab, **{k: meta[lab][k] for k in ("n_slices", "S3D", "R3D", "V_um3", "A_um2")})
  913. row.update(same[lab])
  914. if indep:
  915. row.update(indep[lab])
  916. rows.append(row)
  917. return pd.DataFrame(rows)
  918. def per_slice_table(masks3d, voxel_um=VOXEL_UM, progress=True):
  919. """One row per (label, z) cross-section with the full 2D descriptor set: circularity,
  920. bbox sphericity, ellipse aspect, eccentricity, solidity, convexity, Zheng-Hryciw
  921. corrected roundness and centroid-radii sphericity (r_ins / r_circ)."""
  922. _require_pandas()
  923. xy = voxel_um[1]
  924. bboxes = find_objects(masks3d)
  925. rows = []
  926. it = enumerate(bboxes)
  927. it = tqdm(it, total=len(bboxes), desc="per_slice_table", unit="obj", leave=False) if progress else it
  928. for idx, sl in it:
  929. if sl is None:
  930. continue
  931. lab = idx + 1
  932. for z in range(sl[0].start, sl[0].stop):
  933. blob = masks3d[z, sl[1], sl[2]] == lab
  934. if not blob.any():
  935. continue
  936. s, r = blob_metrics_2d(blob)
  937. p = regionprops(blob.astype(np.uint8))[0]
  938. area, perim = int(p.area), float(p.perimeter)
  939. maj = float(p.axis_major_length)
  940. aspect = float(p.axis_minor_length) / maj if maj > 0 else float("nan")
  941. yy, xx = np.nonzero(blob)
  942. dy, dx = int(yy.max() - yy.min()), int(xx.max() - xx.min())
  943. try:
  944. hull = ConvexHull(np.column_stack((yy, xx)))
  945. convex = perim / hull.area if hull.area > 0 else float("nan") # 2D: hull.area = perimeter
  946. except (QhullError, ValueError):
  947. convex = float("nan")
  948. if math.isnan(aspect) or math.isnan(r):
  949. r_zh = float("nan")
  950. else:
  951. r_zh = r + (0.913 - (-0.528 * aspect * aspect + 1.259 * aspect + 0.270))
  952. by_, bx_ = np.nonzero(blob & ~binary_erosion(blob))
  953. d = np.hypot(by_ - yy.mean(), bx_ - xx.mean())
  954. sph_r = float(d.min() / d.max()) if d.size and d.max() > 0 else float("nan")
  955. rows.append(dict(label=lab, z_index=z, area_px=area, area_um2=area * xy * xy,
  956. perimeter_um=perim * xy, roundness_2d=r, sphericity_2d=s,
  957. aspect_2d=aspect, elong_2d=float(p.eccentricity),
  958. eccbox_2d=dy / dx if dx > 0 else float("nan"),
  959. convex_2d=convex, solid_2d=float(p.solidity),
  960. roundness_zh_2d=r_zh, sphericity_radii_2d=sph_r))
  961. return pd.DataFrame(rows).sort_values(["label", "z_index"]).reset_index(drop=True)
  962. # ============================================================================
  963. # 7. CLI
  964. # ============================================================================
  965. SYNTH_TAG = {"v1": "synth1", "v2": "synth2", "v3": "synth3"}
  966. def _parse(argv):
  967. p = argparse.ArgumentParser(
  968. prog=os.path.basename(__file__), description=__doc__,
  969. formatter_class=argparse.RawDescriptionHelpFormatter)
  970. p.add_argument("--raw-dir", required=True,
  971. help="folder of per-z tiffs (Zeiss ..._b0z<N>c0...tiff) or a single 3D tiff")
  972. p.add_argument("--out-dir", required=True)
  973. p.add_argument("--cellpose-weights", help="custom cpsam model; omit if M_indep.tiff exists")
  974. p.add_argument("--cellpose-mode", choices=list(CELLPOSE_MODES), default="independent",
  975. help="stage 1 front-end (default: independent = current pipeline)")
  976. p.add_argument("--stardist-model-dir", help="StarDist basedir; omit if M_stardist_<v>.tiff exists")
  977. p.add_argument("--stardist-python", help="python.exe of the stardist-gpu env; "
  978. "omit to run stage 3 in this interpreter")
  979. p.add_argument("--synth", choices=["v1", "v2", "v3"], default="v3")
  980. p.add_argument("--variant", choices=list(VARIANTS), default="reduced")
  981. p.add_argument("--split", action="store_true", help="also run the shape-aware split")
  982. p.add_argument("--anchor", action="store_true", help="stage 5: back-match to the 2D masks")
  983. p.add_argument("--hybrid", action="store_true",
  984. help="additionally write the HybridSegment 2D->3D linking of the final masks")
  985. p.add_argument("--metrics", action="store_true", help="stage 6: per_nucleus.csv + per_slice.csv")
  986. p.add_argument("--mips", action="store_true",
  987. help="write XY/XZ/YZ projections of the final masks (cast to uint16, so "
  988. "label ids above 65535 wrap -- for visual QC only)")
  989. p.add_argument("--voxel", default=",".join(str(v) for v in VOXEL_UM), help="Z,Y,X in um")
  990. p.add_argument("--cpu", action="store_true", help="no GPU for cellpose and cap rounding")
  991. p.add_argument("--no-progress", action="store_true")
  992. p.add_argument("--force", action="store_true", help="recompute every stage")
  993. return p.parse_args(argv)
  994. def _stardist_stage(argv):
  995. """Stage 3 in its own interpreter: this file re-invoked with the stardist-gpu python."""
  996. p = argparse.ArgumentParser(prog=f"{os.path.basename(__file__)} --stardist-stage")
  997. p.add_argument("--stardist-stage", action="store_true")
  998. p.add_argument("--input", required=True)
  999. p.add_argument("--output", required=True)
  1000. p.add_argument("--model-dir", required=True)
  1001. p.add_argument("--model", default=STARDIST["model_name"])
  1002. p.add_argument("--tiles", default=",".join(str(t) for t in STARDIST["n_tiles"]))
  1003. p.add_argument("--prob", type=float, default=None)
  1004. p.add_argument("--nms", type=float, default=None)
  1005. a = p.parse_args(argv)
  1006. labels = run_stardist3d(imread_3d(a.input), a.model_dir, a.model,
  1007. tuple(int(t) for t in a.tiles.split(",")), a.prob, a.nms)
  1008. tifffile.imwrite(a.output, labels)
  1009. print(f"segmented {int(labels.max())} nuclei -> {a.output}")
  1010. return 0
  1011. def main(argv=None):
  1012. argv = sys.argv[1:] if argv is None else argv
  1013. if "--stardist-stage" in argv:
  1014. return _stardist_stage(argv)
  1015. a = _parse(argv)
  1016. voxel = tuple(float(v) for v in a.voxel.split(","))
  1017. if len(voxel) != 3:
  1018. raise SystemExit("--voxel needs three comma-separated values (Z,Y,X)")
  1019. progress = not a.no_progress
  1020. os.makedirs(a.out_dir, exist_ok=True)
  1021. P = lambda n: os.path.join(a.out_dir, n)
  1022. fresh = lambda path: a.force or not os.path.exists(path)
  1023. raw = load_raw(a.raw_dir)
  1024. print(f"raw {raw.shape} {raw.dtype}")
  1025. if not os.path.exists(P("raw_dapi.tiff")):
  1026. tifffile.imwrite(P("raw_dapi.tiff"), raw)
  1027. # --- 1. Cellpose 2D, independent -------------------------------------------------
  1028. if fresh(P("M_indep.tiff")):
  1029. if not a.cellpose_weights:
  1030. raise SystemExit("M_indep.tiff missing: pass --cellpose-weights")
  1031. kw = dict(voxel_um=voxel) if a.cellpose_mode == "3d" else {}
  1032. M_indep = CELLPOSE_MODES[a.cellpose_mode](raw, a.cellpose_weights, gpu=not a.cpu, **kw)
  1033. save_labels(P("M_indep.tiff"), M_indep)
  1034. else:
  1035. M_indep = tifffile.imread(P("M_indep.tiff"))
  1036. print(f"M_indep: {n_labels(M_indep)} 2D blobs ({a.cellpose_mode})")
  1037. # --- 2. synthetic DAPI -----------------------------------------------------------
  1038. synth_path = P(f"synth_dapi_{a.synth}.tiff")
  1039. if fresh(synth_path):
  1040. tifffile.imwrite(synth_path, build_synth_dapi(raw, M_indep, a.synth))
  1041. print(f"synth: {synth_path}")
  1042. # --- 3. StarDist 3D --------------------------------------------------------------
  1043. sd_path = P(f"M_stardist_{a.synth}.tiff")
  1044. if fresh(sd_path):
  1045. if not a.stardist_model_dir:
  1046. raise SystemExit(f"{os.path.basename(sd_path)} missing: pass --stardist-model-dir")
  1047. if a.stardist_python:
  1048. subprocess.check_call([a.stardist_python, os.path.abspath(__file__), "--stardist-stage",
  1049. "--input", synth_path, "--output", sd_path,
  1050. "--model-dir", a.stardist_model_dir,
  1051. "--model", STARDIST["model_name"]])
  1052. else:
  1053. labels = run_stardist3d(imread_3d(synth_path), a.stardist_model_dir)
  1054. tifffile.imwrite(sd_path, labels)
  1055. M_sd = tifffile.imread(sd_path)
  1056. print(f"stardist: {n_labels(M_sd)} labels")
  1057. # --- 4 + 5. Yukimatsu variant (+ optional split, + optional anchoring) -----------
  1058. rounding, cleanup = VARIANTS[a.variant]
  1059. tag = (f"M_{SYNTH_TAG[a.synth]}_yuki_{a.variant}"
  1060. + ("_split" if a.split else "") + ("_anchored" if a.anchor else ""))
  1061. out_path = P(tag + ".tiff")
  1062. if fresh(out_path):
  1063. vol, stats = refine(M_sd, voxel, rounding=rounding,
  1064. split=SPLIT_DEFAULT if a.split else None,
  1065. cleanup=cleanup, anchor_to=M_indep if a.anchor else None,
  1066. use_gpu=not a.cpu, progress=progress)
  1067. save_labels(out_path, vol)
  1068. stats.update(synth=a.synth, variant=a.variant, split=a.split, anchor=a.anchor,
  1069. cellpose_mode=a.cellpose_mode, rounding=rounding, cleanup=cleanup,
  1070. voxel_um=voxel)
  1071. with open(P(tag + "_stats.json"), "w") as f:
  1072. json.dump(stats, f, indent=2, default=str)
  1073. print(json.dumps(stats, indent=2, default=str))
  1074. M_out = tifffile.imread(out_path)
  1075. print(f"final: {out_path} ({n_labels(M_out)} nuclei)")
  1076. if a.mips:
  1077. save_mips(P(tag), M_out)
  1078. # --- optional: HybridSegment 2D -> 3D linking of the final volume ----------------
  1079. if a.hybrid:
  1080. hyb_path = P(tag + "_hybrid.tiff")
  1081. if fresh(hyb_path):
  1082. hyb, st = link_2d_to_3d(M_indep, M_out, progress=progress)
  1083. save_labels(hyb_path, hyb)
  1084. print(f"hybrid: {hyb_path} {json.dumps(st)}")
  1085. # --- 6. shape metrics ------------------------------------------------------------
  1086. if a.metrics:
  1087. pn_path, ps_path = P("per_nucleus.csv"), P("per_slice.csv")
  1088. if fresh(pn_path):
  1089. per_nucleus_table(M_out, M_indep, voxel, progress).to_csv(pn_path, index=False)
  1090. print(f"metrics: {pn_path}")
  1091. if fresh(ps_path):
  1092. per_slice_table(M_out, voxel, progress).to_csv(ps_path, index=False)
  1093. print(f"metrics: {ps_path}")
  1094. return 0
  1095. if __name__ == "__main__":
  1096. sys.exit(main())

cbhs3d_onefile.py at commit 063985e, no license · at the source

Overview

  1. Max‐Planck‐Zentrum für Physik und Medizin Erlangen Germany
  2. Max Planck Institute for the Science of Light Erlangen Germany
  3. Faculty of Medicine, Laboratory of Neural Epigenomics, Medical Institute of Biophysics Friedrich‐Alexander Universität Erlangen‐Nürnberg Erlangen Germany
  4. Nuclear Architecture in Neural Plasticity and Aging German Center for Neurodegenerative Diseases (DZNE) Dresden Germany
  5. Department of Physics, Institute of Optics, Information and Photonics FAU Erlangen Germany
  6. Department of Stem Cell Biology University Hospital Erlangen, Friedrich‐Alexander‐Universität Erlangen‐Nürnberg Erlangen Germany
  7. Department of Physics “A. Pontremoli” University of Milan and INFN Milan Italy
Journal: Hippocampus, volume 36, issue 6, article e70139
Dates: received 2 July 2026; accepted 11 September 2026; published online 26 October 2026; in print November 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/hipo.70139 · PMID 42934120 · PMCID PMC13710209 · OpenAlex W7223684705
Open access: hybrid, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: mouse (organism), cellular / molecular (subfield)
Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Preprocessing, Machine learning, fMRI & imaging
Keywords: cellular packing order, hippocampus, quantitative image analysis, topology
MeSH: Hippocampus*, Image Processing, Computer-Assisted*, Neurons*, Animals, Astrocytes, Cell Nucleus, Male, Mice, Mice, Inbred C57BL, Software (* major topic)
Topic: Cell Image Analysis Techniques (Biophysics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (470322152 ‐ TO1347/3‐1, 497658532 ‐ TO1347/4‐1, 507965872 ‐ TO1347/5‐1, 563413271 ‐ TO1347/6‐1, 460333672 (CRC1540 EBM), 455449456); Schram Stifung; European Research Council (ERC‐2018‐STG, 804468 EAGER, ERC‐2023‐COG, 101125034 NEUTIME); Interdisciplinary Center for Clinical Research (IZKF) Erlangen (P162)
Citations: not cited yet (Europe PMC); 32 references in the paper

Abstract

The hippocampus is characterized by a stereotypical macroscopic structure, where the nuclei are densely and heterogeneously packed among different subregions of the hippocampus. Despite the fact that tissue‐specific cellular organization has been implicated in neural function, it has been technically challenging to quantitatively analyze mesoscopic cellular organization in the hippocampus due to its high cellular density. To overcome this technical hurdle, we developed Computational Biophysical Histomorphometry Software (CBHS), an automated image‐analysis pipeline aimed at quantifying nuclear shape and the order of the cellular ensemble in high‐density areas. When applied to the subfields of the hippocampus, we found that denser regions, most notably the dentate gyrus, were the most positionally, but least orientationally ordered. Nuclear shape exhibited a dependence on the local environment in a packing‐dependent manner. This association was cell‐type specific, with neurons, but not astrocytes displaying nuclear shape that varied with neighbor proximity, although astrocytes demonstrated greater intrinsic shape variance. The results reveal the presence of reproducible mesoscale cell packing order in hippocampal tissue, and are consistent with a nucleus‐driven mechanical coupling between neighboring cells. The present study provides a quantitative framework with which to understand mesoscopic tissue organization, thus enabling the formulation of testable hypotheses for future investigation.

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

Repository

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

gitlab.mpcdf.mpg.de/konmi/image-analysis-ui

License: none: the authors keep all their rights
State: the link answers, verified on 1 November 2026
Evidence: files inventoried
Commit: 063985e1c2633ad5ccde4bfe29717844fb44e1e2, 18 September 2026
Languages: Python (7)
Size: 23 files, 7 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, environment (requirements.txt), continuous integration
Not found: license file, CITATION.cff, tests, documentation
Tools: NumPy (6 files), pandas (5 files), tifffile (4 files), scikit-image (3 files), SciPy (3 files), Cellpose (2 files), Matplotlib (2 files), CuPy (1 file), OpenCV (1 file), PyTorch (1 file)
Availability: 1 check, the latest on 1 November 2026: the link answers
  • 1 November 2026: the link answers
8 files

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

Tracing map

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

What the map holds:

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

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

Data

No dataset and no data link were found in the paper.

Data Availability Statement

Images can be made available upon request to the leading author. All original code has been deposited at Gitlab under the permalink https://gitlab.mpcdf.mpg.de/konmi/image‐analysis‐ui/‐/tree/e1eea9740df5e9c3e5e0a185510e0f82145f9623/ (https://gitlab.mpcdf.mpg.de/konmi/image-analysis-ui/-/tree/e1eea9740df5e9c3e5e0a185510e0f82145f9623/) and is publicly available as of the date of publication. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

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

Versions

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

Version 1, 1 November 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 4 keywords, 10 MeSH terms, 4 funders, 28 references.

Cite

This paper

Hein, K., Romero‐Limon, H., Möckel, C., Karasinsky, A., Möllmert, S., Kayser, J., Zaccone, A., Guck, J., & Toda, T. (2026). Quantitative Assessment of Mesoscale Cellular Order and Organization in the Mouse Hippocampus. Hippocampus, 36(6), e70139. https://doi.org/10.1002/hipo.70139

BibTeX

@article{hein2026quantitative,
author = {Hein, Konstantin and Romero‐Limon, Humberto and Möckel, Conrad and Karasinsky, Anne and Möllmert, Stephanie and Kayser, Jona and Zaccone, Alessio and Guck, Jochen and Toda, Tomohisa},
title = {{Quantitative Assessment of Mesoscale Cellular Order and Organization in the Mouse Hippocampus}},
journal = {Hippocampus},
year = {2026},
month = nov,
volume = {36},
number = {6},
pages = {e70139},
publisher = {Wiley},
issn = {1050-9631},
doi = {10.1002/hipo.70139},
url = {https://doi.org/10.1002/hipo.70139},
pmid = {42934120},
pmcid = {PMC13710209}
}

RIS

TY - JOUR
AU - Hein, Konstantin
AU - Romero‐Limon, Humberto
AU - Möckel, Conrad
AU - Karasinsky, Anne
AU - Möllmert, Stephanie
AU - Kayser, Jona
AU - Zaccone, Alessio
AU - Guck, Jochen
AU - Toda, Tomohisa
TI - Quantitative Assessment of Mesoscale Cellular Order and Organization in the Mouse Hippocampus
T2 - Hippocampus
J2 - Hippocampus
PY - 2026
DA - 2026/11/01
VL - 36
IS - 6
SP - e70139
SN - 1050-9631
PB - Wiley
DO - 10.1002/hipo.70139
UR - https://doi.org/10.1002/hipo.70139
LA - en
ER -

CSL-JSON

{
"id": "10.1002/hipo.70139",
"type": "article-journal",
"title": "Quantitative Assessment of Mesoscale Cellular Order and Organization in the Mouse Hippocampus",
"container-title": "Hippocampus",
"author": [
{
"family": "Hein",
"given": "Konstantin"
},
{
"family": "Romero‐Limon",
"given": "Humberto"
},
{
"family": "Möckel",
"given": "Conrad"
},
{
"family": "Karasinsky",
"given": "Anne"
},
{
"family": "Möllmert",
"given": "Stephanie"
},
{
"family": "Kayser",
"given": "Jona"
},
{
"family": "Zaccone",
"given": "Alessio"
},
{
"family": "Guck",
"given": "Jochen"
},
{
"family": "Toda",
"given": "Tomohisa"
}
],
"container-title-short": "Hippocampus",
"volume": "36",
"issue": "6",
"page": "e70139",
"DOI": "10.1002/hipo.70139",
"PMID": "42934120",
"PMCID": "PMC13710209",
"ISSN": "1050-9631",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/hipo.70139",
"language": "en",
"issued": {
"date-parts": [
[
2026,
11,
1
]
]
}
}

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41586-025-09151-3 [code]
Morphodynamics of human early brain organoid development.
Journal: Nature
In common: Cellpose, CuPy, tifffile, 7 other tools, 1 reference
[2] doi:10.1113/ep093267 [code]
Spontaneous Ca&lt;sup&gt;2+&lt;/sup&gt; signals in the developing mammalian cochlea of live mice under different anaesthetic regimes.
Journal: Experimental physiology
In common: Cellpose, CuPy, tifffile, 6 other tools, mouse, cellular / molecular, 1 reference
[3] doi:10.1038/s41467-024-55519-w [code]
In vivo spontaneous Ca&lt;sup&gt;2+&lt;/sup&gt; activity in the pre-hearing mammalian cochlea.
Journal: Nature communications
In common: Cellpose, CuPy, tifffile, 6 other tools, mouse, cellular / molecular, 1 reference
[4] doi:10.1371/journal.pcbi.1014571 [code]
SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: Cellpose, tifffile, OpenCV, 6 other tools, 2 references
[5] doi:10.1038/s41587-022-01455-3 [code]
Scalable in situ single-cell profiling by electrophoretic capture of mRNA using EEL FISH.
Journal: Nature biotechnology
In common: Cellpose, tifffile, OpenCV, 6 other tools, mouse, cellular / molecular, 1 reference
[6] doi:10.1126/science.adq2084 [code]
Deep-tissue transcriptomics and subcellular imaging at high spatial resolution.
Journal: Science (New York, N.Y.)
In common: Cellpose, tifffile, OpenCV, 6 other tools, mouse, 1 reference
[7] doi:10.1016/j.cell.2022.12.010 [code]
Molecular and spatial signatures of mouse brain aging at single-cell resolution.
Journal: Cell
In common: CuPy, tifffile, OpenCV, 6 other tools, mouse, cellular / molecular, 1 reference
[8] doi:10.1016/j.isci.2024.109009 [code]
Spatial and single-nucleus transcriptomics decoding the molecular landscape and cellular organization of avian optic tectum.
Journal: iScience
In common: Cellpose, CuPy, OpenCV, 6 other tools, cellular / molecular
[9] doi:10.1073/pnas.2609132123 [code]
A human lysosomal storage disorder toolkit for decoding proteome landscapes in cortical-like and dopaminergic-like induced neurons.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: Cellpose, tifffile, scikit-image, 5 other tools, cellular / molecular, 1 reference
[10] doi:10.1186/s40478-025-02048-2 [code]
Copper supplementation mitigates Parkinson-like wild-type SOD1 pathology and nigrostriatal degeneration in a novel mouse model.
Journal: Acta neuropathologica communications
In common: Cellpose, tifffile, OpenCV, 6 other tools, mouse

Contribute

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

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

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1002/hipo.70139.svg)](https://openscicode.org/badge/doi/10.1002/hipo.70139)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1002/hipo.70139.svg
   :target: https://openscicode.org/badge/doi/10.1002/hipo.70139
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1002/hipo.70139"><img src="https://openscicode.org/badge/doi/10.1002/hipo.70139.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

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

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

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