From cb7d7d3e4033084962163ee0480f530508213f97 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Thu, 16 Sep 2021 22:21:22 -0400
Subject: [PATCH 01/23] DOC: First pass at SurfaceImage BIAP

---
 doc/source/devel/biaps/biap_0009.rst | 159 +++++++++++++++++++++++++++
 doc/source/devel/biaps/index.rst     |   1 +
 2 files changed, 160 insertions(+)
 create mode 100644 doc/source/devel/biaps/biap_0009.rst

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
new file mode 100644
index 0000000000..33412a07ee
--- /dev/null
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -0,0 +1,159 @@
+.. _biap9:
+
+#############################
+BIAP9 - The Surface 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
+requires 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
+* Patch - a connected subset of vertices (can be the full topology)
+* Data array - an n-dimensional array with one axis corresponding to the the
+  vertices (typical) OR faces (more rare) in a patch
+
+
+Currently supported surface formats
+===================================
+
+* FreeSurfer
+   * Geometry (e.g. ``lh.pial``):
+     :py:func:`~nibabel.freesurfer.io.read_geometry` / 
+     :py:func:`~nibabel.freesurfer.io.write_geometry`
+   * Data
+      * Morphometry:
+        :py:func:`~nibabel.freesurfer.io.read_morph_data` /
+        :py:func:`~nibabel.freesurfer.io.write_morph_data`
+      * Labels: :py:func:`~nibabel.freesurfer.io.read_label`
+      * MGH: :py:class:`~nibabel.freesurfer.mghformat.MGHImage`
+* GIFTI: :py:class:`~nibabel.gifti.gifti.GiftiImage`
+   * Every image contains a collection of data arrays, which may be
+     coordinates, topology, or data (further subdivided by type and intent)
+* CIFTI-2: :py:class:`~nibabel.cifti2.cifti2.Cifti2Image`
+   * Pure data array, with image header containing flexible axes
+   * Geometry referred to by an associated ``wb.spec`` file
+     (no current implementation in NiBabel)
+
+
+*********************************
+Desiderata for a SurfaceImage API
+*********************************
+
+The following are provisional guiding principles 
+
+1. A surface image (data array) should carry a reference to geometric metadata
+   that is easily transferred to a new image.
+2. 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.
+3. All arrays (coordinates, triangles, data arrays) should be proxied to
+   avoid excess memory consumption
+4. Selecting among coordinates (e.g., gray/white boundary, inflated surface)
+   for a single topology should be possible.
+5. Combining multiple brain structures (canonically, left and right hemispheres)
+   in memory should be easy; serializing to file may be format-specific.
+6. Splitting a data array into independent patches that can be separately
+   operated on and serialized should be possible.
+
+
+Prominent use cases
+===================
+
+* 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
+
+These are not necessarily operations that NiBabel must provide, but the
+components needed for each should be readily retrieved.
+
+********
+Proposal
+********
+
+.. code-block:: python
+
+    class SurfaceHeader:
+        @property
+        def ncoords(self):
+            """ Number of coordinates """
+
+        @property
+        def nfaces(self):
+            """ Number of faces """
+
+        def get_coords(self, name=None):
+            """ Nx3 array of coordinates in RAS+ space """
+
+        def get_faces(self, name=None):
+            """ Mx3 array of indices into coordinate table """
+
+        def get_mesh(self, name=None):
+            return self.get_coords(name=name), self.get_faces(name=name)
+
+        def decimate(self, *, ncoords=None, ratio=None):
+            """ Return a SurfaceHeader with a smaller number of vertices that
+            preserves the geometry of the original """
+            # To be overridden when a format provides optimization opportunities
+
+        def load_vertex_data(self):
+            """ Return a SurfaceImage with data corresponding to each vertex """
+
+        def load_face_data(self):
+            """ Return a SurfaceImage with data corresponding to each face """
+
+
+    class SurfaceImage:
+        @property
+        def header(self):
+            """ A SurfaceHeader or None """
+
+        @property
+        def dataobj(self):
+            """ An ndarray or ArrayProxy with one of the following properties:
+            
+            1) self.dataobj.shape[0] == self.header.ncoords
+            2) self.dataobj.shape[0] == self.header.nfaces
+            """
+
+        def load_header(self, pathlike):
+            """ Specify a header to a data-only image """
diff --git a/doc/source/devel/biaps/index.rst b/doc/source/devel/biaps/index.rst
index 1c30b3cd3c..42dcd276e3 100644
--- a/doc/source/devel/biaps/index.rst
+++ b/doc/source/devel/biaps/index.rst
@@ -19,6 +19,7 @@ proposals.
    biap_0006
    biap_0007
    biap_0008
+   biap_0009
 
 .. toctree::
    :hidden:

From ea2a466ae1f7d816eca07981ece77e6f5420fe4a Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Wed, 29 Sep 2021 13:12:24 -0400
Subject: [PATCH 02/23] DOC: Address suggestions

---
 doc/source/devel/biaps/biap_0009.rst | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 33412a07ee..4d4fcde05f 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -43,7 +43,7 @@ For the purposes of this BIAP, the following terms are used:
   typically in the form of a list of triangles
 * Geometry - topology + a specific set of coordinates for a surface
 * Patch - a connected subset of vertices (can be the full topology)
-* Data array - an n-dimensional array with one axis corresponding to the the
+* Data array - an n-dimensional array with one axis corresponding to the
   vertices (typical) OR faces (more rare) in a patch
 
 
@@ -67,6 +67,7 @@ Currently supported surface formats
    * Pure data array, with image header containing flexible axes
    * Geometry 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
 
 
 *********************************
@@ -114,21 +115,26 @@ Proposal
 
     class SurfaceHeader:
         @property
-        def ncoords(self):
+        def n_coords(self):
             """ Number of coordinates """
 
         @property
-        def nfaces(self):
+        def n_triangles(self):
             """ Number of faces """
 
         def get_coords(self, name=None):
             """ Nx3 array of coordinates in RAS+ space """
 
-        def get_faces(self, name=None):
+        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_faces(name=name)
+            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, *, ncoords=None, ratio=None):
             """ Return a SurfaceHeader with a smaller number of vertices that

From 0f4fddcd2bd2c3ef65f7ec0e8f6103570a282423 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Thu, 7 Oct 2021 21:42:27 -0400
Subject: [PATCH 03/23] DOC: Clarify use cases are motivating, not necessarily
 implementation target

---
 doc/source/devel/biaps/biap_0009.rst | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 4d4fcde05f..61680e8311 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -94,6 +94,10 @@ The following are provisional guiding principles
 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
@@ -104,8 +108,9 @@ Prominent use cases
   * Inter-subject resampling by using ``?h.sphere.reg``
 * Interpolation of per-vertex and per-face data arrays
 
-These are not necessarily operations that NiBabel must provide, but the
-components needed for each should be readily retrieved.
+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

From 3b2b575e31240371b1db9a0f3112a4a1abf30f66 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Thu, 7 Oct 2021 21:43:05 -0400
Subject: [PATCH 04/23] DOC: Small updates

---
 doc/source/devel/biaps/biap_0009.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 61680e8311..1d5c3fa086 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -141,15 +141,15 @@ Proposal
             ``get_{coords,triangles,mesh}``
             """
 
