Skip to content

geometry

Morphology geometry and transformation utilities.

toric_spines_sim.geometry

Morphology geometry and transformation utilities.

SinkGeometry(radius=0.5, length=100.0, n_cylinders=1, connector_length=1.0, axis='x') dataclass

Straight cylindrical sink (tag 5 body, optional tip tag).

Attributes:

Name Type Description
radius float

Cylinder radius in the same units as the SWC.

length float

Axial length of the sink (excluding connector_length).

n_cylinders int

Number of frusta along the axis.

connector_length float

Short segment from the snapped neck node to the first sink node.

axis str or sequence of float

'x'/'y'/'z' (optional leading -) or a 3-vector.

SpinyDendriteParams(length, trunk_neck_radius, trunk_tip_radius, n_spines, spine_length, spine_neck_radius, spine_head_radius, spine_neck_length_fraction=0.5, max_spines_per_node=1, distribution='even', seed=None, azimuth0=0.0, axis='z', trunk_tag=3, spine_neck_tag=3, spine_head_tag=3) dataclass

Explicit parameters for building a spiny dendrite (mode A).

ToricSpineMatchParams(spine_length=_DEFAULT_SPINE_LENGTH, spine_neck_radius=_DEFAULT_SPINE_NECK_RADIUS, spine_head_radius=_DEFAULT_SPINE_HEAD_RADIUS, trunk_neck_radius=_DEFAULT_TRUNK_NECK_RADIUS, trunk_tip_radius=_DEFAULT_TRUNK_TIP_RADIUS, spine_neck_length_fraction=0.5, trunk_length_over_spine_length=_DEFAULT_TRUNK_LENGTH_OVER_SPINE_LENGTH, max_spines_per_node=1, distribution='even', seed=None, azimuth0=0.0, axis='z', scale_strategy='relative', trunk_tag=3, spine_neck_tag=3, spine_head_tag=3, sink_tags=(lambda: _DEFAULT_SINK_TAGS)()) dataclass

Template / layout options when matching a toric spine (mode B).

SpinyDendriteMorphology(nodes, az_points, neck_point, trunk_node_ids, spine_head_node_ids, spines_per_attach_node, params, diagnostics=dict()) dataclass

In-memory spiny dendrite subsystem.

SWCNodeRecord(node_id, tag, x, y, z, radius, parent) dataclass

One SWC data line.

MeshToSwcResult(mesh_path, polylines_path, swc_path) dataclass

Outputs from :func:mesh_to_swc.

NeckCandidate(centroid, area, n_faces, planarity, mean_signed, mean_align, ray_frac_inside, dendrite_frac, source='cap', rank_score=0.0) dataclass

One detected neck interface.

NeckpointParams(crop_margin=120.0, align_thr=0.2, align_signed_max=-3.0, signed_thr=-10.0, min_faces=6, planarity_max=0.22, mean_signed_max=-5.0, mean_align_max=0.5, ray_frac_min=0.6, open_score_min=0.75, ray_offsets=(10.0, 25.0, 50.0, 80.0, 120.0), dendrite_radius=150.0, dendrite_far_thr=20.0, dendrite_frac_floor=0.15, dendrite_frac_of_max=0.5, fallback_dist_thr=15.0, fallback_align_thr=0.5, fallback_min_loop_verts=6, attachment_dendrite_dist_max=35.0, attachment_facing_min=0.25, attachment_min_faces=5, attachment_planarity_max=0.3, attachment_ray_offsets=(10.0, 30.0, 60.0, 100.0), max_necks=None) dataclass

Thresholds for TS-vs-cell neck-cap detection.

scale_one_radius_in_segment_tree_by_coordinates(tree, scale_factor, target_xyz, tolerance=1e-06)

Scale the radius of a specific node in an Arbor segment_tree, identified by its (x, y, z) coordinates.

Every proximal or distal endpoint whose position matches target_xyz (within tolerance) has its radius multiplied by scale_factor.

Source code in toric_spines_sim/geometry/rescale.py
def scale_one_radius_in_segment_tree_by_coordinates(
    tree: A.segment_tree,
    scale_factor: float,
    target_xyz: Tuple[float, float, float],
    tolerance: float = 1e-6,
):
    """
    Scale the radius of a specific node in an Arbor segment_tree, identified
    by its (x, y, z) coordinates.

    Every proximal or distal endpoint whose position matches *target_xyz*
    (within *tolerance*) has its radius multiplied by *scale_factor*.
    """
    tx, ty, tz = target_xyz
    tol2 = tolerance * tolerance
    new_tree = A.segment_tree()
    matched = 0
    for parent_idx, seg in zip(tree.parents, tree.segments):
        p, d = seg.prox, seg.dist

        pr = p.radius
        if (p.x - tx) ** 2 + (p.y - ty) ** 2 + (p.z - tz) ** 2 <= tol2:
            pr = p.radius * scale_factor
            matched += 1

        dr = d.radius
        if (d.x - tx) ** 2 + (d.y - ty) ** 2 + (d.z - tz) ** 2 <= tol2:
            dr = d.radius * scale_factor
            matched += 1

        new_prox = A.mpoint(p.x, p.y, p.z, pr)
        new_dist = A.mpoint(d.x, d.y, d.z, dr)
        new_tree.append(parent_idx, new_prox, new_dist, seg.tag)

    if matched == 0:
        logger.warning(
            "scale_radius_at_node: no endpoints matched (%.6f, %.6f, %.6f) "
            "within tolerance %.2e",
            tx,
            ty,
            tz,
            tolerance,
        )
    else:
        logger.debug(
            "scale_radius_at_node: scaled %d endpoint(s) at (%.6f, %.6f, %.6f) by %.4f",
            matched,
            tx,
            ty,
            tz,
            scale_factor,
        )
    return new_tree

scale_radii_in_segment_tree_by_tag(tree, scale_factor, scale_tag)

Scale the radii of all segments in an Arbor segment_tree with a given tag.

The distal radius of a matching segment is always scaled. The proximal radius is only scaled when the parent segment also carries scale_tag, so that boundary nodes (e.g. a neck node shared with a different-tagged region) are not affected.

Source code in toric_spines_sim/geometry/rescale.py
def scale_radii_in_segment_tree_by_tag(
    tree: A.segment_tree, scale_factor: float, scale_tag: int
):
    """
    Scale the radii of all segments in an Arbor segment_tree with a given tag.

    The distal radius of a matching segment is always scaled.  The proximal
    radius is only scaled when the parent segment also carries *scale_tag*,
    so that boundary nodes (e.g. a neck node shared with a different-tagged
    region) are not affected.
    """
    segments = tree.segments
    parents = tree.parents
    new_tree = A.segment_tree()
    for i, (parent_idx, seg) in enumerate(zip(parents, segments)):
        p, d = seg.prox, seg.dist

        if seg.tag == scale_tag:
            # Scale proximal only when the parent segment shares the same tag
            parent_has_same_tag = (
                parent_idx >= 0 and segments[parent_idx].tag == scale_tag
            )
            pr = p.radius * scale_factor if parent_has_same_tag else p.radius
            dr = d.radius * scale_factor
        else:
            pr = p.radius
            dr = d.radius

        new_prox = A.mpoint(p.x, p.y, p.z, pr)
        new_dist = A.mpoint(d.x, d.y, d.z, dr)
        new_tree.append(parent_idx, new_prox, new_dist, seg.tag)
    return new_tree

compute_geodesic_distances(swc_filepath, source_xyz)

Compute geodesic (path-length) distances from a source point to every node.

Parameters:

Name Type Description Default
swc_filepath path - like

Path to an SWC file (with optional CYCLE_BREAK reconnections).

required
source_xyz (x, y, z)

3-D coordinate of the source point. The nearest graph node is used.

required

Returns:

Type Description
dict[int, float]

Mapping from node ID to shortest-path distance (in the same spatial units as the SWC file, typically µm). Unreachable nodes are omitted.

Source code in toric_spines_sim/geometry/graph.py
def compute_geodesic_distances(
    swc_filepath: Union[str, Path],
    source_xyz: Tuple[float, float, float],
) -> Dict[int, float]:
    """Compute geodesic (path-length) distances from a source point to every node.

    Parameters
    ----------
    swc_filepath : path-like
        Path to an SWC file (with optional ``CYCLE_BREAK`` reconnections).
    source_xyz : (x, y, z)
        3-D coordinate of the source point.  The nearest graph node is used.

    Returns
    -------
    dict[int, float]
        Mapping from node ID to shortest-path distance (in the same spatial
        units as the SWC file, typically µm).  Unreachable nodes are omitted.
    """
    graph = _load_undirected_graph(swc_filepath)
    source_node = _nearest_node(graph, source_xyz)
    distances = dict(
        nx.single_source_dijkstra_path_length(graph, source_node, weight="length")
    )
    logger.info(
        "Geodesic distances from node %d: %d reachable nodes, max=%.2f",
        source_node,
        len(distances),
        max(distances.values()) if distances else 0.0,
    )
    return distances

classify_compartments(swc_filepath, source_xyz, target_xyz, mode='path_based', spine_tag=3, sink_tag=5)

Classify every graph node into one of four categories.

Categories

"main_path" Nodes on the shortest path from source to target. "branch" Spine nodes (tag = spine_tag) that are off the main path but structurally related to it (definition depends on mode). "lateral" All other spine nodes. "sink" Nodes with tag = sink_tag.

Parameters:

Name Type Description Default
swc_filepath path - like

Path to an SWC file.

required
source_xyz (x, y, z)

Coordinates of source and target landmarks. Nearest graph nodes are used.

required
target_xyz (x, y, z)

Coordinates of source and target landmarks. Nearest graph nodes are used.

required
mode ``"path_based"`` | ``"distance_based"``

How to distinguish branch from lateral:

  • "path_based": a spine node is branch if its shortest path to the source passes through at least one interior main-path node.
  • "distance_based": a spine node is branch if its geodesic distance to the target is less than the source's geodesic distance to the target (i.e., it lies "closer to the target").
'path_based'
spine_tag int

Tag value identifying spine compartments (default 3).

3
sink_tag int

Tag value identifying sink compartments (default 5).

5

Returns:

Type Description
dict[int, str]

Mapping {node_id: category_label} for every node in the graph.

Source code in toric_spines_sim/geometry/graph.py
def classify_compartments(
    swc_filepath: Union[str, Path],
    source_xyz: Tuple[float, float, float],
    target_xyz: Tuple[float, float, float],
    mode: Literal["path_based", "distance_based"] = "path_based",
    spine_tag: int = 3,
    sink_tag: int = 5,
) -> Dict[int, str]:
    """Classify every graph node into one of four categories.

    Categories
    ----------
    ``"main_path"``
        Nodes on the shortest path from *source* to *target*.
    ``"branch"``
        Spine nodes (tag = *spine_tag*) that are off the main path but
        structurally related to it (definition depends on *mode*).
    ``"lateral"``
        All other spine nodes.
    ``"sink"``
        Nodes with tag = *sink_tag*.

    Parameters
    ----------
    swc_filepath : path-like
        Path to an SWC file.
    source_xyz, target_xyz : (x, y, z)
        Coordinates of source and target landmarks.  Nearest graph nodes are
        used.
    mode : ``"path_based"`` | ``"distance_based"``
        How to distinguish *branch* from *lateral*:

        * ``"path_based"``: a spine node is *branch* if its shortest path to
          the source passes through at least one interior main-path node.
        * ``"distance_based"``: a spine node is *branch* if its geodesic
          distance to the target is less than the source's geodesic distance
          to the target (i.e., it lies "closer to the target").
    spine_tag : int
        Tag value identifying spine compartments (default 3).
    sink_tag : int
        Tag value identifying sink compartments (default 5).

    Returns
    -------
    dict[int, str]
        Mapping ``{node_id: category_label}`` for every node in the graph.
    """
    graph = _load_undirected_graph(swc_filepath)
    source_node = _nearest_node(graph, source_xyz)
    target_node = _nearest_node(graph, target_xyz)

    # Shortest path (list of node IDs) from source to target
    main_path_nodes: list[int] = nx.shortest_path(
        graph, source_node, target_node, weight="length"
    )
    main_path_set: Set[int] = set(main_path_nodes)
    # Interior main-path nodes (excluding source and target themselves)
    interior_main_path: Set[int] = main_path_set - {source_node, target_node}

    logger.info(
        "Main path %d%d: %d nodes",
        source_node,
        target_node,
        len(main_path_nodes),
    )

    # Pre-compute distances needed by both modes
    distances_from_source = dict(
        nx.single_source_dijkstra_path_length(graph, source_node, weight="length")
    )

    if mode == "distance_based":
        distances_from_target = dict(
            nx.single_source_dijkstra_path_length(
                graph, target_node, weight="length"
            )
        )
        source_to_target_distance = distances_from_source.get(
            target_node, float("inf")
        )

    if mode == "path_based":
        # For each non-main-path node, check whether its shortest path to
        # the source passes through an interior main-path node.
        paths_from_source = nx.single_source_dijkstra_path(
            graph, source_node, weight="length"
        )

    classification: Dict[int, str] = {}
    for node_id, attrs in graph.nodes(data=True):
        tag = attrs.get("t", 0)

        # Sink
        if tag == sink_tag:
            classification[node_id] = "sink"
            continue

        # Main path
        if node_id in main_path_set:
            classification[node_id] = "main_path"
            continue

        # Only spine nodes get branch / lateral distinction
        if tag != spine_tag:
            # Nodes with other tags (shouldn't normally happen for TS morphologies)
            classification[node_id] = "lateral"
            continue

        if mode == "path_based":
            # Check if shortest path from source to this node crosses a
            # main-path interior node.
            path_to_node = paths_from_source.get(node_id, [])
            if interior_main_path.intersection(path_to_node):
                classification[node_id] = "branch"
            else:
                classification[node_id] = "lateral"

        elif mode == "distance_based":
            node_to_target = distances_from_target.get(node_id, float("inf"))
            if node_to_target < source_to_target_distance:
                classification[node_id] = "branch"
            else:
                classification[node_id] = "lateral"
        else:
            raise ValueError(
                f"Unknown mode '{mode}'. Use 'path_based' or 'distance_based'."
            )

    # Log summary
    counts = {}
    for label in classification.values():
        counts[label] = counts.get(label, 0) + 1
    logger.info("Compartment classification (mode=%s): %s", mode, counts)

    return classification

