# BIAP9 - The Coordinate Image API¶

- Author:
Chris Markiewicz

- Status:
Draft

- Type:
Standards

- Created:
2021-09-16

## Background¶

### Surface data is generally kept separate from geometric metadata¶

In contrast to volumetric data, whose geometry can be fully encoded in the shape of a data array and a 4x4 affine matrix, data sampled to a surface require the location of each sample to be explicitly represented by a coordinate. In practice, the most common approach is to have a geometry file and a data file.

A geometry file consists of a vertex coordinate array and a triangle array describing the adjacency of vertices, while a data file is an n-dimensional array with one axis corresponding to vertex.

Keeping these files separate is a pragmatic optimization to avoid costly reproductions of geometric data, but presents an administrative burden to direct consumers of the data.

### Terminology¶

For the purposes of this BIAP, the following terms are used:

Coordinate - a triplet of floating point values in RAS+ space

Vertex - an index into a table of coordinates

Triangle (or face) - a triplet of adjacent vertices (A-B-C); the normal vector for the face is (\(\overline{AB}\times\overline{AC}\))

Topology - vertex adjacency data, independent of vertex coordinates, typically in the form of a list of triangles

Geometry - topology + a specific set of coordinates for a surface

Parcel - a subset of vertices; can be the full topology. Special cases include: * Patch - a connected parcel * Decimated mesh - a parcel that has a desired density of vertices

Parcel sequence - an ordered set of parcels

Data array - an n-dimensional array with one axis corresponding to the vertices (typical) OR faces (more rare) in a patch sequence

### Currently supported surface formats¶

- FreeSurfer
Geometry (e.g.

`lh.pial`

):`read_geometry()`

/`write_geometry()`

- Data
Morphometry:

`read_morph_data()`

/`write_morph_data()`

Labels:

`read_label()`

MGH:

`MGHImage`

- GIFTI:
`GiftiImage`

Every image contains a collection of data arrays, which may be coordinates, topology, or data (further subdivided by type and intent)

- GIFTI:
- CIFTI-2:
`Cifti2Image`

Pure data array, with image header containing flexible axes

The

`BrainModelAxis`

is a subspace sequence including patches for each hemisphere (cortex without the medial wall) and subcortical structures defined by indices into three-dimensional array and an affine matrixGeometry referred to by an associated

`wb.spec`

file (no current implementation in NiBabel)Possible to have one with no geometric information, e.g., parcels x time

- CIFTI-2:

### Other relevant formats¶

- MNE’s STC (source time course) format. Contains:
Subject name (resolvable with a FreeSurfer

`SUBJECTS_DIR`

)Index arrays into left and right hemisphere surfaces (subspace sequence)

Data, one of: * ndarray of shape

`(n_verts, n_times)`

* tuple of ndarrays of shapes`(n_verts, n_sensors)`

and`(n_sensors, n_times)`

Time start

Time step

## Desiderata for an API supporting surfaces¶

The following are provisional guiding principles:

A surface image (data array) should carry a reference to geometric metadata that is easily transferred to a new image.

Partial images (data only or geometry only) should be possible. Absence of components should have a well-defined signature, such as a property that is

`None`

or a specific`Exception`

is raised.All arrays (coordinates, triangles, data arrays) should be proxied to avoid excess memory consumption

Selecting among coordinates (e.g., gray/white boundary, inflated surface) for a single topology should be possible.

Combining multiple brain structures (canonically, left and right hemispheres) in memory should be easy; serializing to file may be format-specific.

Splitting a data array into independent patches that can be separately operated on and serialized should be possible.

### Prominent use cases¶

We consider the following use cases for working with surface data. A good API will make retrieving the components needed for each use case straightforward, as well as storing the results in new images.

Arithmetic/modeling - per-vertex mathematical operations

Smoothing - topology/geometry-respecting smoothing

Plotting - paint the data array as a texture on a surface

Decimation - subsampling a topology (possibly a subset, possibly with interpolated vertex locations)

Resampling to a geometrically-aligned surface * Downsampling by decimating, smoothing, resampling * Inter-subject resampling by using

`?h.sphere.reg`

Interpolation of per-vertex and per-face data arrays

When possible, we prefer to expose NumPy `ndarray`

s and
allow use of numpy, scipy, scikit-learn. In some cases, it may
make sense for NiBabel to provide methods.

## Proposal¶

A `CoordinateImage`

is an N-dimensional array, where one axis corresponds
to a sequence of points in one or more parcels.

```
class CoordinateImage:
"""
Attributes
----------
header : a file-specific header
coordaxis : ``CoordinateAxis``
dataobj : array-like
"""
class CoordinateAxis:
"""
Attributes
----------
parcels : list of ``Parcel`` objects
"""
def load_structures(self, mapping):
"""
Associate parcels to ``Pointset`` structures
"""
def __getitem__(self, slicer):
"""
Return a sub-sampled CoordinateAxis containing structures
matching the indices provided.
"""
def get_indices(self, parcel, indices=None):
"""
Return the indices in the full axis that correspond to the
requested parcel. If indices are provided, further subsample
the requested parcel.
"""
class Parcel:
"""
Attributes
----------
name : str
structure : ``Pointset``
indices : object that selects a subset of coordinates in structure
"""
```

