TEST: Read FreeSurfer subject, save/read H5 subject

Author: effigies

nibabel/tests/test_pointset.py

@@ -1,13 +1,17 @@
 from pathlib import Path
+from unittest import skipUnless
 
 import numpy as np
 
 from nibabel import pointset as ps
 from nibabel.arrayproxy import ArrayProxy
 from nibabel.optpkg import optional_package
+from nibabel.tests.nibabel_data import get_nibabel_data
 
 h5, has_h5py, _ = optional_package('h5py')
 
+FS_DATA = Path(get_nibabel_data()) / 'nitest-freesurfer'
+
 
 class H5ArrayProxy:
     def __init__(self, file_like, dataset_name):
@@ -38,7 +42,7 @@ def __array__(self, dtype=None):
         with h5.File(self.file_like, "r") as h5f:
             return np.asanyarray(h5f[self.dataset_name], dtype)
 
-    def __slicer__(self, slicer):
+    def __getitem__(self, slicer):
         with h5.File(self.file_like, "r") as h5f:
             return h5f[self.dataset_name][slicer]
 
@@ -47,19 +51,22 @@ class H5Geometry(ps.TriangularMesh):
     """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 = {}
         with h5.File(pathlike, "r") as h5f:
-            triangles = h5f['topology']
-            for name, coords in h5f['coordinates'].items():
-                meshes[name] = (coords, triangles)
+            triangles = H5ArrayProxy(pathlike, '/topology')
+            for name in h5f['coordinates']:
+                meshes[name] = (H5ArrayProxy(pathlike, f'/coordinates/{name}'), triangles)
         return klass(meshes)
 
     def to_filename(self, pathlike):
         topology = None
         coordinates = {}
-        for name, mesh in self.meshes.items():
+        for name, mesh in self._meshes.items():
             coords, faces = mesh
             if topology is None:
                 topology = faces
@@ -68,9 +75,9 @@ def to_filename(self, pathlike):
             coordinates[name] = coords
 
         with h5.File(pathlike, "w") as h5f:
-            h5f.create_dataset("/topology", topology)
+            h5f.create_dataset("/topology", data=topology)
             for name, coord in coordinates.items():
-                h5f.create_dataset(f"/coordinates/{name}", coord)
+                h5f.create_dataset(f"/coordinates/{name}", data=coord)
 
     def get_coords(self, name=None):
         if name is None:
@@ -137,6 +144,9 @@ def triangles(self):
 
 
 class FreeSurferHemisphere(ps.TriangularMesh):
+    def __init__(self, meshes):
+        self._meshes = meshes
+
     @classmethod
     def from_filename(klass, pathlike):
         path = Path(pathlike)
@@ -167,8 +177,26 @@ def get_triangles(self, name=None):
 
     @property
     def n_coords(self):
-        return self.meshes[self._default].vnum
+        return self._meshes[self._default].vnum
 
     @property
     def n_triangles(self):
-        return self.meshes[self._default].fnum
+        return self._meshes[self._default].fnum
+
+
+def test_FreeSurferHemisphere():
+    lh = FreeSurferHemisphere.from_filename(FS_DATA / 'fsaverage/surf/lh.white')
+    assert lh.n_coords == 163842
+    assert lh.n_triangles == 327680
+
+
+@skipUnless(has_h5py, reason="Test requires h5py")
+def test_make_H5Geometry(tmp_path):
+    lh = FreeSurferHemisphere.from_filename(FS_DATA / 'fsaverage/surf/lh.white')
+    h5geo = H5Geometry({name: lh.get_mesh(name) for name in ('white', 'pial')})
+    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 np.array_equal(lh.get_coords('white'), rt_h5geo.get_coords('white'))
+    assert np.array_equal(lh.get_triangles(), rt_h5geo.get_triangles())