geodesic_distances_from_probe(swc_filepath, record_points, source_probe, target_probes)

Compute geodesic distances from one probe to target probes.

Parameters:

Name Type Description Default
swc_filepath path - like

Path to the SWC morphology.

required
record_points dict

Mapping {probe_label: (x, y, z)} for all simulation record points.

required
source_probe str

Probe label used as source for distance calculation.

required
target_probes iterable[str]

Probe labels to measure from source_probe.

required

Returns:

Type Description
dict[str, float]

Mapping {probe_label: geodesic_distance} in SWC spatial units. Targets that cannot be mapped/reached are omitted.

Source code in toric_spines_sim/geometry/graph.py
def geodesic_distances_from_probe(
    swc_filepath: Union[str, Path],
    record_points: Dict[str, Tuple[float, float, float]],
    source_probe: str,
    target_probes: Iterable[str],
) -> Dict[str, float]:
    """Compute geodesic distances from one probe to target probes.

    Parameters
    ----------
    swc_filepath : path-like
        Path to the SWC morphology.
    record_points : dict
        Mapping ``{probe_label: (x, y, z)}`` for all simulation record points.
    source_probe : str
        Probe label used as source for distance calculation.
    target_probes : iterable[str]
        Probe labels to measure from ``source_probe``.

    Returns
    -------
    dict[str, float]
        Mapping ``{probe_label: geodesic_distance}`` in SWC spatial units.
        Targets that cannot be mapped/reached are omitted.
    """
    if source_probe not in record_points:
        raise KeyError(f"Source probe '{source_probe}' not found in record_points.")

    graph = _load_undirected_graph(swc_filepath)
    probe_to_node = map_probes_to_nodes(swc_filepath, record_points)
    if source_probe not in probe_to_node:
        raise KeyError(f"Source probe '{source_probe}' could not be mapped to an SWC node.")

    source_node = probe_to_node[source_probe]
    node_distances = dict(
        nx.single_source_dijkstra_path_length(graph, source_node, weight="length")
    )

    out: Dict[str, float] = {}
    for probe in target_probes:
        node = probe_to_node.get(probe)
        if node is None:
            continue
        dist = node_distances.get(node)
        if dist is None:
            continue
        out[probe] = float(dist)
    return out

map_probes_to_nodes(swc_filepath, record_points)

Map simulation probe labels to their nearest graph node IDs.

Parameters:

Name Type Description Default
swc_filepath path - like

Path to an SWC file.

required
record_points dict

Mapping {probe_label: (x, y, z)} as returned by :func:~toric_spines_sim.swc.get_center_coordinates_for_all_segments or stored in SimulationResults.record_points.

required

Returns:

Type Description
dict[str, int]

Mapping {probe_label: node_id}.

Source code in toric_spines_sim/geometry/graph.py
def map_probes_to_nodes(
    swc_filepath: Union[str, Path],
    record_points: Dict[str, Tuple[float, float, float]],
) -> Dict[str, int]:
    """Map simulation probe labels to their nearest graph node IDs.

    Parameters
    ----------
    swc_filepath : path-like
        Path to an SWC file.
    record_points : dict
        Mapping ``{probe_label: (x, y, z)}`` as returned by
        :func:`~toric_spines_sim.swc.get_center_coordinates_for_all_segments`
        or stored in ``SimulationResults.record_points``.

    Returns
    -------
    dict[str, int]
        Mapping ``{probe_label: node_id}``.
    """
    graph = _load_undirected_graph(swc_filepath)

    # Build array of node coordinates for vectorised lookup
    node_ids = list(graph.nodes())
    node_coords = np.array(
        [[graph.nodes[n]["x"], graph.nodes[n]["y"], graph.nodes[n]["z"]] for n in node_ids]
    )

    probe_to_node: Dict[str, int] = {}
    for probe_label, (px, py, pz) in record_points.items():
        diffs = node_coords - np.array([px, py, pz])
        dist_sq = np.sum(diffs**2, axis=1)
        nearest_idx = int(np.argmin(dist_sq))
        probe_to_node[probe_label] = node_ids[nearest_idx]

    logger.info("Mapped %d probes to graph nodes", len(probe_to_node))
    return probe_to_node

map_xyz_to_nearest_probes(record_points, xyz_points)

Map named xyz coordinates to nearest recording probe labels.

Parameters:

Name Type Description Default
record_points dict

Mapping {probe_label: (x, y, z)} for available compartments/probes.

required
xyz_points dict

Mapping {name: (x, y, z)} for query points.

required

Returns:

Type Description
dict[str, str]

Mapping {name: nearest_probe_label}.

Source code in toric_spines_sim/geometry/graph.py
def map_xyz_to_nearest_probes(
    record_points: Dict[str, Tuple[float, float, float]],
    xyz_points: Dict[str, Tuple[float, float, float]],
) -> Dict[str, str]:
    """Map named xyz coordinates to nearest recording probe labels.

    Parameters
    ----------
    record_points : dict
        Mapping ``{probe_label: (x, y, z)}`` for available compartments/probes.
    xyz_points : dict
        Mapping ``{name: (x, y, z)}`` for query points.

    Returns
    -------
    dict[str, str]
        Mapping ``{name: nearest_probe_label}``.
    """
    if not record_points:
        raise ValueError("record_points cannot be empty.")

    probe_labels = list(record_points.keys())
    probe_coords = np.array([record_points[label] for label in probe_labels], dtype=float)

    nearest: Dict[str, str] = {}
    for name, (x, y, z) in xyz_points.items():
        diffs = probe_coords - np.array([x, y, z], dtype=float)
        dist_sq = np.sum(diffs**2, axis=1)
        nearest_idx = int(np.argmin(dist_sq))
        nearest[name] = probe_labels[nearest_idx]

    logger.info(
        "Mapped %d xyz points to nearest probes from %d record points",
        len(xyz_points),
        len(record_points),
    )
    return nearest

parse_cycle_breaks(swc_path)

Parse # CYCLE_BREAK reconnect i j annotations from an SWC file.

Returns a list of integer pairs [(i, j), ...].

Raises ValueError if any CYCLE_BREAK directive is malformed (e.g. non-integer IDs).

Source code in toric_spines_sim/geometry/swc.py
def parse_cycle_breaks(swc_path: Path):
    """Parse ``# CYCLE_BREAK reconnect i j`` annotations from an SWC file.

    Returns a list of integer pairs ``[(i, j), ...]``.

    Raises ValueError if any CYCLE_BREAK directive is malformed (e.g. non-integer IDs).
    """
    try:
        result = parse_swc(str(swc_path), validate_reconnections=False)
    except Exception:
        logger.error(
            "Failed to read SWC file for cycle breaks: %s", swc_path, exc_info=True
        )
        raise ValueError(f"Invalid SWC file: {swc_path}")

    # Count CYCLE_BREAK directives in the header to detect malformed ones that
    # parse_swc silently skipped (e.g. non-integer node IDs).
    directive_count = sum(
        1
        for line in result.header
        if "cycle_break" in line.lower() and "reconnect" in line.lower()
    )
    if directive_count != len(result.reconnections):
        raise ValueError(
            f"Found {directive_count} CYCLE_BREAK directive(s) in header but only "
            f"{len(result.reconnections)} parsed successfully; check for malformed lines"
        )

    reconnect_pairs = list(result.reconnections)
    logger.debug(
        "Parsed %d cycle break reconnect pairs from %s", len(reconnect_pairs), swc_path
    )
    return reconnect_pairs

parse_multi_neck_reconnects(swc_path)

Parse # MULTI_NECK reconnect i j annotations from an SWC header.

Extra necks are electrically tied to the sink start by a gap junction between node i (sink start) and node j (colocated sink-start copy).

Source code in toric_spines_sim/geometry/swc.py
def parse_multi_neck_reconnects(swc_path: Path) -> List[Tuple[int, int]]:
    """Parse ``# MULTI_NECK reconnect i j`` annotations from an SWC header.

    Extra necks are electrically tied to the sink start by a gap junction
    between node ``i`` (sink start) and node ``j`` (colocated sink-start copy).
    """
    return _parse_reconnect_kind(swc_path, "MULTI_NECK")

parse_reconnect_pairs(swc_path)

Return CYCLE_BREAK pairs followed by MULTI_NECK pairs (no duplicates).

Source code in toric_spines_sim/geometry/swc.py
def parse_reconnect_pairs(swc_path: Path) -> List[Tuple[int, int]]:
    """Return CYCLE_BREAK pairs followed by MULTI_NECK pairs (no duplicates)."""
    cycle_pairs = parse_cycle_breaks(swc_path)
    extra = parse_multi_neck_reconnects(swc_path)
    seen = set(cycle_pairs)
    merged = list(cycle_pairs)
    for pair in extra:
        if pair not in seen and (pair[1], pair[0]) not in seen:
            merged.append(pair)
            seen.add(pair)
    logger.debug(
        "Parsed %d reconnect pairs (%d cycle-break, %d extra) from %s",
        len(merged),
        len(cycle_pairs),
        len(merged) - len(cycle_pairs),
        swc_path,
    )
    return merged

read_swc_points(swc_path)

Return dict id -> (x, y, z, r) from SWC content (ignores non-data lines).

Source code in toric_spines_sim/geometry/swc.py
def read_swc_points(swc_path: Path):
    """Return dict id -> (x, y, z, r) from SWC content (ignores non-data lines)."""
    try:
        swc_model = SWCModel.from_swc_file(str(swc_path), validate_reconnections=False)
    except Exception:
        logger.error("Failed to read SWC file: %s", swc_path, exc_info=True)
        raise ValueError(f"Invalid SWC file: {swc_path}")
    points_by_id = {}
    for node_id in swc_model.nodes():
        xyz = swc_model.get_node_xyz(node_id)
        radius = swc_model.get_node_radius(node_id)
        points_by_id[node_id] = (*xyz, radius)
    logger.debug("Read %d SWC points from %s", len(points_by_id), swc_path)
    return points_by_id

get_center_coordinates_for_all_segments(swc_filepath, use_radius_weighting=False)

Compute the center coordinate of each SWC segment (parent→child edge).

Returns a dict mapping probe labels like 'probe_seg_0' to (x, y, z) centers.

