From 916bff9a397cd441c80f1b4fcb912bef767a5d39 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 7 Sep 2023 12:15:41 -0400 Subject: [PATCH 01/10] ENH: Add pointset data structures [BIAP9] --- nibabel/pointset.py | 128 +++++++++++++++++++++++++ nibabel/tests/test_pointset.py | 166 +++++++++++++++++++++++++++++++++ 2 files changed, 294 insertions(+) create mode 100644 nibabel/pointset.py create mode 100644 nibabel/tests/test_pointset.py diff --git a/nibabel/pointset.py b/nibabel/pointset.py new file mode 100644 index 000000000..6c2523751 --- /dev/null +++ b/nibabel/pointset.py @@ -0,0 +1,128 @@ +import operator as op +from functools import reduce + +import numpy as np + +from nibabel.affines import apply_affine + + +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 self.get_coords().shape[0] + + def get_coords(self, name=None): + """Nx3 array of coordinates in RAS+ space""" + 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 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 + + ## This method is called for by the BIAP, but it now seems simpler to wait to + ## provide it until there are any proposed implementations + # def decimate(self, *, n_coords=None, ratio=None): + # """ Return a TriangularMesh with a smaller number of vertices that + # preserves the geometry of the original """ + # # To be overridden when a format provides optimization opportunities + # 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 + ---------- + shape : 3-tuple + number of coordinates in each dimension of grid + """ + + def __init__(self, shape, affines): + self.shape = tuple(shape) + try: + self._affines = dict(affines) + except (TypeError, ValueError): + self._affines = {'world': np.array(affines)} + if 'voxels' not in self._affines: + self._affines['voxels'] = np.eye(4, dtype=np.uint8) + + def get_affine(self, name=None): + """4x4 array""" + if name is None: + name = next(iter(self._affines)) + return self._affines[name] + + def get_coords(self, name=None): + if name is None: + name = next(iter(self._affines)) + aff = self.get_affine(name) + dt = np.result_type(*(np.min_scalar_type(dim) for dim in self.shape)) + # This is pretty wasteful; we almost certainly want instead an + # object that will retrieve a coordinate when indexed, but where + # np.array(obj) returns this + ijk_coords = np.array(list(np.ndindex(self.shape)), dtype=dt) + return apply_affine(aff, ijk_coords) + + @property + def n_coords(self): + return reduce(op.mul, self.shape) diff --git a/nibabel/tests/test_pointset.py b/nibabel/tests/test_pointset.py new file mode 100644 index 000000000..efea8bbd7 --- /dev/null +++ b/nibabel/tests/test_pointset.py @@ -0,0 +1,166 @@ +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.onetime import auto_attr +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): + 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()) From 5dceb6490bea99b5d0847459a8840bf58f8c6df9 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 7 Sep 2023 20:56:35 -0400 Subject: [PATCH 02/10] Update nibabel/pointset.py Co-authored-by: Oscar Esteban --- nibabel/pointset.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/nibabel/pointset.py b/nibabel/pointset.py index 6c2523751..91b753140 100644 --- a/nibabel/pointset.py +++ b/nibabel/pointset.py @@ -19,7 +19,15 @@ def n_coords(self): return self.get_coords().shape[0] def get_coords(self, name=None): - """Nx3 array of coordinates in RAS+ space""" + """Nx3 array of coordinates. + + Parameters + ---------- + name : :obj:`str` + Select a particular coordinate system if more than one may exist. + By default, `None` is equivalent to `"world"` and corresponds to + an RAS+ coordinate system. + """ return self._coords From 2b765df29e6338968b49807c99093d8b7280ded4 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 18 Sep 2023 15:08:17 -0400 Subject: [PATCH 03/10] MNT: Update pre-commit hooks --- .pre-commit-config.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1fc7efd0b..137aa4946 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,7 @@ exclude: ".*/data/.*" repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.1.0 + rev: v4.4.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer @@ -21,12 +21,12 @@ repos: hooks: - id: isort - repo: https://github.com/pycqa/flake8 - rev: 6.0.0 + rev: 6.1.0 hooks: - id: flake8 exclude: "^(doc|nisext|tools)/" - repo: https://github.com/pre-commit/mirrors-mypy - rev: v0.991 + rev: v1.5.1 hooks: - id: mypy # Sync with project.optional-dependencies.typing From 3f9b623f448305ce6c79c521201978dfbd315f19 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 18 Sep 2023 15:09:04 -0400 Subject: [PATCH 04/10] RF: Recast Pointset as a dataclass with associated affines --- nibabel/pointset.py | 226 +++++++++++++++++++++++++++++++++----------- 1 file changed, 173 insertions(+), 53 deletions(-) diff --git a/nibabel/pointset.py b/nibabel/pointset.py index 91b753140..c131b8131 100644 --- a/nibabel/pointset.py +++ b/nibabel/pointset.py @@ -1,34 +1,151 @@ -import operator as op -from functools import reduce +"""Point-set structures + +Imaging data are sampled at points in space, and these points +can be described by coordinates. +These structures are designed to enable operations on sets of +points, as opposed to the data sampled at those points. + +Abstractly, a point set is any collection of points, but there are +two types that warrant special consideration in the neuroimaging +context: grids and meshes. + +A *grid* is a collection of regularly-spaced points. The canonical +examples of grids are the indices of voxels and their affine +projection into a reference space. + +A *mesh* is a collection of points and some structure that enables +adjacent points to be identified. A *triangular mesh* in particular +uses triplets of adjacent vertices to describe faces. +""" +from __future__ import annotations + +import math +import typing as ty +from dataclasses import dataclass, replace import numpy as np -from nibabel.affines import apply_affine +from nibabel.casting import able_int_type +from nibabel.fileslice import strided_scalar +from nibabel.spatialimages import SpatialImage + +if ty.TYPE_CHECKING: # pragma: no cover + from typing_extensions import Self + _DType = ty.TypeVar('_DType', bound=np.dtype[ty.Any]) + +class CoordinateArray(ty.Protocol): + ndim: int + shape: tuple[int, int] + + @ty.overload + def __array__(self, dtype: None = ..., /) -> np.ndarray[ty.Any, np.dtype[ty.Any]]: + ... # pragma: no cover + + @ty.overload + def __array__(self, dtype: _DType, /) -> np.ndarray[ty.Any, _DType]: + ... # pragma: no cover + + +@dataclass class Pointset: - def __init__(self, coords): - self._coords = coords + """A collection of points described by coordinates. + + Parameters + ---------- + coords : array-like + 2-dimensional array with coordinates as rows + affine : :class:`numpy.ndarray` + Affine transform to be applied to coordinates array + homogeneous : :class:`bool` + Indicate whether the provided coordinates are homogeneous, + i.e., homogeneous 3D coordinates have the form ``(x, y, z, 1)`` + """ + + coordinates: CoordinateArray + affine: np.ndarray + homogeneous: bool = False + ndim = 2 + __array_priority__ = 99 + + def __init__( + self, + coordinates: CoordinateArray, + affine: np.ndarray | None = None, + homogeneous: bool = False, + ): + self.coordinates = coordinates + self.homogeneous = homogeneous + + if affine is None: + self.affine = np.eye(self.dim + 1) + else: + self.affine = np.asanyarray(affine) + + if self.affine.shape != (self.dim + 1,) * 2: + raise ValueError(f'Invalid affine for {self.dim}D coordinates:\n{self.affine}') + if np.any(self.affine[-1, :-1] != 0) or self.affine[-1, -1] != 1: + raise ValueError(f'Invalid affine matrix:\n{self.affine}') + + @property + def shape(self) -> tuple[int, int]: + """The shape of the coordinate array""" + return self.coordinates.shape @property - def n_coords(self): + def n_coords(self) -> int: """Number of coordinates Subclasses should override with more efficient implementations. """ - return self.get_coords().shape[0] + return self.coordinates.shape[0] + + @property + def dim(self) -> int: + """The dimensionality of the space the coordinates are in""" + return self.coordinates.shape[1] - self.homogeneous + + def __rmatmul__(self, affine: np.ndarray) -> Self: + """Apply an affine transformation to the pointset + + This will return a new pointset with an updated affine matrix only. + """ + return replace(self, affine=np.asanyarray(affine) @ self.affine) + + def _homogeneous_coords(self): + if self.homogeneous: + return np.asanyarray(self.coordinates) + + ones = strided_scalar( + shape=(self.coordinates.shape[0], 1), + scalar=np.array(1, dtype=self.coordinates.dtype), + ) + return np.hstack((self.coordinates, ones)) + + def get_coords(self, *, as_homogeneous: bool = False): + """Retrieve the coordinates - def get_coords(self, name=None): - """Nx3 array of coordinates. - Parameters ---------- - name : :obj:`str` + as_homogeneous : :class:`bool` + Return homogeneous coordinates if ``True``, or Cartesian + coordiantes if ``False``. + + name : :class:`str` Select a particular coordinate system if more than one may exist. By default, `None` is equivalent to `"world"` and corresponds to an RAS+ coordinate system. """ - return self._coords + ident = np.allclose(self.affine, np.eye(self.affine.shape[0])) + if self.homogeneous == as_homogeneous and ident: + return np.asanyarray(self.coordinates) + coords = self._homogeneous_coords() + if not ident: + coords = (self.affine @ coords.T).T + if not as_homogeneous: + coords = coords[:, :-1] + return coords class TriangularMesh(Pointset): @@ -65,14 +182,6 @@ def get_names(self): """ raise NotImplementedError - ## This method is called for by the BIAP, but it now seems simpler to wait to - ## provide it until there are any proposed implementations - # def decimate(self, *, n_coords=None, ratio=None): - # """ Return a TriangularMesh with a smaller number of vertices that - # preserves the geometry of the original """ - # # To be overridden when a format provides optimization opportunities - # raise NotImplementedError - class TriMeshFamily(TriangularMesh): def __init__(self, mapping, default=None): @@ -97,40 +206,51 @@ def get_coords(self, name=None): return self._coords[name] -class NdGrid(Pointset): - """ - Attributes - ---------- - shape : 3-tuple - number of coordinates in each dimension of grid +class Grid(Pointset): + r"""A regularly-spaced collection of coordinates + + This class provides factory methods for generating Pointsets from + :class:`~nibabel.spatialimages.SpatialImage`\s and generating masks + from coordinate sets. """ - def __init__(self, shape, affines): - self.shape = tuple(shape) - try: - self._affines = dict(affines) - except (TypeError, ValueError): - self._affines = {'world': np.array(affines)} - if 'voxels' not in self._affines: - self._affines['voxels'] = np.eye(4, dtype=np.uint8) - - def get_affine(self, name=None): - """4x4 array""" - if name is None: - name = next(iter(self._affines)) - return self._affines[name] + @classmethod + def from_image(cls, spatialimage: SpatialImage) -> Self: + return cls(coordinates=GridIndices(spatialimage.shape[:3]), affine=spatialimage.affine) - def get_coords(self, name=None): - if name is None: - name = next(iter(self._affines)) - aff = self.get_affine(name) - dt = np.result_type(*(np.min_scalar_type(dim) for dim in self.shape)) - # This is pretty wasteful; we almost certainly want instead an - # object that will retrieve a coordinate when indexed, but where - # np.array(obj) returns this - ijk_coords = np.array(list(np.ndindex(self.shape)), dtype=dt) - return apply_affine(aff, ijk_coords) + @classmethod + def from_mask(cls, mask: SpatialImage) -> Self: + mask_arr = np.bool_(mask.dataobj) + return cls( + coordinates=np.c_[np.nonzero(mask_arr)].astype(able_int_type(mask.shape)), + affine=mask.affine, + ) - @property - def n_coords(self): - return reduce(op.mul, self.shape) + def to_mask(self, shape=None) -> SpatialImage: + if shape is None: + shape = tuple(np.max(self.coordinates, axis=1)[: self.dim]) + mask_arr = np.zeros(shape, dtype='bool') + mask_arr[np.asanyarray(self.coordinates)[:, : self.dim]] = True + return SpatialImage(mask_arr, self.affine) + + +class GridIndices: + """Class for generating indices just-in-time""" + + __slots__ = ('gridshape', 'dtype', 'shape') + ndim = 2 + + def __init__(self, shape, dtype=None): + self.gridshape = shape + self.dtype = dtype or able_int_type(shape) + self.shape = (math.prod(self.gridshape), len(self.gridshape)) + + def __repr__(self): + return f'<{self.__class__.__name__}{self.gridshape}>' + + def __array__(self, dtype=None): + if dtype is None: + dtype = self.dtype + + axes = [np.arange(s, dtype=dtype) for s in self.gridshape] + return np.reshape(np.meshgrid(*axes, copy=False, indexing='ij'), (len(axes), -1)).T From f19ef3348b5a3a1dfe48eac59ec5f8e4d9ff0054 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 18 Sep 2023 15:37:20 -0400 Subject: [PATCH 05/10] TEST: Test Pointset and GridIndices classes --- nibabel/tests/test_pointset.py | 122 +++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) diff --git a/nibabel/tests/test_pointset.py b/nibabel/tests/test_pointset.py index efea8bbd7..88001b401 100644 --- a/nibabel/tests/test_pointset.py +++ b/nibabel/tests/test_pointset.py @@ -2,8 +2,10 @@ from unittest import skipUnless import numpy as np +import pytest from nibabel import pointset as ps +from nibabel.affines import apply_affine from nibabel.arrayproxy import ArrayProxy from nibabel.onetime import auto_attr from nibabel.optpkg import optional_package @@ -14,6 +16,126 @@ FS_DATA = Path(get_nibabel_data()) / 'nitest-freesurfer' +class TestPointsets: + rng = np.random.default_rng() + + @pytest.mark.parametrize('shape', [(5, 2), (5, 3), (5, 4)]) + @pytest.mark.parametrize('homogeneous', [True, False]) + def test_init(self, shape, homogeneous): + coords = self.rng.random(shape) + + if homogeneous: + coords = np.column_stack([coords, np.ones(shape[0])]) + + expected_shape = (shape[0], shape[1] + homogeneous) + + points = ps.Pointset(coords, homogeneous=homogeneous) + assert points.shape == expected_shape + assert np.allclose(points.affine, np.eye(shape[1] + 1)) + assert points.homogeneous is homogeneous + assert points.ndim == 2 + assert points.n_coords == shape[0] + assert points.dim == shape[1] + + points = ps.Pointset(coords, affine=np.diag([2] * shape[1] + [1]), homogeneous=homogeneous) + assert points.shape == expected_shape + assert np.allclose(points.affine, np.diag([2] * shape[1] + [1])) + assert points.homogeneous is homogeneous + assert points.ndim == 2 + assert points.n_coords == shape[0] + assert points.dim == shape[1] + + # Badly shaped affine + with pytest.raises(ValueError): + ps.Pointset(coords, affine=[0, 1]) + + # Badly valued affine + with pytest.raises(ValueError): + ps.Pointset(coords, affine=np.ones((shape[1] + 1, shape[1] + 1))) + + @pytest.mark.parametrize('shape', [(5, 2), (5, 3), (5, 4)]) + @pytest.mark.parametrize('homogeneous', [True, False]) + def test_affines(self, shape, homogeneous): + orig_coords = coords = self.rng.random(shape) + + if homogeneous: + coords = np.column_stack([coords, np.ones(shape[0])]) + + points = ps.Pointset(coords, homogeneous=homogeneous) + assert np.allclose(points.get_coords(), orig_coords) + + # Apply affines + scaler = np.diag([2] * shape[1] + [1]) + scaled = scaler @ points + assert np.array_equal(scaled.coordinates, points.coordinates) + assert np.array_equal(scaled.affine, scaler) + assert np.allclose(scaled.get_coords(), 2 * orig_coords) + + flipper = np.eye(shape[1] + 1) + # [[1, 0, 0], [0, 1, 0], [0, 0, 1]] becomes [[0, 1, 0], [1, 0, 0], [0, 0, 1]] + flipper[:-1] = flipper[-2::-1] + flipped = flipper @ points + assert np.array_equal(flipped.coordinates, points.coordinates) + assert np.array_equal(flipped.affine, flipper) + assert np.allclose(flipped.get_coords(), orig_coords[:, ::-1]) + + # Concatenate affines, with any associativity + for doubledup in [(scaler @ flipper) @ points, scaler @ (flipper @ points)]: + assert np.array_equal(doubledup.coordinates, points.coordinates) + assert np.allclose(doubledup.affine, scaler @ flipper) + assert np.allclose(doubledup.get_coords(), 2 * orig_coords[:, ::-1]) + + def test_homogeneous_coordinates(self): + ccoords = self.rng.random((5, 3)) + hcoords = np.column_stack([ccoords, np.ones(5)]) + + cartesian = ps.Pointset(ccoords) + homogeneous = ps.Pointset(hcoords, homogeneous=True) + + for points in (cartesian, homogeneous): + assert np.array_equal(points.get_coords(), ccoords) + assert np.array_equal(points.get_coords(as_homogeneous=True), hcoords) + + affine = np.diag([2, 3, 4, 1]) + cart2 = affine @ cartesian + homo2 = affine @ homogeneous + + exp_c = apply_affine(affine, ccoords) + exp_h = (affine @ hcoords.T).T + for points in (cart2, homo2): + assert np.array_equal(points.get_coords(), exp_c) + assert np.array_equal(points.get_coords(as_homogeneous=True), exp_h) + + +def test_GridIndices(): + # 2D case + shape = (2, 3) + gi = ps.GridIndices(shape) + + assert gi.dtype == np.dtype('u1') + assert gi.shape == (6, 2) + assert repr(gi) == '' + + gi_arr = np.asanyarray(gi) + assert gi_arr.dtype == np.dtype('u1') + assert gi_arr.shape == (6, 2) + # Tractable to write out + assert np.array_equal(gi_arr, [[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2]]) + + shape = (2, 3, 4) + gi = ps.GridIndices(shape) + + assert gi.dtype == np.dtype('u1') + assert gi.shape == (24, 3) + assert repr(gi) == '' + + gi_arr = np.asanyarray(gi) + assert gi_arr.dtype == np.dtype('u1') + assert gi_arr.shape == (24, 3) + # Separate implementation + assert np.array_equal(gi_arr, np.mgrid[:2, :3, :4].reshape(3, -1).T) + + class H5ArrayProxy: def __init__(self, file_like, dataset_name): self.file_like = file_like From 1e246154a4f72224a626b5dcce5f244dc4034c1f Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Mon, 18 Sep 2023 16:34:10 -0400 Subject: [PATCH 06/10] TEST: Test Grid methods --- nibabel/tests/test_pointset.py | 54 ++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/nibabel/tests/test_pointset.py b/nibabel/tests/test_pointset.py index 88001b401..35d47428e 100644 --- a/nibabel/tests/test_pointset.py +++ b/nibabel/tests/test_pointset.py @@ -1,3 +1,4 @@ +from math import prod from pathlib import Path from unittest import skipUnless @@ -7,8 +8,10 @@ from nibabel import pointset as ps from nibabel.affines import apply_affine from nibabel.arrayproxy import ArrayProxy +from nibabel.fileslice import strided_scalar from nibabel.onetime import auto_attr from nibabel.optpkg import optional_package +from nibabel.spatialimages import SpatialImage from nibabel.tests.nibabel_data import get_nibabel_data h5, has_h5py, _ = optional_package('h5py') @@ -136,6 +139,57 @@ def test_GridIndices(): assert np.array_equal(gi_arr, np.mgrid[:2, :3, :4].reshape(3, -1).T) +class TestGrids(TestPointsets): + @pytest.mark.parametrize('shape', [(5, 5, 5), (5, 5, 5, 5), (5, 5, 5, 5, 5)]) + def test_from_image(self, shape): + # Check image is generates voxel coordinates + affine = np.diag([2, 3, 4, 1]) + img = SpatialImage(strided_scalar(shape), affine) + grid = ps.Grid.from_image(img) + grid_coords = grid.get_coords() + + assert grid.shape == (prod(shape[:3]), 3) + assert np.allclose(grid.affine, affine) + + assert np.allclose(grid_coords[0], [0, 0, 0]) + # Final index is [4, 4, 4], scaled by affine + assert np.allclose(grid_coords[-1], [8, 12, 16]) + + def test_from_mask(self): + affine = np.diag([2, 3, 4, 1]) + mask = np.zeros((3, 3, 3)) + mask[1, 1, 1] = 1 + img = SpatialImage(mask, affine) + + grid = ps.Grid.from_mask(img) + grid_coords = grid.get_coords() + + assert grid.shape == (1, 3) + assert np.array_equal(grid_coords, [[2, 3, 4]]) + + def test_to_mask(self): + coords = np.array([[1, 1, 1]]) + + grid = ps.Grid(coords) + + mask_img = grid.to_mask() + assert mask_img.shape == (2, 2, 2) + assert np.array_equal(mask_img.get_fdata(), [[[0, 0], [0, 0]], [[0, 0], [0, 1]]]) + assert np.array_equal(mask_img.affine, np.eye(4)) + + mask_img = grid.to_mask(shape=(3, 3, 3)) + assert mask_img.shape == (3, 3, 3) + assert np.array_equal( + mask_img.get_fdata(), + [ + [[0, 0, 0], [0, 0, 0], [0, 0, 0]], + [[0, 0, 0], [0, 1, 0], [0, 0, 0]], + [[0, 0, 0], [0, 0, 0], [0, 0, 0]], + ], + ) + assert np.array_equal(mask_img.affine, np.eye(4)) + + class H5ArrayProxy: def __init__(self, file_like, dataset_name): self.file_like = file_like From 8ac45d061847a3305fdfba044c01f2d008de6cb2 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 19 Sep 2023 09:37:01 -0400 Subject: [PATCH 07/10] Apply suggestions from code review Co-authored-by: Oscar Esteban --- nibabel/pointset.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/pointset.py b/nibabel/pointset.py index c131b8131..cdb08c8cc 100644 --- a/nibabel/pointset.py +++ b/nibabel/pointset.py @@ -55,7 +55,7 @@ class Pointset: Parameters ---------- coords : array-like - 2-dimensional array with coordinates as rows + (*N*, *n*) array with *N* being points and columns their *n*-dimensional coordinates affine : :class:`numpy.ndarray` Affine transform to be applied to coordinates array homogeneous : :class:`bool` From f55c2868670676c35d7bc66a5fca6ea34c5057aa Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 19 Sep 2023 12:57:37 -0400 Subject: [PATCH 08/10] RF: Drop ndim/shape attributes, make explicit comment on __array_priority__ --- nibabel/pointset.py | 8 ++------ nibabel/tests/test_pointset.py | 18 ++++++------------ 2 files changed, 8 insertions(+), 18 deletions(-) diff --git a/nibabel/pointset.py b/nibabel/pointset.py index cdb08c8cc..162466d90 100644 --- a/nibabel/pointset.py +++ b/nibabel/pointset.py @@ -66,7 +66,8 @@ class Pointset: coordinates: CoordinateArray affine: np.ndarray homogeneous: bool = False - ndim = 2 + + # Force use of __rmatmul__ with numpy arrays __array_priority__ = 99 def __init__( @@ -88,11 +89,6 @@ def __init__( if np.any(self.affine[-1, :-1] != 0) or self.affine[-1, -1] != 1: raise ValueError(f'Invalid affine matrix:\n{self.affine}') - @property - def shape(self) -> tuple[int, int]: - """The shape of the coordinate array""" - return self.coordinates.shape - @property def n_coords(self) -> int: """Number of coordinates diff --git a/nibabel/tests/test_pointset.py b/nibabel/tests/test_pointset.py index 35d47428e..49c51251c 100644 --- a/nibabel/tests/test_pointset.py +++ b/nibabel/tests/test_pointset.py @@ -30,23 +30,15 @@ def test_init(self, shape, homogeneous): if homogeneous: coords = np.column_stack([coords, np.ones(shape[0])]) - expected_shape = (shape[0], shape[1] + homogeneous) - points = ps.Pointset(coords, homogeneous=homogeneous) - assert points.shape == expected_shape assert np.allclose(points.affine, np.eye(shape[1] + 1)) assert points.homogeneous is homogeneous - assert points.ndim == 2 - assert points.n_coords == shape[0] - assert points.dim == shape[1] + assert (points.n_coords, points.dim) == shape points = ps.Pointset(coords, affine=np.diag([2] * shape[1] + [1]), homogeneous=homogeneous) - assert points.shape == expected_shape assert np.allclose(points.affine, np.diag([2] * shape[1] + [1])) assert points.homogeneous is homogeneous - assert points.ndim == 2 - assert points.n_coords == shape[0] - assert points.dim == shape[1] + assert (points.n_coords, points.dim) == shape # Badly shaped affine with pytest.raises(ValueError): @@ -148,7 +140,8 @@ def test_from_image(self, shape): grid = ps.Grid.from_image(img) grid_coords = grid.get_coords() - assert grid.shape == (prod(shape[:3]), 3) + assert grid.n_coords == prod(shape[:3]) + assert grid.dim == 3 assert np.allclose(grid.affine, affine) assert np.allclose(grid_coords[0], [0, 0, 0]) @@ -164,7 +157,8 @@ def test_from_mask(self): grid = ps.Grid.from_mask(img) grid_coords = grid.get_coords() - assert grid.shape == (1, 3) + assert grid.n_coords == 1 + assert grid.dim == 3 assert np.array_equal(grid_coords, [[2, 3, 4]]) def test_to_mask(self): From c3ba28d558086b35baf2a9009d8ae099f397dcb5 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 19 Sep 2023 12:59:47 -0400 Subject: [PATCH 09/10] FIX: to_mask() implementation --- nibabel/pointset.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nibabel/pointset.py b/nibabel/pointset.py index 162466d90..324b76d36 100644 --- a/nibabel/pointset.py +++ b/nibabel/pointset.py @@ -224,9 +224,9 @@ def from_mask(cls, mask: SpatialImage) -> Self: def to_mask(self, shape=None) -> SpatialImage: if shape is None: - shape = tuple(np.max(self.coordinates, axis=1)[: self.dim]) + shape = tuple(np.max(self.coordinates, axis=0)[: self.dim] + 1) mask_arr = np.zeros(shape, dtype='bool') - mask_arr[np.asanyarray(self.coordinates)[:, : self.dim]] = True + mask_arr[tuple(np.asanyarray(self.coordinates)[:, : self.dim].T)] = True return SpatialImage(mask_arr, self.affine) From 5ded8517cb230712c8ecd70c9f0c510d2533874a Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 19 Sep 2023 14:33:27 -0400 Subject: [PATCH 10/10] RF: Drop triangular meshes for now --- nibabel/pointset.py | 58 ------------- nibabel/tests/test_pointset.py | 152 --------------------------------- 2 files changed, 210 deletions(-) diff --git a/nibabel/pointset.py b/nibabel/pointset.py index 324b76d36..b40449801 100644 --- a/nibabel/pointset.py +++ b/nibabel/pointset.py @@ -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 diff --git a/nibabel/tests/test_pointset.py b/nibabel/tests/test_pointset.py index 49c51251c..fb9a7c5c8 100644 --- a/nibabel/tests/test_pointset.py +++ b/nibabel/tests/test_pointset.py @@ -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())