-        def decimate(self, *, ncoords=None, ratio=None):
+        def decimate(self, *, n_coords=None, ratio=None):
             """ Return a SurfaceHeader with a smaller number of vertices that
             preserves the geometry of the original """
             # To be overridden when a format provides optimization opportunities
 
-        def load_vertex_data(self):
+        def load_vertex_data(self, filename):
             """ Return a SurfaceImage with data corresponding to each vertex """
 
-        def load_face_data(self):
+        def load_face_data(self, filename):
             """ Return a SurfaceImage with data corresponding to each face """
 
 

From 9ddbcc79fea7c9d20f57fcffbcf003f75f0c6255 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 8 Oct 2021 08:58:56 -0400
Subject: [PATCH 05/23] DOC: Update BIAP with a couple examples leading to
 further questions

---
 doc/source/devel/biaps/biap_0009.rst | 56 ++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 1d5c3fa086..0530652489 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -116,6 +116,8 @@ make sense for NiBabel to provide methods.
 Proposal
 ********
 
+The basic API is as follows:
+
 .. code-block:: python
 
     class SurfaceHeader:
@@ -168,3 +170,57 @@ Proposal
 
         def load_header(self, pathlike):
             """ Specify a header to a data-only image """
+
+
+Here we present common use cases:
+
+
+Modeling
+--------
+
+.. code-block:: python
+
+    from nilearn.glm.first_level import make_first_level_design_matrix, run_glm
+
+    bold = SurfaceImage.from_filename("/data/func/hemi-L_bold.func.gii")
+    dm = make_first_level_design_matrix(...)
+    labels, results = run_glm(bold.get_data(), dm)
+    betas = bold.__class__(results["betas"], bold.header)
+    betas.to_filename("/data/stats/hemi-L_betas.mgz")
+
+Data images have their own metadata apart from geometry that needs handling in a header:
+
+* Axis labels (time series, labels, tensors)
+* Data dtype
+
+
+Plotting
+--------
+
+Nilearn currently provides a
+`plot_surf <https://nilearn.github.io/modules/generated/nilearn.plotting.plot_surf.html>`_ function.
+With the proposed API, we could interface as follows:
+
+.. code-block:: python
+
+    def plot_surf_img(img, surface="inflated"):
+        from nilearn.plotting import plot_surf
+        coords, triangles = img.header.get_mesh(name=surface)
+
+        data = tstats.get_data()
+
+        return plot_surf((triangles, coords), data)
+
+    tstats = SurfaceImage.from_filename("/data/stats/hemi-L_contrast-taskVsBase_tstat.mgz")
+    tstats.load_header("/data/subjects/fsaverage5")
+    plot_surf_img(tstats)
+
+This assumes that ``load_header()`` will be able to identify a FreeSurfer subject directory.
+
+******************
+* Open questions *
+******************
+
+1) Can/should image and header objects be promiscuous? If I load a GIFTI header,
+   should ``header.load_vertex_data()`` accept a FreeSurfer curvature file, or
+   should there be a bit more ecosystem-local logic?

From 27fd6f4e9216db1dc6f05712417ca5e422a21e5d Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 8 Oct 2021 10:10:16 -0400
Subject: [PATCH 06/23] DOC: Separate Geometry and Header objects

---
 doc/source/devel/biaps/biap_0009.rst | 28 ++++++++++++++++++----------
 1 file changed, 18 insertions(+), 10 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 0530652489..ed1982a063 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -120,7 +120,7 @@ The basic API is as follows:
 
 .. code-block:: python
 
-    class SurfaceHeader:
+    class SurfaceGeometry:
         @property
         def n_coords(self):
             """ Number of coordinates """
@@ -155,10 +155,18 @@ The basic API is as follows:
             """ Return a SurfaceImage with data corresponding to each face """
 
 
-    class SurfaceImage:
+    class SurfaceHeader(FileBaseeHeader):
+        ...
+
+
+    class SurfaceImage(FileBasedImage):
         @property
         def header(self):
-            """ A SurfaceHeader or None """
+            """ A SurfaceHeader """
+
+        @property
+        def geometry(self):
+            """ A SurfaceGeometry or None """
 
         @property
         def dataobj(self):