Source code in toric_spines_sim/geometry/swc.py
def get_center_coordinates_for_all_segments(
    swc_filepath: Union[str, Path],
    use_radius_weighting: bool = False,
) -> Dict[str, Tuple[float, float, float]]:
    """Compute the center coordinate of each SWC segment (parent→child edge).

    Returns a dict mapping probe labels like 'probe_seg_0' to (x, y, z) centers.
    """
    swc_filepath = Path(swc_filepath)
    logger.debug(
        "Preparing record points for all segments from %s (radius_weighting=%s)",
        swc_filepath,
        use_radius_weighting,
    )
    swc_model = SWCModel.from_swc_file(str(swc_filepath), validate_reconnections=False)
    node_ids = sorted(swc_model.nodes())
    if not node_ids:
        logger.warning("SWCModel has no nodes for %s", swc_filepath)
        return {}

    logger.debug("Loaded SWCModel with %d nodes from %s", len(node_ids), swc_filepath)

    record_points: Dict[str, Tuple[float, float, float]] = {}
    segment_index = 0
    skipped_root = 0
    for node_id in node_ids:
        parent_id = swc_model.parent_of(node_id)
        if parent_id is None:
            skipped_root += 1
            continue

        x0, y0, z0 = swc_model.get_node_xyz(parent_id)
        x1, y1, z1 = swc_model.get_node_xyz(node_id)

        if use_radius_weighting:
            r0 = swc_model.get_node_radius(parent_id)
            r1 = swc_model.get_node_radius(node_id)
            w0 = r0 * r0
            w1 = r1 * r1
            denom = w0 + w1
            if denom == 0.0:
                cx = 0.5 * (x0 + x1)
                cy = 0.5 * (y0 + y1)
                cz = 0.5 * (z0 + z1)
            else:
                cx = (w0 * x0 + w1 * x1) / denom
                cy = (w0 * y0 + w1 * y1) / denom
                cz = (w0 * z0 + w1 * z1) / denom
        else:
            cx = 0.5 * (x0 + x1)
            cy = 0.5 * (y0 + y1)
            cz = 0.5 * (z0 + z1)

        record_points[f"probe_seg_{segment_index}"] = (cx, cy, cz)
        if segment_index < 5:
            logger.debug(
                "probe_seg_%d (node=%s parent=%s) -> (%.6f, %.6f, %.6f)",
                segment_index,
                node_id,
                parent_id,
                cx,
                cy,
                cz,
            )
        segment_index += 1

    logger.info(
        "Prepared %d record points from %s (skipped roots=%d)",
        len(record_points),
        swc_filepath,
        skipped_root,
    )
    return record_points

arbor_locations_for_swc_nodes(swc_path, morphology, segment_tree, node_ids)

Map SWC sample IDs to Arbor (branch, pos) at that sample's endpoint.

A.load_swc_arbor makes one segment per non-root SWC sample, in file order. The distal end of segment k is the k-th non-root sample. Using this mapping (not place_pwlin.closest on XYZ) keeps colocated cycle-break / multi-neck copies on different branches.

Parameters:

Name Type Description Default
swc_path path - like

SWC used to build segment_tree.

required
morphology morphology

Morphology corresponding to segment_tree.

required
segment_tree segment_tree

Tree produced by arbor.load_swc_arbor.

required
node_ids iterable of int

SWC sample IDs to map.

required

Returns:

Type Description
dict[int, location]

Distal (branch, pos) for each requested sample.

Raises:

Type Description
ValueError

If sample/segment counts disagree or a node ID is missing.

Examples:

>>> locs = arbor_locations_for_swc_nodes(path, morph, tree, [12, 13])
>>> locs[12].branch != locs[13].branch  # colocated CYCLE_BREAK copies
True
Source code in toric_spines_sim/geometry/swc.py
def arbor_locations_for_swc_nodes(
    swc_path: Path,
    morphology: A.morphology,
    segment_tree: A.segment_tree,
    node_ids: Iterable[int],
) -> Dict[int, A.location]:
    """Map SWC sample IDs to Arbor ``(branch, pos)`` at that sample's endpoint.

    ``A.load_swc_arbor`` makes one segment per non-root SWC sample, in file
    order. The distal end of segment *k* is the *k*-th non-root sample. Using
    this mapping (not ``place_pwlin.closest`` on XYZ) keeps colocated cycle-break
    / multi-neck copies on different branches.

    Parameters
    ----------
    swc_path : path-like
        SWC used to build ``segment_tree``.
    morphology : arbor.morphology
        Morphology corresponding to ``segment_tree``.
    segment_tree : arbor.segment_tree
        Tree produced by ``arbor.load_swc_arbor``.
    node_ids : iterable of int
        SWC sample IDs to map.

    Returns
    -------
    dict[int, arbor.location]
        Distal ``(branch, pos)`` for each requested sample.

    Raises
    ------
    ValueError
        If sample/segment counts disagree or a node ID is missing.

    Examples
    --------
    >>> locs = arbor_locations_for_swc_nodes(path, morph, tree, [12, 13])  # doctest: +SKIP
    >>> locs[12].branch != locs[13].branch  # colocated CYCLE_BREAK copies
    True
    """
    wanted = set(int(i) for i in node_ids)
    samples = iter_swc_samples(swc_path)
    non_root = [s for s in samples if s[6] != -1]
    if len(non_root) != len(segment_tree.segments):
        raise ValueError(
            f"SWC {swc_path} has {len(non_root)} non-root samples but Arbor "
            f"segment tree has {len(segment_tree.segments)} segments"
        )

    node_to_seg = {sample[0]: i for i, sample in enumerate(non_root)}
    root_ids = {sample[0] for sample in samples if sample[6] == -1}
    branches = segment_tree_branches(segment_tree)
    seg_to_branch_local: Dict[int, Tuple[int, int]] = {}
    for bid, segs in enumerate(branches):
        for local, seg_idx in enumerate(segs):
            seg_to_branch_local[seg_idx] = (bid, local)

    def distal_location(seg_idx: int) -> A.location:
        bid, local = seg_to_branch_local[seg_idx]
        msegs = morphology.branch_segments(bid)
        lengths = [_segment_length(seg) for seg in msegs]
        total = float(sum(lengths))
        if total <= 0.0:
            pos = 1.0
        else:
            pos = float(sum(lengths[: local + 1])) / total
        return A.location(bid, pos)

    locations: Dict[int, A.location] = {}
    for node_id in wanted:
        if node_id in node_to_seg:
            locations[node_id] = distal_location(node_to_seg[node_id])
        elif node_id in root_ids:
            locations[node_id] = A.location(0, 0.0)
        else:
            raise ValueError(f"SWC node {node_id} not found in {swc_path}")
    return locations

sink_endpoint_location_from_swc_file(filepath)

Return XYZ of the distal sink sample referenced by end= in the header.

Expects a header line of the form::

# SINK: start=15, end=22, ...

Parameters:

Name Type Description Default
filepath path - like

SWC file with a # SINK: header.

required

Returns:

Type Description
tuple of float

(x, y, z) of the node whose id is end.

Raises:

Type Description
ValueError

If the # SINK: header or end= field is missing.

Examples:

>>> xyz = sink_endpoint_location_from_swc_file("TS1_wsink_r10um.swc")
>>> len(xyz)
3
Source code in toric_spines_sim/geometry/sink.py
def sink_endpoint_location_from_swc_file(
    filepath: Union[str, Path],
) -> Tuple[float, float, float]:
    """Return XYZ of the distal sink sample referenced by ``end=`` in the header.

    Expects a header line of the form::

        # SINK: start=15, end=22, ...

    Parameters
    ----------
    filepath : path-like
        SWC file with a ``# SINK:`` header.

    Returns
    -------
    tuple of float
        ``(x, y, z)`` of the node whose id is ``end``.

    Raises
    ------
    ValueError
        If the ``# SINK:`` header or ``end=`` field is missing.

    Examples
    --------
    >>> xyz = sink_endpoint_location_from_swc_file("TS1_wsink_r10um.swc")
    >>> len(xyz)
    3
    """
    fields = _parse_sink_header_fields(filepath)
    if "end" not in fields:
        raise ValueError(f"No end= field in SINK header of {filepath}")
    end_idx = int(fields["end"])

    from swctools import SWCModel

    swc_model = SWCModel.from_swc_file(filepath)
    return (
        swc_model.nodes[end_idx]["x"],
        swc_model.nodes[end_idx]["y"],
        swc_model.nodes[end_idx]["z"],
    )

neck_point_from_swc_file(filepath)

Read the neck point coordinates from the SINK header in an SWC file.

Expects a header line of the form::

# SINK: ..., neck_xyz=<x> <y> <z>

Parameters:

Name Type Description Default
filepath path - like

SWC file with a # SINK: header.

required

Returns:

Type Description
tuple of float

(x, y, z) of the neck attachment.

Raises:

Type Description
ValueError

If no # SINK: header or neck_xyz field is found.

Examples:

>>> xyz = neck_point_from_swc_file("TS1_wsink_r10um.swc")
>>> xyz
(12.3, 4.5, 6.7)
Source code in toric_spines_sim/geometry/sink.py
def neck_point_from_swc_file(filepath: Union[str, Path]) -> Tuple[float, float, float]:
    """Read the neck point coordinates from the SINK header in an SWC file.

    Expects a header line of the form::

        # SINK: ..., neck_xyz=<x> <y> <z>

    Parameters
    ----------
    filepath : path-like
        SWC file with a ``# SINK:`` header.

    Returns
    -------
    tuple of float
        ``(x, y, z)`` of the neck attachment.

    Raises
    ------
    ValueError
        If no ``# SINK:`` header or ``neck_xyz`` field is found.

    Examples
    --------
    >>> xyz = neck_point_from_swc_file("TS1_wsink_r10um.swc")
    >>> xyz  # doctest: +SKIP
    (12.3, 4.5, 6.7)
    """
    fields = _parse_sink_header_fields(filepath)
    if "neck_xyz" not in fields:
        raise ValueError(f"No neck_xyz found in SINK header of {filepath}")
    parts = fields["neck_xyz"].split()
    if len(parts) < 3:
        raise ValueError(
            f"neck_xyz field has fewer than 3 values: {fields['neck_xyz']}"
        )
    return (float(parts[0]), float(parts[1]), float(parts[2]))

optimal_sink_direction(neck_coords, swc, average_multiple=True)

Compute optimal sink direction(s) pointing away from the morphology.

Args: neck_coords: Either a single (x, y, z) tuple, or a path to a file containing one or more neck points (one per line: x y z). swc: SWC file path, SWCModel instance, or dict of node_id -> (x, y, z, r). average_multiple: If True and neck_coords is a file with multiple points, return the average direction. If False, return a list of directions (one per neck point).

Returns: Single (dx, dy, dz) direction tuple if average_multiple=True or single neck point. List of direction tuples if average_multiple=False and multiple neck points.

Source code in toric_spines_sim/geometry/sink.py
def optimal_sink_direction(
    neck_coords: Union[Tuple[float, float, float], str, Path],
    swc: Union[str, Path, dict],
    average_multiple: bool = True,
) -> Union[Tuple[float, float, float], List[Tuple[float, float, float]]]:
    """
    Compute optimal sink direction(s) pointing away from the morphology.

    Args:
        neck_coords: Either a single (x, y, z) tuple, or a path to a file containing
                  one or more neck points (one per line: x y z).
        swc: SWC file path, SWCModel instance, or dict of node_id -> (x, y, z, r).
        average_multiple: If True and neck_coords is a file with multiple points,
                         return the average direction. If False, return a list
                         of directions (one per neck point).

    Returns:
        Single (dx, dy, dz) direction tuple if average_multiple=True or single neck point.
        List of direction tuples if average_multiple=False and multiple neck points.
    """
    from swctools import SWCModel

    # Load SWC model points
    if isinstance(swc, (str, Path)):
        pts_dict = read_swc_points(Path(swc))
    elif isinstance(swc, SWCModel):
        # SWCModel inherits from networkx.DiGraph
        # Iterate over nodes using the NetworkX API
        pts_dict = {}
        for node_id in swc.nodes():
            node_data = swc.nodes[node_id]
            x = float(node_data["x"])
            y = float(node_data["y"])
            z = float(node_data["z"])
            r = float(node_data.get("r", node_data.get("radius", 0.0)))
            pts_dict[node_id] = (x, y, z, r)
    else:
        # Assume it's already a dict
        pts_dict = swc

    if not pts_dict:
        raise ValueError("SWC model has no points")

    # Load neck point(s)
    if isinstance(neck_coords, (str, Path)):
        # Load from file
        neck_points = load_xyz_points(neck_coords)
    else:
        # Single point provided as tuple
        neck_points = [neck_coords]

    if not neck_points:
        raise ValueError("No neck points provided")

    # Compute direction for each neck point
    directions = []
    for neck_xyz_single in neck_points:
        neck = np.asarray(neck_xyz_single, dtype=float)
        vectors = []
        for _nid, (x, y, z, _r) in pts_dict.items():
            vectors.append(np.array([x, y, z], dtype=float) - neck)
        mean_vec = np.mean(np.stack(vectors, axis=0), axis=0)
        away = -mean_vec
        n = float(np.linalg.norm(away))
        if n == 0.0:
            away = np.array([1.0, 0.0, 0.0], dtype=float)
        else:
            away = away / n
        directions.append((float(away[0]), float(away[1]), float(away[2])))

    # Return based on number of points and averaging preference
    if len(directions) == 1:
        return directions[0]

    if average_multiple:
        # Average all directions and normalize
        avg_direction = np.mean(np.array(directions), axis=0)
        n = float(np.linalg.norm(avg_direction))
        if n == 0.0:
            avg_direction = np.array([1.0, 0.0, 0.0], dtype=float)
        else:
            avg_direction = avg_direction / n
        return (
            float(avg_direction[0]),
            float(avg_direction[1]),
            float(avg_direction[2]),
        )
    else:
        return directions

