RF: Add TriMeshFamily to formalize 1-mesh-to-many-coords structure

Author: effigies

nibabel/pointset.py

@@ -7,34 +7,49 @@
 
 
 class Pointset:
+    def __init__(self, coords):
+        self._coords = coords
+
     @property
     def n_coords(self):
         """ Number of coordinates
 
         Subclasses should override with more efficient implementations.
         """
-        return len(self.get_coords())
+        return self.get_coords().shape[0]
 
     def get_coords(self, name=None):
         """ Nx3 array of coordinates in RAS+ space """
-        raise NotImplementedError
+        return self._coords
 
 
 class TriangularMesh(Pointset):
+    def __init__(self, mesh):
+        if isinstance(mesh, tuple) and len(mesh) == 2:
+            coords, self._triangles = mesh
+        elif hasattr(mesh, 'coords') and hasattr(mesh, 'triangles'):
+            coords = mesh.coords
+            self._triangles = mesh.triangles
+        elif hasattr(mesh, 'get_mesh'):
+            coords, self._triangles = mesh.get_mesh()
+        else:
+            raise ValueError("Cannot interpret input as triangular mesh")
+        super().__init__(coords)
+
     @property
     def n_triangles(self):
         """ Number of faces
 
         Subclasses should override with more efficient implementations.
         """
-        return len(self.get_triangles())
+        return self._triangles.shape[0]
 
-    def get_triangles(self, name=None):
+    def get_triangles(self):
         """ Mx3 array of indices into coordinate table """
-        raise NotImplementedError
+        return self._triangles
 
     def get_mesh(self, name=None):
-        return self.get_coords(name=name), self.get_triangles(name=name)
+        return self.get_coords(name=name), self.get_triangles()
 
     def get_names(self):
         """ List of surface names that can be passed to
@@ -51,6 +66,29 @@ def get_names(self):
     #     raise NotImplementedError
 
 
+class TriMeshFamily(TriangularMesh):
+    def __init__(self, mapping, default=None):
+        self._triangles = None
+        self._coords = {}
+        for name, mesh in dict(mapping).items():
+            coords, triangles = TriangularMesh(mesh).get_mesh()
+            if self._triangles is None:
+                self._triangles = triangles
+            self._coords[name] = coords
+
+        if default is None:
+            default = next(iter(self._coords))
+        self._default = default
+
+    def get_names(self):
+        return list(self._coords)
+
+    def get_coords(self, name=None):
+        if name is None:
+            name = self._default
+        return self._coords[name]
+
+
 class NdGrid(Pointset):
     """
     Attributes

nibabel/tests/test_pointset.py

@@ -47,13 +47,10 @@ def __getitem__(self, slicer):
             return h5f[self.dataset_name][slicer]
 
 
-class H5Geometry(ps.TriangularMesh):
+class H5Geometry(ps.TriMeshFamily):
     """Simple Geometry file structure that combines a single topology
     with one or more coordinate sets
     """
-    def __init__(self, meshes):
-        self._meshes = meshes
-
     @classmethod
     def from_filename(klass, pathlike):
         meshes = {}
@@ -64,33 +61,11 @@ def from_filename(klass, pathlike):
         return klass(meshes)
 
     def to_filename(self, pathlike):
-        topology = None
-        coordinates = {}
-        for name, mesh in self._meshes.items():
-            coords, faces = mesh
-            if topology is None:
-                topology = faces
-            elif not np.array_equal(faces, topology):
-                raise ValueError("Inconsistent topology")
-            coordinates[name] = coords
-
         with h5.File(pathlike, "w") as h5f:
-            h5f.create_dataset("/topology", data=topology)
-            for name, coord in coordinates.items():
+            h5f.create_dataset("/topology", data=self.get_triangles())
+            for name, coord in self._coords.items():
                 h5f.create_dataset(f"/coordinates/{name}", data=coord)
 
-    def get_coords(self, name=None):
-        if name is None:
-            name = next(iter(self._meshes))
-        coords, _ = self._meshes[name]
-        return coords
-
-    def get_triangles(self, name=None):
-        if name is None:
-            name = next(iter(self._meshes))
-        _, triangles = self._meshes[name]
-        return triangles
-
 
 class FSGeometryProxy:
     def __init__(self, pathlike):
@@ -143,10 +118,7 @@ def triangles(self):
         return ap
 
 
-class FreeSurferHemisphere(ps.TriangularMesh):
-    def __init__(self, meshes):
-        self._meshes = meshes
-
+class FreeSurferHemisphere(ps.TriMeshFamily):
     @classmethod
     def from_filename(klass, pathlike):
         path = Path(pathlike)
@@ -165,24 +137,6 @@ def from_filename(klass, pathlike):
         hemi._default = default
         return hemi
 
-    def get_coords(self, name=None):
-        if name is None:
-            name = self._default
-        return self._meshes[name].coords
-
-    def get_triangles(self, name=None):
-        if name is None:
-            name = self._default
-        return self._meshes[name].triangles
-
-    @property
-    def n_coords(self):
-        return self._meshes[self._default].vnum
-
-    @property
-    def n_triangles(self):
-        return self._meshes[self._default].fnum
-
 
 def test_FreeSurferHemisphere():
     lh = FreeSurferHemisphere.from_filename(FS_DATA / 'fsaverage/surf/lh.white')
@@ -197,6 +151,6 @@ def test_make_H5Geometry(tmp_path):
     h5geo.to_filename(tmp_path / "geometry.h5")
 
     rt_h5geo = H5Geometry.from_filename(tmp_path / "geometry.h5")
-    assert set(h5geo._meshes) == set(rt_h5geo._meshes)
+    assert set(h5geo._coords) == set(rt_h5geo._coords)
     assert np.array_equal(lh.get_coords('white'), rt_h5geo.get_coords('white'))
     assert np.array_equal(lh.get_triangles(), rt_h5geo.get_triangles())