Dev radii

.gitignore

@@ -37,4 +37,6 @@ dev/
 *copy*
 *draft*
 
-.python-version
+.python-version
+
+*.obj

README.md

@@ -1,29 +1,33 @@
-# skeliner 
+# skeliner
 
 [![PyPI version](https://badge.fury.io/py/skeliner.svg)](https://badge.fury.io/py/skeliner)
 
-A lightweight skeletonizer that converts neuron meshes into biophysical‑modelling‑ready morphologies. It extracts an acyclic center‑line skeleton, estimates per‑node radii, detects the soma, and bridges small gaps. It also provides robust contact‑site mapping between pairs of cells using both skeletons and meshes.
+A lightweight skeletonizer that converts neuron meshes into biophysical‑modelling‑ready morphologies. It extracts an
+acyclic center‑line skeleton, estimates per‑node radii, detects the soma, and bridges small gaps. It also provides
+robust contact‑site mapping between pairs of cells using both skeletons and meshes.
 
 ![](.github/banner.png)
 
 ## Features
 
-- Mesh → SWC skeletons 
-  - Center‑line, acyclic skeleton.
-  - Per‑node radii with multiple estimators (median/mean/trim) and an automatic recommendation.
-  - Ellipsoidal soma detection near the soma centroid.
-  - Bridges disconnected mesh components back to the soma (gap‑closing) and prunes tiny peri‑soma artefacts.
-  - SWC and NPZ export with metadata;
+- Mesh → SWC skeletons
+    - Center‑line, acyclic skeleton.
+    - Per‑node radii with multiple estimators (median/mean/trim) and an automatic recommendation.
+    - Ellipsoidal soma detection near the soma centroid.
+    - Bridges disconnected mesh components back to the soma (gap‑closing) and prunes tiny peri‑soma artefacts.
+    - SWC and NPZ export with metadata;
 
 - Pairwise contact sites (two cells)
-  - Skeleton‑based seeding: fast KD‑tree search finds node↔node proximity pairs with tolerance; returns contact midpoints and per‑node loci.
-  - Mesh‑based patches: around seeds, projects to both meshes and extracts touching triangle patches; returns per‑site areas on A/B and their mean, plus compact AABBs for downstream use.
-  - Optional skeleton‑only approximation of contact sites/areas when meshes are unavailable.
+    - Skeleton‑based seeding: fast KD‑tree search finds node↔node proximity pairs with tolerance; returns contact
+      midpoints and per‑node loci.
+    - Mesh‑based patches: around seeds, projects to both meshes and extracts touching triangle patches; returns per‑site
+      areas on A/B and their mean, plus compact AABBs for downstream use.
+    - Optional skeleton‑only approximation of contact sites/areas when meshes are unavailable.
 
 - I/O, visualization, diagnostics and post-processing
-  - Read/write SWC and compact NPZ; load meshes via `trimesh`.
-  - 2D projections and 3D viewers for meshes, skeletons, and contact patches.
-  - Various diagnostics and post-processing tools.
+    - Read/write SWC and compact NPZ; load meshes via `trimesh`.
+    - 2D projections and 3D viewers for meshes, skeletons, and contact patches.
+    - Various diagnostics and post-processing tools.
 
 ## Installation
 
@@ -62,6 +66,43 @@ skel.to_swc("cellA.swc", scale=1e-3)
 skel.to_npz("cellA.npz")
 ```
 
+### Calibrate skeleton radii
+
+The quick skeletonization process often yields biased radius estimates, for two reasons:
+
+1. The raw radius estimators compute the distance to the node center, which is biased too be too low for cylindrical
+   segments.
+2. The mesh may have internal strucutre, typically from mitochondria or other organelles, which reduces the apparent
+   radius.
+
+Problem (1) is adressed by using the distance to the midline rather than the center point.  
+Problem (2) is adressed using ray‑casting to remove inner vertices of the mesh.
+This is computationally expensive, so per default it is only done for nodes that have many vertices.
+This can be controlled via the `min_verts_q_outer` parameter.
+
+You can use sk.post.calibrate_radii to improve it:
+
+```python
+import skeliner as sk
+
+mesh = sk.io.load_mesh("cellA.obj")  # or trimesh.load_mesh directly
+skel = sk.skeletonize(mesh, unit="nm", verbose=True)
+
+# Calibrate radii against the mesh (takes a few minutes)
+sk.post.calibrate_radii(
+    skel=skel,
+    mesh=mesh,
+    min_n_outer=20,
+    min_frac_outer=0.33,
+    min_verts_q_outer=80.,  # Set to 0 to check all nodes for internal structure, set to 100 to skip this step
+    rays_num_outer=30,
+    rays_thresh_outer=0.2,
+    verbose=True,
+    aggregate='trim',
+    store_key="calibrated",
+)
+```
+
 ### Find contact sites between two cells (skeleton + mesh)
 
 ```python

notebooks/data/ball_and_stick_neuron_v1.swc

@@ -0,0 +1,41 @@
+1 1 0.0 0.0 0 10.0 -1
+2 3 14.5 0.0 0 1.05 1
+3 3 19.0 0.0 0 1.1 2
+4 3 23.5 0.0 0 1.15 3
+5 3 28.0 0.0 0 1.2 4
+6 3 32.5 0.0 0 1.25 5
+7 3 37.0 0.0 0 1.3 6
+8 3 41.5 0.0 0 2.7 7
+9 3 46.0 0.0 0 1.4 8
+10 3 50.5 0.0 0 1.45 9
+11 3 55.0 0.0 0 1.5 10
+12 3 59.5 0.0 0 1.55 11
+13 3 64.0 0.0 0 1.6 12
+14 3 68.5 0.0 0 1.65 13
+15 3 73.0 0.0 0 1.7000000000000002 14
+16 3 77.5 0.0 0 1.75 15
+17 3 82.0 0.0 0 1.8 16
+18 3 86.5 0.0 0 1.85 17
+19 3 91.0 0.0 0 1.9 18
+20 3 95.5 0.0 0 1.9500000000000002 19
+21 3 100.0 0.0 0 4.0 20
+22 3 8.878689293818311e-16 14.5 0 2.1 1
+23 3 1.1634144591899855e-15 19.0 0 1.1 22
+24 3 1.43895998899814e-15 23.5 0 1.15 23
+25 3 1.7145055188062944e-15 28.0 0 1.2 24
+26 3 1.990051048614449e-15 32.5 0 1.25 25
+27 3 2.2655965784226034e-15 37.0 0 1.3 26
+28 3 2.541142108230758e-15 41.5 0 1.35 27
+29 3 2.8166876380389124e-15 46.0 0 1.4 28
+30 3 3.092233167847067e-15 50.5 0 1.45 29
+31 3 3.3677786976552213e-15 55.0 0 1.5 30
+32 3 3.643324227463376e-15 59.5 0 1.55 31
+33 3 3.91886975727153e-15 64.0 0 1.6 32
+34 3 4.194415287079684e-15 68.5 0 1.65 33
+35 3 4.469960816887839e-15 73.0 0 1.7000000000000002 34
+36 3 4.745506346695994e-15 77.5 0 1.75 35
+37 3 5.021051876504148e-15 82.0 0 1.8 36
+38 3 5.296597406312302e-15 86.5 0 1.85 37
+39 3 5.572142936120457e-15 91.0 0 1.9 38
+40 3 5.847688465928612e-15 95.5 0 1.9500000000000002 39
+41 3 6.123233995736766e-15 100.0 0 2.0 40

pyproject.toml

@@ -15,6 +15,11 @@ dependencies = [
     "trimesh==4.7.4",
     "matplotlib>=3.10.3",
     "rtree>=1.4.1",
+    "pandas>=2.3.3",
+    "jupyterlab>=4.4.10",
+    "notebook>=7.4.7",
+    "jupyter>=1.1.1",
+    "pytest>=8.3.5",
 ]
 license = "GPL-3.0-or-later"
 readme = {file = "README.md", content-type = "text/markdown"}

skeliner/_core.py

@@ -107,8 +107,8 @@ def _estimate_radius(
     d: np.ndarray,
     *,
     method: str = "median",
-    trim_fraction: float = 0.05,
-    q: float = 0.90,
+    trim_fraction: float = 0.15,
+    q: float = 80.,
 ) -> float:
     """Return one scalar radius according to *method*."""
     if method == "median":
@@ -119,6 +119,8 @@ def _estimate_radius(
         return float(d.max())
     if method == "min":
         return float(d.min())
+    if method == "percentile":
+        return float(np.percentile(d, q=q))
     if method == "trim":
         lo, hi = np.quantile(d, [trim_fraction, 1.0 - trim_fraction])
         mask = (d >= lo) & (d <= hi)

skeliner/dx.py

@@ -161,10 +161,39 @@ def distance(
     k_nearest: int = 4,
     radius_metric: str | None = None,
     mode: str = "surface",
+    allowed_nodes: Sequence[int] | None = None,
+    allowed_edges: Sequence[tuple[int, int]] | None = None,
 ) -> float | np.ndarray:
     """
     Distance from an arbitrary point (or collection of points) to the skeleton.
 
+    Parameters
+    ----------
+    skel
+        :class:`skeliner.Skeleton` instance.
+    point
+        3-vector or array of shape (M, 3) giving query locations.
+    point_unit
+        Unit of the input coordinates and the returned distance. If ``None`` or
+        identical to ``skel.meta['unit']``, no conversion is performed.
+    k_nearest
+        Number of nearest skeleton nodes considered when refining the distance
+        against neighbouring edges (≥ 1). Ignored when a whitelist is provided.
+    radius_metric
+        Which column of ``skel.radii`` to use. Defaults to the recommended estimator.
+        Only consulted when *mode* is ``'surface'``.
+    mode
+        ``'surface'`` (default) returns the distance to the radius-aware capsule
+        envelope (values inside clamp to ``0``). ``'centerline'`` measures distance
+        to the centreline alone.
+    allowed_nodes
+        Optional per-call whitelist of node IDs. When provided, distances are
+        restricted to these node centres and to edges incident to any of these nodes.
+    allowed_edges
+        Optional per-call whitelist of edges (u,v). When provided, distances are
+        refined only against these edges (u and v are also considered as centres).
+        Edges are treated as undirected; order is ignored.
+
     Parameters
     ----------
     skel
@@ -239,25 +268,75 @@ def distance(
     max_k = min(int(k_nearest), len(skel.nodes))
     nodes = skel.nodes
 
+    # Prepare whitelist sets if provided
+    use_whitelist = (allowed_nodes is not None) or (allowed_edges is not None)
+    allowed_nodes_set: Set[int] | None = None
+    allowed_edges_set: Set[tuple[int, int]] | None = None
+    if allowed_nodes is not None:
+        allowed_nodes_set = {int(n) for n in allowed_nodes if 0 <= int(n) < len(nodes)}
+    if allowed_edges is not None:
+        allowed_edges_set = set()
+        for (u, v) in allowed_edges:
+            u2 = int(u)
+            v2 = int(v)
+            if u2 == v2:
+                continue
+            if not (0 <= u2 < len(nodes) and 0 <= v2 < len(nodes)):
+                continue
+            a, b = (u2, v2) if u2 < v2 else (v2, u2)
+            allowed_edges_set.add((a, b))
+
     for i, p in enumerate(pts):
         p_skel = p * scale_in
 
-        nn_dist, nn_idx = tree.query(p_skel, k=max_k)
-        nn_idx_arr = np.atleast_1d(nn_idx).astype(np.int64, copy=False)
-        nn_dist_arr = np.atleast_1d(nn_dist)
-        if surface:
-            best = float(np.min(nn_dist_arr - radii[nn_idx_arr]))
-        else:
-            best = float(nn_dist_arr.min())
-
-        # Collect unique edges incident to the nearest nodes
-        candidates: set[tuple[int, int]] = set()
-        for nid in nn_idx_arr:
-            for nb in neighbours[nid]:
-                if nid < nb:
-                    candidates.add((nid, nb))
+        # Initialize best distance from node centres
+        if use_whitelist:
+            # Collect centres to consider: explicit allowed nodes and endpoints of allowed edges
+            centres: Set[int] = set()
+            if allowed_nodes_set is not None:
+                centres.update(allowed_nodes_set)
+            if allowed_edges_set is not None:
+                for a, b in allowed_edges_set:
+                    centres.add(a)
+                    centres.add(b)
+
+            if centres:
+                # compute min distance to allowed centres
+                centres_list = list(centres)
+                diffs = nodes[centres_list] - p_skel
+                nn_dist_arr = np.linalg.norm(diffs, axis=1)
+                if surface:
+                    rad = (np.asarray([radii[c] for c in centres_list], dtype=np.float64) if radii is not None else 0.0)
+                    best = float(np.min(nn_dist_arr - rad))
                 else:
-                    candidates.add((nb, nid))
+                    best = float(np.min(nn_dist_arr))
+            else:
+                best = float("inf")
+
+            # Candidate edges: explicit allowed_edges plus edges incident to allowed_nodes
+            candidates: Set[tuple[int, int]] = set()
+            if allowed_edges_set is not None:
+                candidates.update(allowed_edges_set)
+            if allowed_nodes_set is not None:
+                for nid in allowed_nodes_set:
+                    for nb in neighbours[nid]:
+                        a, b = (nid, nb) if nid < nb else (nb, nid)
+                        candidates.add((a, b))
+        else:
+            # default global behaviour via KD-tree + incident edges
+            nn_dist, nn_idx = tree.query(p_skel, k=max_k)
+            nn_idx_arr = np.atleast_1d(nn_idx).astype(np.int64, copy=False)
+            nn_dist_arr = np.atleast_1d(nn_dist)
+            if surface:
+                best = float(np.min(nn_dist_arr - radii[nn_idx_arr]))
+            else:
+                best = float(nn_dist_arr.min())
+
+            candidates: set[tuple[int, int]] = set()
+            for nid in nn_idx_arr:
+                for nb in neighbours[nid]:
+                    a, b = (nid, nb) if nid < nb else (nb, nid)
+                    candidates.add((a, b))
 
         if candidates:
             for a_idx, b_idx in candidates: