diff --git a/cortex/freesurfer.py b/cortex/freesurfer.py index 52f8a9fe..095e89c0 100644 --- a/cortex/freesurfer.py +++ b/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, diff --git a/examples/fsaverage/README.txt b/examples/fsaverage/README.txt new file mode 100644 index 00000000..6e31f972 --- /dev/null +++ b/examples/fsaverage/README.txt @@ -0,0 +1,4 @@ +Examples with fsaverage +------------------------------ + +Examples showing how to use PyCortex to visualize data on fsaverage. \ No newline at end of file diff --git a/examples/fsaverage/upsample_to_fsaverage.py b/examples/fsaverage/upsample_to_fsaverage.py new file mode 100644 index 00000000..36593725 --- /dev/null +++ b/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()