@@ -168,8 +176,8 @@ The basic API is as follows:
             2) self.dataobj.shape[0] == self.header.nfaces
             """
 
-        def load_header(self, pathlike):
-            """ Specify a header to a data-only image """
+        def load_geometry(self, pathlike):
+            """ Specify a geometry for a data-only image """
 
 
 Here we present common use cases:
@@ -205,22 +213,22 @@ 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.header.get_mesh(name=surface)
+        coords, triangles = img.geometry.get_mesh(name=surface)
 
         data = tstats.get_data()
 
         return plot_surf((triangles, coords), data)
 
     tstats = SurfaceImage.from_filename("/data/stats/hemi-L_contrast-taskVsBase_tstat.mgz")
-    tstats.load_header("/data/subjects/fsaverage5")
+    tstats.load_geometry("/data/subjects/fsaverage5")
     plot_surf_img(tstats)
 
-This assumes that ``load_header()`` will be able to identify a FreeSurfer subject directory.
+This assumes that ``load_geometry()`` will be able to identify a FreeSurfer subject directory.
 
 ******************
 * Open questions *
 ******************
 
-1) Can/should image and header objects be promiscuous? If I load a GIFTI header,
-   should ``header.load_vertex_data()`` accept a FreeSurfer curvature file, or
+1) Can/should image and geometry objects be promiscuous? If I load a GIFTI geometry,
+   should ``geometry.load_vertex_data()`` accept a FreeSurfer curvature file, or
    should there be a bit more ecosystem-local logic?

From 3e89a507205cad2e1714f3e8f728ffa16a357e24 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Wed, 20 Oct 2021 16:07:34 -0400
Subject: [PATCH 07/23] DOC: Smoothing example, typo

---
 doc/source/devel/biaps/biap_0009.rst | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index ed1982a063..870d4667f7 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -155,7 +155,7 @@ The basic API is as follows:
             """ Return a SurfaceImage with data corresponding to each face """
 
 
-    class SurfaceHeader(FileBaseeHeader):
+    class SurfaceHeader(FileBasedHeader):
         ...
 
 
@@ -202,6 +202,21 @@ Data images have their own metadata apart from geometry that needs handling in a
 * Data dtype
 
 
+Smoothing
+---------
+
+.. code-block:: python
+
+    bold = SurfaceImage.from_filename("/data/func/hemi-L_bold.func.gii")
+    bold.load_geometry("/data/anat/hemi-L_midthickness.surf.gii")
+    # Not implementing networkx weighted graph here, so assume we have a method
+    graph = bold.geometry.get_graph()
+    distances = distance_matrix(graph)  # n_coords x n_coords matrix
+    weights = normalize(gaussian(distances, sigma))
+    smoothed = bold.__class__(bold.get_fdata() @ weights, bold.header)
+    smoothed.to_filename(f"/data/func/hemi-L_smooth-{sigma}_bold.func.gii")
+
+
 Plotting
 --------
 

From b51e7a081c25036eedd4d395779c08ba92261010 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Wed, 20 Oct 2021 16:18:23 -0400
Subject: [PATCH 08/23] ENH: First pass at surfaceimage template classes

---
 nibabel/surfaceimages.py | 94 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 94 insertions(+)
 create mode 100644 nibabel/surfaceimages.py

diff --git a/nibabel/surfaceimages.py b/nibabel/surfaceimages.py
new file mode 100644
index 0000000000..88b404dad8
--- /dev/null
+++ b/nibabel/surfaceimages.py
@@ -0,0 +1,94 @@
+from nibabel.filebasedimages import FileBasedHeader
+from nibabel.dataobj_images import DataobjImage
+
+
+class SurfaceGeometry:
+    def __init__(self, meshes=None):
+        """ Surface header objects have access to an internal
+        ``_meshes`` dictionary that has keys that are mesh names.
+        The values may be any structure that permits the class to
+        provide coordinates and triangles on request.
+        """
+        if meshes is None:
+            meshes = {}
+        self._meshes = meshes
+        super().__init__(*args, **kwargs)
+
+    @property
+    def n_coords(self):
+        """ Number of coordinates (vertices)
+
+        The default implementation loads coordinates. Subclasses may
+        override with more efficient implementations.
+        """
+        return self.get_coords().shape[0]
+
+    @property
+    def n_triangles(self):
+        """ Number of triangles (faces)
+
+        The default implementation loads triangles. Subclasses may
+        override with more efficient implementations.
+        """
+        return self.get_triangles().shape[0]
+
+    def get_coords(self, name=None):
+        """ Nx3 array of coordinates in RAS+ space """
+        raise NotImplementedError
+
+    def get_triangles(self, name=None):
+        """ Mx3 array of indices into coordinate table """
+        raise NotImplementedError
+
+    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}``
+        """
+        return list(self._meshes.keys())
+
+    def decimate(self, *, ncoords=None, ratio=None):
+        """ Return a SurfaceHeader with a smaller number of vertices that
+        preserves the geometry of the original.
+
+        Please contribute a generic decimation algorithm at
+        https://github.com/nipy/nibabel
+        """
+        raise NotImplementedError
+
+    def load_vertex_data(self, pathlike):
+        """ Return a SurfaceImage with data corresponding to each vertex """
+        raise NotImplementedError
+
+    def load_face_data(self, pathlike):
+        """ Return a SurfaceImage with data corresponding to each face """
+        raise NotImplementedError
+
+
+class SurfaceHeader(FileBasedHeader):
+    """ Template class to implement SurfaceHeader protocol """
+    def get_geometry(self):
+        """ Generate ``SurfaceGeometry`` object from header object
+
+        If no default geometry can be provided, returns ``None``.
+        """
+        return None
+
+
+class SurfaceImage(DataobjImage):
+    header_class = SurfaceHeader
+
+    def __init__(self, dataobj, header=None, geometry=None, extra=None, file_map=None):
+        """ Initialize image
+
+        The image is a combination of 
+        """
+        super().__init__(dataobj, header=header, extra=extra, file_map=file_map)
+        if geometry is None:
+            geometry = self.header.get_geometry()
+        self._geometry = geometry
+
+    def load_geometry(self, pathlike):
+        """ Specify a header to a data-only image """

From 0f720556a9f25d8c16db297e471eeb38a0d672d4 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Thu, 21 Oct 2021 10:51:42 -0400
Subject: [PATCH 09/23] TEST: Build HDF5/Numpy-based surface classes

---
 nibabel/tests/test_surfaceimages.py | 110 ++++++++++++++++++++++++++++
 1 file changed, 110 insertions(+)
 create mode 100644 nibabel/tests/test_surfaceimages.py

