add geodesic path functionality

src/cedalion/dataclasses/geometry.py

@@ -38,6 +38,7 @@ def __lt__(self, other):
 @dataclass
 class Surface(ABC):
     """Abstract base class for surfaces."""
+
     mesh: Any
     crs: str
     units: pint.Unit
@@ -97,14 +98,15 @@ def snap(self, points: cdt.LabeledPointCloud):
 
 
 @dataclass
-class Voxels():
+class Voxels:
     """3D voxels represented by a np.array.
 
     Attributes:
         voxels (np.ndarray): The voxels.
         crs (str): The coordinate reference system of the voxels.
         units (pint.Unit): The units of the voxels.
     """
+
     voxels: np.ndarray
     crs: str
     units: pint.Unit
@@ -128,12 +130,12 @@ def nvertices(self) -> int:
     def apply_transform(self, transform: cdt.AffineTransform) -> "Voxels":
         # convert to homogeneous coordinates
         num, dim = self.voxels.shape
-        hom = np.ones((num,dim+1))
-        hom[:,:3] = self.voxels
+        hom = np.ones((num, dim + 1))
+        hom[:, :3] = self.voxels
         # apply transformation
         hom = (transform.pint.dequantify().values.dot(hom.T)).T
         # backtransformation
-        transformed = np.array([hom[i,:3] / hom[i,3] for i in range(hom.shape[0])])
+        transformed = np.array([hom[i, :3] / hom[i, 3] for i in range(hom.shape[0])])
 
         new_units = self.units * transform.pint.units
         new_crs = transform.dims[0]
@@ -162,6 +164,7 @@ class TrimeshSurface(Surface):
         crs (str): The coordinate reference system of the surface.
         units (pint.Unit): The units of the surface.
     """
+
     mesh: trimesh.Trimesh
 
     @property
@@ -571,7 +574,6 @@ def laplace_operator(self):
         B = (Be1 + Be1.T + Be2 + Be2.T + Be3 + Be3.T) / 12 + dBd
         return B, D, W, V
 
-
     @property
     def avg_edge_length(self):
         """Average length of all edges in the surface."""
@@ -582,7 +584,6 @@ def avg_edge_length(self):
         )
         return edgelens.mean()
 
-
     def surface_gradient(self, scalars, at_verts=True):
         """Gradient of a function with values `scalars` at each vertex on the surface.
 
@@ -608,12 +609,10 @@ def surface_gradient(self, scalars, at_verts=True):
 
         gradu = np.nan_to_num(((fe12 * pu3 + fe23 * pu1 + fe31 * pu2) / (2 * fa)).T)
 
-
         if at_verts:
             return (self.connected.dot(gradu).T / self.connected.sum(1).A.squeeze()).T
         return gradu
 
-
     @property
     def _facenorm_cross_edge(self):
         ppts = self.ppts
@@ -623,9 +622,8 @@ def _facenorm_cross_edge(self):
         fe31 = np.cross(fnorms, ppts[:, 0] - ppts[:, 2])
         return fe12, fe23, fe31
 
-
     def geodesic_distance(self, verts, m=1.0, fem=False):
-        """Calcualte the inimum mesh geodesic distance (in mm).
+        """Calcualte the minimum mesh geodesic distance (in mm).
 
         The geodesic distance is calculated from each vertex in surface to any vertex in
         the collection `verts`.
@@ -680,9 +678,7 @@ def geodesic_distance(self, verts, m=1.0, fem=False):
             self._rlfac_solvers[m] = sparse.linalg.factorized(
                 lfac[goodrows][:, goodrows]
             )
-            self._nLC_solvers[m] = sparse.linalg.factorized(
-                nLC[goodrows][:, goodrows]
-            )
+            self._nLC_solvers[m] = sparse.linalg.factorized(nLC[goodrows][:, goodrows])
 
         # I. "Integrate the heat flow ̇u = ∆u for some fixed time t"
         # ---------------------------------------------------------
@@ -701,7 +697,7 @@ def geodesic_distance(self, verts, m=1.0, fem=False):
         gradu = self.surface_gradient(u, at_verts=False)
 
         # Compute X (normalized grad u)
-        gusum = np.sum(gradu ** 2, axis=1)
+        gusum = np.sum(gradu**2, axis=1)
         X = np.nan_to_num((-gradu.T / np.sqrt(gusum)).T)
 
         # III. "Solve the Poisson equation ∆φ = ∇·X"
@@ -726,7 +722,6 @@ def geodesic_distance(self, verts, m=1.0, fem=False):
 
         return phi
 
-
     def geodesic_path(self, a, b, max_len=1000, d=None, **kwargs):
         """Finds the shortest path between two points `a` and `b`.
 
@@ -772,13 +767,29 @@ def geodesic_path(self, a, b, max_len=1000, d=None, **kwargs):
                 return path
         return path
 
+    @property
+    def graph(self):
+        """NetworkX undirected graph representing this Surface."""
+        import networkx as nx
+
+        graph = nx.Graph()
+        graph.add_edges_from(self.iter_surfedges)
+        return graph
+
+    @property
+    def iter_surfedges(self):
+        for a, b, c in self.mesh.polys:
+            yield a, b
+            yield b, c
+            yield a, c
+
     @property
     def _cot_edge(self):
         ppts = self.ppts
         cots1, cots2, cots3 = self.cotangent_weights
-        c3 = cots3[:,np.newaxis] * (ppts[:,1] - ppts[:,0])
-        c2 = cots2[:,np.newaxis] * (ppts[:,0] - ppts[:,2])
-        c1 = cots1[:,np.newaxis] * (ppts[:,2] - ppts[:,1])
+        c3 = cots3[:, np.newaxis] * (ppts[:, 1] - ppts[:, 0])
+        c2 = cots2[:, np.newaxis] * (ppts[:, 0] - ppts[:, 2])
+        c1 = cots1[:, np.newaxis] * (ppts[:, 2] - ppts[:, 1])
         c32 = c3 - c2
         c13 = c1 - c3
         c21 = c2 - c1
@@ -790,9 +801,15 @@ def _polyconn(self):
         npoly = len(self.mesh.polys)
         o = np.ones((npoly,))
 
-        c1 = sparse.coo_matrix((o, (self.mesh.polys[:,0], range(npoly))), (npt, npoly)).tocsr() # noqa: E501
-        c2 = sparse.coo_matrix((o, (self.mesh.polys[:,1], range(npoly))), (npt, npoly)).tocsr() # noqa: E501
-        c3 = sparse.coo_matrix((o, (self.mesh.polys[:,2], range(npoly))), (npt, npoly)).tocsr() # noqa: E501
+        c1 = sparse.coo_matrix(
+            (o, (self.mesh.polys[:, 0], range(npoly))), (npt, npoly)
+        ).tocsr()  # noqa: E501
+        c2 = sparse.coo_matrix(
+            (o, (self.mesh.polys[:, 1], range(npoly))), (npt, npoly)
+        ).tocsr()  # noqa: E501
+        c3 = sparse.coo_matrix(
+            (o, (self.mesh.polys[:, 2], range(npoly))), (npt, npoly)
+        ).tocsr()  # noqa: E501
 
         return c1, c2, c3