To describe coordinate geometry, the following structures are proposed:

```
class Pointset:
@property
def n_coords(self):
""" Number of coordinates """
def get_coords(self, name=None):
""" Nx3 array of coordinates in RAS+ space """
class TriangularMesh(Pointset):
@property
def n_triangles(self):
""" Number of faces """
def get_triangles(self, name=None):
""" Mx3 array of indices into coordinate table """
def get_mesh(self, name=None):
return self.get_coords(name=name), self.get_triangles(name=name)
def get_names(self):
""" List of surface names that can be passed to
``get_{coords,triangles,mesh}``
"""
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
class NdGrid(Pointset):
"""
Attributes
----------
shape : 3-tuple
number of coordinates in each dimension of grid
"""
def get_affine(self, name=None):
""" 4x4 array """
```

The `NdGrid`

class allows raveled volumetric data to be treated the same as
triangular mesh or other coordinate data.

Finally, a structure for containing a collection of related geometric files is defined:

```
class GeometryCollection:
"""
Attributes
----------
structures : dict
Mapping from structure names to ``Pointset``
"""
@classmethod
def from_spec(klass, pathlike):
""" Load a collection of geometries from a specification. """
```

The canonical example of a geometry collection is a left hemisphere mesh, right hemisphere mesh.

Here we present common use cases:

### Modeling¶

```
from nilearn.glm.first_level import make_first_level_design_matrix, run_glm
bold = CoordinateImage.from_filename("/data/func/hemi-L_bold.func.gii")
dm = make_first_level_design_matrix(...)
labels, results = run_glm(bold.get_fdata(), dm)
betas = CoordinateImage(results["betas"], bold.coordaxis, bold.header)
betas.to_filename("/data/stats/hemi-L_betas.mgz")
```

In this case, no reference to the surface structure is needed, as the operations occur on a per-vertex basis. The coordinate axis and header are preserved to ensure that any metadata is not lost.

Here we assume that `CoordinateImage`

is able to make the appropriate
translations between formats (GIFTI, MGH). This is not guaranteed in the final
API.

### Smoothing¶

```
bold = CoordinateImage.from_filename("/data/func/hemi-L_bold.func.gii")
bold.coordaxis.load_structures({"lh": "/data/anat/hemi-L_midthickness.surf.gii"})
# Not implementing networkx weighted graph here, so assume we have a function
# that retrieves a graph for each structure
graphs = get_graphs(bold.coordaxis)
distances = distance_matrix(graphs['lh']) # n_coords x n_coords matrix
weights = normalize(gaussian(distances, sigma))
# Wildly inefficient smoothing algorithm
smoothed = CoordinateImage(weights @ bold.get_fdata(), bold.coordaxis, bold.header)
smoothed.to_filename(f"/data/func/hemi-L_smooth-{sigma}_bold.func.gii")
```

### Plotting¶

Nilearn currently provides a plot_surf function. With the proposed API, we could interface as follows:

```
def plot_surf_img(img, surface="inflated"):
from nilearn.plotting import plot_surf
coords, triangles = img.coordaxis.parcels[0].get_mesh(name=surface)
data = img.get_fdata()
return plot_surf((triangles, coords), data)
tstats = CoordinateImage.from_filename("/data/stats/hemi-L_contrast-taskVsBase_tstat.mgz")
# Assume a GeometryCollection that reads a FreeSurfer subject directory
fs_subject = FreeSurferSubject.from_spec("/data/subjects/fsaverage5")
tstats.coordaxis.load_structures(fs_subject.get_structure("lh"))
plot_surf_img(tstats)
```

### Subsampling CIFTI-2¶

```
img = nb.load("sub-01_task-rest_bold.dtseries.nii") # Assume CIFTI CoordinateImage
parcel = nb.load("sub-fsLR_hemi-L_label-DLPFC_mask.label.gii") # GiftiImage
structure = parcel.meta.metadata['AnatomicalStructurePrimary'] # "CortexLeft"
vtx_idcs = np.where(parcel.agg_data())[0]
dlpfc_idcs = img.coordaxis.get_indices(parcel=structure, indices=vtx_idcs)
# Subsampled coordinate axes will override any duplicate information from header
dlpfc_img = CoordinateImage(img.dataobj[dlpfc_idcs], img.coordaxis[dlpfc_idcs], img.header)
# Now load geometry so we can plot
wbspec = CaretSpec("fsLR.wb.spec")
dlpfc_img.coordaxis.load_structures(wbspec)
...
```