append_sink_to_swc(swc_in, swc_out, neck_coords, geom, tag=5, last_segment_tag=None)

Append a cylindrical sink as a new tree and write # SINK: metadata.

The sink starts at the SWC node nearest neck_coords and extends along geom.axis. A connector frustum of length geom.connector_length links the neck to the first sink node. When last_segment_tag is set, the distal tip uses that tag (typically 6 for HH) instead of tag.

Parameters:

Name Type Description Default
swc_in path - like

Input SWC and destination path.

required
swc_out path - like

Input SWC and destination path.

required
neck_coords path-like or sequence of float

Neck XYZ, or a file with x y z [r].

required
geom SinkGeometry

Cylinder radius, length, axis, and segmentation.

required
tag int

SWC tag for sink body nodes (default 5).

5
last_segment_tag int

Tag for the distal tip node.

None

Returns:

Type Description
Path

swc_out.

Examples:

>>> geom = SinkGeometry(radius=0.5, length=100.0, axis="x")
>>> append_sink_to_swc("TS1.swc", "TS1_wsink.swc", (0, 0, 0), geom)
Source code in toric_spines_sim/geometry/sink.py
def append_sink_to_swc(
    swc_in: Union[str, Path],
    swc_out: Union[str, Path],
    neck_coords: Union[str, Path, Tuple[float, float, float], List[float]],
    geom: SinkGeometry,
    tag: int = 5,
    last_segment_tag: Optional[int] = None,
) -> Path:
    """Append a cylindrical sink as a new tree and write ``# SINK:`` metadata.

    The sink starts at the SWC node nearest ``neck_coords`` and extends along
    ``geom.axis``. A connector frustum of length ``geom.connector_length``
    links the neck to the first sink node. When ``last_segment_tag`` is set,
    the distal tip uses that tag (typically 6 for HH) instead of ``tag``.

    Parameters
    ----------
    swc_in, swc_out : path-like
        Input SWC and destination path.
    neck_coords : path-like or sequence of float
        Neck XYZ, or a file with ``x y z [r]``.
    geom : SinkGeometry
        Cylinder radius, length, axis, and segmentation.
    tag : int
        SWC tag for sink body nodes (default 5).
    last_segment_tag : int, optional
        Tag for the distal tip node.

    Returns
    -------
    pathlib.Path
        ``swc_out``.

    Examples
    --------
    >>> geom = SinkGeometry(radius=0.5, length=100.0, axis="x")
    >>> append_sink_to_swc("TS1.swc", "TS1_wsink.swc", (0, 0, 0), geom)  # doctest: +SKIP
    """
    swc_in = Path(swc_in)
    swc_out = Path(swc_out)
    if not swc_in.exists():
        raise FileNotFoundError(f"Input SWC not found: {swc_in}")

    neck_x, neck_y, neck_z, neck_r = _as_xyzr(neck_coords)
    neck_xyz = (neck_x, neck_y, neck_z)
    pts_dict = read_swc_points(swc_in)
    max_id = max(pts_dict.keys()) if pts_dict else 0

    # Snap to nearest existing SWC node so the sink connects exactly.
    snapped_xyz = neck_xyz
    neck_id = None
    if pts_dict:
        neck_id, snapped_xyz, _snapped_r = _snap_to_nearest_node(neck_xyz, pts_dict)
        logger.debug("snapped neck point to node %d at %s", neck_id, snapped_xyz)

    # Build sink nodes (n_cylinders + 1 nodes defining n_cylinders frusta)
    pts = _gen_sink_points(snapped_xyz, geom)
    radius_default = float(geom.radius)

    # Prepare new SWC lines
    logger.debug(
        "building %d sink nodes (defining %d cylinders):", len(pts), geom.n_cylinders
    )
    new_lines = []
    nid = max_id
    sink_start_id = max_id + 1
    logger.debug("sink_start_id: %d", sink_start_id)
    parent = neck_id
    for idx, p in enumerate(pts):
        nid += 1
        x, y, z = p
        node_tag = (
            last_segment_tag
            if last_segment_tag is not None and idx == len(pts) - 1
            else tag
        )
        new_lines.append(
            f"{nid} {node_tag} {x:.6f} {y:.6f} {z:.6f} {radius_default:.6f} {parent}\n"
        )
        parent = nid
        logger.debug("  nid=%d: %s", nid, p)
    sink_end_id = nid
    logger.debug("sink_end_id: %d", sink_end_id)

    header_tag_fields = f"tag={tag}"
    if last_segment_tag is not None:
        header_tag_fields += f", last_segment_tag={last_segment_tag}"
    header = [
        f"# SINK: start={sink_start_id}, end={sink_end_id}, axis={geom.axis}, segments={geom.n_cylinders+1}, nodes={len(pts)}, length={geom.length}, radius={geom.radius}, {header_tag_fields}, neck_xyz={snapped_xyz[0]:.6f} {snapped_xyz[1]:.6f} {snapped_xyz[2]:.6f}\n",
    ]

    _write_swc_with_header(swc_in, swc_out, header, new_lines)
    return swc_out

append_sink_to_swc_multi_neck_points(swc_in, swc_out, neck_points, geom, tag=5, last_segment_tag=None)

Append one sink, then extra-neck copies of the sink start node.

neck_points is a file of x y z rows. The first row is the primary neck; each later row gets a copy of the sink-start sample parented at that neck, recorded as # MULTI_NECK reconnect i j.

Source code in toric_spines_sim/geometry/sink.py
def append_sink_to_swc_multi_neck_points(
    swc_in: Union[str, Path],
    swc_out: Union[str, Path],
    neck_points: Union[str, Path],
    geom: SinkGeometry,
    tag: int = 5,
    last_segment_tag: Optional[int] = None,
) -> Path:
    """Append one sink, then extra-neck copies of the sink start node.

    ``neck_points`` is a file of ``x y z`` rows. The first row is the primary
    neck; each later row gets a copy of the sink-start sample parented at that
    neck, recorded as ``# MULTI_NECK reconnect i j``.
    """
    swc_in = Path(swc_in)
    swc_out = Path(swc_out)
    if not swc_in.exists():
        raise FileNotFoundError(f"Input SWC not found: {swc_in}")

    neck_xyzs = load_xyz_points(neck_points)
    if len(neck_xyzs) == 0:
        raise ValueError("neck point file contained no points")

    pts_dict = read_swc_points(swc_in)
    max_id = max(pts_dict.keys()) if pts_dict else 0

    # Snap primary neck point
    primary_neck_xyz = neck_xyzs[0]
    snapped_xyz = primary_neck_xyz
    primary_neck_id = None
    if pts_dict:
        primary_neck_id, snapped_xyz, _snapped_r = _snap_to_nearest_node(
            primary_neck_xyz, pts_dict
        )
        logger.debug(
            "snapped primary neck point to node %d at %s", primary_neck_id, snapped_xyz
        )

    pts = _gen_sink_points(snapped_xyz, geom)
    radius_default = float(geom.radius)
    sink_start_xyz = pts[0]

    new_lines: List[str] = []
    nid = max_id
    sink_start_id = max_id + 1
    parent = primary_neck_id if primary_neck_id is not None else -1
    for idx, p in enumerate(pts):
        nid += 1
        x, y, z = p
        node_tag = (
            last_segment_tag
            if last_segment_tag is not None and idx == len(pts) - 1
            else tag
        )
        new_lines.append(
            f"{nid} {node_tag} {x:.6f} {y:.6f} {z:.6f} {radius_default:.6f} {parent}\n"
        )
        parent = nid
    sink_end_id = nid

    # For each additional neck point, connect the (snapped) neck node to a *copy* of the
    # sink start node. We then annotate the SWC header with a directive that downstream
    # code can interpret as a direct connection (e.g. a gap junction) between the true
    # sink start node and its copies.
    reconnect_pairs: List[Tuple[int, int]] = []
    for neck_xyz in neck_xyzs[1:]:
        branch_neck_id = None
        branch_xyz = neck_xyz
        if pts_dict:
            branch_neck_id, branch_xyz, _branch_r = _snap_to_nearest_node(
                neck_xyz, pts_dict
            )

        if (
            abs(branch_xyz[0] - snapped_xyz[0]) < 1e-12
            and abs(branch_xyz[1] - snapped_xyz[1]) < 1e-12
            and abs(branch_xyz[2] - snapped_xyz[2]) < 1e-12
        ):
            continue

        # Copy of the sink start node: colocated with sink start, but parented to this neck.
        nid += 1
        x, y, z = sink_start_xyz
        sink_copy_parent = branch_neck_id if branch_neck_id is not None else -1
        new_lines.append(
            f"{nid} {tag} {x:.6f} {y:.6f} {z:.6f} {radius_default:.6f} {sink_copy_parent}\n"
        )
        reconnect_pairs.append((sink_start_id, nid))

    header_tag_fields = f"tag={tag}"
    if last_segment_tag is not None:
        header_tag_fields += f", last_segment_tag={last_segment_tag}"
    header = [
        f"# SINK: start={sink_start_id}, end={sink_end_id}, axis={geom.axis}, segments={geom.n_cylinders+1}, nodes={len(pts)}, length={geom.length}, radius={geom.radius}, {header_tag_fields}, n_necks={len(neck_xyzs)}, neck_xyz={snapped_xyz[0]:.6f} {snapped_xyz[1]:.6f} {snapped_xyz[2]:.6f}\n",
    ]
    for id_a, id_b in reconnect_pairs:
        header.append(f"# MULTI_NECK reconnect {id_a} {id_b}\n")

    _write_swc_with_header(swc_in, swc_out, header, new_lines)
    return swc_out

build_spiny_dendrite(params)

Build an in-memory spiny dendrite from explicit parameters (mode A).

Source code in toric_spines_sim/geometry/dendrite.py
def build_spiny_dendrite(params: SpinyDendriteParams) -> SpinyDendriteMorphology:
    """Build an in-memory spiny dendrite from explicit parameters (mode A)."""
    params.validate()
    direction = _axis_unit_vector(params.axis)
    u_hat, v_hat = _perpendicular_basis(direction)

    n_attach = _n_attach_nodes(params.n_spines, params.max_spines_per_node)
    # Trunk: neck (no spines) + n_attach attachment nodes (includes tip).
    n_trunk = 1 + n_attach
    spines_per_attach = distribute_spine_counts(
        params.n_spines,
        n_attach,
        params.max_spines_per_node,
        distribution=params.distribution,
        seed=params.seed,
    )

    nodes: List[SWCNodeRecord] = []
    trunk_node_ids: List[int] = []
    next_id = 1

    # Proximal neck at origin.
    neck_xyz = (0.0, 0.0, 0.0)
    nodes.append(
        SWCNodeRecord(
            node_id=next_id,
            tag=params.trunk_tag,
            x=neck_xyz[0],
            y=neck_xyz[1],
            z=neck_xyz[2],
            radius=params.trunk_neck_radius,
            parent=-1,
        )
    )
    trunk_node_ids.append(next_id)
    next_id += 1

    # Remaining trunk nodes evenly spaced along the axis.
    for i in range(1, n_trunk):
        t = i / (n_trunk - 1)
        pos = t * params.length * direction
        radius = (1.0 - t) * params.trunk_neck_radius + t * params.trunk_tip_radius
        parent_id = trunk_node_ids[-1]
        nodes.append(
            SWCNodeRecord(
                node_id=next_id,
                tag=params.trunk_tag,
                x=float(pos[0]),
                y=float(pos[1]),
                z=float(pos[2]),
                radius=float(radius),
                parent=parent_id,
            )
        )
        trunk_node_ids.append(next_id)
        next_id += 1

    attach_trunk_ids = trunk_node_ids[1:]  # skip neck
    az_points: List[Tuple[float, float, float]] = []
    spine_head_node_ids: List[int] = []

    neck_seg_len = params.spine_length * params.spine_neck_length_fraction
    for trunk_id, m_spines in zip(attach_trunk_ids, spines_per_attach):
        if m_spines == 0:
            continue
        trunk_node = next(n for n in nodes if n.node_id == trunk_id)
        trunk_pos = np.array([trunk_node.x, trunk_node.y, trunk_node.z], dtype=float)
        for k in range(m_spines):
            theta = params.azimuth0 + (2.0 * math.pi * k) / m_spines
            spine_dir = math.cos(theta) * u_hat + math.sin(theta) * v_hat
            neck_pos = trunk_pos + neck_seg_len * spine_dir
            head_pos = trunk_pos + params.spine_length * spine_dir

            neck_id = next_id
            nodes.append(
                SWCNodeRecord(
                    node_id=neck_id,
                    tag=params.spine_neck_tag,
                    x=float(neck_pos[0]),
                    y=float(neck_pos[1]),
                    z=float(neck_pos[2]),
                    radius=params.spine_neck_radius,
                    parent=trunk_id,
                )
            )
            next_id += 1

            head_id = next_id
            head_xyz = (
                float(head_pos[0]),
                float(head_pos[1]),
                float(head_pos[2]),
            )
            nodes.append(
                SWCNodeRecord(
                    node_id=head_id,
                    tag=params.spine_head_tag,
                    x=head_xyz[0],
                    y=head_xyz[1],
                    z=head_xyz[2],
                    radius=params.spine_head_radius,
                    parent=neck_id,
                )
            )
            next_id += 1
            spine_head_node_ids.append(head_id)
            az_points.append(head_xyz)

    morph = SpinyDendriteMorphology(
        nodes=nodes,
        az_points=az_points,
        neck_point=neck_xyz,
        trunk_node_ids=trunk_node_ids,
        spine_head_node_ids=spine_head_node_ids,
        spines_per_attach_node=spines_per_attach,
        params=params,
    )
    morph.diagnostics = {
        "surface_area": morph.surface_area(),
        "volume": morph.volume(),
        "n_trunk_nodes": float(n_trunk),
        "n_attach_nodes": float(n_attach),
        "n_spines": float(params.n_spines),
    }
    logger.debug(
        "Built spiny dendrite: %d trunk nodes, %d spines, SA=%.6f",
        n_trunk,
        params.n_spines,
        morph.diagnostics["surface_area"],
    )
    return morph

