Skeletons (connectomekg.skeletons)

connectomekg.skeletons

SWC skeleton reading, writing and simplification -- NumPy only.

FlyWire's Codex download carries one navis-written .swc file per neuron (fafb_v783/sk_lod1_783_healed/<root_id>.swc): a header of # comment lines, then one row per traced point -- PointNo Label X Y Z Radius Parent, coordinates in nanometres, sharing the graph's own x/y/z frame. Parent names another row's PointNo, or -1 for a root; PointNo itself may have gaps and is not the row index. Files mix CRLF and LF line endings. Label values: 0 undefined, 1 soma, 5 fork point, 6 end point.

This module has no PyVista or Qt import anywhere, so it stays importable with no extra installed -- :mod:connectomekg.scene builds on it lazily.

Skeleton(root_id, points, radius, labels, parent) dataclass

One neuron's traced skeleton.

Parameters:

Name	Type	Description	Default
root_id	int	Neuron root id, from the SWC file's stem.	required
points	ndarray	(n, 3) float64 point coordinates, nanometres.	required
radius	ndarray	(n,) float64 radius per point, nanometres.	required
labels	ndarray	(n,) int SWC label per point (0 undefined, 1 soma, 5 fork, 6 end).	required
parent	ndarray	(n,) int, each point's parent as a row index into these arrays, or -1 for a root.	required

load_skeletons(data_dir, root_ids)

Load whichever of root_ids have a skeleton file, reporting the rest.

Parameters:

Name	Type	Description	Default
data_dir	str | Path	Root of the skeleton download.	required
root_ids	Iterable[int]	Neuron root ids to load.	required

Returns:

Type	Description
tuple[dict[int, Skeleton], list[int]]	(skeletons, missing) -- loaded skeletons keyed by root id, and the root ids with no file at :func:skeleton_path.

Source code in src/connectomekg/skeletons.py

def load_skeletons(
    data_dir: str | Path, root_ids: Iterable[int]
) -> tuple[dict[int, Skeleton], list[int]]:
    """Load whichever of ``root_ids`` have a skeleton file, reporting the rest.

    :param data_dir: Root of the skeleton download.
    :param root_ids: Neuron root ids to load.
    :return: ``(skeletons, missing)`` -- loaded skeletons keyed by root id,
        and the root ids with no file at :func:`skeleton_path`.
    """
    loaded: dict[int, Skeleton] = {}
    missing: list[int] = []
    for root_id in root_ids:
        path = skeleton_path(data_dir, root_id)
        if path.exists():
            loaded[int(root_id)] = read_swc(path)
        else:
            missing.append(int(root_id))
    return loaded, missing

polylines(skeleton, *, step=1)

Maximal unbranched chains, as index arrays into the skeleton's points.

:func:segments returns each child-to-parent edge on its own, which a tube filter extrudes into one short cylinder per edge: consecutive cylinders meet at an angle with nothing joining them, and the drawn neuron is a heap of faceted stubs. A chain is the run of points between one branch point (or root) and the next branch point (or tip), so tubing it produces a single continuous surface with mitred joins. It is also cheaper: a chain's interior points are shared rather than duplicated once per edge.

Every edge belongs to exactly one chain, so the drawn geometry is the same as :func:segments would give -- only joined up.

Parameters:

Name	Type	Description	Default
skeleton	Skeleton	The skeleton to draw.	required
step	int	Simplification stride, as :func:segments takes it. Chains are found over the simplified topology, not the full one.	1

Returns:

Type	Description
list[ndarray]	One (k,) int array per chain, k >= 2, indexing skeleton.points and skeleton.radius.

Source code in src/connectomekg/skeletons.py

