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]
 
@@ -48,19 +52,22 @@ class H5Geometry(ps.TriangularMesh):
     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
@@ -69,9 +76,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:
@@ -139,6 +146,9 @@ def triangles(self):
 
 
 class FreeSurferHemisphere(ps.TriangularMesh):
+    def __init__(self, meshes):
+        self._meshes = meshes
+
     @classmethod
     def from_filename(klass, pathlike):
         path = Path(pathlike)
@@ -176,8 +186,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())