ENH: Add pointset data structures [BIAP9]

[effigies]

Breaking this section out of #1090, as these will be useful for the transformation effort, as well as surfaces.

We will want factory functions or classmethods to generate these types from some common structures, such as SpatialImage, GiftiDataArray or FreeSurfer geometry arrays (which are just numpy arrays).

We may also consider pulling Parcel out of the CoordinateImage PR, for the sake of more efficiently handling subsets of coordinates.

nibabel/nibabel/coordimage.py

Lines 147 to 168 in 42e712e

	class Parcel:
	"""
	Attributes
	----------
	name : str
	structure : ``Pointset``
	indices : object that selects a subset of coordinates in structure
	"""
	def __init__(self, name, structure, indices):
	self.name = name
	self.structure = structure
	self.indices = indices
	def __repr__(self):
	return f'<Parcel {self.name}({len(self.indices)})>'
	def __len__(self):
	return len(self.indices)
	def __getitem__(self, slicer):
	return self.__class__(self.name, self.structure, self.indices[slicer])

This API is up for debate. I'm not sure I like the inclusion of named spaces and default spaces anymore. And we may want to expand these from such bare bones, though I think we should be careful of relying too much on an object.

Mapping/comparison to nitransforms concepts:

• Pointset = SampledSpatialData
  • SampledSpatialData also includes the dimensionality of the points (generally 3)
• NDGrid = ImageGrid
  • ImageGrid().ndindex corresponds to NDGrid().get_coords('voxels')
  • ImageGrid().ndcoords corresponds to NDGrid().get_coords()
  • ImageGrid preserves an original image header and pre-calculates an inverse affine
  • ImageGrid().ras() and ImageGrid().index() apply the affine (or inverse) to passed indices/coordinates

[codecov]

Codecov Report

Patch coverage: 98.78% and project coverage change: +0.04% 🎉

Comparison is base (5f37398) 92.16% compared to head (5ded851) 92.21%.
Report is 7 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1251      +/-   ##
==========================================
+ Coverage   92.16%   92.21%   +0.04%     
==========================================
  Files          98       99       +1     
  Lines       12371    12453      +82     
  Branches     2542     2559      +17     
==========================================
+ Hits        11402    11483      +81     
  Misses        646      646              
- Partials      323      324       +1     

Files Changed	Coverage Δ
nibabel/pointset.py	98.78% <98.78%> (ø)

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[oesteban]

A few questions, but it looks good to me overall

[effigies]

@oesteban What do you think about adding an affine directly to Pointset?

class Pointset:
    def __init__(self, coordinates, affine, is_homogeneous): ...

    def __rmatmul__(self, affine):
        """Transform coordinates without evaluating"""
        assert is_affine(affine)
        return replace(self, affine @ self.affine)

    def as_homogeneous(self):
        if self.is_homogeneous:
            return self.coordinates
        return np.hstack(self.coordinates, np.ones(...))

    def get_coords(self, homogeneous=False):
        """Evaluate coordinates"""
        coords = self.affine @ self.as_homogeneous()
        if not homogeneous:
            coords = coords[:, :-1]
        return coords

Then transforming can simply be:

transformed = xfm @ pointset

(It takes a little bit more to get this actually to work.)

from dataclasses import dataclass, replace
import numpy as np

@dataclass
class Pointset:
    coordinates: np.ndarray
    affine: np.ndarray
    is_homogeneous: bool

    def __rmatmul__(self, affine):
        """Transform coordinates without evaluating"""
        assert is_affine(affine)
        return replace(self, affine=affine @ self.affine)

    def as_homogeneous(self):
        if self.is_homogeneous:
            return self.coordinates
        return np.hstack(self.coordinates, np.ones((self.coordinates.shape[0], 1)))

    def get_coords(self, homogeneous=False):
        """Evaluate coordinates"""
        coords = self.affine @ self.as_homogeneous()
        if not homogeneous:
            coords = coords[:, :-1]
        return coords

    # Trick numpy into using our __rmatmul__
    ndim = 2
    __array_priority__ = 99


def is_affine(arr):
    return arr.shape == (4, 4) and np.array_equal(arr[3], [0, 0, 0, 1])

This would work for clouds, meshes or grids and would let us add format-defined affines without applying them prematurely. And for grids, it would allow us to wait to generate the voxel indices.

I would not attempt to apply warps with __matmul__, so a final call to get_coords() would be necessary.

[effigies]

I've gone ahead and pushed forward on bundling affines with Pointsets, defaulting to identity. Grid is now a simple wrapper that provides factory functions for creating a Pointset from a spatial image (all coordinates) or a mask (selected coordinates only).

Updated comparisons to ImageGrid:

• ImageGrid().ndindex corresponds to NDGrid().coordinates (untransformed)
• ImageGrid().ndcoords corresponds to NDGrid().get_coords() (transformed)

LMK what you think. It can be reverted, if this doesn't make sense to anybody else.

[oesteban]

Looks great. I left a few thoughts, but nothing stood out in my fast pass

[effigies]

I've removed the triangular meshes from this for now. The existing architecture depended too much on the previous Pointset implementation, and I don't want to hold this up.