New API: PathGraph dataclass (#246)

The Skeleton class is a bit of a hodgepodge of data, generated data,
things that should be properties, manually cached properties, and
functions that should not be part of the class at all. Thanks to the
helpful prodding of @kevinyamauchi, this PR attempts to simplify the
concepts and data structures in skan into a simple class that can be
created without an image — since the path graph, not the image, is the
core of the computational abilities of skan.

This is still a work in progress but should already be useful: if you
have a graph as a `scipy.sparse.csr_array` (note: `csr_array`, not the
deprecated `csr_matrix` used by skan so far) generated through your own
means, you can make a "Skeleton" as follows:

```python
from skan.csr import PathGraph, Skeleton, summarize

g = PathGraph.from_graph(
        node_coordinates=coordinates,  # (n, ndim) NumPy array
        graph=graph,  # scipy.sparse.csr_array
        )

s = Skeleton.from_path_graph(g)
summary = summarize(s, separator='_')
```

Still to do in this PR:

- allow making PathGraphs from Skeletons
- allow `summarize` to take in PathGraphs directly

In the future, we might want to deprecate Skeleton altogether, but I'm
happy to do this over a long time.

@kevinyamauchi, I'm curious what you think about `paths` being a data
attribute, even though it is generated. I think it's expensive enough to
compute that we want it to be data and serializable. But it kinda breaks
the dataclass paradigm a little bit.

@DragaDoncila, you should be able to pull down this branch and use it on
your networkx tracking graphs after exporting them to csr arrays (I'm
pretty sure that's built-in to nx).

Author: jni

doc/misc/release-notes/release-v0.13.md

@@ -0,0 +1,29 @@
+# skan v0.13.0
+
+This is a minor step forward from v0.12.x which adds a new API:
+`skan.csr.PathGraph` is a more abstract version of a Skeleton, which only
+needs a pixel adjacency matrix to work, rather than a full skeleton image.
+The motivations for PathGraph are numerous:
+
+- a simpler dataclass object that doesn't require complex instantiation
+  logic. See e.g. Glyph's [Stop Writing `__init__`
+  Methods](https://blog.glyph.im/2025/04/stop-writing-init-methods.html)
+- making it easier to compute the pixel adjacency matrix separately, for
+  example [using dask when the images don't fit in
+  memory](https://blog.dask.org/2021/05/07/skeleton-analysis), and having
+  an API for which you can provide this matrix (rather than having to
+  modify a Skeleton instance in-place).
+- allowing more flexible use cases, for example to use skan to measure
+  tracks, as in
+  [live-image-tracking-tools/traccuracy#251](https://github.com/live-image-tracking-tools/traccuracy/pull/251).
+
+Due to some urgent need to use this code in the wild, this release doesn't
+provide any documentation examples. Indeed, the new class may see some changes
+in upcoming releases based on user feedback. See the discussion in
+[jni/skan#246](https://github.com/jni/skan/pull/246) for details. Look for
+further refinement of this idea in the 0.13.x releases!
+
+## New features
+
+- [#246](https://github.com/jni/skan/pull/246): New API: PathGraph dataclass
+

src/skan/csr.py

@@ -1,8 +1,12 @@
 from __future__ import annotations
 
+from dataclasses import dataclass
+from functools import cached_property
+
 import networkx as nx
 import numpy as np
 import pandas as pd
+import scipy
 from scipy import sparse, ndimage as ndi
 from scipy.sparse import csgraph
 from scipy.spatial import distance_matrix
@@ -203,8 +207,8 @@ def csr_to_nbgraph(csr, node_props=None):
         node_props = np.broadcast_to(1., csr.shape[0])
         node_props.flags.writeable = True
     return NBGraph(
-            csr.indptr,
-            csr.indices,
+            csr.indptr.astype(np.int32, copy=False),
+            csr.indices.astype(np.int32, copy=False),
             csr.data,
             np.array(csr.shape, dtype=np.int32),
             node_props.astype(np.float64),
@@ -410,8 +414,9 @@ def _build_skeleton_path_graph(graph):
     degrees = np.diff(graph.indptr)
     visited_data = np.zeros(graph.data.shape, dtype=bool)
     visited = NBGraphBool(
-            graph.indptr, graph.indices, visited_data, graph.shape,
-            np.broadcast_to(1.0, graph.shape[0])
+            graph.indptr.astype(np.int32, copy=False),
+            graph.indices.astype(np.int32, copy=False), visited_data,
+            graph.shape, np.broadcast_to(1.0, graph.shape[0])
             )
     endpoints = (degrees != 2)
     endpoint_degrees = degrees[endpoints]
@@ -440,6 +445,122 @@ def _build_skeleton_path_graph(graph):
     return paths
 
 
+@dataclass
+class PathGraph:
+    """Generalization of a skeleton.
+
+    A morphological skeleton is a collection of single pixel wide paths in
+    a binary (or, more generally, quantitative) image. Skan was created to
+    make measurements of those paths in image data. It turns out, though,
+    that the neighborhood topology of those paths (a graph containing long
+    strings of nodes of degree 2) is more generally useful, and indeed makes
+    sense without a source image or without representing pixels. (One example:
+    tracklets in cell tracking data.)
+
+    In the text below, we use the following notation:
+
+    - N: the number of points in the pixel skeleton,
+    - ndim: the dimensionality of the skeleton
+    - P: the number of paths in the skeleton (also the number of links in the
+      junction graph).
+    - J: the number of junction nodes
+    - Sd: the sum of the degrees of all the junction nodes
+    - [Nt], [Np], Nr, Nc: the dimensions of the source image
+    """
+    adj: scipy.sparse.csr_array  # pixel/node-id neighbor adjacency matrix
+    node_coordinates: np.ndarray | None
+    node_values: np.ndarray | None
+    paths: scipy.sparse.csr_array  # paths[i, j] = 1 iff coord j is in path i
+    spacing: float | tuple[float, ...] = 1  # spatial scale between coordinates
+
+    @classmethod
+    def from_graph(cls, *, adj, node_coordinates, node_values=None, spacing=1):
+        """Build a PathGraph from an adjacency matrix and node coordinates.
+
+        Parameters
+        ----------
+        adj : scipy.sparse.csr_array
+            An adjacency matrix where adj[i, j] is nonzero iff there is an
+            edge between node i and node j.
+        node_coordinates : np.ndarray, shape (N, ndim)
+            The coordinates of the nodes. node_coordinates[i] is the
+            coordinate of node i. The indices of these coordinates must match
+            the indices of adj.
+        node_values : np.ndarray, shape (N,)
+            Values of the nodes. Could be image intensity, height, or some
+            other quantity of which you want to compute statistics along the
+            path.
+        spacing : float or tuple of float, shape (ndim,)
+            The pixel/voxel spacing along each axis coordinate.
+        """
+        nbgraph = csr_to_nbgraph(adj, node_values)
+        paths = _build_skeleton_path_graph(nbgraph)
+        return cls(adj, node_coordinates, node_values, paths, spacing)
+
+    @classmethod
+    def from_image(cls, skeleton_image, *, spacing=1, value_is_height=False):
+        """Build a PathGraph from a skeleton image.
+
+        This is just a convenience meant to mirror Skeleton.__init__.
+        """
+        graph, coords = skeleton_to_csgraph(
+                skeleton_image,
+                spacing=spacing,
+                value_is_height=value_is_height,
+                )
+        values = _extract_values(skeleton_image, coords)
+        return cls.from_graph(
+                adj=graph,
+                node_coordinates=np.transpose(coords),
+                node_values=values,
+                spacing=spacing,
+                )
+
+    @cached_property
+    def nbgraph(self):
+        return csr_to_nbgraph(self.adj, self.node_values)
+
+    @cached_property
+    def distances(self):
+        """The path distances.
+
+        Returns
+        -------
+        distances : np.ndarray of float, shape (P,)
+            distances[i] contains the distance of path i.
+        """
+        distances = np.empty(self.n_paths, dtype=float)
+        _compute_distances(
+                self.nbgraph, self.paths.indptr, self.paths.indices, distances
+                )
+        return distances
+
+    @property
+    def n_paths(self):
+        return self.paths.shape[0]
+
+    @cached_property
+    def degrees(self):
+        """The degree (number of neighbors) of each node/pixel.
+
+        Returns
+        -------
+        degrees : np.ndarray of int, shape (N,)
+        """
+        return np.diff(self.adj.indptr)
+
+
+def _extract_values(image, coords):
+    if image.dtype == np.bool_:
+        return None
+    values = image[coords]
+    output_dtype = (
+            np.float64 if np.issubdtype(image.dtype, np.integer) else
+            image.dtype
+            )
+    return values.astype(output_dtype, copy=False)
+
+
 class Skeleton:
     """Object to group together all the properties of a skeleton.
 
@@ -522,12 +643,7 @@ def __init__(
                 spacing=spacing,
                 value_is_height=value_is_height,
                 )
-        if np.issubdtype(skeleton_image.dtype, np.floating):
-            self.pixel_values = skeleton_image[coords]
-        elif np.issubdtype(skeleton_image.dtype, np.integer):
-            self.pixel_values = skeleton_image.astype(np.float64)[coords]
-        else:
-            self.pixel_values = None
+        self.pixel_values = _extract_values(skeleton_image, coords)
         self.graph = graph
         self.nbgraph = csr_to_nbgraph(graph, self.pixel_values)
         self.coordinates = np.transpose(coords)
@@ -549,6 +665,26 @@ def __init__(
             self.skeleton_image = skeleton_image
             self.source_image = source_image
 
+    @classmethod
+    def from_path_graph(cls, pg: PathGraph):
+        dtype = bool if pg.node_values is None else pg.node_values.dtype
+        ndim = pg.node_coordinates.shape[1]
+        dummy_image = np.zeros((4,) * ndim, dtype=dtype)
+        dummy_image[([1, 2],) * ndim] = 1  # make a single diagonal branch
+        obj = cls(dummy_image, spacing=pg.spacing, keep_images=False)
+        obj.pixel_values = pg.node_values
+        obj.graph = pg.adj
+        obj.nbgraph = pg.nbgraph
+        obj.coordinates = pg.node_coordinates
+        obj.paths = pg.paths
+        obj.n_paths = pg.n_paths
+        obj.distances = pg.distances
+        obj.skeleton_shape = np.max(obj.coordinates, axis=0) + 1
+        obj.skeleton_dtype = dtype
+        obj.degrees = pg.degrees
+        obj.spacing = pg.spacing
+        return obj
+
     def path(self, index):
         """Return the pixel indices of path number `index`.
 
@@ -721,7 +857,7 @@ def __array__(self, dtype=None):
 
 
 def summarize(
-        skel: Skeleton,
+        skel: Skeleton | PathGraph,
         *,
         value_is_height: bool = False,
         find_main_branch: bool = False,
@@ -754,6 +890,8 @@ def summarize(
         A summary of the branches including branch length, mean branch value,
         branch euclidean distance, etc.
     """
+    if isinstance(skel, PathGraph):
+        skel = Skeleton.from_path_graph(skel)
     if separator is None:
         warnings.warn(
                 "separator in column name will change to _ in version 0.13; "

src/skan/test/test_csr.py

@@ -384,6 +384,39 @@ def test_skeletonlabel():
     assert stats['mean-pixel-value'].max() > 1
 
 
+@pytest.mark.parametrize(('np_skeleton', 'spacing', 'dtype'),
+                         product(
+                                 [
+                                         tinycycle, tinyline, skeleton0,
+                                         skeleton1, skeleton2, skeleton3d
+                                         ],
+                                 [1.0, 2.0, (5.0, 2.5, 2.5)],
+                                 [bool, np.float32, np.float64],
+                                 ))
+def test_pathgraph_skeleton_equiv(np_skeleton, spacing, dtype):
+    if dtype is not bool:
+        np_skeleton = (
+                np_skeleton
+                * np.random.random(np_skeleton.shape).astype(dtype)
+                )
+    if isinstance(spacing, tuple):
+        spacing = spacing[:np_skeleton.ndim]  # truncate spacing to image ndim
+    s = csr.Skeleton(np_skeleton, spacing=spacing)
+    p = csr.PathGraph.from_image(np_skeleton, spacing=spacing)
+    p2 = csr.PathGraph.from_graph(
+            adj=p.adj,
+            node_coordinates=p.node_coordinates,
+            node_values=p.node_values,
+            spacing=spacing,
+            )
+    ss = csr.summarize(s)
+    sp = csr.summarize(p)
+    sp2 = csr.summarize(p2)
+
+    np.testing.assert_allclose(ss.to_numpy(), sp.to_numpy())
+    np.testing.assert_allclose(sp.to_numpy(), sp2.to_numpy())
+
+
 @pytest.mark.parametrize(
         ('np_skeleton', 'summary', 'nodes', 'edges'),
         [