build_from_toric_spine(swc_path, az_path, match=None, neck_path=None)

Build a comparable spiny dendrite matched to a toric spine subsystem.

Hard constraints: n_spines == n_AZ and lateral surface area ≈ TS SA (excluding sink-tagged segments).

Scale strategies (match.scale_strategy):

  • relative (default): scale a typical-proportion template uniformly so SA matches (preserves relative geometry).
  • absolute_spines: keep literature absolute spine sizes; solve trunk length so SA matches.

Parameters:

Name Type Description Default
swc_path path - like

Toric-spine SWC (sink tags excluded from SA).

required
az_path path - like

Active-zone XYZ file; one spine per point.

required
match ToricSpineMatchParams

Template proportions and scale strategy.

None
neck_path path - like

Recorded in diagnostics only (geometry uses a single proximal neck).

None

Returns:

Type Description
SpinyDendriteMorphology

Synthetic morphology plus SA/volume diagnostics.

Examples:

>>> morph = build_from_toric_spine("TS1.swc", "TS1_AZ.txt")
>>> morph.diagnostics["surface_area_rel_error"] < 0.05
True
Source code in toric_spines_sim/geometry/dendrite.py
def build_from_toric_spine(
    swc_path: Union[str, Path],
    az_path: Union[str, Path],
    match: Optional[ToricSpineMatchParams] = None,
    neck_path: Optional[Union[str, Path]] = None,
) -> SpinyDendriteMorphology:
    """Build a comparable spiny dendrite matched to a toric spine subsystem.

    Hard constraints: ``n_spines == n_AZ`` and lateral surface area ≈ TS SA
    (excluding sink-tagged segments).

    Scale strategies (``match.scale_strategy``):

    - ``relative`` (default): scale a typical-proportion template uniformly so
      SA matches (preserves relative geometry).
    - ``absolute_spines``: keep literature absolute spine sizes; solve trunk
      length so SA matches.

    Parameters
    ----------
    swc_path : path-like
        Toric-spine SWC (sink tags excluded from SA).
    az_path : path-like
        Active-zone XYZ file; one spine per point.
    match : ToricSpineMatchParams, optional
        Template proportions and scale strategy.
    neck_path : path-like, optional
        Recorded in diagnostics only (geometry uses a single proximal neck).

    Returns
    -------
    SpinyDendriteMorphology
        Synthetic morphology plus SA/volume diagnostics.

    Examples
    --------
    >>> morph = build_from_toric_spine("TS1.swc", "TS1_AZ.txt")  # doctest: +SKIP
    >>> morph.diagnostics["surface_area_rel_error"] < 0.05
    True
    """
    match = match or ToricSpineMatchParams()
    swc_path = Path(swc_path)
    az_path = Path(az_path)

    n_spines = count_points_file(az_path)
    target_sa = swc_subsystem_surface_area(swc_path, sink_tags=match.sink_tags)
    target_volume = swc_subsystem_volume(swc_path, sink_tags=match.sink_tags)

    if target_sa <= 0:
        raise ValueError(f"Target surface area from {swc_path} is non-positive")

    if match.scale_strategy == "relative":
        template_params = _params_from_template(match, n_spines, scale=1.0)
        template_morph = build_spiny_dendrite(template_params)
        template_sa = template_morph.surface_area()
        if template_sa <= 0:
            raise ValueError("Template morphology has non-positive surface area")
        scale = math.sqrt(target_sa / template_sa)
        params = _params_from_template(match, n_spines, scale=scale)
        morph = build_spiny_dendrite(params)
        scale_factor = scale
    elif match.scale_strategy == "absolute_spines":
        params, morph = _solve_trunk_length_for_sa(match, n_spines, target_sa)
        scale_factor = 1.0
    else:
        raise ValueError(f"Unknown scale_strategy: {match.scale_strategy!r}")

    achieved_sa = morph.surface_area()
    achieved_vol = morph.volume()
    morph.diagnostics.update(
        {
            "target_surface_area": target_sa,
            "achieved_surface_area": achieved_sa,
            "surface_area_rel_error": abs(achieved_sa - target_sa) / target_sa,
            "target_volume": target_volume,
            "achieved_volume": achieved_vol,
            "scale_factor": scale_factor,
            "spine_length": morph.params.spine_length,
            "trunk_length": morph.params.length,
            "literature_spine_length": _DEFAULT_SPINE_LENGTH,
            "spine_length_over_literature": morph.params.spine_length
            / _DEFAULT_SPINE_LENGTH,
        }
    )
    if neck_path is not None:
        # Record only; geometry uses a single proximal neck at the origin.
        morph.diagnostics["source_n_neck_points"] = float(count_points_file(neck_path))

    logger.info(
        "Matched toric spine %s: n_spines=%d, SA target=%.6f achieved=%.6f "
        "(rel err=%.3e), strategy=%s",
        swc_path,
        n_spines,
        target_sa,
        achieved_sa,
        match.scale_strategy,
    )
    return morph

write_subsystem(morph, swc_path, az_path, neck_path)

Write SWC + AZ + neckpoint files for a subsystem.

Source code in toric_spines_sim/geometry/dendrite.py
def write_subsystem(
    morph: SpinyDendriteMorphology,
    swc_path: Union[str, Path],
    az_path: Union[str, Path],
    neck_path: Union[str, Path],
) -> Tuple[Path, Path, Path]:
    """Write SWC + AZ + neckpoint files for a subsystem."""
    swc_out = write_swc(swc_path, morph)
    az_out = write_xyz_points(az_path, morph.az_points)
    neck_out = write_xyz_points(neck_path, [morph.neck_point])
    logger.info(
        "Wrote subsystem: swc=%s az=%s neck=%s (%d synapses)",
        swc_out,
        az_out,
        neck_out,
        morph.n_spines,
    )
    return swc_out, az_out, neck_out

write_swc(path, morph, extra_header=None)

Write morphology nodes as an SWC file.

Source code in toric_spines_sim/geometry/dendrite.py
def write_swc(
    path: Union[str, Path],
    morph: SpinyDendriteMorphology,
    extra_header: Optional[Sequence[str]] = None,
) -> Path:
    """Write morphology nodes as an SWC file."""
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    p = morph.params
    header = [
        "# generated by toric_spines_sim.geometry.dendrite\n",
        (
            f"# SPINY_DENDRITE: n_spines={p.n_spines}, "
            f"length={p.length}, trunk_neck_radius={p.trunk_neck_radius}, "
            f"trunk_tip_radius={p.trunk_tip_radius}, spine_length={p.spine_length}, "
            f"max_spines_per_node={p.max_spines_per_node}, "
            f"distribution={p.distribution}\n"
        ),
    ]
    if extra_header:
        for line in extra_header:
            header.append(line if line.endswith("\n") else line + "\n")

    body = [
        (
            f"{n.node_id} {n.tag} {n.x:.6f} {n.y:.6f} {n.z:.6f} "
            f"{n.radius:.6f} {n.parent}\n"
        )
        for n in morph.nodes
    ]
    path.write_text("".join(header + body), encoding="utf-8")
    return path

write_xyz_points(path, points)

Write whitespace-delimited XYZ points (%.6f), one per line.

Source code in toric_spines_sim/geometry/dendrite.py
def write_xyz_points(
    path: Union[str, Path],
    points: Sequence[Tuple[float, float, float]],
) -> Path:
    """Write whitespace-delimited XYZ points (``%.6f``), one per line."""
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    lines = [f"{x:.6f} {y:.6f} {z:.6f}\n" for x, y, z in points]
    path.write_text("".join(lines), encoding="utf-8")
    return path

distribute_spine_counts(n_spines, n_attach, max_spines_per_node, distribution='even', seed=None)

Distribute n_spines across n_attach trunk nodes (each ≤ max).

Source code in toric_spines_sim/geometry/dendrite.py
def distribute_spine_counts(
    n_spines: int,
    n_attach: int,
    max_spines_per_node: int,
    distribution: DistributionMode = "even",
    seed: Optional[int] = None,
) -> List[int]:
    """Distribute *n_spines* across *n_attach* trunk nodes (each ≤ max)."""
    if n_attach < 1:
        raise ValueError("n_attach must be >= 1")
    if n_spines < 0:
        raise ValueError("n_spines must be >= 0")
    if n_spines > n_attach * max_spines_per_node:
        raise ValueError(
            f"Cannot place {n_spines} spines on {n_attach} nodes with "
            f"max_spines_per_node={max_spines_per_node}"
        )
    if n_spines == 0:
        return [0] * n_attach

    if distribution == "even":
        base = n_spines // n_attach
        rem = n_spines % n_attach
        counts = [base + (1 if i < rem else 0) for i in range(n_attach)]
        if any(c > max_spines_per_node for c in counts):
            raise ValueError(
                "even distribution exceeds max_spines_per_node; "
                "increase n_attach or max_spines_per_node"
            )
        return counts

    if distribution == "random":
        rng = np.random.default_rng(seed)
        counts = [0] * n_attach
        for _ in range(n_spines):
            candidates = [i for i, c in enumerate(counts) if c < max_spines_per_node]
            if not candidates:
                raise ValueError("no capacity left while assigning spines")
            choice = int(rng.choice(candidates))
            counts[choice] += 1
        return counts

    raise ValueError(f"Unknown distribution: {distribution!r}")

frustum_lateral_area(r1, r2, length)

Lateral surface area of a conical frustum (excluding end caps).

Source code in toric_spines_sim/geometry/dendrite.py
def frustum_lateral_area(r1: float, r2: float, length: float) -> float:
    """Lateral surface area of a conical frustum (excluding end caps)."""
    slant = math.sqrt((r1 - r2) ** 2 + length**2)
    return math.pi * (r1 + r2) * slant

morphology_surface_area(nodes)

Sum lateral frustum areas over parent→child edges.

Source code in toric_spines_sim/geometry/dendrite.py
def morphology_surface_area(nodes: Sequence[SWCNodeRecord]) -> float:
    """Sum lateral frustum areas over parent→child edges."""
    by_id = {n.node_id: n for n in nodes}
    total = 0.0
    for node in nodes:
        if node.parent < 0:
            continue
        parent = by_id[node.parent]
        length = math.dist(
            (parent.x, parent.y, parent.z), (node.x, node.y, node.z)
        )
        total += frustum_lateral_area(parent.radius, node.radius, length)
    return total

morphology_volume(nodes)

Sum frustum volumes over parent→child edges.

Source code in toric_spines_sim/geometry/dendrite.py
def morphology_volume(nodes: Sequence[SWCNodeRecord]) -> float:
    """Sum frustum volumes over parent→child edges."""
    by_id = {n.node_id: n for n in nodes}
    total = 0.0
    for node in nodes:
        if node.parent < 0:
            continue
        parent = by_id[node.parent]
        length = math.dist(
            (parent.x, parent.y, parent.z), (node.x, node.y, node.z)
        )
        total += frustum_volume(parent.radius, node.radius, length)
    return total

swc_subsystem_surface_area(swc_path, sink_tags=None)

