Skip to content

Commit

Permalink
RF: Drop triangular meshes for now
Browse files Browse the repository at this point in the history
  • Loading branch information
@effigies
effigies committed Sep 19, 2023
1 parent c3ba28d commit 5ded851
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 210 deletions.
58 changes: 0 additions & 58 deletions nibabel/pointset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,64 +144,6 @@ def get_coords(self, *, as_homogeneous: bool = False):
return 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 self._triangles.shape[0]

def get_triangles(self):
"""Mx3 array of indices into coordinate table"""
return self._triangles

def get_mesh(self, name=None):
return self.get_coords(name=name), self.get_triangles()

def get_names(self):
"""List of surface names that can be passed to
``get_{coords,triangles,mesh}``
"""
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 Grid(Pointset):
r"""A regularly-spaced collection of coordinates
Expand Down
152 changes: 0 additions & 152 deletions nibabel/tests/test_pointset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,155 +182,3 @@ def test_to_mask(self):
],
)
assert np.array_equal(mask_img.affine, np.eye(4))


class H5ArrayProxy:
def __init__(self, file_like, dataset_name):
self.file_like = file_like
self.dataset_name = dataset_name
with h5.File(file_like, 'r') as h5f:
arr = h5f[dataset_name]
self._shape = arr.shape
self._dtype = arr.dtype

@property
def is_proxy(self):
return True

@property
def shape(self):
return self._shape

@property
def ndim(self):
return len(self.shape)

@property
def dtype(self):
return self._dtype

def __array__(self, dtype=None):
with h5.File(self.file_like, 'r') as h5f:
return np.asanyarray(h5f[self.dataset_name], dtype)

def __getitem__(self, slicer):
with h5.File(self.file_like, 'r') as h5f:
return h5f[self.dataset_name][slicer]


class H5Geometry(ps.TriMeshFamily):
"""Simple Geometry file structure that combines a single topology
with one or more coordinate sets
"""

@classmethod
def from_filename(klass, pathlike):
meshes = {}
with h5.File(pathlike, 'r') as h5f:
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):
with h5.File(pathlike, 'w') as h5f:
h5f.create_dataset('/topology', data=self.get_triangles())
for name, coord in self._coords.items():
h5f.create_dataset(f'/coordinates/{name}', data=coord)


class FSGeometryProxy:
def __init__(self, pathlike):
self._file_like = str(Path(pathlike))
self._offset = None
self._vnum = None
self._fnum = None

def _peek(self):
from nibabel.freesurfer.io import _fread3

with open(self._file_like, 'rb') as fobj:
magic = _fread3(fobj)
if magic != 16777214:
raise NotImplementedError('Triangle files only!')
fobj.readline()
fobj.readline()
self._vnum = np.fromfile(fobj, '>i4', 1)[0]
self._fnum = np.fromfile(fobj, '>i4', 1)[0]
self._offset = fobj.tell()

@property
def vnum(self):
if self._vnum is None:
self._peek()
return self._vnum

@property
def fnum(self):
if self._fnum is None:
self._peek()
return self._fnum

@property
def offset(self):
if self._offset is None:
self._peek()
return self._offset

@auto_attr
def coords(self):
ap = ArrayProxy(self._file_like, ((self.vnum, 3), '>f4', self.offset))
ap.order = 'C'
return ap

@auto_attr
def triangles(self):
offset = self.offset + 12 * self.vnum
ap = ArrayProxy(self._file_like, ((self.fnum, 3), '>i4', offset))
ap.order = 'C'
return ap


class FreeSurferHemisphere(ps.TriMeshFamily):
@classmethod
def from_filename(klass, pathlike):
path = Path(pathlike)
hemi, default = path.name.split('.')
mesh_names = (
'orig',
'white',
'smoothwm',
'pial',
'inflated',
'sphere',
'midthickness',
'graymid',
) # Often created
if default not in mesh_names:
mesh_names.append(default)
meshes = {}
for mesh in mesh_names:
fpath = path.parent / f'{hemi}.{mesh}'
if fpath.exists():
meshes[mesh] = FSGeometryProxy(fpath)
hemi = klass(meshes)
hemi._default = default
return hemi


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._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())

0 comments on commit 5ded851

Please sign in to comment.