diff --git a/nibabel/tests/test_surfaceimages.py b/nibabel/tests/test_surfaceimages.py
new file mode 100644
index 0000000000..d1eab3e3f8
--- /dev/null
+++ b/nibabel/tests/test_surfaceimages.py
@@ -0,0 +1,110 @@
+from nibabel.surfaceimages import SurfaceGeometry, SurfaceHeader, SurfaceImage
+from nibabel.optpkg import optional_package
+
+from nibabel.tests.test_filebasedimages import FBNumpyImage
+
+import numpy as np
+from pathlib import Path
+
+h5, has_h5py, _ = optional_package('h5py')
+
+
+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 __slicer__(self, slicer):
+        with h5.File(self.file_like, "r") as h5f:
+            return h5f[self.dataset_name][slicer]
+
+
+class H5Geometry(SurfaceGeometry):
+    """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 = h5f['topology']
+            for name, coords in h5f['coordinates'].items():
+                meshes[name] = (coords, triangles)
+        return klass(meshes)
+
+    
+    def to_filename(self, pathlike):
+        topology = None
+        coordinates = {}
+        for name, mesh in self.meshes.items():
+            coords, faces = mesh
+            if topology is None:
+                topology = faces
+            elif not np.array_equal(faces, topology):
+                raise ValueError("Inconsistent topology")
+            coordinates[name] = coords
+
+        with h5.File(pathlike, "w") as h5f:
+            h5f.create_dataset("/topology", topology)
+            for name, coord in coordinates.items():
+                h5f.create_dataset(f"/coordinates/{name}", coord)
+
+
+    def get_coords(self, name=None):
+        if name is None:
+            name = next(iter(self._meshes))
+        coords, _ = self._meshes[name]
+        return coords
+
+
+    def get_triangles(self, name=None):
+        if name is None:
+            name = next(iter(self._meshes))
+        _, triangles = self._meshes[name]
+        return triangles
+
+
+class NPSurfaceImage(SurfaceImage):
+    valid_exts = ('.npy',)
+    files_types = (('image', '.npy'),)
+
+    @classmethod
+    def from_file_map(klass, file_map):
+        with file_map['image'].get_prepare_fileobj('rb') as fobj:
+            arr = np.load(fobj)
+        return klass(arr)
+
+    def to_file_map(self, file_map=None):
+        file_map = self.file_map if file_map is None else file_map
+        with file_map['image'].get_prepare_fileobj('wb') as fobj:
+            np.save(fobj, self.arr)
+
+    def get_data_dtype(self):
+        return self.dataobj.dtype
+
+    def set_data_dtype(self, dtype):
+        self.dataobj = self.dataobj.astype(dtype)

From 408a2278e0cb81ef47c4cea839ee13e660e37599 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 08:18:09 -0400
Subject: [PATCH 10/23] DOC: Add VolumeGeometry stub

---
 doc/source/devel/biaps/biap_0009.rst | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 870d4667f7..392ccf2b2d 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -180,6 +180,22 @@ The basic API is as follows:
             """ Specify a geometry for a data-only image """
 
 
+To enable a similar interface to raveled voxel data:
+
+.. code-block:: python
+
+    class VolumeGeometry:
+        _affine  # Or _affines, if we want multiple, e.g. qform, sform
+        _shape
+        _ijk_coords
+
+        def n_coords(self):
+            """ Number of coordinates """
+
+        def get_coords(self):
+            """ Nx3 array of coordinates in RAS+ space """
+
+
 Here we present common use cases:
 
 

From 79a801e8f2017c0cf3b97f2c9bc62ee253a9c766 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 08:18:45 -0400
Subject: [PATCH 11/23] DOC: Fix header formatting

---
 doc/source/devel/biaps/biap_0009.rst | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 392ccf2b2d..f879fd4b54 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -200,7 +200,7 @@ Here we present common use cases:
 
 
 Modeling
---------
+========
 
 .. code-block:: python
 
@@ -219,7 +219,7 @@ Data images have their own metadata apart from geometry that needs handling in a
 
 
 Smoothing
----------
+=========
 
 .. code-block:: python
 
@@ -234,7 +234,7 @@ Smoothing
 
 
 Plotting
---------
+========
 
 Nilearn currently provides a
 `plot_surf <https://nilearn.github.io/modules/generated/nilearn.plotting.plot_surf.html>`_ function.
@@ -256,9 +256,10 @@ With the proposed API, we could interface as follows:
 
 This assumes that ``load_geometry()`` will be able to identify a FreeSurfer subject directory.
 
-******************
-* Open questions *
-******************
+
+**************
+Open questions
+**************
 
 1) Can/should image and geometry objects be promiscuous? If I load a GIFTI geometry,
    should ``geometry.load_vertex_data()`` accept a FreeSurfer curvature file, or

From eabc77fcf09bd174830eeee8de2622c36dbcc054 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 08:48:20 -0400
Subject: [PATCH 12/23] TEST: Example FreeSurfer subject (not fitting into
 class hierarchy)

---
 nibabel/tests/test_surfaceimages.py | 104 +++++++++++++++++++++++++++-
 1 file changed, 103 insertions(+), 1 deletion(-)

diff --git a/nibabel/tests/test_surfaceimages.py b/nibabel/tests/test_surfaceimages.py
index d1eab3e3f8..6f31f56e49 100644
--- a/nibabel/tests/test_surfaceimages.py
+++ b/nibabel/tests/test_surfaceimages.py
@@ -1,4 +1,6 @@
-from nibabel.surfaceimages import SurfaceGeometry, SurfaceHeader, SurfaceImage
+import os
+from nibabel.surfaceimages import *
+from nibabel.arrayproxy import ArrayProxy
 from nibabel.optpkg import optional_package
 
 from nibabel.tests.test_filebasedimages import FBNumpyImage
@@ -108,3 +110,103 @@ def get_data_dtype(self):
 
     def set_data_dtype(self, dtype):
         self.dataobj = self.dataobj.astype(dtype)
+
+
+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
+
+    @property
+    def coords(self):
+        ap = ArrayProxy(self._file_like, ((self.vnum, 3), ">f4", self.offset))
+        ap.order = 'C'
+        return ap
+
+    @property
+    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(SurfaceGeometry):
+    @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 get_coords(self, name=None):
+        if name is None:
+            name = self._default
+        return self._meshes[name].coords
+
+    def get_triangles(self, name=None):
+        if name is None:
+            name = self._default
+        return self._meshes[name].triangles
+
+    @property
+    def n_coords(self):
+        return self.meshes[self._default].vnum
+
+    @property
+    def n_triangles(self):
+        return self.meshes[self._default].fnum
+
+
+class FreeSurferSubject(Geometry):
+    def __init__(self, pathlike):
+        self._subject_dir = Path(pathlike)
+        surfs = self._subject_dir / "surf"
+        self._structures = {
+            "lh": FreeSurferHemisphere.from_filename(surfs / "lh.white"),
+            "rh": FreeSurferHemisphere.from_filename(surfs / "rh.white"),
+            }
+        super().__init__()

From 199547f927ead4f9ce0d2d0a7ae912cf20bbd0f3 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 11:39:13 -0400
Subject: [PATCH 13/23] DOC: Add concatenable structure proposal

---
 doc/source/devel/biaps/biap_0009.rst | 53 +++++++++++++++++++++++++---
 1 file changed, 48 insertions(+), 5 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index f879fd4b54..2317d6fcdb 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -120,18 +120,20 @@ The basic API is as follows:
 
 .. code-block:: python
 
-    class SurfaceGeometry:
+    class Geometry:
         @property
         def n_coords(self):
             """ Number of coordinates """
 
+        def get_coords(self, name=None):
+            """ Nx3 array of coordinates in RAS+ space """
+
+
+    class SurfaceGeometry(Geometry):
         @property
         def n_triangles(self):
             """ Number of faces """
 
-        def get_coords(self, name=None):
-            """ Nx3 array of coordinates in RAS+ space """
-
         def get_triangles(self, name=None):
             """ Mx3 array of indices into coordinate table """
 
@@ -184,7 +186,7 @@ To enable a similar interface to raveled voxel data:
 
 .. code-block:: python
 
-    class VolumeGeometry:
+    class VolumeGeometry(Geometry):
         _affine  # Or _affines, if we want multiple, e.g. qform, sform
         _shape
         _ijk_coords
@@ -264,3 +266,44 @@ Open questions
 1) Can/should image and geometry objects be promiscuous? If I load a GIFTI geometry,
    should ``geometry.load_vertex_data()`` accept a FreeSurfer curvature file, or
    should there be a bit more ecosystem-local logic?
+
+
+********************************
+Concatenable structures proposal
+********************************
+
+A brain is typically described with multiple structures,
+such as left and right hemispheres and a volumetric subcortical region.
+We will consider two reference cases:
+
+1) FreeSurfer, in which a subject directory can be used to retrieve
+   multiple surfaces or volumetric regions. Data arrays are specific to a
+   single structure.
+2) CIFTI-2, which permits an axis to have type ``CIFTI_INDEX_TYPE_BRAIN_MODELS``,
+   indicating that each index along an axis corresponds to a vertex or voxel in
+   an associated surface or volume.
+   The data array thus may span multiple structures.
+
+In the former case, the structures may be considered an unordered collection;
+in the latter, the sequence of structures directly corresponds to segments of
+the data array.
+
+.. code-block:: python
+
+    class GeometryCollection:
+        def get_geometry(self, name):
+            """ Return a geometry by name """
+
+        @property
+        def names(self):
+            """ A set of structure names """
+
+        @classmethod
+        def from_spec(klass, pathlike):
+        """ Load a collection of geometries from a specification, broadly construed. """
+
+
+    class GeometrySequence(GeometryCollection, Geometry):
+        # Inherit and override get_coords, n_coords from Geometry
+        def get_indices(self, name):
+            """ Return indices for data array corresponding to a named geometry """

From efef0273ce722f1ff9ed4e4ec82d765252de6a9e Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 11:40:44 -0400
Subject: [PATCH 14/23] ENH: Add structure collection API

---
 nibabel/surfaceimages.py | 54 ++++++++++++++++++++++++++++++++++++++--
 1 file changed, 52 insertions(+), 2 deletions(-)

diff --git a/nibabel/surfaceimages.py b/nibabel/surfaceimages.py
index 88b404dad8..e01a70b3a2 100644
--- a/nibabel/surfaceimages.py
+++ b/nibabel/surfaceimages.py
@@ -2,7 +2,13 @@
 from nibabel.dataobj_images import DataobjImage
 
 
-class SurfaceGeometry:
+class Geometry:
+    def get_coords(self, name=None):
+        """ Nx3 array of coordinates in RAS+ space """
+        raise NotImplementedError
+
+
+class SurfaceGeometry(Geometry):
     def __init__(self, meshes=None):
         """ Surface header objects have access to an internal
         ``_meshes`` dictionary that has keys that are mesh names.