Lateral SA of an SWC, excluding segments whose child tag is a sink tag.

Source code in toric_spines_sim/geometry/dendrite.py
def swc_subsystem_surface_area(
    swc_path: Union[str, Path],
    sink_tags: Optional[Sequence[int]] = None,
) -> float:
    """Lateral SA of an SWC, excluding segments whose child tag is a sink tag."""
    sink = set(sink_tags) if sink_tags is not None else set(_DEFAULT_SINK_TAGS)
    total = 0.0
    for prox, dist, child_tag in _parse_swc_segments_with_tags(Path(swc_path)):
        if child_tag in sink:
            continue
        length = math.dist(prox[:3], dist[:3])
        total += frustum_lateral_area(prox[3], dist[3], length)
    return total

fit_swc(mesh_path, polylines_path, swc_path, *, max_edge_length_frac=0.08, radius_strategy='equivalent_area', scale_radii=True, basis_optimize=False, basis_optimizer_options=None, scale_metric='surface_area')

Fit a cable SWC to a mesh + skeleton with mascaf.

Parameters:

Name Type Description Default
mesh_path PathLike

Mesh file path, or a bare filename under data/mesh/.

required
polylines_path PathLike

Skeleton polylines text file.

required
swc_path PathLike

Destination SWC path.

required
max_edge_length_frac float

FitOptions.max_edge_length as a fraction of the mesh bounding-box diagonal (default 0.08).

0.08
radius_strategy str

Radius estimation strategy (default "equivalent_area").

'equivalent_area'
scale_radii bool

If True, call scale_radii_to_match_mesh before export.

True
basis_optimize bool

If True, enable mascaf BasisOptimizer via FitOptions.

False
basis_optimizer_options Optional[dict[str, Any]]

Optional kwargs for BasisOptimizerOptions when basis_optimize is True.

None
scale_metric str

Metric for radius scaling (default "surface_area").

'surface_area'

Returns:

Type Description
Path

Path to the written SWC file.

Source code in toric_spines_sim/geometry/mesh_pipeline.py
def fit_swc(
    mesh_path: PathLike,
    polylines_path: PathLike,
    swc_path: PathLike,
    *,
    max_edge_length_frac: float = 0.08,
    radius_strategy: str = "equivalent_area",
    scale_radii: bool = True,
    basis_optimize: bool = False,
    basis_optimizer_options: Optional[dict[str, Any]] = None,
    scale_metric: str = "surface_area",
) -> Path:
    """Fit a cable SWC to a mesh + skeleton with mascaf.

    Parameters
    ----------
    mesh_path
        Mesh file path, or a bare filename under ``data/mesh/``.
    polylines_path
        Skeleton polylines text file.
    swc_path
        Destination SWC path.
    max_edge_length_frac
        ``FitOptions.max_edge_length`` as a fraction of the mesh bounding-box
        diagonal (default ``0.08``).
    radius_strategy
        Radius estimation strategy (default ``\"equivalent_area\"``).
    scale_radii
        If True, call ``scale_radii_to_match_mesh`` before export.
    basis_optimize
        If True, enable mascaf ``BasisOptimizer`` via ``FitOptions``.
    basis_optimizer_options
        Optional kwargs for ``BasisOptimizerOptions`` when ``basis_optimize``
        is True.
    scale_metric
        Metric for radius scaling (default ``\"surface_area\"``).

    Returns
    -------
    Path
        Path to the written SWC file.
    """
    (
        BasisOptimizerOptions,
        CableFitter,
        FitOptions,
        MeshManager,
        SkeletonGraph,
    ) = _require_mascaf()

    mesh_path = resolve_mesh_path(mesh_path)
    polylines_path = Path(polylines_path)
    swc_path = Path(swc_path)
    if not polylines_path.is_file():
        raise FileNotFoundError(f"Polylines file not found: {polylines_path}")

    swc_path.parent.mkdir(parents=True, exist_ok=True)

    optimizer = None
    if basis_optimize:
        opts = dict(basis_optimizer_options or {})
        optimizer = BasisOptimizerOptions(**opts)

    mesh_mgr = MeshManager(mesh_path=str(mesh_path))
    bbox_diag = float(mesh_mgr.bounding_box_diagonal())
    max_edge_length = float(max_edge_length_frac) * bbox_diag

    fit_options = FitOptions(
        max_edge_length=max_edge_length,
        radius_strategy=radius_strategy,
        basis_optimizer_options=optimizer,
    )

    logger.info(
        "Fitting SWC from mesh %s and skeleton %s "
        "(max_edge_length_frac=%s → max_edge_length=%s, bbox_diag=%s, "
        "radius_strategy=%s, scale_radii=%s, basis_optimize=%s)",
        mesh_path,
        polylines_path,
        max_edge_length_frac,
        max_edge_length,
        bbox_diag,
        radius_strategy,
        scale_radii,
        basis_optimize,
    )
    skeleton = SkeletonGraph.from_txt(str(polylines_path))
    morphology = CableFitter(fit_options).fit(mesh_mgr, skeleton)

    if scale_radii:
        morphology.scale_radii_to_match_mesh(mesh_mgr, metric=scale_metric)

    morphology.to_swc_file(str(swc_path))
    logger.info("Wrote SWC to %s", swc_path)
    return swc_path.resolve()

list_ts_meshes()

Return sorted TS*.obj paths under data/mesh/.

Source code in toric_spines_sim/geometry/mesh_pipeline.py
def list_ts_meshes() -> list[Path]:
    """Return sorted ``TS*.obj`` paths under ``data/mesh/``."""
    return sorted(MESH_DIR.glob("TS*.obj"))

mesh_to_swc(mesh_path, *, polylines_path=None, swc_path=None, skip_skeletonize=False, polylines_only=False, profile='auto', max_edge_length_frac=0.08, radius_strategy='equivalent_area', scale_radii=True, basis_optimize=False, basis_optimizer_options=None, scale_metric='surface_area', **skeletonize_kwargs)

Run mesh → polylines → SWC (or a subset of those steps).

Default outputs are data/skeletons/<spine_id>.polylines.txt and data/swc/pixels/<spine_id>.swc.

Parameters:

Name Type Description Default
mesh_path path - like

Mesh file or spine id (TS1).

required
polylines_path path - like

Override default output paths.

None
swc_path path - like

Override default output paths.

None
skip_skeletonize bool

If True, reuse existing polylines (must already exist).

False
polylines_only bool

Stop after skeletonization; swc_path in the result is None.

False
profile str

pymcfs skeletonize profile (auto or a named profile).

'auto'
max_edge_length_frac float

Forwarded to fit_swc.

0.08
radius_strategy float

Forwarded to fit_swc.

0.08
scale_radii float

Forwarded to fit_swc.

0.08
basis_optimize float

Forwarded to fit_swc.

0.08
skeletonize_kwargs Any

Extra kwargs for skeletonize_mesh.

{}

Returns:

Type Description
MeshToSwcResult

Resolved mesh / polylines / SWC paths.

Examples:

>>> mesh_to_swc("TS1", skip_skeletonize=True)
>>> mesh_to_swc("TS1.obj", polylines_only=True)
Source code in toric_spines_sim/geometry/mesh_pipeline.py
def mesh_to_swc(
    mesh_path: PathLike,
    *,
    polylines_path: Optional[PathLike] = None,
    swc_path: Optional[PathLike] = None,
    skip_skeletonize: bool = False,
    polylines_only: bool = False,
    profile: str = "auto",
    max_edge_length_frac: float = 0.08,
    radius_strategy: str = "equivalent_area",
    scale_radii: bool = True,
    basis_optimize: bool = False,
    basis_optimizer_options: Optional[dict[str, Any]] = None,
    scale_metric: str = "surface_area",
    **skeletonize_kwargs: Any,
) -> MeshToSwcResult:
    """Run mesh → polylines → SWC (or a subset of those steps).

    Default outputs are ``data/skeletons/<spine_id>.polylines.txt`` and
    ``data/swc/pixels/<spine_id>.swc``.

    Parameters
    ----------
    mesh_path : path-like
        Mesh file or spine id (``TS1``).
    polylines_path, swc_path : path-like, optional
        Override default output paths.
    skip_skeletonize : bool
        If True, reuse existing polylines (must already exist).
    polylines_only : bool
        Stop after skeletonization; ``swc_path`` in the result is ``None``.
    profile : str
        pymcfs skeletonize profile (``auto`` or a named profile).
    max_edge_length_frac, radius_strategy, scale_radii, basis_optimize
        Forwarded to ``fit_swc``.
    skeletonize_kwargs
        Extra kwargs for ``skeletonize_mesh``.

    Returns
    -------
    MeshToSwcResult
        Resolved mesh / polylines / SWC paths.

    Examples
    --------
    >>> mesh_to_swc("TS1", skip_skeletonize=True)  # doctest: +SKIP
    >>> mesh_to_swc("TS1.obj", polylines_only=True)  # doctest: +SKIP
    """
    mesh_path = resolve_mesh_path(mesh_path)
    output_polylines_path = (
        Path(polylines_path)
        if polylines_path is not None
        else default_polylines_path(mesh_path)
    )
    output_swc_path = Path(swc_path) if swc_path is not None else default_swc_path(mesh_path)

    SKELETONS_DIR.mkdir(parents=True, exist_ok=True)
    SWC_PIXELS_DIR.mkdir(parents=True, exist_ok=True)

    if not skip_skeletonize:
        skeletonize_mesh(
            mesh_path,
            output_polylines_path,
            profile=profile,
            **skeletonize_kwargs,
        )
    elif not output_polylines_path.is_file():
        raise FileNotFoundError(
            f"skip_skeletonize=True but polylines not found: {output_polylines_path}"
        )

    if polylines_only:
        return MeshToSwcResult(
            mesh_path=mesh_path,
            polylines_path=output_polylines_path.resolve(),
            swc_path=None,
        )

    written_swc = fit_swc(
        mesh_path,
        output_polylines_path,
        output_swc_path,
        max_edge_length_frac=max_edge_length_frac,
        radius_strategy=radius_strategy,
        scale_radii=scale_radii,
        basis_optimize=basis_optimize,
        basis_optimizer_options=basis_optimizer_options,
        scale_metric=scale_metric,
    )
    return MeshToSwcResult(
        mesh_path=mesh_path,
        polylines_path=output_polylines_path.resolve(),
        swc_path=written_swc,
    )

resolve_mesh_targets(meshes=None, *, all_meshes=False)

Resolve CLI mesh targets.

Parameters:

Name Type Description Default
meshes Optional[Sequence[str]]

Explicit mesh names or paths. Ignored when all_meshes is True.

None
all_meshes bool

If True, return every TS*.obj under data/mesh/.

False

Raises:

Type Description
ValueError

If neither --all nor any mesh arguments were provided, or if --all is combined with explicit mesh arguments.

FileNotFoundError

If a requested mesh cannot be resolved, or if --all finds none.

Source code in toric_spines_sim/geometry/mesh_pipeline.py
def resolve_mesh_targets(
    meshes: Optional[Sequence[str]] = None,
    *,
    all_meshes: bool = False,
) -> list[Path]:
    """Resolve CLI mesh targets.

    Parameters
    ----------
    meshes
        Explicit mesh names or paths. Ignored when ``all_meshes`` is True.
    all_meshes
        If True, return every ``TS*.obj`` under ``data/mesh/``.

    Raises
    ------
    ValueError
        If neither ``--all`` nor any mesh arguments were provided, or if
        ``--all`` is combined with explicit mesh arguments.
    FileNotFoundError
        If a requested mesh cannot be resolved, or if ``--all`` finds none.
    """
    if all_meshes and meshes:
        raise ValueError("Pass either --all or explicit mesh arguments, not both")
    if all_meshes:
        targets = list_ts_meshes()
        if not targets:
            raise FileNotFoundError(f"No TS*.obj meshes found under {MESH_DIR}")
        return [p.resolve() for p in targets]
    if not meshes:
        raise ValueError("Provide one or more mesh arguments, or pass --all")
    return [resolve_mesh_path(m) for m in meshes]

skeletonize_mesh(mesh_path, polylines_path, *, profile='auto', branching='sparse', **skeletonize_kwargs)

Skeletonize a closed triangle mesh with pymcfs and write polylines.

Defaults match pymcfs toric-spine batch settings (profile="auto", branching="sparse", tip extension on). See :data:TORIC_SPINES_SKELETONIZE_DEFAULTS.

Parameters:

Name Type Description Default
mesh_path PathLike

Mesh file path, or a bare filename under data/mesh/.

required
polylines_path PathLike

Destination .polylines.txt path.

required
profile str

pymcfs skeletonization profile (default "auto" for TS meshes).

'auto'
branching str

Branching preference when profile="auto" (default "sparse").