def polylines(skeleton: Skeleton, *, step: int = 1) -> list[np.ndarray]:
    """Maximal unbranched chains, as index arrays into the skeleton's points.

    :func:`segments` returns each child-to-parent edge on its own, which a
    tube filter extrudes into one short cylinder per edge: consecutive
    cylinders meet at an angle with nothing joining them, and the drawn
    neuron is a heap of faceted stubs. A chain is the run of points between
    one branch point (or root) and the next branch point (or tip), so tubing
    it produces a single continuous surface with mitred joins. It is also
    cheaper: a chain's interior points are shared rather than duplicated once
    per edge.

    Every edge belongs to exactly one chain, so the drawn geometry is the same
    as :func:`segments` would give -- only joined up.

    :param skeleton: The skeleton to draw.
    :param step: Simplification stride, as :func:`segments` takes it. Chains
        are found over the simplified topology, not the full one.
    :return: One ``(k,)`` int array per chain, ``k >= 2``, indexing
        ``skeleton.points`` and ``skeleton.radius``.
    """
    parent = skeleton.parent
    if step <= 1:
        keep = np.ones(len(parent), dtype=bool)
        anc = parent
    else:
        keep = _simplify_keep_mask(parent, skeleton.labels, step)
        anc = _nearest_kept_ancestor(parent, keep)

    kept = np.nonzero(keep)[0]
    child = kept[anc[kept] != -1]
    if child.size == 0:
        return []
    par = anc[child]
    n = len(parent)
    n_children = np.bincount(par, minlength=n)
    # The one child of each single-child point, for walking a chain onward.
    only_child = np.full(n, -1, dtype=np.int64)
    single = n_children[par] == 1
    only_child[par[single]] = child[single]
    # A chain starts at an edge whose parent forks, or is itself a root: any
    # other parent has exactly one edge in and one out, so the chain runs on
    # through it.
    starts = child[(n_children[par] != 1) | (anc[par] == -1)]

    chains: list[np.ndarray] = []
    for start in starts:
        chain = [int(anc[start]), int(start)]
        cur = int(start)
        while (nxt := only_child[cur]) != -1:
            cur = int(nxt)
            chain.append(cur)
        chains.append(np.asarray(chain, dtype=np.int64))
    return chains

read_swc(path)

Parse one SWC file into a :class:Skeleton.

Parameters:

Name	Type	Description	Default
path	str | Path	Path to the .swc file; root_id comes from its stem.	required

Returns:

Type	Description
Skeleton	The parsed skeleton.

Raises:

Type	Description
ValueError	If a row has too few fields, a Parent names no PointNo in the file, or the file has no point rows at all.

Source code in src/connectomekg/skeletons.py

def read_swc(path: str | Path) -> Skeleton:
    """Parse one SWC file into a :class:`Skeleton`.

    :param path: Path to the ``.swc`` file; ``root_id`` comes from its stem.
    :return: The parsed skeleton.
    :raises ValueError: If a row has too few fields, a ``Parent`` names no
        ``PointNo`` in the file, or the file has no point rows at all.
    """
    path = Path(path)
    # str.splitlines() treats CRLF, LF and lone CR uniformly -- the files mix
    # line endings, and this sidesteps that rather than special-casing it.
    raw_lines = path.read_text(encoding="utf-8", errors="replace").splitlines()

    point_ids: list[int] = []
    labels: list[int] = []
    coords: list[tuple[float, float, float]] = []
    radii: list[float] = []
    parent_ids: list[int] = []
    line_numbers: list[int] = []

    for lineno, raw in enumerate(raw_lines, start=1):
        line = raw.strip()
        if not line or line.startswith("#"):
            continue
        parts = line.split()
        if len(parts) < 7:
            raise ValueError(f"{path}: line {lineno}: expected 7 fields, got {len(parts)}")
        pid, label, x, y, z, r, parent = parts[:7]
        point_ids.append(int(pid))
        labels.append(int(label))
        coords.append((float(x), float(y), float(z)))
        radii.append(float(r))
        parent_ids.append(int(parent))
        line_numbers.append(lineno)

    if not point_ids:
        raise ValueError(f"{path}: no SWC point rows found")

    index_of = {pid: i for i, pid in enumerate(point_ids)}
    parent = np.empty(len(point_ids), dtype=np.int64)
    for i, (praw, lineno) in enumerate(zip(parent_ids, line_numbers, strict=True)):
        if praw == -1:
            parent[i] = -1
        elif praw in index_of:
            parent[i] = index_of[praw]
        else:
            raise ValueError(f"{path}: line {lineno}: parent {praw} names no PointNo in this file")

    try:
        root_id = int(path.stem)
    except ValueError as exc:
        raise ValueError(f"{path}: file stem {path.stem!r} is not a root id") from exc

    return Skeleton(
        root_id=root_id,
        points=np.asarray(coords, dtype=np.float64),
        radius=np.asarray(radii, dtype=np.float64),
        labels=np.asarray(labels, dtype=np.int64),
        parent=parent,
    )

segments(skeleton, *, step=1)

Child-to-parent line segments, optionally simplified.

Parameters:

Name	Type	Description	Default
skeleton	Skeleton	The skeleton to draw.	required
step	int	1 (default) draws every point-to-parent edge. Above 1, simplifies first via :func:_simplify_keep_mask and connects each kept point straight to its nearest kept ancestor, so the simplified skeleton's topology (which points connect to which) matches the full one.	1

Returns:

Type	Description
ndarray	(m, 2, 3) array of (start, end) point pairs, nanometres.

Source code in src/connectomekg/skeletons.py