@@ -12,7 +18,7 @@ def __init__(self, meshes=None):
         if meshes is None:
             meshes = {}
         self._meshes = meshes
-        super().__init__(*args, **kwargs)
+        super().__init__()
 
     @property
     def n_coords(self):
@@ -92,3 +98,47 @@ def __init__(self, dataobj, header=None, geometry=None, extra=None, file_map=Non
 
     def load_geometry(self, pathlike):
         """ Specify a header to a data-only image """
+
+
+class GeometryCollection:
+    def __init__(self, structures=()):
+        self._structures = dict(structures)
+
+    def get_structure(self, name):
+        return self._structures[name]
+
+    @property
+    def names(self):
+        return list(self._structures)
+
+    @classmethod
+    def from_spec(klass, pathlike):
+        """ Load a collection of geometries from a specification, broadly construed. """
+        raise NotImplementedError
+
+
+class GeometrySequence(GeometryCollection, Geometry):
+    def __init__(self, structures=()):
+        super().__init__(structures)
+        self._indices = {}
+        next_index = 0
+        for name, struct in self._structures.items():
+            end = next_index + struct.n_coords
+            self._indices[name] = slice(next_index, end)
+            next_index = end + 1
+
+    def get_indices(self, *names):
+        if len(names) == 1:
+            return self._indices[name]
+        return [self._indices[name] for name in names]
+
+    # def get_structures(self, *, names=None, indices=None):
+    #     """ We probably want some way to get a subset of structures """
+
+    def get_coords(self, name=None):
+        return np.vstack([struct.get_coords(name=name)
+                          for struct in self._structures.values()])
+
+    @property
+    def n_coords(self):
+        return sum(struct.n_coords for struct in self._structures.values())

From d88c7c60133aca2db6751381372f002d5e193207 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 11:41:11 -0400
Subject: [PATCH 15/23] TEST: Rewrite FreeSurferSubject as GeometryCollection

---
 nibabel/tests/test_surfaceimages.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/nibabel/tests/test_surfaceimages.py b/nibabel/tests/test_surfaceimages.py
index 6f31f56e49..3429f7d25b 100644
--- a/nibabel/tests/test_surfaceimages.py
+++ b/nibabel/tests/test_surfaceimages.py
@@ -201,8 +201,17 @@ def n_triangles(self):
         return self.meshes[self._default].fnum
 
 
-class FreeSurferSubject(Geometry):
-    def __init__(self, pathlike):
+class FreeSurferSubject(GeometryCollection):
+    @classmethod
+    def from_subject(klass, subject_id, subjects_dir=None):
+        """ Load a FreeSurfer subject by ID """
+        if subjects_dir is None:
+            subjects_dir = os.environ["SUBJECTS_DIR"]
+        return klass.from_directory(Path(subjects_dir) / subject_id)
+
+    @classmethod
+    def from_spec(klass, pathlike):
+        """ Load a FreeSurfer subject from its directory structure """
         self._subject_dir = Path(pathlike)
         surfs = self._subject_dir / "surf"
         self._structures = {

From 25e6d523c7e686ca8889bd175973ed2482478184 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 5 Nov 2021 11:59:18 -0400
Subject: [PATCH 16/23] ENH: Possible VolumeGeometry

---
 nibabel/surfaceimages.py | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/nibabel/surfaceimages.py b/nibabel/surfaceimages.py
index e01a70b3a2..b04e521d9a 100644
--- a/nibabel/surfaceimages.py
+++ b/nibabel/surfaceimages.py
@@ -1,3 +1,4 @@
+from nibabel.affines import apply_affine
 from nibabel.filebasedimages import FileBasedHeader
 from nibabel.dataobj_images import DataobjImage
 
@@ -100,6 +101,29 @@ def load_geometry(self, pathlike):
         """ Specify a header to a data-only image """
 
 
+class VolumeGeometry(Geometry):
+    def __init__(self, affines, *, indices=None):
+        try:
+            self._affines = dict(affines)
+        except TypeError:
+            self._affines = {"affine": np.array(affines)}
+        self._default = next(iter(self._affines))
+
+        self._indices = indices
+
+    def get_coords(self, name=None):
+        if name is None:
+            name = self._default
+        return apply_affine(self._affines[name], self._indices)
+
+    @property
+    def n_coords(self):
+        return self._indices.shape[0]
+
+    def get_indices(self):
+        return self._indices
+
+
 class GeometryCollection:
     def __init__(self, structures=()):
         self._structures = dict(structures)

From e6be497ff02a482a68baca5f036fb19c92e124c8 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Sun, 7 Nov 2021 19:28:28 -0500
Subject: [PATCH 17/23] STY: Geometry -> Pointset

---
 nibabel/surfaceimages.py | 52 ++++++++++++++++++++++++++--------------
 1 file changed, 34 insertions(+), 18 deletions(-)

diff --git a/nibabel/surfaceimages.py b/nibabel/surfaceimages.py
index b04e521d9a..8120510d6b 100644
--- a/nibabel/surfaceimages.py
+++ b/nibabel/surfaceimages.py
@@ -3,13 +3,42 @@
 from nibabel.dataobj_images import DataobjImage
 
 
-class Geometry:
+class Pointset:
+    """ A collection of related sets of coordinates in 3D space.
+
+    Coordinates are in RAS+ orientation, i.e., $(x, y, z)$ refers to
+    a point $x\mathrm{mm}$ right, $y\mathrm{mm}$ anterior and $z\mathrm{mm}$
+    superior of the origin.
+
+    The order of coordinates should be assumed to be significant.
+    """
+
     def get_coords(self, name=None):
-        """ Nx3 array of coordinates in RAS+ space """
+        """ Get coordinate data in RAS+ space
+
+        Parameters
+        ----------
+        name : str
+            Name of one of a family of related coordinates.
+            For example, ``"white"`` and ``"pial"`` surfaces
+
+        Returns
+        -------
+        coords : (N, 3) array-like
+            Coordinates in RAS+ space
+        """
         raise NotImplementedError
 
+    @property
+    def n_coords(self):
+        """ Number of coordinates
+
+        The default implementation loads coordinates. Subclasses may
+        override with more efficient implementations.
+        """
+        return self.get_coords().shape[0]
 
-class SurfaceGeometry(Geometry):
+class SurfaceGeometry(Pointset):
     def __init__(self, meshes=None):
         """ Surface header objects have access to an internal
         ``_meshes`` dictionary that has keys that are mesh names.
@@ -21,15 +50,6 @@ def __init__(self, meshes=None):
         self._meshes = meshes
         super().__init__()
 
-    @property
-    def n_coords(self):
-        """ Number of coordinates (vertices)
-
-        The default implementation loads coordinates. Subclasses may
-        override with more efficient implementations.
-        """
-        return self.get_coords().shape[0]
-
     @property
     def n_triangles(self):
         """ Number of triangles (faces)
@@ -39,10 +59,6 @@ def n_triangles(self):
         """
         return self.get_triangles().shape[0]
 
-    def get_coords(self, name=None):
-        """ Nx3 array of coordinates in RAS+ space """
-        raise NotImplementedError
-
     def get_triangles(self, name=None):
         """ Mx3 array of indices into coordinate table """
         raise NotImplementedError
@@ -101,7 +117,7 @@ def load_geometry(self, pathlike):
         """ Specify a header to a data-only image """
 
 
-class VolumeGeometry(Geometry):
+class VolumeGeometry(Pointset):
     def __init__(self, affines, *, indices=None):
         try:
             self._affines = dict(affines)
@@ -141,7 +157,7 @@ def from_spec(klass, pathlike):
         raise NotImplementedError
 
 
-class GeometrySequence(GeometryCollection, Geometry):
+class GeometrySequence(GeometryCollection, Pointset):
     def __init__(self, structures=()):
         super().__init__(structures)
         self._indices = {}

From 3f62a80250fa2f6363bd678a6c16d71e002cdaed Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 19 Nov 2021 08:19:27 -0500
Subject: [PATCH 18/23] Rename SurfaceGeometry to TriangularMesh

---
 doc/source/devel/biaps/biap_0009.rst |  4 ++--
 nibabel/surfaceimages.py             | 13 ++++++++-----
 nibabel/tests/test_surfaceimages.py  |  7 ++-----
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index 2317d6fcdb..fa2c1cb528 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -129,7 +129,7 @@ The basic API is as follows:
             """ Nx3 array of coordinates in RAS+ space """
 
 
-    class SurfaceGeometry(Geometry):
+    class TriangularMesh(Geometry):
         @property
         def n_triangles(self):
             """ Number of faces """
@@ -146,7 +146,7 @@ The basic API is as follows:
             """
 
         def decimate(self, *, n_coords=None, ratio=None):
-            """ Return a SurfaceHeader with a smaller number of vertices that
+            """ 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
 
diff --git a/nibabel/surfaceimages.py b/nibabel/surfaceimages.py
index 8120510d6b..2195c77373 100644
--- a/nibabel/surfaceimages.py
+++ b/nibabel/surfaceimages.py
@@ -38,9 +38,12 @@ def n_coords(self):
         """
         return self.get_coords().shape[0]
 
-class SurfaceGeometry(Pointset):
+
+class TriangularMesh(Pointset):
+    """ A triangular mesh is a description of a surface tesselated into triangles """
+
     def __init__(self, meshes=None):
-        """ Surface header objects have access to an internal
+        """ Triangular mesh objects have access to an internal
         ``_meshes`` dictionary that has keys that are mesh names.
         The values may be any structure that permits the class to
         provide coordinates and triangles on request.
@@ -73,7 +76,7 @@ def get_names(self):
         return list(self._meshes.keys())
 
     def decimate(self, *, ncoords=None, ratio=None):
-        """ Return a SurfaceHeader with a smaller number of vertices that
+        """ Return a TriangularMesh with a smaller number of vertices that
         preserves the geometry of the original.
 
         Please contribute a generic decimation algorithm at
@@ -93,7 +96,7 @@ def load_face_data(self, pathlike):
 class SurfaceHeader(FileBasedHeader):
     """ Template class to implement SurfaceHeader protocol """
     def get_geometry(self):
-        """ Generate ``SurfaceGeometry`` object from header object
+        """ Generate ``TriangularMesh`` object from header object
 
         If no default geometry can be provided, returns ``None``.
         """
@@ -157,7 +160,7 @@ def from_spec(klass, pathlike):
         raise NotImplementedError
 
 
-class GeometrySequence(GeometryCollection, Pointset):
+class PointsetSequence(GeometryCollection, Pointset):
     def __init__(self, structures=()):
         super().__init__(structures)
         self._indices = {}
diff --git a/nibabel/tests/test_surfaceimages.py b/nibabel/tests/test_surfaceimages.py
index 3429f7d25b..59589b2f0d 100644
--- a/nibabel/tests/test_surfaceimages.py
+++ b/nibabel/tests/test_surfaceimages.py
@@ -45,7 +45,7 @@ def __slicer__(self, slicer):
             return h5f[self.dataset_name][slicer]
 
 
-class H5Geometry(SurfaceGeometry):
+class H5Geometry(TriangularMesh):
     """Simple Geometry file structure that combines a single topology
     with one or more coordinate sets
     """
@@ -57,7 +57,6 @@ def from_filename(klass, pathlike):
             for name, coords in h5f['coordinates'].items():
                 meshes[name] = (coords, triangles)
         return klass(meshes)
-
     
     def to_filename(self, pathlike):
         topology = None
@@ -75,14 +74,12 @@ def to_filename(self, pathlike):
             for name, coord in coordinates.items():
                 h5f.create_dataset(f"/coordinates/{name}", coord)
 
-
     def get_coords(self, name=None):
         if name is None:
             name = next(iter(self._meshes))
         coords, _ = self._meshes[name]
         return coords
 
-
     def get_triangles(self, name=None):
         if name is None:
             name = next(iter(self._meshes))
@@ -163,7 +160,7 @@ def triangles(self):
         return ap
 
 
-class FreeSurferHemisphere(SurfaceGeometry):
+class FreeSurferHemisphere(TriangularMesh):
     @classmethod
     def from_filename(klass, pathlike):
         path = Path(pathlike)

From dcd2050b21c27cab37b5b77487858b21e5ba44a3 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 19 Nov 2021 08:20:08 -0500
Subject: [PATCH 19/23] FIX: FreeSurfer example implementation

---
 nibabel/tests/test_surfaceimages.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/nibabel/tests/test_surfaceimages.py b/nibabel/tests/test_surfaceimages.py
index 59589b2f0d..a822cfcb22 100644
--- a/nibabel/tests/test_surfaceimages.py
+++ b/nibabel/tests/test_surfaceimages.py
@@ -204,15 +204,17 @@ def from_subject(klass, subject_id, subjects_dir=None):
         """ Load a FreeSurfer subject by ID """
         if subjects_dir is None:
             subjects_dir = os.environ["SUBJECTS_DIR"]
-        return klass.from_directory(Path(subjects_dir) / subject_id)
+        return klass.from_spec(Path(subjects_dir) / subject_id)
 
     @classmethod
     def from_spec(klass, pathlike):
         """ Load a FreeSurfer subject from its directory structure """
-        self._subject_dir = Path(pathlike)
-        surfs = self._subject_dir / "surf"
-        self._structures = {
+        subject_dir = Path(pathlike)
+        surfs = subject_dir / "surf"
+        structures = {
             "lh": FreeSurferHemisphere.from_filename(surfs / "lh.white"),
             "rh": FreeSurferHemisphere.from_filename(surfs / "rh.white"),
             }
-        super().__init__()
+        subject = super().__init__(structures)
+        subject._subject_dir = subject_dir
+        return subject

From ff19edc1b1079c912565eac4f129b0c4cf4405f1 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 19 Nov 2021 08:32:48 -0500
Subject: [PATCH 20/23] ENH: Flesh out VolumeGeometry

---
 nibabel/surfaceimages.py | 43 ++++++++++++++++++++++++++++++++++++++--
 1 file changed, 41 insertions(+), 2 deletions(-)

diff --git a/nibabel/surfaceimages.py b/nibabel/surfaceimages.py
index 2195c77373..719bae9047 100644
--- a/nibabel/surfaceimages.py
+++ b/nibabel/surfaceimages.py
@@ -121,14 +121,38 @@ def load_geometry(self, pathlike):
 
 
 class VolumeGeometry(Pointset):
-    def __init__(self, affines, *, indices=None):
+    """
+
+    .. testsetup:
+
+        >>> from nibabel.tests.nibabel_data import needs_nibabel_data
+        >>> needs_nibabel_data('nitest-cifti2')()
+
+    >>> import nibabel as nb
+    >>> from nibabel.tests.nibabel_data import get_nibabel_data
+    >>> from pathlib import Path
+    >>> dscalar = nb.load(Path(get_nibabel_data()) / 'nitest-cifti2' / "ones.dscalar.nii")
+    >>> brainmodel_axis = dscalar.header.get_axis(1)
+    >>> volgeom = VolumeGeometry(
+    ...     brainmodel_axis.affine,
+    ...     brainmodel_axis.voxel[brainmodel_axis.volume_mask],
+    ...     brainmodel_axis.volume_shape)
+    """
+
+    def __init__(self, affines, indices, shape=None):
         try:
             self._affines = dict(affines)
         except TypeError:
             self._affines = {"affine": np.array(affines)}
         self._default = next(iter(self._affines))
 
-        self._indices = indices
+        self._indices = np.array(indices)
+        if self._indices.ndim != 2 or self._indices.shape[1] != 3:
+            raise ValueError("Indices must be an Nx3 matrix")
+
+        if shape is None:
+            shape = self._indices.max(axis=0)
+        self._shape = shape
 
     def get_coords(self, name=None):
         if name is None:
@@ -142,6 +166,21 @@ def n_coords(self):
     def get_indices(self):
         return self._indices
 
+    def to_mask(self, name=None, shape=None):
+        if name is None:
+            name = self._default
+        if shape is None:
+            shape = self._shape
+        dataobj = np.zeros(shape, dtype=np.bool_)
+        dataobj[self._indices] = True
+        return SpatialImage(dataobj, self._affines[name])
+
+    @classmethod
+    def from_mask(klass, img):
+        affine = img.affine
+        indices = np.vstack(np.where(img.dataobj)).T
+        return klass(affine, indices=indices)
+
 
 class GeometryCollection:
     def __init__(self, structures=()):

From 4af4897fe365c970438499a483bdf99d53c7d7c5 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Fri, 19 Nov 2021 08:33:54 -0500
Subject: [PATCH 21/23] BIAP: Add SurfaceHeader.get_geometry() method

---
 doc/source/devel/biaps/biap_0009.rst | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index fa2c1cb528..af1a212c07 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -158,7 +158,8 @@ The basic API is as follows:
 
 
     class SurfaceHeader(FileBasedHeader):
-        ...
+        def get_geometry(self):
+            """ Construct a geometry from the header if possible, or return None """
 
 
     class SurfaceImage(FileBasedImage):

From a3adfe74298d7f72c865030fa7927355a60e47b2 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Thu, 13 Jan 2022 20:30:27 -0500
Subject: [PATCH 22/23] Rename Geometry -> Pointset

---
 doc/source/devel/biaps/biap_0009.rst | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index af1a212c07..a26922b127 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -120,7 +120,7 @@ The basic API is as follows:
 
 .. code-block:: python
 
-    class Geometry:
+    class Pointset:
         @property
         def n_coords(self):
             """ Number of coordinates """
@@ -129,7 +129,7 @@ The basic API is as follows:
             """ Nx3 array of coordinates in RAS+ space """
 
 
-    class TriangularMesh(Geometry):
+    class TriangularMesh(Pointset):
         @property
         def n_triangles(self):
             """ Number of faces """
@@ -169,7 +169,7 @@ The basic API is as follows:
 
         @property
         def geometry(self):
-            """ A SurfaceGeometry or None """
+            """ A Pointset or None """
 
         @property
         def dataobj(self):
@@ -187,7 +187,7 @@ To enable a similar interface to raveled voxel data:
 
 .. code-block:: python
 
-    class VolumeGeometry(Geometry):
+    class VolumeGeometry(Pointset):
         _affine  # Or _affines, if we want multiple, e.g. qform, sform
         _shape
         _ijk_coords
@@ -304,7 +304,7 @@ the data array.
         """ Load a collection of geometries from a specification, broadly construed. """
 
 
-    class GeometrySequence(GeometryCollection, Geometry):
-        # Inherit and override get_coords, n_coords from Geometry
+    class PointsetSequence(GeometryCollection, Pointset):
+        # Inherit and override get_coords, n_coords from Pointset
         def get_indices(self, name):
             """ Return indices for data array corresponding to a named geometry """

From e278981f60ee9324c3ddf78f175cf90d26918222 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Thu, 10 Feb 2022 21:10:36 -0500
Subject: [PATCH 23/23] DOC: Commit current thinking on BIAP0009

---
 doc/source/devel/biaps/biap_0009.rst | 251 +++++++++++++++------------
 1 file changed, 138 insertions(+), 113 deletions(-)

diff --git a/doc/source/devel/biaps/biap_0009.rst b/doc/source/devel/biaps/biap_0009.rst
index a26922b127..897bd1c3ba 100644
--- a/doc/source/devel/biaps/biap_0009.rst
+++ b/doc/source/devel/biaps/biap_0009.rst
@@ -1,8 +1,8 @@
 .. _biap9:
 
-#############################
-BIAP9 - The Surface Image API
-#############################
+################################
+BIAP9 - The Coordinate Image API
+################################
 
 :Author: Chris Markiewicz
 :Status: Draft
@@ -18,7 +18,7 @@ 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
-requires the location of each sample to be explicitly represented by a
+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.
 
@@ -42,10 +42,12 @@ For the purposes of this BIAP, the following terms are used:
 * 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
-* Patch - a connected subset of vertices (can be the full topology)
+* Parcel - a subset of vertices; can be the full topology. Special cases include:
+  * Patch - a connected subspace
+  * Decimated mesh - a subspace that has a desired density of vertices
+* Subspace sequence - an ordered set of subspaces
 * Data array - an n-dimensional array with one axis corresponding to the
-  vertices (typical) OR faces (more rare) in a patch
-
+  vertices (typical) OR faces (more rare) in a patch sequence
 
 Currently supported surface formats
 ===================================
@@ -65,14 +67,29 @@ Currently supported surface formats
      coordinates, topology, or data (further subdivided by type and intent)
 * CIFTI-2: :py:class:`~nibabel.cifti2.cifti2.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 matrix
    * Geometry 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
 
+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 a SurfaceImage API
-*********************************
+*****************************************
+Desiderata for an API supporting surfaces
+*****************************************
 
 The following are provisional guiding principles 
 
@@ -116,7 +133,55 @@ make sense for NiBabel to provide methods.
 Proposal
 ********
 
-The basic API is as follows:
+A ``CoordinateImage`` is an N-dimensional array, where one axis corresponds
+to a sequence of points in one or more parcels.
+
+.. code-block:: python
+
+    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:
 
 .. code-block:: python
 
@@ -150,54 +215,40 @@ The basic API is as follows:
             preserves the geometry of the original """
             # To be overridden when a format provides optimization opportunities
 
-        def load_vertex_data(self, filename):
-            """ Return a SurfaceImage with data corresponding to each vertex """
-
-        def load_face_data(self, filename):
-            """ Return a SurfaceImage with data corresponding to each face """
-
-
-    class SurfaceHeader(FileBasedHeader):
-        def get_geometry(self):
-            """ Construct a geometry from the header if possible, or return None """
 
+    class NdGrid(Pointset):
+        """
+        Attributes
+        ----------
+        shape : 3-tuple
+            number of coordinates in each dimension of grid
+        """
+        def get_affine(self, name=None):
+            """ 4x4 array """
 
-    class SurfaceImage(FileBasedImage):
-        @property
-        def header(self):
-            """ A SurfaceHeader """
-
-        @property
-        def geometry(self):
-            """ A Pointset or None """
-
-        @property
-        def dataobj(self):
-            """ An ndarray or ArrayProxy with one of the following properties:
-            
-            1) self.dataobj.shape[0] == self.header.ncoords
-            2) self.dataobj.shape[0] == self.header.nfaces
-            """
-
-        def load_geometry(self, pathlike):
-            """ Specify a geometry for a data-only image """
 
+The ``NdGrid`` class allows raveled volumetric data to be treated the same as
+triangular mesh or other coordinate data.
 
-To enable a similar interface to raveled voxel data:
+Finally, a structure for containing a collection of related geometric files is
+defined:
 
 .. code-block:: python
 
-    class VolumeGeometry(Pointset):
-        _affine  # Or _affines, if we want multiple, e.g. qform, sform
-        _shape
-        _ijk_coords
-
-        def n_coords(self):
-            """ Number of coordinates """
+    class GeometryCollection:
+        """
+        Attributes
+        ----------
+        structures : dict
+            Mapping from structure names to ``Pointset``
+        """
 
-        def get_coords(self):
-            """ Nx3 array of coordinates in RAS+ space """
+        @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:
 
@@ -209,30 +260,35 @@ Modeling
 
     from nilearn.glm.first_level import make_first_level_design_matrix, run_glm
 
-    bold = SurfaceImage.from_filename("/data/func/hemi-L_bold.func.gii")
+    bold = CoordinateImage.from_filename("/data/func/hemi-L_bold.func.gii")
     dm = make_first_level_design_matrix(...)
-    labels, results = run_glm(bold.get_data(), dm)
-    betas = bold.__class__(results["betas"], bold.header)
+    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")
 
-Data images have their own metadata apart from geometry that needs handling in a header:
-
-* Axis labels (time series, labels, tensors)
-* Data dtype
+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
 =========
 
 .. code-block:: python
 
-    bold = SurfaceImage.from_filename("/data/func/hemi-L_bold.func.gii")
-    bold.load_geometry("/data/anat/hemi-L_midthickness.surf.gii")
-    # Not implementing networkx weighted graph here, so assume we have a method
-    graph = bold.geometry.get_graph()
-    distances = distance_matrix(graph)  # n_coords x n_coords matrix
+    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))
-    smoothed = bold.__class__(bold.get_fdata() @ weights, bold.header)
+    # 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")
 
 
@@ -247,64 +303,33 @@ 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.geometry.get_mesh(name=surface)
+        coords, triangles = img.coordaxis.parcels[0].get_mesh(name=surface)
 
-        data = tstats.get_data()
+        data = img.get_fdata()
 
         return plot_surf((triangles, coords), data)
 
-    tstats = SurfaceImage.from_filename("/data/stats/hemi-L_contrast-taskVsBase_tstat.mgz")
-    tstats.load_geometry("/data/subjects/fsaverage5")
+    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)
 
-This assumes that ``load_geometry()`` will be able to identify a FreeSurfer subject directory.
-
-
-**************
-Open questions
-**************
-
-1) Can/should image and geometry objects be promiscuous? If I load a GIFTI geometry,
-   should ``geometry.load_vertex_data()`` accept a FreeSurfer curvature file, or
-   should there be a bit more ecosystem-local logic?
-
-
-********************************
-Concatenable structures proposal
-********************************
-
-A brain is typically described with multiple structures,
-such as left and right hemispheres and a volumetric subcortical region.
-We will consider two reference cases:
-
-1) FreeSurfer, in which a subject directory can be used to retrieve
-   multiple surfaces or volumetric regions. Data arrays are specific to a
-   single structure.
-2) CIFTI-2, which permits an axis to have type ``CIFTI_INDEX_TYPE_BRAIN_MODELS``,
-   indicating that each index along an axis corresponds to a vertex or voxel in
-   an associated surface or volume.
-   The data array thus may span multiple structures.
-
-In the former case, the structures may be considered an unordered collection;
-in the latter, the sequence of structures directly corresponds to segments of
-the data array.
+Subsampling CIFTI-2
+===================
 
 .. code-block:: python
 
-    class GeometryCollection:
-        def get_geometry(self, name):
-            """ Return a geometry by name """
-
-        @property
-        def names(self):
-            """ A set of structure names """
-
-        @classmethod
-        def from_spec(klass, pathlike):
-        """ Load a collection of geometries from a specification, broadly construed. """
+    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)
 
-    class PointsetSequence(GeometryCollection, Pointset):
-        # Inherit and override get_coords, n_coords from Pointset
-        def get_indices(self, name):
-            """ Return indices for data array corresponding to a named geometry """
+    # Now load geometry so we can plot
+    wbspec = CaretSpec("fsLR.wb.spec")
+    dlpfc_img.coordaxis.load_structures(wbspec)
+    ...