'sparse'
**skeletonize_kwargs Any

Forwarded to pymcfs.skeletonize (overrides toric-spine defaults).

{}

Returns:

Type Description
Path

Path to the written polylines file.

Source code in toric_spines_sim/geometry/mesh_pipeline.py
def skeletonize_mesh(
    mesh_path: PathLike,
    polylines_path: PathLike,
    *,
    profile: str = "auto",
    branching: str = "sparse",
    **skeletonize_kwargs: Any,
) -> Path:
    """Skeletonize a closed triangle mesh with pymcfs and write polylines.

    Defaults match pymcfs toric-spine batch settings
    (``profile=\"auto\"``, ``branching=\"sparse\"``, tip extension on).
    See :data:`TORIC_SPINES_SKELETONIZE_DEFAULTS`.

    Parameters
    ----------
    mesh_path
        Mesh file path, or a bare filename under ``data/mesh/``.
    polylines_path
        Destination ``.polylines.txt`` path.
    profile
        pymcfs skeletonization profile (default ``\"auto\"`` for TS meshes).
    branching
        Branching preference when ``profile=\"auto\"`` (default ``\"sparse\"``).
    **skeletonize_kwargs
        Forwarded to ``pymcfs.skeletonize`` (overrides toric-spine defaults).

    Returns
    -------
    Path
        Path to the written polylines file.
    """
    load_and_repair, skeletonize = _require_pymcfs()

    mesh_path = resolve_mesh_path(mesh_path)
    polylines_path = Path(polylines_path)
    polylines_path.parent.mkdir(parents=True, exist_ok=True)

    options = dict(TORIC_SPINES_SKELETONIZE_DEFAULTS)
    options["profile"] = profile
    options["branching"] = branching
    options.update(skeletonize_kwargs)

    logger.info(
        "Skeletonizing mesh %s (profile=%s branching=%s extend_tips=%s)",
        mesh_path,
        options.get("profile"),
        options.get("branching"),
        options.get("extend_tips"),
    )
    mesh = load_and_repair(str(mesh_path))
    skeleton = skeletonize(mesh, **options)
    skeleton.write_polylines(str(polylines_path))
    logger.info("Wrote polylines to %s", polylines_path)
    return polylines_path.resolve()

compute_neck_candidates(ts_mesh, cell_mesh, params=None)

Detect neck interfaces between a spine mesh and the parent cell mesh.

Tries planar cap detection first, then attachment-face and boundary-loop fallbacks. params.max_necks truncates the ranked list.

Parameters:

Name Type Description Default
ts_mesh

Trimesh objects (spine and local cell).

required
cell_mesh

Trimesh objects (spine and local cell).

required
params NeckpointParams

Thresholds, crop margin, and max_necks.

None

Returns:

Type Description
list of NeckCandidate

Largest-first after truncation.

Examples:

>>> cands = compute_neck_candidates(ts, cell, NeckpointParams(max_necks=2))
>>> cands[0].centroid.shape
(3,)
Source code in toric_spines_sim/geometry/neckpoint.py
def compute_neck_candidates(
    ts_mesh,
    cell_mesh,
    params: Optional[NeckpointParams] = None,
) -> list[NeckCandidate]:
    """Detect neck interfaces between a spine mesh and the parent cell mesh.

    Tries planar cap detection first, then attachment-face and boundary-loop
    fallbacks. ``params.max_necks`` truncates the ranked list.

    Parameters
    ----------
    ts_mesh, cell_mesh
        Trimesh objects (spine and local cell).
    params : NeckpointParams, optional
        Thresholds, crop margin, and ``max_necks``.

    Returns
    -------
    list of NeckCandidate
        Largest-first after truncation.

    Examples
    --------
    >>> cands = compute_neck_candidates(ts, cell, NeckpointParams(max_necks=2))  # doctest: +SKIP
    >>> cands[0].centroid.shape
    (3,)
    """
    params = params or NeckpointParams()
    local = _crop_local_cell(cell_mesh, ts_mesh, params.crop_margin)
    candidates = _cap_candidates(ts_mesh, cell_mesh, local, params)
    if not candidates:
        logger.info("No neck caps passed filters; trying attachment fallback")
        candidates = _fallback_attachment_candidates(
            ts_mesh, cell_mesh, local, params
        )
    if not candidates:
        logger.info("Attachment fallback empty; trying boundary-loop fallback")
        candidates = _fallback_boundary_loop_candidates(ts_mesh, local, params)
    return _apply_max_necks(candidates, params.max_necks)

compute_neck_points(ts_mesh, cell_mesh, params=None)

Return neckpoint XYZ arrays (pixel space), largest-first after max_necks.

Source code in toric_spines_sim/geometry/neckpoint.py
def compute_neck_points(
    ts_mesh,
    cell_mesh,
    params: Optional[NeckpointParams] = None,
) -> list[np.ndarray]:
    """Return neckpoint XYZ arrays (pixel space), largest-first after ``max_necks``."""
    return [c.centroid.copy() for c in compute_neck_candidates(ts_mesh, cell_mesh, params)]

compute_neck_points_all(*, meshes=None, all_meshes=False, cell_mesh_path=None, params=None, write=True, overwrite=True)

Compute neckpoints for one or more TS meshes.

Source code in toric_spines_sim/geometry/neckpoint.py
def compute_neck_points_all(
    *,
    meshes: Optional[Sequence[str]] = None,
    all_meshes: bool = False,
    cell_mesh_path: Optional[PathLike] = None,
    params: Optional[NeckpointParams] = None,
    write: bool = True,
    overwrite: bool = True,
) -> dict[str, list[tuple[float, float, float]]]:
    """Compute neckpoints for one or more TS meshes."""
    targets = resolve_mesh_targets(meshes, all_meshes=all_meshes)
    results: dict[str, list[tuple[float, float, float]]] = {}
    for mesh_path in targets:
        results[mesh_path.stem] = compute_neck_points_for_spine(
            mesh_path,
            cell_mesh_path=cell_mesh_path,
            params=params,
            write=write,
            overwrite=overwrite,
        )
    return results

compute_neck_points_for_spine(spine_id, *, cell_mesh_path=None, params=None, write=True, overwrite=True, output_path=None)

Compute neckpoints for one TS mesh and optionally write the pixel file.

Parameters:

Name Type Description Default
spine_id PathLike

Mesh spine id, filename, or path (e.g. TS1, TS1.obj).

required
cell_mesh_path Optional[PathLike]

Full cell mesh; defaults to cell_wrapped_simplified.obj.

None
params Optional[NeckpointParams]

Detection parameters including max_necks.

None
write bool

If True, write data/pointsets/pixels/<spine_id>_neckpoint.txt.

True
overwrite bool

If False and the output exists, skip writing and return existing points only when write would be skipped after a successful compute — still recomputes unless you check existence first in the CLI.

True
output_path Optional[PathLike]

Override output path.

None
Source code in toric_spines_sim/geometry/neckpoint.py
def compute_neck_points_for_spine(
    spine_id: PathLike,
    *,
    cell_mesh_path: Optional[PathLike] = None,
    params: Optional[NeckpointParams] = None,
    write: bool = True,
    overwrite: bool = True,
    output_path: Optional[PathLike] = None,
) -> list[tuple[float, float, float]]:
    """Compute neckpoints for one TS mesh and optionally write the pixel file.

    Parameters
    ----------
    spine_id
        Mesh spine id, filename, or path (e.g. ``TS1``, ``TS1.obj``).
    cell_mesh_path
        Full cell mesh; defaults to ``cell_wrapped_simplified.obj``.
    params
        Detection parameters including ``max_necks``.
    write
        If True, write ``data/pointsets/pixels/<spine_id>_neckpoint.txt``.
    overwrite
        If False and the output exists, skip writing and return existing points
        only when write would be skipped after a successful compute — still
        recomputes unless you check existence first in the CLI.
    output_path
        Override output path.
    """
    params = params or NeckpointParams()
    mesh_path = resolve_mesh_path(spine_id)
    resolved_spine_id = mesh_path.stem
    cell_path = (
        Path(cell_mesh_path)
        if cell_mesh_path is not None
        else get_cell_mesh_path()
    )
    if not cell_path.is_file():
        # Allow bare name under data/mesh/
        resolved_cell_path = get_mesh_path(Path(cell_path).name)
        if resolved_cell_path.is_file():
            cell_path = resolved_cell_path
        else:
            raise FileNotFoundError(f"Cell mesh not found: {cell_path}")

    neckpoint_output_path = (
        Path(output_path)
        if output_path is not None
        else default_neckpoint_path(resolved_spine_id)
    )
    if write and neckpoint_output_path.exists() and not overwrite:
        logger.info("Skipping existing neckpoint file: %s", neckpoint_output_path)
        from toric_spines_sim.utils import load_xyz_points

        return load_xyz_points(neckpoint_output_path)

    ts = _load_trimesh(mesh_path)
    cell = _load_trimesh(cell_path)
    candidates = compute_neck_candidates(ts, cell, params)
    points = [
        (float(c.centroid[0]), float(c.centroid[1]), float(c.centroid[2]))
        for c in candidates
    ]
    for i, cand in enumerate(candidates):
        logger.info(
            "%s neck[%d] source=%s area=%.1f n=%d plan=%.3f signed=%.1f "
            "ray=%.2f dend=%.2f xyz=(%.3f, %.3f, %.3f)",
            resolved_spine_id,
            i,
            cand.source,
            cand.area,
            cand.n_faces,
            cand.planarity,
            cand.mean_signed,
            cand.ray_frac_inside,
            cand.dendrite_frac,
            cand.centroid[0],
            cand.centroid[1],
            cand.centroid[2],
        )

    if not points:
        logger.warning("No neckpoints found for %s", resolved_spine_id)
        return []

    if write:
        write_xyz_points(neckpoint_output_path, points)
        logger.info("Wrote %d neckpoint(s) to %s", len(points), neckpoint_output_path)
    return points

default_neckpoint_path(spine_id)

Pixel-space neckpoint path for a mesh spine id (e.g. TS1).

Source code in toric_spines_sim/geometry/neckpoint.py
def default_neckpoint_path(spine_id: str) -> Path:
    """Pixel-space neckpoint path for a mesh spine id (e.g. ``TS1``)."""
    return get_pointset_path(f"{spine_id}_neckpoint.txt", units="pixels")

append_sink_write(swc_path_pixels, *, neckpoint_path_pixels=None, radius_um=DEFAULT_SINK_RADIUS_UM, connector_length_um=DEFAULT_SINK_CONNECTOR_LENGTH_UM, n_cylinders=DEFAULT_SINK_N_CYLINDERS, um_per_px, output_swc_path_pixels=None, output_swc_path_microns=None, output_neckpoint_path_microns=None, tag=DEFAULT_SINK_TAG, last_segment_tag=DEFAULT_SINK_TIP_TAG)

Append a sink in pixel space, then write both pixel and micron SWCs.

Sink dimensions are specified in microns and converted to pixels for attachment. Outputs:

  • data/swc/pixels/<spine_id>_wsink_r<R>um.swc
  • data/swc/microns/<spine_id>_wsink_r<R>um.swc (scaled copy with SINK header fixed)
  • data/pointsets/microns/<spine_id>_neckpoint.txt

Returns (output_swc_path_pixels, output_swc_path_microns).

