ENH,EXA add function to upsample from fsaverageX to fsaverage (#519)

* ENH,EXA add function to upsample from fsaverageX to fsaverage

* EXA add more comments

cortex/freesurfer.py

@@ -1,27 +1,26 @@
 """Contains functions for interfacing with freesurfer
 """
 from __future__ import print_function
-import os
+
 import copy
+import os
+import shlex
 import shutil
 import struct
+import subprocess as sp
 import tempfile
 import warnings
-import shlex
-import subprocess as sp
 from builtins import input
+from tempfile import NamedTemporaryFile
 
-import numpy as np
 import nibabel
+import numpy as np
 from nibabel import gifti
-from tempfile import NamedTemporaryFile
-from scipy.spatial import KDTree
 from scipy.linalg import lstsq
 from scipy.sparse import coo_matrix
+from scipy.spatial import KDTree
 
-
-from . import database
-from . import anat
+from . import anat, database
 
 
 def get_paths(fs_subject, hemi, type="patch", freesurfer_subject_dir=None):
@@ -895,7 +894,7 @@ def read_dot(fname, pts):
 def write_decimated(path, pts, polys):
     """
     """
-    from .polyutils import decimate, boundary_edges
+    from .polyutils import boundary_edges, decimate
     dpts, dpolys = decimate(pts, polys)
     write_surf(path+'.smoothwm', dpts, dpolys)
     edges = boundary_edges(dpolys)
@@ -1012,6 +1011,86 @@ def stretch_mwall(pts, polys, mwall):
     pts[mwall, 2] = radius * np.sin(angles) + center[2]
     return SpringLayout(pts, polys, inflated, pins=mwall)
 
+
+def upsample_to_fsaverage(
+        data, data_space="fsaverage6", freesurfer_subjects_dir=None
+):
+    """Project data from fsaverage6 (or other fsaverage surface) to fsaverage to
+    visualize it in pycortex.
+
+    Parameters
+    ----------
+    data : array (n_samples, n_vertices)
+        Data in space `space`. The first n_vertices/2 vertices correspond to the left
+        hemisphere, and the last n_vertices/2 vertices correspond to the right
+        hemisphere.
+    data_space : str
+        One of fsaverage[1-6], corresponding to the source template space of `data`.
+    freesurfer_subjects_dir : str or None
+        Path to Freesurfer subjects directory. If None, defaults to the value of the
+        environment variable $SUBJECTS_DIR.
+
+    Returns
+    -------
+    projected_data : array (n_samples, 327684)
+        Data projected to fsaverage(7).
+
+    Notes
+    -----
+    Data in the lower resolution fsaverage template is upsampled to the full resolution
+    fsaverage template by nearest-neighbor interpolation. To project the data from a 
+    lower resolution version of fsaverage, this code exploits the structure of fsaverage 
+    surfaces. (That is, each hemisphere in fsaverage6 corresponds to the first 
+    40,962 vertices of fsaverage; fsaverage5 corresponds to the first 10,242 vertices of 
+    fsaverage, etc.)
+    """
+
+
+    def get_n_vertices_ico(icoorder):
+        return 4 ** icoorder * 10 + 2
+
+    ico_order = int(data_space[-1])
+    n_ico_vertices = get_n_vertices_ico(ico_order)
+    ndim = data.ndim
+    data = np.atleast_2d(data)
+    _, n_vertices = data.shape
+    if n_vertices != 2 * n_ico_vertices:
+        raise ValueError(
+            f"data has {n_vertices} vertices, but {2 * n_ico_vertices} "
+            f"are expected for both hemispheres in {data_space}"
+        )
+
+    if freesurfer_subjects_dir is None:
+        freesurfer_subjects_dir = os.environ.get("SUBJECTS_DIR", None)
+    if freesurfer_subjects_dir is None:
+        raise ValueError(
+            "freesurfer_subjects_dir must be specified or $SUBJECTS_DIR must be set"
+        )
+
+    data_hemi = np.split(data, 2, axis=-1)
+    hemis = ["lh", "rh"]
+    projected_data = []
+    for i, (hemi, dt) in enumerate(zip(hemis, data_hemi)):
+        # Load fsaverage sphere for this hemisphere
+        pts, faces = nibabel.freesurfer.read_geometry(
+            os.path.join(
+                freesurfer_subjects_dir, "fsaverage", "surf", f"{hemi}.sphere.reg"
+            )
+        )
+        # build kdtree using only vertices in reduced fsaverage surface
+        kdtree = KDTree(pts[:n_ico_vertices])
+        # figure out neighbors in reduced version for all other vertices in fsaverage
+        _, neighbors = kdtree.query(pts[n_ico_vertices:], k=1)
+        # now simply fill remaining vertices with original values
+        projected_data.append(
+            np.concatenate([dt, dt[:, neighbors]], axis=-1)
+        )
+    projected_data = np.hstack(projected_data)
+    if ndim == 1:
+        projected_data = projected_data[0]
+    return projected_data
+
+
 # aseg partition labels (up to 256 only)
 fs_aseg_dict = {'Unknown': 0,
                 'Left-Cerebral-Exterior': 1,

examples/fsaverage/README.txt

@@ -0,0 +1,4 @@
+Examples with fsaverage
+------------------------------
+
+Examples showing how to use PyCortex to visualize data on fsaverage.

examples/fsaverage/upsample_to_fsaverage.py

@@ -0,0 +1,36 @@
+"""
+===================
+Upsample data from a lower resolution fsaverage template to fsaverage for visualization
+===================
+
+This example shows how data in a lower resolution fsaverage template 
+(e.g., fsaverage5 or fsaverage6) can be upsampled to the high resolution fsaverage 
+template for visualization.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import cortex
+
+subject = "fsaverage"
+
+# First we check if the fsaverage template is already in the pycortex filestore. If not,
+# we download the template from the web and add it to the filestore.
+if subject not in cortex.db.subjects:
+    cortex.download_subject(subject)
+
+# Next we create some data on fsaverage5. Each hemisphere has 10242 vertices.
+n_vertices_fsaverage5 = 10242
+data_fs5 = np.arange(1, n_vertices_fsaverage5 + 1)
+# We concatenate the data to itself to create a vector of length 20484, corresponding to
+# the two hemispheres together.
+data_fs5 = np.concatenate((data_fs5, data_fs5))
+# Finally, we upsample the data to fsaverage.
+data_fs7 = cortex.freesurfer.upsample_to_fsaverage(data_fs5, "fsaverage5")
+
+# Now that the data is in the fsaverage template, we can visualize it in PyCortex as any
+# other vertex dataset.
+vtx = cortex.Vertex(data_fs7, subject, vmin=0, vmax=n_vertices_fsaverage5, cmap="turbo")
+cortex.quickshow(vtx, with_curvature=False, with_colorbar=False)
+plt.show()