def segments(skeleton: Skeleton, *, step: int = 1) -> np.ndarray:
    """Child-to-parent line segments, optionally simplified.

    :param skeleton: The skeleton to draw.
    :param step: ``1`` (default) draws every point-to-parent edge. Above
        ``1``, simplifies first via :func:`_simplify_keep_mask` and connects
        each kept point straight to its nearest kept ancestor, so the
        simplified skeleton's topology (which points connect to which)
        matches the full one.
    :return: ``(m, 2, 3)`` array of ``(start, end)`` point pairs, nanometres.
    """
    parent = skeleton.parent
    if step <= 1:
        keep = np.ones(len(parent), dtype=bool)
        anc = parent
    else:
        keep = _simplify_keep_mask(parent, skeleton.labels, step)
        anc = _nearest_kept_ancestor(parent, keep)
    idx = np.nonzero(keep & (anc != -1))[0]
    if idx.size == 0:
        return np.empty((0, 2, 3), dtype=np.float64)
    starts = skeleton.points[idx]
    ends = skeleton.points[anc[idx]]
    return np.stack([starts, ends], axis=1)

skeleton_path(data_dir, root_id)

The SWC file path for one neuron in the Codex skeleton download.

Parameters:

Name	Type	Description	Default
data_dir	str | Path	Root of the download (default layout: fafb_v783).	required
root_id	int	Neuron root id.	required

Returns:

Type	Description
Path	<data_dir>/sk_lod1_783_healed/<root_id>.swc.

Source code in src/connectomekg/skeletons.py

def skeleton_path(data_dir: str | Path, root_id: int) -> Path:
    """The SWC file path for one neuron in the Codex skeleton download.

    :param data_dir: Root of the download (default layout: ``fafb_v783``).
    :param root_id: Neuron root id.
    :return: ``<data_dir>/sk_lod1_783_healed/<root_id>.swc``.
    """
    return Path(data_dir) / SKELETON_SUBDIR / f"{int(root_id)}.swc"

soma(skeleton)

The skeleton's soma point, or a fallback.

Parameters:

Name	Type	Description	Default
skeleton	Skeleton	The skeleton to inspect.	required

Returns:

Type	Description
tuple[ndarray, bool]	(point, is_soma) -- the Label 1 point's coordinates and True when one exists; otherwise the first root point's coordinates and False.

Raises:

Type	Description
ValueError	If the skeleton has neither a Label 1 row nor a root (Parent == -1) row.

Source code in src/connectomekg/skeletons.py

def soma(skeleton: Skeleton) -> tuple[np.ndarray, bool]:
    """The skeleton's soma point, or a fallback.

    :param skeleton: The skeleton to inspect.
    :return: ``(point, is_soma)`` -- the ``Label 1`` point's coordinates and
        ``True`` when one exists; otherwise the first root point's
        coordinates and ``False``.
    :raises ValueError: If the skeleton has neither a ``Label 1`` row nor a
        root (``Parent == -1``) row.
    """
    labeled = np.nonzero(skeleton.labels == _LABEL_SOMA)[0]
    if labeled.size:
        return skeleton.points[labeled[0]], True
    roots = np.nonzero(skeleton.parent == -1)[0]
    if roots.size == 0:
        raise ValueError("skeleton has no soma row and no root point")
    return skeleton.points[roots[0]], False

write_swc(skeleton, path)

Write a :class:Skeleton back out as an SWC file, for fixtures and round trips.

PointNo is the 1-based row index; this is only ever read back by :func:read_swc, which does not assume that, so the choice is arbitrary but self-consistent.

Parameters:

Name	Type	Description	Default
skeleton	Skeleton	The skeleton to write.	required
path	str | Path	Destination path.	required

Source code in src/connectomekg/skeletons.py

def write_swc(skeleton: Skeleton, path: str | Path) -> None:
    """Write a :class:`Skeleton` back out as an SWC file, for fixtures and round trips.

    ``PointNo`` is the 1-based row index; this is only ever read back by
    :func:`read_swc`, which does not assume that, so the choice is arbitrary
    but self-consistent.

    :param skeleton: The skeleton to write.
    :param path: Destination path.
    """
    path = Path(path)
    lines = [
        "# SWC format file",
        f'# Meta: {{"id": "{skeleton.root_id}"}}',
        "# PointNo Label X Y Z Radius Parent",
    ]
    for i in range(len(skeleton)):
        point_no = i + 1
        parent_no = -1 if skeleton.parent[i] == -1 else int(skeleton.parent[i]) + 1
        x, y, z = (float(v) for v in skeleton.points[i])
        lines.append(
            f"{point_no} {int(skeleton.labels[i])} {x} {y} {z} "
            f"{float(skeleton.radius[i])} {parent_no}"
        )
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")