Source code in toric_spines_sim/geometry/prepare.py
def append_sink_write(
    swc_path_pixels: PathLike,
    *,
    neckpoint_path_pixels: Optional[PathLike] = None,
    radius_um: float = DEFAULT_SINK_RADIUS_UM,
    connector_length_um: float = DEFAULT_SINK_CONNECTOR_LENGTH_UM,
    n_cylinders: int = DEFAULT_SINK_N_CYLINDERS,
    um_per_px: float,
    output_swc_path_pixels: Optional[PathLike] = None,
    output_swc_path_microns: Optional[PathLike] = None,
    output_neckpoint_path_microns: Optional[PathLike] = None,
    tag: int = DEFAULT_SINK_TAG,
    last_segment_tag: int = DEFAULT_SINK_TIP_TAG,
) -> Tuple[Path, Path]:
    """Append a sink in pixel space, then write both pixel and micron SWCs.

    Sink dimensions are specified in microns and converted to pixels for
    attachment. Outputs:

    - ``data/swc/pixels/<spine_id>_wsink_r<R>um.swc``
    - ``data/swc/microns/<spine_id>_wsink_r<R>um.swc`` (scaled copy with SINK header fixed)
    - ``data/pointsets/microns/<spine_id>_neckpoint.txt``

    Returns ``(output_swc_path_pixels, output_swc_path_microns)``.
    """
    swc_path_pixels = Path(swc_path_pixels)
    spine_id = swc_path_pixels.stem
    if output_swc_path_pixels is None:
        output_swc_path_pixels = get_swc_path(
            f"{spine_id}_wsink_r{radius_um:g}um.swc", units="pixels"
        )
    else:
        output_swc_path_pixels = Path(output_swc_path_pixels)
    if output_swc_path_microns is None:
        output_swc_path_microns = get_swc_path(
            f"{spine_id}_wsink_r{radius_um:g}um.swc", units="microns"
        )
    else:
        output_swc_path_microns = Path(output_swc_path_microns)
    output_swc_path_pixels.parent.mkdir(parents=True, exist_ok=True)
    output_swc_path_microns.parent.mkdir(parents=True, exist_ok=True)

    neck_points_pixels, neck_source = resolve_neck_points_px(
        swc_path_pixels, neckpoint_path_pixels
    )
    neck_points_microns = [
        (x * um_per_px, y * um_per_px, z * um_per_px) for x, y, z in neck_points_pixels
    ]

    if output_neckpoint_path_microns is None:
        output_neckpoint_path_microns = get_pointset_path(
            f"{spine_id}_neckpoint.txt", units="microns"
        )
    write_xyz_points(output_neckpoint_path_microns, neck_points_microns)

    # Ensure a pixel neckpoint file exists for multi-neck append / direction.
    resolved_neckpoint_path_pixels = (
        Path(neckpoint_path_pixels)
        if neckpoint_path_pixels is not None
        else (POINTSETS_PIXELS_DIR / f"{spine_id}_neckpoint.txt")
    )
    if not resolved_neckpoint_path_pixels.exists():
        write_xyz_points(resolved_neckpoint_path_pixels, neck_points_pixels)

    radius_px = float(radius_um) / float(um_per_px)
    connector_length_px = float(connector_length_um) / float(um_per_px)
    direction = optimal_sink_direction(resolved_neckpoint_path_pixels, swc_path_pixels)

    geom = SinkGeometry(
        radius=radius_px,
        length=2.0 * radius_px,
        n_cylinders=int(n_cylinders),
        connector_length=connector_length_px,
        axis=direction,
    )
    logger.info(
        "%s sink axis=%s neck_source=%s radius=%.3g µm (%.3g px) necks=%d",
        spine_id,
        direction,
        neck_source,
        radius_um,
        radius_px,
        len(neck_points_pixels),
    )

    if len(neck_points_pixels) > 1:
        written_px = append_sink_to_swc_multi_neck_points(
            swc_in=swc_path_pixels,
            swc_out=output_swc_path_pixels,
            neck_points=resolved_neckpoint_path_pixels,
            geom=geom,
            tag=tag,
            last_segment_tag=last_segment_tag,
        )
    else:
        written_px = append_sink_to_swc(
            swc_in=swc_path_pixels,
            swc_out=output_swc_path_pixels,
            neck_coords=neck_points_pixels[0],
            geom=geom,
            tag=tag,
            last_segment_tag=last_segment_tag,
        )
    written_px = Path(written_px).resolve()
    written_um = scale_swc_file(written_px, output_swc_path_microns, um_per_px)
    logger.info("Wrote spine+sink (px) to %s", written_px)
    logger.info("Wrote spine+sink (µm) to %s", written_um)
    return written_px, written_um

append_sink_write_microns(swc_path_pixels, *, neckpoint_path_pixels=None, radius_um=DEFAULT_SINK_RADIUS_UM, connector_length_um=DEFAULT_SINK_CONNECTOR_LENGTH_UM, n_cylinders=DEFAULT_SINK_N_CYLINDERS, um_per_px, output_swc_path_microns=None, output_neckpoint_path_microns=None, tag=DEFAULT_SINK_TAG, last_segment_tag=DEFAULT_SINK_TIP_TAG, output_swc_path_pixels=None)

Append a sink and write micron (and pixel) SWCs; return the micron path.

Prefer :func:append_sink_write when both output paths are needed.

Source code in toric_spines_sim/geometry/prepare.py
def append_sink_write_microns(
    swc_path_pixels: PathLike,
    *,
    neckpoint_path_pixels: Optional[PathLike] = None,
    radius_um: float = DEFAULT_SINK_RADIUS_UM,
    connector_length_um: float = DEFAULT_SINK_CONNECTOR_LENGTH_UM,
    n_cylinders: int = DEFAULT_SINK_N_CYLINDERS,
    um_per_px: float,
    output_swc_path_microns: Optional[PathLike] = None,
    output_neckpoint_path_microns: Optional[PathLike] = None,
    tag: int = DEFAULT_SINK_TAG,
    last_segment_tag: int = DEFAULT_SINK_TIP_TAG,
    output_swc_path_pixels: Optional[PathLike] = None,
) -> Path:
    """Append a sink and write micron (and pixel) SWCs; return the micron path.

    Prefer :func:`append_sink_write` when both output paths are needed.
    """
    _px, um = append_sink_write(
        swc_path_pixels,
        neckpoint_path_pixels=neckpoint_path_pixels,
        radius_um=radius_um,
        connector_length_um=connector_length_um,
        n_cylinders=n_cylinders,
        um_per_px=um_per_px,
        output_swc_path_pixels=output_swc_path_pixels,
        output_swc_path_microns=output_swc_path_microns,
        output_neckpoint_path_microns=output_neckpoint_path_microns,
        tag=tag,
        last_segment_tag=last_segment_tag,
    )
    return um

convert_all_nff_active_zones()

Convert every data/nff/*.nff file that has s points to a pixel AZ file.

Source code in toric_spines_sim/geometry/prepare.py
def convert_all_nff_active_zones() -> list[Path]:
    """Convert every ``data/nff/*.nff`` file that has ``s`` points to a pixel AZ file."""
    written: list[Path] = []
    paths = sorted(NFF_DIR.glob("*.nff"))
    if not paths:
        raise FileNotFoundError(f"No .nff files found under {NFF_DIR}")
    for nff_path in paths:
        output_path = convert_nff_active_zone(nff_path)
        if output_path.exists() and output_path.stat().st_size > 0:
            written.append(output_path)
            logger.info("Wrote AZ pointset %s from %s", output_path, nff_path.name)
    return written

convert_nff_active_zone(nff_path, output_path=None)

Write NFF s points to data/pointsets/pixels/<spine_id>.txt.

If the NFF has no s points, writes an empty file and logs a warning.

Source code in toric_spines_sim/geometry/prepare.py
def convert_nff_active_zone(nff_path: PathLike, output_path: Optional[PathLike] = None) -> Path:
    """Write NFF ``s`` points to ``data/pointsets/pixels/<spine_id>.txt``.

    If the NFF has no ``s`` points, writes an empty file and logs a warning.
    """
    nff_path = Path(nff_path)
    if output_path is None:
        POINTSETS_PIXELS_DIR.mkdir(parents=True, exist_ok=True)
        output_path = POINTSETS_PIXELS_DIR / f"{nff_path.stem}.txt"
    else:
        output_path = Path(output_path)
        output_path.parent.mkdir(parents=True, exist_ok=True)
    points = read_nff_s_points(nff_path)
    if not points:
        logger.warning("No s-points in %s; skipping AZ file", nff_path.name)
        if Path(output_path).exists() and Path(output_path).stat().st_size == 0:
            Path(output_path).unlink()
        return Path(output_path)
    read_nff_points_and_write_txt_file(nff_path, output_path)
    return Path(output_path)

list_ts_spine_swcs()

Return sorted TS{n}.swc paths under data/swc/pixels (no _wsink_).

Source code in toric_spines_sim/geometry/prepare.py
def list_ts_spine_swcs() -> list[Path]:
    """Return sorted ``TS{n}.swc`` paths under ``data/swc/pixels`` (no ``_wsink_``)."""
    return sorted(
        p for p in SWC_PIXELS_DIR.glob("TS*.swc") if _SPINE_ID_PATTERN.match(p.stem)
    )

resolve_swc_targets(swcs=None, *, all_swcs=False)

Resolve CLI SWC targets under data/swc/pixels.

Source code in toric_spines_sim/geometry/prepare.py
def resolve_swc_targets(
    swcs: Optional[Sequence[str]] = None,
    *,
    all_swcs: bool = False,
) -> list[Path]:
    """Resolve CLI SWC targets under ``data/swc/pixels``."""
    if all_swcs and swcs:
        raise ValueError("Pass either --all or explicit SWC arguments, not both")
    if all_swcs:
        targets = list_ts_spine_swcs()
        if not targets:
            raise FileNotFoundError(f"No TS{{n}}.swc files found under {SWC_PIXELS_DIR}")
        return [p.resolve() for p in targets]
    if not swcs:
        raise ValueError("Provide one or more SWC arguments, or pass --all")
    return [resolve_swc_path(s) for s in swcs]

scale_swc_file(swc_in, swc_out, scale)

Scale SWC coordinates and radii, preserving CYCLE_BREAK / SINK headers.

SWCModel.scale keeps header text verbatim, so # SINK: length, radius, and neck_xyz are rewritten here to match the scaled geometry.

Source code in toric_spines_sim/geometry/prepare.py
def scale_swc_file(swc_in: PathLike, swc_out: PathLike, scale: float) -> Path:
    """Scale SWC coordinates and radii, preserving CYCLE_BREAK / SINK headers.

    ``SWCModel.scale`` keeps header text verbatim, so ``# SINK:`` length, radius,
    and ``neck_xyz`` are rewritten here to match the scaled geometry.
    """
    swc_in = Path(swc_in)
    swc_out = Path(swc_out)
    swc_out.parent.mkdir(parents=True, exist_ok=True)
    model = SWCModel.from_swc_file(str(swc_in), validate_reconnections=False)
    scaled = model.scale(scale)
    scaled.to_swc_file(str(swc_out))

    lines = swc_out.read_text(encoding="utf-8").splitlines(keepends=True)
    rewritten = [
        _scale_sink_header_line(line, scale) if line.startswith("# SINK:") else line
        for line in lines
    ]
    swc_out.write_text("".join(rewritten), encoding="utf-8")
    return swc_out.resolve()

write_synpts_microns(swc_path_pixels, az_path_pixels, synpts_path_microns=None, *, um_per_px)

Project pixel AZ points onto the SWC, scale to microns, and write synpts.

Parameters:

Name Type Description Default
swc_path_pixels path - like

Pixel-space SWC.

required
az_path_pixels path - like

Pixel-space active-zone points.

required
synpts_path_microns path - like

Output path. Default data/pointsets/microns/<spine_id>_synpts.txt.

None
um_per_px float

Scale factor (typically 0.005).

required

Returns:

Type Description
Path

Written micron synpts file.

Examples:

>>> write_synpts_microns("TS1.swc", "TS1_AZ.txt", um_per_px=0.005)
Source code in toric_spines_sim/geometry/prepare.py
def write_synpts_microns(
    swc_path_pixels: PathLike,
    az_path_pixels: PathLike,
    synpts_path_microns: Optional[PathLike] = None,
    *,
    um_per_px: float,
) -> Path:
    """Project pixel AZ points onto the SWC, scale to microns, and write synpts.

    Parameters
    ----------
    swc_path_pixels : path-like
        Pixel-space SWC.
    az_path_pixels : path-like
        Pixel-space active-zone points.
    synpts_path_microns : path-like, optional
        Output path. Default ``data/pointsets/microns/<spine_id>_synpts.txt``.
    um_per_px : float
        Scale factor (typically 0.005).

    Returns
    -------
    pathlib.Path
        Written micron synpts file.

    Examples
    --------
    >>> write_synpts_microns("TS1.swc", "TS1_AZ.txt", um_per_px=0.005)  # doctest: +SKIP
    """
    swc_path_pixels = Path(swc_path_pixels)
    az_path_pixels = Path(az_path_pixels)
    if synpts_path_microns is None:
        synpts_path_microns = get_pointset_path(
            f"{swc_path_pixels.stem}_synpts.txt", units="microns"
        )
    else:
        synpts_path_microns = Path(synpts_path_microns)
    synpts_path_microns.parent.mkdir(parents=True, exist_ok=True)

    swc_model = SWCModel.from_swc_file(str(swc_path_pixels), validate_reconnections=False)
    frusta = FrustaSet.from_swc_model(swc_model, sides=20, end_caps=False)
    az_pointset = PointSet.from_txt_file(str(az_path_pixels))
    if az_pointset is None or not az_pointset.points:
        raise ValueError(f"No AZ points in {az_path_pixels}")
    projected = az_pointset.project_onto_frusta(frusta)
    synpts = projected.scale(um_per_px)
    synpts.to_txt_file(str(synpts_path_microns))
    logger.info(
        "Wrote %d synpts (µm) to %s",
        len(synpts.points),
        synpts_path_microns,
    )
    return synpts_path_microns.resolve()