diff --git a/cortex/__init__.py b/cortex/__init__.py index 0c0bc11c..28c3dc4e 100644 --- a/cortex/__init__.py +++ b/cortex/__init__.py @@ -1,6 +1,15 @@ # emacs: -*- coding: utf-8; mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set fileencoding=utf-8 ft=python sts=4 ts=4 sw=4 et: -from cortex.dataset import Dataset, Volume, Vertex, VolumeRGB, VertexRGB, Volume2D, Vertex2D, Colors +from cortex.dataset import ( + Dataset, + Volume, + Vertex, + VolumeRGB, + VertexRGB, + Volume2D, + Vertex2D, + Colors, +) from cortex import align, volume, quickflat, webgl, segment, options from cortex.database import db from cortex.utils import * @@ -12,9 +21,11 @@ try: from cortex import formats except ImportError: - raise ImportError("Either are running pycortex from the source directory, or the build is broken. " - "If your current working directory is 'cortex', where pycortex is installed, then change this. " - "If your current working directory is somewhere else, then you may have to rebuild pycortex.") + raise ImportError( + "Either are running pycortex from the source directory, or the build is broken. " + "If your current working directory is 'cortex', where pycortex is installed, then change this. " + "If your current working directory is somewhere else, then you may have to rebuild pycortex." + ) load = Dataset.from_file @@ -29,20 +40,29 @@ except ImportError: pass +try: + from cortex import hcp +except ImportError: + pass + + # Create deprecated interface for database class dep(object): def __getattr__(self, name): warnings.warn("cortex.surfs is deprecated, use cortex.db instead", Warning) return getattr(db, name) + def __dir__(self): warnings.warn("cortex.surfs is deprecated, use cortex.db instead", Warning) return db.__dir__() + surfs = dep() import sys + if sys.version_info < (3,): stdout = sys.stdout reload(sys) - sys.setdefaultencoding('utf8') + sys.setdefaultencoding("utf8") sys.stdout = stdout diff --git a/cortex/data/upsample_fsaverage_neighbors.npz b/cortex/data/upsample_fsaverage_neighbors.npz new file mode 100644 index 00000000..664cd196 Binary files /dev/null and b/cortex/data/upsample_fsaverage_neighbors.npz differ diff --git a/cortex/freesurfer.py b/cortex/freesurfer.py index a76a1700..5fb3441c 100644 --- a/cortex/freesurfer.py +++ b/cortex/freesurfer.py @@ -1,14 +1,12 @@ -"""Contains functions for interfacing with freesurfer -""" +"""Contains functions for interfacing with freesurfer""" + from __future__ import print_function import copy import os import shlex -import shutil import struct import subprocess as sp -import tempfile import warnings from builtins import input from tempfile import NamedTemporaryFile @@ -19,7 +17,7 @@ from scipy.sparse import coo_matrix from scipy.spatial import KDTree -from . import anat, database +from . import appdirs, database def get_paths(fs_subject, hemi, type="patch", freesurfer_subject_dir=None): @@ -39,16 +37,16 @@ def get_paths(fs_subject, hemi, type="patch", freesurfer_subject_dir=None): by freesurfer) """ if freesurfer_subject_dir is None: - freesurfer_subject_dir = os.environ['SUBJECTS_DIR'] + freesurfer_subject_dir = os.environ["SUBJECTS_DIR"] base = os.path.join(freesurfer_subject_dir, fs_subject) if type == "patch": - return os.path.join(base, "surf", hemi+".{name}.patch.3d") + return os.path.join(base, "surf", hemi + ".{name}.patch.3d") elif type == "surf": - return os.path.join(base, "surf", hemi+".{name}") + return os.path.join(base, "surf", hemi + ".{name}") elif type == "curv": - return os.path.join(base, "surf", hemi+".curv{name}") + return os.path.join(base, "surf", hemi + ".curv{name}") elif type == "slim": - return os.path.join(base, "surf", hemi+".{name}_slim.obj") + return os.path.join(base, "surf", hemi + ".{name}_slim.obj") def autorecon(fs_subject, type="all", parallel=True, n_cores=None): @@ -63,32 +61,31 @@ def autorecon(fs_subject, type="all", parallel=True, n_cores=None): """ types = { - 'all': 'autorecon-all', - '1': "autorecon1", - '2': "autorecon2", - '3': "autorecon3", - 'cp': "autorecon2-cp", - 'wm': "autorecon2-wm", - 'pia': "autorecon-pial"} - - times = { - 'all': "12 hours", - '2': "6 hours", - 'cp': "8 hours", - 'wm': "4 hours" - } + "all": "autorecon-all", + "1": "autorecon1", + "2": "autorecon2", + "3": "autorecon3", + "cp": "autorecon2-cp", + "wm": "autorecon2-wm", + "pia": "autorecon-pial", + } + + times = {"all": "12 hours", "2": "6 hours", "cp": "8 hours", "wm": "4 hours"} if str(type) in times: - resp = input("recon-all will take approximately %s to run! Continue? "%times[str(type)]) + resp = input( + "recon-all will take approximately %s to run! Continue? " % times[str(type)] + ) if resp.lower() not in ("yes", "y"): return cmd = "recon-all -s {subj} -{cmd}".format(subj=fs_subject, cmd=types[str(type)]) - if parallel and type in ('2', 'wm', 'all'): + if parallel and type in ("2", "wm", "all"): # Parallelization only works for autorecon2 or autorecon2-wm if n_cores is None: import multiprocessing as mp + n_cores = mp.cpu_count() - cmd += ' -parallel -openmp {n_cores:d}'.format(n_cores=n_cores) + cmd += " -parallel -openmp {n_cores:d}".format(n_cores=n_cores) print("Calling:\n{cmd}".format(cmd=cmd)) sp.check_call(shlex.split(cmd)) @@ -123,15 +120,21 @@ def flatten(fs_subject, hemi, patch, freesurfer_subject_dir=None, save_every=Non shows that a flattening iteration has completed) https://stackoverflow.com/questions/4417546/constantly-print-subprocess-output-while-process-is-running """ - resp = input('Flattening takes approximately 2 hours! Continue? ') - if resp.lower() in ('y', 'yes'): - inpath = get_paths(fs_subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir).format(name=patch) - outpath = get_paths(fs_subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir).format(name=patch+".flat") + resp = input("Flattening takes approximately 2 hours! Continue? ") + if resp.lower() in ("y", "yes"): + inpath = get_paths( + fs_subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir + ).format(name=patch) + outpath = get_paths( + fs_subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir + ).format(name=patch + ".flat") if save_every is None: - save_every_str = '' + save_every_str = "" else: - save_every_str = ' -w %d'%save_every - cmd = "mris_flatten -O fiducial{save_every_str} {inpath} {outpath}".format(inpath=inpath, outpath=outpath, save_every_str=save_every_str) + save_every_str = " -w %d" % save_every + cmd = "mris_flatten -O fiducial{save_every_str} {inpath} {outpath}".format( + inpath=inpath, outpath=outpath, save_every_str=save_every_str + ) print("Calling: ") print(cmd) sp.check_call(shlex.split(cmd)) @@ -149,28 +152,28 @@ def import_subj( ): """Imports a subject from freesurfer - This will overwrite (after giving a warning and an option to continue) the + This will overwrite (after giving a warning and an option to continue) the pre-existing subject, including all blender cuts, masks, transforms, etc., and - re-generate surface info files (curvature, sulcal depth, thickness) stored in - the surfinfo/ folder for the subject. All cached files for the subject will be - deleted. + re-generate surface info files (curvature, sulcal depth, thickness) stored in + the surfinfo/ folder for the subject. All cached files for the subject will be + deleted. Parameters ---------- freesurfer_subject : str Freesurfer subject name pycortex_subject : str, optional - Pycortex subject name. By default it uses the freesurfer subject name. - It is advised to stick to that convention, if possible + Pycortex subject name. By default it uses the freesurfer subject name. + It is advised to stick to that convention, if possible (your life will go more smoothly.) freesurfer_subject_dir : str, optional Freesurfer subject directory to pull data from. By default uses the directory given by the environment variable $SUBJECTS_DIR. whitematter_surf : str, optional - Which whitematter surface to import as 'wm'. By default uses 'smoothwm', but - that surface is smoothed and may not be appropriate. + Which whitematter surface to import as 'wm'. By default uses 'smoothwm', but + that surface is smoothed and may not be appropriate. A good alternative is 'white'. - + Notes ----- This function uses command line functions from freesurfer, so you should make sure @@ -178,7 +181,7 @@ def import_subj( This function will also generate the fiducial surfaces for the subject, which are halfway between the white matter and pial surfaces. The surfaces will be stored - in the freesurfer subject's directory. These fiducial surfaces are used for + in the freesurfer subject's directory. These fiducial surfaces are used for cutting and flattening. """ # Check if freesurfer is sourced or if subjects dir is passed @@ -235,14 +238,14 @@ def import_subj( # system as the volume data, rather than the TKR coordinate system, which # has the center set to FOV/2. # NOTE: the resulting gifti surfaces will look misaligned with respect to - # the anatomical volumes when visualized in freeview, because freeview + # the anatomical volumes when visualized in freeview, because freeview # expects the surfaces to be in TKR coordinates (with center set to FOV/2). # But the surfaces stored in the pycortex database are only to be used by # pycortex, so that's fine. cmd = f"mris_convert --to-scanner {in_surface} {out_surface}" sp.check_output(shlex.split(cmd)) - # Dictionary mapping for curvature and extra info to be imported + # Dictionary mapping for curvature and extra info to be imported info_fs2pycortex = { "sulc": "sulcaldepth", "thickness": "thickness", @@ -254,23 +257,26 @@ def import_subj( fs_surf_template.format(hemi=hemi, name=fsname) for hemi in ["lh", "rh"] ] lh, rh = [parse_curv(in_info) for in_info in in_info_lhrh] - np.savez( - surfinfo_template.format(name=name), - left=-lh, - right=-rh - ) + np.savez(surfinfo_template.format(name=name), left=-lh, right=-rh) # Finally update the database by re-initializing it database.db = database.Database() -def import_flat(fs_subject, patch, hemis=['lh', 'rh'], cx_subject=None, - flat_type='freesurfer', auto_overwrite=False, - freesurfer_subject_dir=None, clean=True): +def import_flat( + fs_subject, + patch, + hemis=["lh", "rh"], + cx_subject=None, + flat_type="freesurfer", + auto_overwrite=False, + freesurfer_subject_dir=None, + clean=True, +): """Imports a flat brain from freesurfer NOTE: This will delete the overlays.svg file for this subject, since THE - FLATMAPS WILL CHANGE, as well as all cached information (e.g. old flatmap - boundaries, roi svg intermediate renders, etc). + FLATMAPS WILL CHANGE, as well as all cached information (e.g. old flatmap + boundaries, roi svg intermediate renders, etc). Parameters ---------- @@ -295,34 +301,56 @@ def import_flat(fs_subject, patch, hemis=['lh', 'rh'], cx_subject=None, ------- """ if not auto_overwrite: - proceed = input(('Warning: This is intended to over-write .gii files storing\n' - 'flatmap vertex locations for this subject, and will result\n' - 'in deletion of the overlays.svg file and all cached info\n' - 'for this subject (because flatmaps will fundamentally change).\n' - 'Proceed? [y]/n: ')) - if proceed.lower() not in ['y', 'yes', '']: + proceed = input( + ( + "Warning: This is intended to over-write .gii files storing\n" + "flatmap vertex locations for this subject, and will result\n" + "in deletion of the overlays.svg file and all cached info\n" + "for this subject (because flatmaps will fundamentally change).\n" + "Proceed? [y]/n: " + ) + ) + if proceed.lower() not in ["y", "yes", ""]: print(">>> Elected to quit rather than delete & overwrite files.") return if cx_subject is None: cx_subject = fs_subject - surfs = os.path.join(database.default_filestore, cx_subject, "surfaces", "flat_{hemi}.gii") + surfs = os.path.join( + database.default_filestore, cx_subject, "surfaces", "flat_{hemi}.gii" + ) from . import formats + for hemi in hemis: - if flat_type in ['freesurfer', 'blender']: - surf_path = (patch + ".flat") if (flat_type == 'freesurfer') else (patch + ".flat.blender") - pts, polys, _ = get_surf(fs_subject, hemi, "patch", surf_path, freesurfer_subject_dir=freesurfer_subject_dir) + if flat_type in ["freesurfer", "blender"]: + surf_path = ( + (patch + ".flat") + if (flat_type == "freesurfer") + else (patch + ".flat.blender") + ) + pts, polys, _ = get_surf( + fs_subject, + hemi, + "patch", + surf_path, + freesurfer_subject_dir=freesurfer_subject_dir, + ) - if flat_type == 'freesurfer': + if flat_type == "freesurfer": # Reorder axes: X, Y, Z instead of Y, X, Z flat = pts[:, [1, 0, 2]] # Flip Y axis upside down flat[:, 1] = -flat[:, 1] else: flat = pts - elif flat_type == 'slim': - flat_file = get_paths(fs_subject, hemi, type='slim', freesurfer_subject_dir=freesurfer_subject_dir) + elif flat_type == "slim": + flat_file = get_paths( + fs_subject, + hemi, + type="slim", + freesurfer_subject_dir=freesurfer_subject_dir, + ) flat_file = flat_file.format(name=patch + ".flat") flat, polys = formats.read_obj(flat_file) @@ -331,25 +359,25 @@ def import_flat(fs_subject, patch, hemis=['lh', 'rh'], cx_subject=None, flat = _move_disconnect_points_to_zero(flat, polys) fname = surfs.format(hemi=hemi) - print("saving to %s"%fname) + print("saving to %s" % fname) formats.write_gii(fname, pts=flat, polys=polys) # clear the cache, per #81 database.db.clear_cache(cx_subject) # Remove overlays.svg file (FLATMAPS HAVE CHANGED) - overlays_file = database.db.get_paths(cx_subject)['overlays'] + overlays_file = database.db.get_paths(cx_subject)["overlays"] if os.path.exists(overlays_file): os.unlink(overlays_file) - # Regenerate it? + # Regenerate it? def _remove_disconnected_polys(polys): """Remove polygons that are not in the main connected component. - + This function creates a sparse graph based on edges in the input. Then it computes the connected components, and returns only the polygons that are in the largest component. - + This filtering is useful to remove disconnected vertices resulting from a poor surface cut. """ @@ -357,18 +385,15 @@ def _remove_disconnected_polys(polys): import scipy.sparse as sp # create the sparse graph - row = np.concatenate([ - polys[:, 0], polys[:, 1], polys[:, 0], - polys[:, 2], polys[:, 1], polys[:, 2] - ]) - col = np.concatenate([ - polys[:, 1], polys[:, 0], polys[:, 2], - polys[:, 0], polys[:, 2], polys[:, 1] - ]) + row = np.concatenate( + [polys[:, 0], polys[:, 1], polys[:, 0], polys[:, 2], polys[:, 1], polys[:, 2]] + ) + col = np.concatenate( + [polys[:, 1], polys[:, 0], polys[:, 2], polys[:, 0], polys[:, 2], polys[:, 1]] + ) data = np.ones(len(col), dtype=bool) - graph = sp.coo_matrix((data, (row, col)), shape=(n_points, n_points), - dtype=bool) - + graph = sp.coo_matrix((data, (row, col)), shape=(n_points, n_points), dtype=bool) + # compute connected components n_components, labels = sp.csgraph.connected_components(graph) unique_labels, counts = np.unique(labels, return_counts=True) @@ -384,7 +409,7 @@ def _remove_disconnected_polys(polys): def _move_disconnect_points_to_zero(pts, polys): """Change coordinates of points not in polygons to zero. - + This cleaning step is useful after _remove_disconnected_polys, to avoid using this points in boundaries computations (through pts.max(axis=0) here and there). @@ -396,41 +421,44 @@ def _move_disconnect_points_to_zero(pts, polys): def make_fiducial(fs_subject, freesurfer_subject_dir=None): - """Make fiducial surface (halfway between white matter and pial surfaces) - """ - for hemi in ['lh', 'rh']: - spts, polys, _ = get_surf(fs_subject, hemi, "smoothwm", freesurfer_subject_dir=freesurfer_subject_dir) - ppts, _, _ = get_surf(fs_subject, hemi, "pial", freesurfer_subject_dir=freesurfer_subject_dir) - fname = get_paths(fs_subject, hemi, "surf", freesurfer_subject_dir=freesurfer_subject_dir).format(name="fiducial") + """Make fiducial surface (halfway between white matter and pial surfaces)""" + for hemi in ["lh", "rh"]: + spts, polys, _ = get_surf( + fs_subject, hemi, "smoothwm", freesurfer_subject_dir=freesurfer_subject_dir + ) + ppts, _, _ = get_surf( + fs_subject, hemi, "pial", freesurfer_subject_dir=freesurfer_subject_dir + ) + fname = get_paths( + fs_subject, hemi, "surf", freesurfer_subject_dir=freesurfer_subject_dir + ).format(name="fiducial") write_surf(fname, (spts + ppts) / 2, polys) def parse_surf(filename): - """ - """ - with open(filename, 'rb') as fp: - #skip magic + """ """ + with open(filename, "rb") as fp: + # skip magic fp.seek(3) comment = fp.readline() fp.readline() print(comment) - verts, faces = struct.unpack('>2I', fp.read(8)) - pts = np.frombuffer(fp.read(4*3*verts), dtype='f4').byteswap() - polys = np.frombuffer(fp.read(4*3*faces), dtype='i4').byteswap() + verts, faces = struct.unpack(">2I", fp.read(8)) + pts = np.frombuffer(fp.read(4 * 3 * verts), dtype="f4").byteswap() + polys = np.frombuffer(fp.read(4 * 3 * faces), dtype="i4").byteswap() return pts.reshape(-1, 3), polys.reshape(-1, 3) -def write_surf(filename, pts, polys, comment=''): - """Write freesurfer surface file - """ - with open(filename, 'wb') as fp: - fp.write(b'\xff\xff\xfe') - fp.write((comment+'\n\n').encode()) - fp.write(struct.pack('>2I', len(pts), len(polys))) +def write_surf(filename, pts, polys, comment=""): + """Write freesurfer surface file""" + with open(filename, "wb") as fp: + fp.write(b"\xff\xff\xfe") + fp.write((comment + "\n\n").encode()) + fp.write(struct.pack(">2I", len(pts), len(polys))) fp.write(pts.astype(np.float32).byteswap().tobytes()) fp.write(polys.astype(np.uint32).byteswap().tobytes()) - fp.write(b'\n') + fp.write(b"\n") def write_patch(filename, pts, edges=None): @@ -455,53 +483,63 @@ def write_patch(filename, pts, edges=None): if edges is None: edges = set() - with open(filename, 'wb') as fp: - fp.write(struct.pack('>2i', -1, len(pts))) + with open(filename, "wb") as fp: + fp.write(struct.pack(">2i", -1, len(pts))) for i, pt in pts: if i in edges: - fp.write(struct.pack('>i3f', -i-1, *pt)) + fp.write(struct.pack(">i3f", -i - 1, *pt)) else: - fp.write(struct.pack('>i3f', i+1, *pt)) + fp.write(struct.pack(">i3f", i + 1, *pt)) def parse_curv(filename): - """ - """ - with open(filename, 'rb') as fp: + """ """ + with open(filename, "rb") as fp: fp.seek(15) - return np.frombuffer(fp.read(), dtype='>f4').byteswap().view(np.dtype('>f4').newbyteorder('=')) + return ( + np.frombuffer(fp.read(), dtype=">f4") + .byteswap() + .view(np.dtype(">f4").newbyteorder("=")) + ) def parse_patch(filename): - """ - """ - with open(filename, 'rb') as fp: - header, = struct.unpack('>i', fp.read(4)) - nverts, = struct.unpack('>i', fp.read(4)) - data = np.frombuffer(fp.read(), dtype=[('vert', '>i4'), ('x', '>f4'), - ('y', '>f4'), ('z', '>f4')]) + """ """ + with open(filename, "rb") as fp: + (header,) = struct.unpack(">i", fp.read(4)) + (nverts,) = struct.unpack(">i", fp.read(4)) + data = np.frombuffer( + fp.read(), dtype=[("vert", ">i4"), ("x", ">f4"), ("y", ">f4"), ("z", ">f4")] + ) assert len(data) == nverts return data -def get_surf(subject, hemi, type, patch=None, flatten_step=None, freesurfer_subject_dir=None): - """Read freesurfer surface file - """ +def get_surf( + subject, hemi, type, patch=None, flatten_step=None, freesurfer_subject_dir=None +): + """Read freesurfer surface file""" if type == "patch": assert patch is not None - surf_file = get_paths(subject, hemi, 'surf', freesurfer_subject_dir=freesurfer_subject_dir).format(name='smoothwm') + surf_file = get_paths( + subject, hemi, "surf", freesurfer_subject_dir=freesurfer_subject_dir + ).format(name="smoothwm") else: - surf_file = get_paths(subject, hemi, 'surf', freesurfer_subject_dir=freesurfer_subject_dir).format(name=type) + surf_file = get_paths( + subject, hemi, "surf", freesurfer_subject_dir=freesurfer_subject_dir + ).format(name=type) pts, polys = parse_surf(surf_file) if patch is not None: - patch_file = get_paths(subject, hemi, 'patch', freesurfer_subject_dir=freesurfer_subject_dir).format(name=patch) + patch_file = get_paths( + subject, hemi, "patch", freesurfer_subject_dir=freesurfer_subject_dir + ).format(name=patch) if flatten_step is not None: - patch_file += '%04d'%flatten_step + patch_file += "%04d" % flatten_step patch = parse_patch(patch_file) - verts = patch[patch['vert'] > 0]['vert'] - 1 - edges = -patch[patch['vert'] < 0]['vert'] - 1 + verts = patch[patch["vert"] > 0]["vert"] - 1 + edges = -patch[patch["vert"] < 0]["vert"] - 1 idx = np.zeros((len(pts),), dtype=bool) idx[verts] = True @@ -513,42 +551,64 @@ def get_surf(subject, hemi, type, patch=None, flatten_step=None, freesurfer_subj idx[edges] = -1 if type == "patch": - for i, x in enumerate(['x', 'y', 'z']): - pts[verts, i] = patch[patch['vert'] > 0][x] - pts[edges, i] = patch[patch['vert'] < 0][x] + for i, x in enumerate(["x", "y", "z"]): + pts[verts, i] = patch[patch["vert"] > 0][x] + pts[edges, i] = patch[patch["vert"] < 0][x] return pts, polys, idx - return pts, polys, get_curv(subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir) + return ( + pts, + polys, + get_curv(subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir), + ) -def _move_labels(subject, label, hemisphere=('lh','rh'), fs_dir=None, src_subject='fsaverage'): +def _move_labels( + subject, label, hemisphere=("lh", "rh"), fs_dir=None, src_subject="fsaverage" +): """subject is a freesurfer subject""" if fs_dir is None: - fs_dir = os.environ['SUBJECTS_DIR'] + fs_dir = os.environ["SUBJECTS_DIR"] for hemi in hemisphere: - srclabel = os.path.join(fs_dir, src_subject, 'label', - '{hemi}.{label}.label'.format(hemi=hemi, label=label)) - trglabel = os.path.join(fs_dir, subject, 'label', - '{hemi}.{label}.label'.format(hemi=hemi, label=label)) + srclabel = os.path.join( + fs_dir, + src_subject, + "label", + "{hemi}.{label}.label".format(hemi=hemi, label=label), + ) + trglabel = os.path.join( + fs_dir, + subject, + "label", + "{hemi}.{label}.label".format(hemi=hemi, label=label), + ) if not os.path.exists(srclabel): raise ValueError("Label {} doesn't exist!".format(srclabel)) - fs_sub_dir = os.path.join(fs_dir, subject, 'label') + fs_sub_dir = os.path.join(fs_dir, subject, "label") if not os.path.exists(fs_sub_dir): - raise ValueError("Freesurfer subject directory for subject ({}) does not exist!".format(fs_sub_dir)) - cmd = ("mri_label2label --srcsubject {src_subject} --trgsubject {subject} " - "--srclabel {srclabel} --trglabel {trglabel} " - "--regmethod surface --hemi {hemi}") - cmd_f = cmd.format(hemi=hemi, subject=subject, src_subject=src_subject, - srclabel=srclabel, trglabel=trglabel) + raise ValueError( + "Freesurfer subject directory for subject ({}) does not exist!".format( + fs_sub_dir + ) + ) + cmd = ( + "mri_label2label --srcsubject {src_subject} --trgsubject {subject} " + "--srclabel {srclabel} --trglabel {trglabel} " + "--regmethod surface --hemi {hemi}" + ) + cmd_f = cmd.format( + hemi=hemi, + subject=subject, + src_subject=src_subject, + srclabel=srclabel, + trglabel=trglabel, + ) print("Calling: ") print(cmd_f) to_call = shlex.split(cmd_f) - proc = sp.Popen(to_call, - stdin=sp.PIPE, - stdout=sp.PIPE, - stderr=sp.PIPE) + proc = sp.Popen(to_call, stdin=sp.PIPE, stdout=sp.PIPE, stderr=sp.PIPE) stdout, stderr = proc.communicate() - if stderr not in ('', b''): + if stderr not in ("", b""): raise Exception("Error in freesurfer function call:\n{}".format(stderr)) print("Labels transferred") @@ -567,22 +627,34 @@ def _parse_labels(label_files, cx_subject): label_files = [label_files] verts = [] values = [] - lh_surf, _ = database.db.get_surf(cx_subject, 'fiducial', 'left') + lh_surf, _ = database.db.get_surf(cx_subject, "fiducial", "left") for fname in label_files: with open(fname) as fid: lines = fid.readlines() - lines = [[float(xx.strip()) for xx in x.split(' ') if xx.strip()] for x in lines[2:]] + lines = [ + [float(xx.strip()) for xx in x.split(" ") if xx.strip()] + for x in lines[2:] + ] vals = np.array(lines) - if '/lh.' in fname: - verts.append(vals[:,0]) - elif '/rh.' in fname: - verts.append(vals[:,0] + lh_surf.shape[0]) - values.append(vals[:,-1]) + if "/lh." in fname: + verts.append(vals[:, 0]) + elif "/rh." in fname: + verts.append(vals[:, 0] + lh_surf.shape[0]) + values.append(vals[:, -1]) verts = np.hstack(verts) values = np.hstack(values) return verts, values -def get_label(cx_subject, label, fs_subject=None, fs_dir=None, src_subject='fsaverage', hemisphere=('lh', 'rh'), **kwargs): + +def get_label( + cx_subject, + label, + fs_subject=None, + fs_dir=None, + src_subject="fsaverage", + hemisphere=("lh", "rh"), + **kwargs, +): """Get data from a label file for fsaverage subject Parameters @@ -601,18 +673,27 @@ def get_label(cx_subject, label, fs_subject=None, fs_dir=None, src_subject='fsav """ if fs_dir is None: - fs_dir = os.environ['SUBJECTS_DIR'] + fs_dir = os.environ["SUBJECTS_DIR"] else: - os.environ['SUBJECTS_DIR'] = fs_dir + os.environ["SUBJECTS_DIR"] = fs_dir if fs_subject is None: fs_subject = cx_subject - label_files = [os.path.join(fs_dir, fs_subject, 'label', '{}.{}.label'.format(h, label)) for h in hemisphere] - if cx_subject not in ['fsaverage', 'MNI']: + label_files = [ + os.path.join(fs_dir, fs_subject, "label", "{}.{}.label".format(h, label)) + for h in hemisphere + ] + if cx_subject not in ["fsaverage", "MNI"]: # If label file doesn't exist, try to move it there - print('looking for {}'.format(label_files)) + print("looking for {}".format(label_files)) if not all([os.path.exists(f) for f in label_files]): print("Transforming label file to subject's freesurfer directory...") - _move_labels(fs_subject, label, hemisphere=hemisphere, fs_dir=fs_dir, src_subject=src_subject) + _move_labels( + fs_subject, + label, + hemisphere=hemisphere, + fs_dir=fs_dir, + src_subject=src_subject, + ) verts, values = _parse_labels(label_files, cx_subject) idx = verts.astype(int) return idx, values @@ -623,12 +704,19 @@ def _mri_surf2surf_command(src_subj, trg_subj, input_file, output_file, hemi): # --trgsubject --trgsurfval --hemi - cmd = ["mri_surf2surf", "--srcsubject", src_subj, - "--sval", input_file, - "--trgsubject", trg_subj, - "--tval", output_file, - "--hemi", hemi, - ] + cmd = [ + "mri_surf2surf", + "--srcsubject", + src_subj, + "--sval", + input_file, + "--trgsubject", + trg_subj, + "--tval", + output_file, + "--hemi", + hemi, + ] return cmd @@ -645,21 +733,21 @@ def _check_datatype(data): def mri_surf2surf(data, source_subj, target_subj, hemi, subjects_dir=None): """Uses freesurfer mri_surf2surf to transfer vertex data between two freesurfer subjects - + Parameters ========== data: ndarray, shape=(n_imgs, n_verts) data arrays representing vertex data - + source_subj: str freesurfer subject name of source subject - + target_subj: str freesurfer subject name of target subject - + hemi: str in ("lh", "rh") string indicating hemisphere. - + Notes ===== Requires path to mri_surf2surf or freesurfer environment to be active. @@ -671,31 +759,38 @@ def mri_surf2surf(data, source_subj, target_subj, hemi, subjects_dir=None): tf_in = NamedTemporaryFile(suffix=".gii") nibabel.save(gifti_image, tf_in.name) - tf_out = NamedTemporaryFile(suffix='.gii') - cmd = _mri_surf2surf_command(source_subj, target_subj, - tf_in.name, tf_out.name, hemi) + tf_out = NamedTemporaryFile(suffix=".gii") + cmd = _mri_surf2surf_command( + source_subj, target_subj, tf_in.name, tf_out.name, hemi + ) if subjects_dir is not None: env = os.environ.copy() - env['SUBJECTS_DIR'] = subjects_dir + env["SUBJECTS_DIR"] = subjects_dir else: env = None - print('Calling:') - print(' '.join(cmd)) + print("Calling:") + print(" ".join(cmd)) p = sp.Popen(cmd, env=env) exit_code = p.wait() if exit_code != 0: if exit_code == 255: - raise Exception(("Missing file (see above). " - "If lh.sphere.reg is missing,\n" - "you likely need to run the 3rd " - "stage of freesurfer autorecon\n" - "(sphere registration) for this subject:\n" - ">>> cortex.freesurfer.autorecon('{fs_subject}', type='3')" - ).format(fs_subject=source_subj)) - #from subprocess import CalledProcessError # handle with this, maybe? - raise Exception(("Exit code {exit_code} means that " - "mri_surf2surf failed").format(exit_code=exit_code)) + raise Exception( + ( + "Missing file (see above). " + "If lh.sphere.reg is missing,\n" + "you likely need to run the 3rd " + "stage of freesurfer autorecon\n" + "(sphere registration) for this subject:\n" + ">>> cortex.freesurfer.autorecon('{fs_subject}', type='3')" + ).format(fs_subject=source_subj) + ) + # from subprocess import CalledProcessError # handle with this, maybe? + raise Exception( + ("Exit code {exit_code} means that mri_surf2surf failed").format( + exit_code=exit_code + ) + ) tf_in.close() output_img = nibabel.load(tf_out.name) @@ -710,9 +805,9 @@ def _read_sphere_reg(subject, hemi, subjects_dir=None): These are the coordinates on which freesurfer's spherical registration defines the cross-subject vertex correspondence used by ``mri_surf2surf``. """ - surf_file = get_paths(subject, hemi, 'surf', - freesurfer_subject_dir=subjects_dir).format( - name='sphere.reg') + surf_file = get_paths( + subject, hemi, "surf", freesurfer_subject_dir=subjects_dir + ).format(name="sphere.reg") pts, _ = parse_surf(surf_file) return pts @@ -782,8 +877,9 @@ def _surf2surf_nnfr_matrix(src_sphere, trg_sphere): rows = np.concatenate([np.arange(n_trg), rev_trg[orphan]]) cols = np.concatenate([fwd_src, orphan]) - matrix = coo_matrix((np.ones(len(rows)), (rows, cols)), - shape=(n_trg, n_src)).tocsr() + matrix = coo_matrix( + (np.ones(len(rows)), (rows, cols)), shape=(n_trg, n_src) + ).tocsr() # A target/source pair can be added by both passes; collapse duplicates. matrix.sum_duplicates() @@ -800,13 +896,18 @@ def _surf2surf_nnfr_matrix(src_sphere, trg_sphere): # implementation of ``get_mri_surf2surf_matrix``; kept only so existing # callers do not break. _SURF2SURF_LEGACY_KWARGS = frozenset( - ('n_neighbors', 'random_state', 'n_test_images', 'coef_threshold', - 'renormalize')) + ("n_neighbors", "random_state", "n_test_images", "coef_threshold", "renormalize") +) -def get_mri_surf2surf_matrix(source_subj, hemi, surface_type=None, - target_subj='fsaverage', subjects_dir=None, - **kwargs): +def get_mri_surf2surf_matrix( + source_subj, + hemi, + surface_type=None, + target_subj="fsaverage", + subjects_dir=None, + **kwargs, +): """Create a sparse matrix implementing freesurfer's ``mri_surf2surf``. A surface-to-surface resampling is a linear, highly localized transform @@ -858,19 +959,23 @@ def get_mri_surf2surf_matrix(source_subj, hemi, surface_type=None, "get_mri_surf2surf_matrix no longer uses the regression-based " "parameters {}; they are ignored. The matrix is now built " "directly from the spherical registration.".format(sorted(legacy)), - DeprecationWarning, stacklevel=2) + DeprecationWarning, + stacklevel=2, + ) unexpected = set(kwargs) - _SURF2SURF_LEGACY_KWARGS if unexpected: raise TypeError( - "get_mri_surf2surf_matrix got unexpected keyword argument(s) " - "{}".format(sorted(unexpected))) + "get_mri_surf2surf_matrix got unexpected keyword argument(s) {}".format( + sorted(unexpected) + ) + ) src_sphere = _read_sphere_reg(source_subj, hemi, subjects_dir) trg_sphere = _read_sphere_reg(target_subj, hemi, subjects_dir) return _surf2surf_nnfr_matrix(src_sphere, trg_sphere) -def get_curv(fs_subject, hemi, type='wm', freesurfer_subject_dir=None): +def get_curv(fs_subject, hemi, type="wm", freesurfer_subject_dir=None): """Load freesurfer curv file for a freesurfer subject Parameters @@ -882,13 +987,17 @@ def get_curv(fs_subject, hemi, type='wm', freesurfer_subject_dir=None): type : str 'wm' or other type of surface (e.g. 'fiducial' or 'pial') freesurfer_subject_dir : str - directory for Freesurfer subjects (defaults to value for the + directory for Freesurfer subjects (defaults to value for the environment variable $SUBJECTS_DIR if None) """ if type == "wm": - curv_file = get_paths(fs_subject, hemi, 'curv', freesurfer_subject_dir=freesurfer_subject_dir).format(name='') + curv_file = get_paths( + fs_subject, hemi, "curv", freesurfer_subject_dir=freesurfer_subject_dir + ).format(name="") else: - curv_file = get_paths(fs_subject, hemi, 'curv', freesurfer_subject_dir=freesurfer_subject_dir).format(name='.'+type) + curv_file = get_paths( + fs_subject, hemi, "curv", freesurfer_subject_dir=freesurfer_subject_dir + ).format(name="." + type) return parse_curv(curv_file) @@ -910,74 +1019,91 @@ def show_surf(subject, hemi, type, patch=None, curv=True, freesurfer_subject_dir freesurfer_subject_dir : """ - warnings.warn(('Deprecated and probably broken! Try `cortex.segment.show_surf()`\n' - 'which uses a third-party program (meshlab, available for linux & mac\n' - 'instead of buggy mayavi code!')) + warnings.warn( + ( + "Deprecated and probably broken! Try `cortex.segment.show_surf()`\n" + "which uses a third-party program (meshlab, available for linux & mac\n" + "instead of buggy mayavi code!" + ) + ) from mayavi import mlab from tvtk.api import tvtk - pts, polys, idx = get_surf(subject, hemi, type, patch, freesurfer_subject_dir=freesurfer_subject_dir) + pts, polys, idx = get_surf( + subject, hemi, type, patch, freesurfer_subject_dir=freesurfer_subject_dir + ) if curv: curv = get_curv(subject, hemi, freesurfer_subject_dir=freesurfer_subject_dir) else: curv = idx fig = mlab.figure() - src = mlab.pipeline.triangular_mesh_source(pts[:,0], pts[:,1], pts[:,2], polys, scalars=curv, figure=fig) + src = mlab.pipeline.triangular_mesh_source( + pts[:, 0], pts[:, 1], pts[:, 2], polys, scalars=curv, figure=fig + ) norms = mlab.pipeline.poly_data_normals(src, figure=fig) norms.filter.splitting = False surf = mlab.pipeline.surface(norms, figure=fig) - surf.parent.scalar_lut_manager.set(lut_mode='RdBu', data_range=[-1,1], use_default_range=False) + surf.parent.scalar_lut_manager.set( + lut_mode="RdBu", data_range=[-1, 1], use_default_range=False + ) cursors = mlab.pipeline.scalar_scatter([0], [0], [0]) glyphs = mlab.pipeline.glyph(cursors, figure=fig) - glyphs.glyph.glyph_source.glyph_source = glyphs.glyph.glyph_source.glyph_dict['axes'] + glyphs.glyph.glyph_source.glyph_source = glyphs.glyph.glyph_source.glyph_dict[ + "axes" + ] - fig.scene.background = (0,0,0) + fig.scene.background = (0, 0, 0) fig.scene.interactor.interactor_style = tvtk.InteractorStyleTerrain() - path = os.path.join(os.environ['SUBJECTS_DIR'], subject) + path = os.path.join(os.environ["SUBJECTS_DIR"], subject) + def picker_callback(picker): if picker.actor in surf.actor.actors: - npts = np.append(cursors.data.points.to_array(), [pts[picker.point_id]], axis=0) + npts = np.append( + cursors.data.points.to_array(), [pts[picker.point_id]], axis=0 + ) cursors.data.points = npts print(picker.point_id) x, y, z = pts[picker.point_id] - with open(os.path.join(path, 'tmp', 'edit.dat'), 'w') as fp: - fp.write('%f %f %f\n'%(x, y, z)) + with open(os.path.join(path, "tmp", "edit.dat"), "w") as fp: + fp.write("%f %f %f\n" % (x, y, z)) picker = fig.on_mouse_pick(picker_callback) picker.tolerance = 0.01 mlab.show() return fig, surf + def write_dot(fname, pts, polys, name="test"): - """ - """ + """ """ import networkx as nx + def iter_surfedges(tris): - for a,b,c in tris: - yield a,b - yield b,c - yield a,c + for a, b, c in tris: + yield a, b + yield b, c + yield a, c + graph = nx.Graph() graph.add_edges_from(iter_surfedges(polys)) lengths = [] with open(fname, "w") as fp: - fp.write("graph %s {\n"%name) + fp.write("graph %s {\n" % name) fp.write('node [shape=point,label=""];\n') for a, b in graph.edges_iter(): - l = np.sqrt(((pts[a] - pts[b])**2).sum(-1)) + l = np.sqrt(((pts[a] - pts[b]) ** 2).sum(-1)) lengths.append(l) - fp.write("%s -- %s [len=%f];\n"%(a, b, l)) - fp.write("maxiter=1000000;\n"); + fp.write("%s -- %s [len=%f];\n" % (a, b, l)) + fp.write("maxiter=1000000;\n") fp.write("}") def read_dot(fname, pts): - """ - """ + """ """ import re + parse = re.compile(r'\s(\d+)\s\[label="", pos="([\d\.]+),([\d\.]+)".*];') data = np.zeros((len(pts), 2)) with open(fname) as fp: @@ -985,36 +1111,38 @@ def read_dot(fname, pts): fp.readline() fp.readline() fp.readline() - el = fp.readline().split(' ') - while el[1] != '--': - x, y = el[2][5:-2].split(',') + el = fp.readline().split(" ") + while el[1] != "--": + x, y = el[2][5:-2].split(",") data[int(el[0][1:])] = float(x), float(y) - el = fp.readline().split(' ') + el = fp.readline().split(" ") return data def write_decimated(path, pts, polys): - """ - """ + """ """ from .polyutils import boundary_edges, decimate + dpts, dpolys = decimate(pts, polys) - write_surf(path+'.smoothwm', dpts, dpolys) + write_surf(path + ".smoothwm", dpts, dpolys) edges = boundary_edges(dpolys) - data = np.zeros((len(dpts),), dtype=[('vert', '>i4'), ('x', '>f4'), ('y', '>f4'), ('z', '>f4')]) - data['vert'] = np.arange(len(dpts))+1 - data['vert'][edges] *= -1 - data['x'] = dpts[:,0] - data['y'] = dpts[:,1] - data['z'] = dpts[:,2] - with open(path+'.full.patch.3d', 'w') as fp: - fp.write(struct.pack('>i', -1)) - fp.write(struct.pack('>i', len(dpts))) + data = np.zeros( + (len(dpts),), dtype=[("vert", ">i4"), ("x", ">f4"), ("y", ">f4"), ("z", ">f4")] + ) + data["vert"] = np.arange(len(dpts)) + 1 + data["vert"][edges] *= -1 + data["x"] = dpts[:, 0] + data["y"] = dpts[:, 1] + data["z"] = dpts[:, 2] + with open(path + ".full.patch.3d", "w") as fp: + fp.write(struct.pack(">i", -1)) + fp.write(struct.pack(">i", len(dpts))) fp.write(data.tobytes()) class SpringLayout(object): - """ - """ + """ """ + def __init__(self, pts, polys, dpts=None, pins=None, stepsize=1, neighborhood=0): self.pts = pts self.polys = polys @@ -1046,20 +1174,20 @@ def __init__(self, pts, polys, dpts=None, pins=None, stepsize=1, neighborhood=0) if dpts is None: dpts = pts - #self.kdt = cKDTree(self.pts) + # self.kdt = cKDTree(self.pts) self._next = self.pts.copy() width = max(len(n) for n in self.neighbors) self._mask = np.zeros((len(pts), width), dtype=bool) self._move = np.zeros((len(pts), width, 3)) - #self._mean = np.zeros((len(pts), width)) + # self._mean = np.zeros((len(pts), width)) self._num = np.zeros((len(pts),)) self._dists = [] self._idx = [] for i, n in enumerate(self.neighbors): - self._mask[i, :len(n)] = True - self._dists.append(np.sqrt(((dpts[n] - dpts[i])**2).sum(-1))) - self._idx.append(np.ones((len(n),))*i) + self._mask[i, : len(n)] = True + self._dists.append(np.sqrt(((dpts[n] - dpts[i]) ** 2).sum(-1))) + self._idx.append(np.ones((len(n),)) * i) self._num[i] = len(n) self._dists = np.hstack(self._dists) self._idx = np.hstack(self._idx).astype(np.uint) @@ -1069,8 +1197,8 @@ def __init__(self, pts, polys, dpts=None, pins=None, stepsize=1, neighborhood=0) def _spring(self): svec = self.pts[self._neigh] - self.pts[self._idx] slen = np.sqrt((svec**2).sum(-1)) - force = (slen - self._dists) # / self._dists - svec /= slen[:,np.newaxis] + force = slen - self._dists # / self._dists + svec /= slen[:, np.newaxis] fvec = force[:, np.newaxis] * svec self._move[self._mask] = self.stepsize * fvec return self._move.sum(1) / self._num[:, np.newaxis] @@ -1084,10 +1212,10 @@ def _estatic(self, idx): def step(self): move = self._spring()[~self.pins] - self._next[~self.pins] += move #+ self._estatic(i) + self._next[~self.pins] += move # + self._estatic(i) self.pts = self._next.copy() - return dict(x=self.pts[:,0],y=self.pts[:, 1], z=self.pts[:,2]), move - #self.kdt = cKDTree(self.pts) + return dict(x=self.pts[:, 0], y=self.pts[:, 1], z=self.pts[:, 2]), move + # self.kdt = cKDTree(self.pts) def run(self, n=1000): for _ in range(n): @@ -1096,27 +1224,106 @@ def run(self, n=1000): def view_step(self): from mayavi import mlab + if self.figure is None: - self.figure = mlab.triangular_mesh(self.pts[:,0], self.pts[:,1], self.pts[:,2], self.polys, representation='wireframe') + self.figure = mlab.triangular_mesh( + self.pts[:, 0], + self.pts[:, 1], + self.pts[:, 2], + self.polys, + representation="wireframe", + ) self.step() - self.figure.mlab_source.set(x=self.pts[:,0], y=self.pts[:,1], z=self.pts[:,2]) + self.figure.mlab_source.set( + x=self.pts[:, 0], y=self.pts[:, 1], z=self.pts[:, 2] + ) + def stretch_mwall(pts, polys, mwall): - """ - """ + """ """ inflated = pts.copy() center = pts[mwall].mean(0) radius = max((pts.max(0) - pts.min(0))[1:]) - angles = np.arctan2(pts[mwall][:,2], pts[mwall][:,1]) + angles = np.arctan2(pts[mwall][:, 2], pts[mwall][:, 1]) pts[mwall, 0] = center[0] pts[mwall, 1] = radius * np.cos(angles) + center[1] 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 -): +def _n_vertices_ico(icoorder): + return 4**icoorder * 10 + 2 + + +# Precomputed upsampling neighbor arrays bundled with pycortex (see +# ``_get_upsample_neighbors``). Lazily loaded ``NpzFile`` keyed by +# ``"{hemi}_ico{order}"``. +_UPSAMPLE_NEIGHBORS_FILE = os.path.join( + os.path.dirname(__file__), "data", "upsample_fsaverage_neighbors.npz" +) +_upsample_neighbors_bundle = None + + +def _get_upsample_neighbors(hemi, ico_order, freesurfer_subjects_dir=None): + """Nearest low-res neighbor index for each extra fsaverage vertex. + + Maps every full-``fsaverage`` vertex beyond the first ``n_ico`` (i.e. the + ones absent from ``fsaverage{ico_order}``) to its nearest ``fsaverage`` + -sphere neighbor among those first ``n_ico`` vertices. This index array is a + fixed property of the standard fsaverage tessellation -- it depends only on + ``(hemi, ico_order)``, never on the data -- so it is precomputed. + + Resolution order: + + 1. the arrays bundled with pycortex (covers the common fsaverage5/6 case, so + no FreeSurfer install or ``$SUBJECTS_DIR`` is required); + 2. a previously cached array in the pycortex user cache dir; + 3. otherwise it is computed from the fsaverage ``?h.sphere.reg`` (which needs + ``freesurfer_subjects_dir``/``$SUBJECTS_DIR``) and cached for next time. + """ + global _upsample_neighbors_bundle + key = f"{hemi}_ico{ico_order}" + + # 1. Bundled with the package. + if _upsample_neighbors_bundle is None: + if os.path.exists(_UPSAMPLE_NEIGHBORS_FILE): + _upsample_neighbors_bundle = np.load(_UPSAMPLE_NEIGHBORS_FILE) + else: + _upsample_neighbors_bundle = {} + bundle_keys = getattr( + _upsample_neighbors_bundle, "files", list(_upsample_neighbors_bundle) + ) + if key in bundle_keys: + return np.asarray(_upsample_neighbors_bundle[key]) + + # 2. User cache. + cache_dir = os.path.join(appdirs.user_cache_dir("pycortex"), "upsample_fsaverage") + cache_path = os.path.join(cache_dir, key + ".npy") + if os.path.exists(cache_path): + return np.load(cache_path) + + # 3. Compute from the fsaverage sphere and cache. + if freesurfer_subjects_dir is None: + freesurfer_subjects_dir = os.environ.get("SUBJECTS_DIR", None) + if freesurfer_subjects_dir is None: + raise ValueError( + f"Upsampling neighbors for fsaverage{ico_order} are not bundled " + "with pycortex and have not been cached yet. Set $SUBJECTS_DIR (or " + "pass freesurfer_subjects_dir) so they can be computed from the " + "fsaverage ?h.sphere.reg." + ) + n_ico = _n_vertices_ico(ico_order) + pts, _ = nibabel.freesurfer.read_geometry( + os.path.join(freesurfer_subjects_dir, "fsaverage", "surf", f"{hemi}.sphere.reg") + ) + _, neighbors = KDTree(pts[:n_ico]).query(pts[n_ico:], k=1) + neighbors = neighbors.astype(np.uint16 if neighbors.max() < 65536 else np.int64) + os.makedirs(cache_dir, exist_ok=True) + np.save(cache_path, neighbors) + return neighbors + + +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. @@ -1129,8 +1336,10 @@ def upsample_to_fsaverage( 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. + Path to Freesurfer subjects directory. Only needed for source spaces + whose neighbor arrays are neither bundled with pycortex (fsaverage5 and + fsaverage6 are) nor already cached; in that case it defaults to the + ``$SUBJECTS_DIR`` environment variable. Returns ------- @@ -1140,19 +1349,19 @@ def upsample_to_fsaverage( 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 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 + The per-vertex nearest-neighbor mapping is a fixed property of the fsaverage + tessellation, so it is precomputed (see :func:`_get_upsample_neighbors`). + For the common fsaverage5/fsaverage6 sources the arrays ship with pycortex, + so neither FreeSurfer nor ``$SUBJECTS_DIR`` is required. + """ ico_order = int(data_space[-1]) - n_ico_vertices = get_n_vertices_ico(ico_order) + n_ico_vertices = _n_vertices_ico(ico_order) ndim = data.ndim data = np.atleast_2d(data) _, n_vertices = data.shape @@ -1162,31 +1371,13 @@ def get_n_vertices_ico(icoorder): 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) - ) + for hemi, dt in zip(hemis, data_hemi): + neighbors = _get_upsample_neighbors(hemi, ico_order, freesurfer_subjects_dir) + # Fill the extra fsaverage vertices from their nearest low-res neighbor. + projected_data.append(np.concatenate([dt, dt[:, neighbors]], axis=-1)) projected_data = np.hstack(projected_data) if ndim == 1: projected_data = projected_data[0] @@ -1194,206 +1385,209 @@ def get_n_vertices_ico(icoorder): # aseg partition labels (up to 256 only) -fs_aseg_dict = {'Unknown': 0, - 'Left-Cerebral-Exterior': 1, - 'Left-Cerebral-White-Matter': 2, - 'Left-Cerebral-Cortex': 3, - 'Left-Lateral-Ventricle': 4, - 'Left-Inf-Lat-Vent': 5, - 'Left-Cerebellum-Exterior': 6, - 'Left-Cerebellum-White-Matter': 7, - 'Left-Cerebellum-Cortex': 8, - 'Left-Thalamus': 9, - 'Left-Thalamus-Proper': 10, - 'Left-Caudate': 11, - 'Left-Putamen': 12, - 'Left-Pallidum': 13, - '3rd-Ventricle': 14, - '4th-Ventricle': 15, - 'Brain-Stem': 16, - 'Left-Hippocampus': 17, - 'Left-Amygdala': 18, - 'Left-Insula': 19, - 'Left-Operculum': 20, - 'Line-1': 21, - 'Line-2': 22, - 'Line-3': 23, - 'CSF': 24, - 'Left-Lesion': 25, - 'Left-Accumbens-area': 26, - 'Left-Substancia-Nigra': 27, - 'Left-VentralDC': 28, - 'Left-undetermined': 29, - 'Left-vessel': 30, - 'Left-choroid-plexus': 31, - 'Left-F3orb': 32, - 'Left-lOg': 33, - 'Left-aOg': 34, - 'Left-mOg': 35, - 'Left-pOg': 36, - 'Left-Stellate': 37, - 'Left-Porg': 38, - 'Left-Aorg': 39, - 'Right-Cerebral-Exterior': 40, - 'Right-Cerebral-White-Matter': 41, - 'Right-Cerebral-Cortex': 42, - 'Right-Lateral-Ventricle': 43, - 'Right-Inf-Lat-Vent': 44, - 'Right-Cerebellum-Exterior': 45, - 'Right-Cerebellum-White-Matter': 46, - 'Right-Cerebellum-Cortex': 47, - 'Right-Thalamus': 48, - 'Right-Thalamus-Proper': 49, - 'Right-Caudate': 50, - 'Right-Putamen': 51, - 'Right-Pallidum': 52, - 'Right-Hippocampus': 53, - 'Right-Amygdala': 54, - 'Right-Insula': 55, - 'Right-Operculum': 56, - 'Right-Lesion': 57, - 'Right-Accumbens-area': 58, - 'Right-Substancia-Nigra': 59, - 'Right-VentralDC': 60, - 'Right-undetermined': 61, - 'Right-vessel': 62, - 'Right-choroid-plexus': 63, - 'Right-F3orb': 64, - 'Right-lOg': 65, - 'Right-aOg': 66, - 'Right-mOg': 67, - 'Right-pOg': 68, - 'Right-Stellate': 69, - 'Right-Porg': 70, - 'Right-Aorg': 71, - '5th-Ventricle': 72, - 'Left-Interior': 73, - 'Right-Interior': 74, - 'Left-Lateral-Ventricles': 75, - 'Right-Lateral-Ventricles': 76, - 'WM-hypointensities': 77, - 'Left-WM-hypointensities': 78, - 'Right-WM-hypointensities': 79, - 'non-WM-hypointensities': 80, - 'Left-non-WM-hypointensities': 81, - 'Right-non-WM-hypointensities': 82, - 'Left-F1': 83, - 'Right-F1': 84, - 'Optic-Chiasm': 85, - 'Corpus_Callosum': 86, - 'Left-Amygdala-Anterior': 96, - 'Right-Amygdala-Anterior': 97, - 'Dura': 98, - 'Left-wm-intensity-abnormality': 100, - 'Left-caudate-intensity-abnormality': 101, - 'Left-putamen-intensity-abnormality': 102, - 'Left-accumbens-intensity-abnormality': 103, - 'Left-pallidum-intensity-abnormality': 104, - 'Left-amygdala-intensity-abnormality': 105, - 'Left-hippocampus-intensity-abnormality': 106, - 'Left-thalamus-intensity-abnormality': 107, - 'Left-VDC-intensity-abnormality': 108, - 'Right-wm-intensity-abnormality': 109, - 'Right-caudate-intensity-abnormality': 110, - 'Right-putamen-intensity-abnormality': 111, - 'Right-accumbens-intensity-abnormality': 112, - 'Right-pallidum-intensity-abnormality': 113, - 'Right-amygdala-intensity-abnormality': 114, - 'Right-hippocampus-intensity-abnormality': 115, - 'Right-thalamus-intensity-abnormality': 116, - 'Right-VDC-intensity-abnormality': 117, - 'Epidermis': 118, - 'Conn-Tissue': 119, - 'SC-Fat/Muscle': 120, - 'Cranium': 121, - 'CSF-SA': 122, - 'Muscle': 123, - 'Ear': 124, - 'Adipose': 125, - 'Spinal-Cord': 126, - 'Soft-Tissue': 127, - 'Nerve': 128, - 'Bone': 129, - 'Air': 130, - 'Orbital-Fat': 131, - 'Tongue': 132, - 'Nasal-Structures': 133, - 'Globe': 134, - 'Teeth': 135, - 'Left-Caudate/Putamen': 136, - 'Right-Caudate/Putamen': 137, - 'Left-Claustrum': 138, - 'Right-Claustrum': 139, - 'Cornea': 140, - 'Diploe': 142, - 'Vitreous-Humor': 143, - 'Lens': 144, - 'Aqueous-Humor': 145, - 'Outer-Table': 146, - 'Inner-Table': 147, - 'Periosteum': 148, - 'Endosteum': 149, - 'R/C/S': 150, - 'Iris': 151, - 'SC-Adipose/Muscle': 152, - 'SC-Tissue': 153, - 'Orbital-Adipose': 154, - 'Left-IntCapsule-Ant': 155, - 'Right-IntCapsule-Ant': 156, - 'Left-IntCapsule-Pos': 157, - 'Right-IntCapsule-Pos': 158, - 'Left-Cerebral-WM-unmyelinated': 159, - 'Right-Cerebral-WM-unmyelinated': 160, - 'Left-Cerebral-WM-myelinated': 161, - 'Right-Cerebral-WM-myelinated': 162, - 'Left-Subcortical-Gray-Matter': 163, - 'Right-Subcortical-Gray-Matter': 164, - 'Skull': 165, - 'Posterior-fossa': 166, - 'Scalp': 167, - 'Hematoma': 168, - 'Left-Cortical-Dysplasia': 180, - 'Right-Cortical-Dysplasia': 181, - 'Left-hippocampal_fissure': 193, - 'Left-CADG-head': 194, - 'Left-subiculum': 195, - 'Left-fimbria': 196, - 'Right-hippocampal_fissure': 197, - 'Right-CADG-head': 198, - 'Right-subiculum': 199, - 'Right-fimbria': 200, - 'alveus': 201, - 'perforant_pathway': 202, - 'parasubiculum': 203, - 'presubiculum': 204, - 'subiculum': 205, - 'CA1': 206, - 'CA2': 207, - 'CA3': 208, - 'CA4': 209, - 'GC-DG': 210, - 'HATA': 211, - 'fimbria': 212, - 'lateral_ventricle': 213, - 'molecular_layer_HP': 214, - 'hippocampal_fissure': 215, - 'entorhinal_cortex': 216, - 'molecular_layer_subiculum': 217, - 'Amygdala': 218, - 'Cerebral_White_Matter': 219, - 'Cerebral_Cortex': 220, - 'Inf_Lat_Vent': 221, - 'Perirhinal': 222, - 'Cerebral_White_Matter_Edge': 223, - 'Background': 224, - 'Ectorhinal': 225, - 'Fornix': 250, - 'CC_Posterior': 251, - 'CC_Mid_Posterior': 252, - 'CC_Central': 253, - 'CC_Mid_Anterior': 254, - 'CC_Anterior': 255} +fs_aseg_dict = { + "Unknown": 0, + "Left-Cerebral-Exterior": 1, + "Left-Cerebral-White-Matter": 2, + "Left-Cerebral-Cortex": 3, + "Left-Lateral-Ventricle": 4, + "Left-Inf-Lat-Vent": 5, + "Left-Cerebellum-Exterior": 6, + "Left-Cerebellum-White-Matter": 7, + "Left-Cerebellum-Cortex": 8, + "Left-Thalamus": 9, + "Left-Thalamus-Proper": 10, + "Left-Caudate": 11, + "Left-Putamen": 12, + "Left-Pallidum": 13, + "3rd-Ventricle": 14, + "4th-Ventricle": 15, + "Brain-Stem": 16, + "Left-Hippocampus": 17, + "Left-Amygdala": 18, + "Left-Insula": 19, + "Left-Operculum": 20, + "Line-1": 21, + "Line-2": 22, + "Line-3": 23, + "CSF": 24, + "Left-Lesion": 25, + "Left-Accumbens-area": 26, + "Left-Substancia-Nigra": 27, + "Left-VentralDC": 28, + "Left-undetermined": 29, + "Left-vessel": 30, + "Left-choroid-plexus": 31, + "Left-F3orb": 32, + "Left-lOg": 33, + "Left-aOg": 34, + "Left-mOg": 35, + "Left-pOg": 36, + "Left-Stellate": 37, + "Left-Porg": 38, + "Left-Aorg": 39, + "Right-Cerebral-Exterior": 40, + "Right-Cerebral-White-Matter": 41, + "Right-Cerebral-Cortex": 42, + "Right-Lateral-Ventricle": 43, + "Right-Inf-Lat-Vent": 44, + "Right-Cerebellum-Exterior": 45, + "Right-Cerebellum-White-Matter": 46, + "Right-Cerebellum-Cortex": 47, + "Right-Thalamus": 48, + "Right-Thalamus-Proper": 49, + "Right-Caudate": 50, + "Right-Putamen": 51, + "Right-Pallidum": 52, + "Right-Hippocampus": 53, + "Right-Amygdala": 54, + "Right-Insula": 55, + "Right-Operculum": 56, + "Right-Lesion": 57, + "Right-Accumbens-area": 58, + "Right-Substancia-Nigra": 59, + "Right-VentralDC": 60, + "Right-undetermined": 61, + "Right-vessel": 62, + "Right-choroid-plexus": 63, + "Right-F3orb": 64, + "Right-lOg": 65, + "Right-aOg": 66, + "Right-mOg": 67, + "Right-pOg": 68, + "Right-Stellate": 69, + "Right-Porg": 70, + "Right-Aorg": 71, + "5th-Ventricle": 72, + "Left-Interior": 73, + "Right-Interior": 74, + "Left-Lateral-Ventricles": 75, + "Right-Lateral-Ventricles": 76, + "WM-hypointensities": 77, + "Left-WM-hypointensities": 78, + "Right-WM-hypointensities": 79, + "non-WM-hypointensities": 80, + "Left-non-WM-hypointensities": 81, + "Right-non-WM-hypointensities": 82, + "Left-F1": 83, + "Right-F1": 84, + "Optic-Chiasm": 85, + "Corpus_Callosum": 86, + "Left-Amygdala-Anterior": 96, + "Right-Amygdala-Anterior": 97, + "Dura": 98, + "Left-wm-intensity-abnormality": 100, + "Left-caudate-intensity-abnormality": 101, + "Left-putamen-intensity-abnormality": 102, + "Left-accumbens-intensity-abnormality": 103, + "Left-pallidum-intensity-abnormality": 104, + "Left-amygdala-intensity-abnormality": 105, + "Left-hippocampus-intensity-abnormality": 106, + "Left-thalamus-intensity-abnormality": 107, + "Left-VDC-intensity-abnormality": 108, + "Right-wm-intensity-abnormality": 109, + "Right-caudate-intensity-abnormality": 110, + "Right-putamen-intensity-abnormality": 111, + "Right-accumbens-intensity-abnormality": 112, + "Right-pallidum-intensity-abnormality": 113, + "Right-amygdala-intensity-abnormality": 114, + "Right-hippocampus-intensity-abnormality": 115, + "Right-thalamus-intensity-abnormality": 116, + "Right-VDC-intensity-abnormality": 117, + "Epidermis": 118, + "Conn-Tissue": 119, + "SC-Fat/Muscle": 120, + "Cranium": 121, + "CSF-SA": 122, + "Muscle": 123, + "Ear": 124, + "Adipose": 125, + "Spinal-Cord": 126, + "Soft-Tissue": 127, + "Nerve": 128, + "Bone": 129, + "Air": 130, + "Orbital-Fat": 131, + "Tongue": 132, + "Nasal-Structures": 133, + "Globe": 134, + "Teeth": 135, + "Left-Caudate/Putamen": 136, + "Right-Caudate/Putamen": 137, + "Left-Claustrum": 138, + "Right-Claustrum": 139, + "Cornea": 140, + "Diploe": 142, + "Vitreous-Humor": 143, + "Lens": 144, + "Aqueous-Humor": 145, + "Outer-Table": 146, + "Inner-Table": 147, + "Periosteum": 148, + "Endosteum": 149, + "R/C/S": 150, + "Iris": 151, + "SC-Adipose/Muscle": 152, + "SC-Tissue": 153, + "Orbital-Adipose": 154, + "Left-IntCapsule-Ant": 155, + "Right-IntCapsule-Ant": 156, + "Left-IntCapsule-Pos": 157, + "Right-IntCapsule-Pos": 158, + "Left-Cerebral-WM-unmyelinated": 159, + "Right-Cerebral-WM-unmyelinated": 160, + "Left-Cerebral-WM-myelinated": 161, + "Right-Cerebral-WM-myelinated": 162, + "Left-Subcortical-Gray-Matter": 163, + "Right-Subcortical-Gray-Matter": 164, + "Skull": 165, + "Posterior-fossa": 166, + "Scalp": 167, + "Hematoma": 168, + "Left-Cortical-Dysplasia": 180, + "Right-Cortical-Dysplasia": 181, + "Left-hippocampal_fissure": 193, + "Left-CADG-head": 194, + "Left-subiculum": 195, + "Left-fimbria": 196, + "Right-hippocampal_fissure": 197, + "Right-CADG-head": 198, + "Right-subiculum": 199, + "Right-fimbria": 200, + "alveus": 201, + "perforant_pathway": 202, + "parasubiculum": 203, + "presubiculum": 204, + "subiculum": 205, + "CA1": 206, + "CA2": 207, + "CA3": 208, + "CA4": 209, + "GC-DG": 210, + "HATA": 211, + "fimbria": 212, + "lateral_ventricle": 213, + "molecular_layer_HP": 214, + "hippocampal_fissure": 215, + "entorhinal_cortex": 216, + "molecular_layer_subiculum": 217, + "Amygdala": 218, + "Cerebral_White_Matter": 219, + "Cerebral_Cortex": 220, + "Inf_Lat_Vent": 221, + "Perirhinal": 222, + "Cerebral_White_Matter_Edge": 223, + "Background": 224, + "Ectorhinal": 225, + "Fornix": 250, + "CC_Posterior": 251, + "CC_Mid_Posterior": 252, + "CC_Central": 253, + "CC_Mid_Anterior": 254, + "CC_Anterior": 255, +} if __name__ == "__main__": import sys + show_surf(sys.argv[1], sys.argv[2], sys.argv[3]) diff --git a/cortex/hcp.py b/cortex/hcp.py new file mode 100644 index 00000000..f7093c07 --- /dev/null +++ b/cortex/hcp.py @@ -0,0 +1,621 @@ +# emacs: -*- coding: utf-8; mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set fileencoding=utf-8 ft=python sts=4 ts=4 sw=4 et: +"""Interfacing with HCP fs_LR 32k data. + +This module lets you visualize HCP data with pycortex, both natively on the +fs_LR 32k mesh (via the ``32k_fs_LR`` pycortex subject) and projected to +``fsaverage``. + +It provides three things: + +1. :func:`download_fs_lr` -- fetch the prebuilt ``32k_fs_LR`` pycortex subject + into the filestore (a thin wrapper around :func:`cortex.download_subject`), + so HCP data can be rendered with the usual + ``cortex.Vertex(data, "32k_fs_LR")`` + ``cortex.quickshow`` path. +2. :func:`cifti_to_surface` -- expand CIFTI grayordinate data (59412 cortical + vertices, no medial wall) to the full 64984-vertex fs_LR 32k surface, filling + the medial wall with NaN. +3. Resampling fs_LR 32k data to ``fsaverage`` (:func:`project_fslr_to_fsaverage`, + :func:`to_fsaverage`). The projection matrix is built directly from the HCP + standard-mesh spheres using spherical barycentric interpolation -- the same + interpolation ``wb_command -metric-resample ... BARYCENTRIC`` performs -- so + no Connectome Workbench installation is required at runtime. + +Notes +----- +The ``fs_LR-deformed_to-fsaverage`` source sphere and the ``fsaverageN`` +standard spheres live in the same spherical-registration space, so a target +vertex can be located inside a source-mesh triangle directly (see +:func:`_barycentric_resample_matrix`). This mirrors, and reuses the caching +conventions of, :func:`cortex.freesurfer.get_mri_surf2surf_matrix`. +""" + +import logging +from pathlib import Path +from urllib.error import URLError +from urllib.request import urlretrieve + +import numpy as np +import nibabel as nib +from scipy.sparse import coo_matrix, load_npz, save_npz +from scipy.spatial import KDTree + +from cortex import appdirs +from cortex.freesurfer import upsample_to_fsaverage +from cortex.utils import download_subject + +logger = logging.getLogger(__name__) + +# --------------------------------------------------------------------------- +# Constants +# --------------------------------------------------------------------------- + +#: Name of the pycortex subject shipping the HCP fs_LR 32k surfaces. +SUBJECT = "32k_fs_LR" + +#: Vertices per hemisphere on the fs_LR 32k surface. +N_VERTICES_FS_LR_32K_HEM = 32492 +#: Vertices on both fs_LR 32k hemispheres (full surface, includes medial wall). +N_VERTICES_FS_LR_32K = 2 * N_VERTICES_FS_LR_32K_HEM # 64984 +#: Cortical grayordinates in a standard HCP CIFTI file (no medial wall). +N_VERTICES_CIFTI_CORTEX = 59412 + +#: Vertices per hemisphere on the target fsaverage meshes. +N_VERTICES_TARGET_HEM = { + "fsaverage5": 10242, + "fsaverage6": 40962, +} + +# --------------------------------------------------------------------------- +# Sphere file download helpers +# --------------------------------------------------------------------------- + +_SPHERE_FILES_URL = ( + "https://raw.githubusercontent.com/Washington-University/HCPpipelines/" + "master/global/templates/standard_mesh_atlases/resample_fsaverage/" +) + +#: Standard-mesh sphere GIFTIs, keyed by space then hemisphere ("L"/"R"). +_SPHERE_FILES = { + "fs_LR_32k": { + "L": "fs_LR-deformed_to-fsaverage.L.sphere.32k_fs_LR.surf.gii", + "R": "fs_LR-deformed_to-fsaverage.R.sphere.32k_fs_LR.surf.gii", + }, + "fsaverage5": { + "L": "fsaverage5_std_sphere.L.10k_fsavg_L.surf.gii", + "R": "fsaverage5_std_sphere.R.10k_fsavg_R.surf.gii", + }, + "fsaverage6": { + "L": "fsaverage6_std_sphere.L.41k_fsavg_L.surf.gii", + "R": "fsaverage6_std_sphere.R.41k_fsavg_R.surf.gii", + }, +} + + +def _default_cache_dir(): + """Directory for downloaded spheres and cached projection matrices.""" + return Path(appdirs.user_cache_dir("pycortex")) / "hcp" + + +def ensure_sphere_files(space, hemi, cache_dir=None, download=True): + """Return the path to a standard-mesh sphere GIFTI, downloading if needed. + + Parameters + ---------- + space : {"fs_LR_32k", "fsaverage5", "fsaverage6"} + Surface space of the sphere. + hemi : {"L", "R"} + Hemisphere. + cache_dir : str or Path or None + Directory to store sphere files. Defaults to a pycortex user cache dir. + download : bool + If True (default), download the file from the HCPpipelines repository + when it is missing. + + Returns + ------- + path : pathlib.Path + Path to the sphere ``.surf.gii`` file. + """ + if space not in _SPHERE_FILES: + raise ValueError( + f"Unknown space {space!r}; choose from {sorted(_SPHERE_FILES)}" + ) + if hemi not in ("L", "R"): + raise ValueError(f"hemi must be 'L' or 'R', got {hemi!r}") + + if cache_dir is None: + cache_dir = _default_cache_dir() + atlas_dir = Path(cache_dir) / "standard_mesh_atlases" + atlas_dir.mkdir(parents=True, exist_ok=True) + + fname = _SPHERE_FILES[space][hemi] + fpath = atlas_dir / fname + if fpath.exists(): + return fpath + if not download: + raise FileNotFoundError( + f"Sphere file not found: {fpath}. Set download=True to fetch it." + ) + + url = _SPHERE_FILES_URL + fname + logger.info("Downloading %s -> %s", url, fpath) + tmp_path = fpath.with_suffix(".tmp") + try: + urlretrieve(url, tmp_path) + except (URLError, OSError) as exc: + tmp_path.unlink(missing_ok=True) + raise RuntimeError( + f"Failed to download sphere file from {url}. Check your internet " + f"connection or download it manually to {fpath}." + ) from exc + if tmp_path.stat().st_size == 0: + tmp_path.unlink() + raise RuntimeError(f"Downloaded file is empty: {url}") + tmp_path.rename(fpath) + return fpath + + +# --------------------------------------------------------------------------- +# Subject download +# --------------------------------------------------------------------------- + + +def download_fs_lr(pycortex_store=None, download_again=False): + """Download the ``32k_fs_LR`` pycortex subject into the filestore. + + Thin wrapper around :func:`cortex.download_subject` for the HCP fs_LR 32k + surfaces. Once present, HCP data can be rendered natively with + ``cortex.Vertex(data, "32k_fs_LR")`` and ``cortex.quickshow``. + + Parameters + ---------- + pycortex_store : str or None + Directory to place the subject folder. If None, uses the current + filestore (``cortex.db.filestore``). + download_again : bool + Re-download even if the subject is already present. + + Notes + ----- + The ``32k_fs_LR`` subject is derived from HCP S1200 group-average Open + Access surfaces and is redistributed under the WU-Minn HCP Consortium Open + Access Data Use Terms; use of it requires the standard HCP acknowledgment. + """ + download_subject( + subject_id=SUBJECT, pycortex_store=pycortex_store, download_again=download_again + ) + + +# --------------------------------------------------------------------------- +# CIFTI -> fs_LR 32k surface expansion +# --------------------------------------------------------------------------- + + +def get_cifti_vertex_indices(cifti): + """Left/right surface-vertex indices for a CIFTI file's cortical models. + + HCP CIFTI files store only non-medial-wall cortical grayordinates. This + reads, from the CIFTI brain-model axis, which fs_LR 32k surface vertices + each hemisphere's grayordinates correspond to. + + Parameters + ---------- + cifti : str or Path or nibabel.Cifti2Image + A CIFTI-2 file (``.dscalar.nii``/``.dtseries.nii``/...) or loaded image. + + Returns + ------- + left_indices : ndarray, shape (n_left,) + Surface vertex indices (into a 32492-vertex hemisphere) for the left + cortex grayordinates. + right_indices : ndarray, shape (n_right,) + Same for the right cortex. + """ + if isinstance(cifti, (str, Path)): + cifti = nib.load(str(cifti)) + if not isinstance(cifti, nib.Cifti2Image): + raise TypeError("cifti must be a Cifti2Image or a path to a CIFTI-2 file.") + + # Find the brain-model axis (the one carrying the grayordinate structures). + bm_axis = None + for i in range(cifti.ndim): + axis = cifti.header.get_axis(i) + if isinstance(axis, nib.cifti2.BrainModelAxis): + bm_axis = axis + break + if bm_axis is None: + raise ValueError("CIFTI file has no BrainModelAxis (grayordinate axis).") + + indices = {} + for name, data_slice, model in bm_axis.iter_structures(): + if name == "CIFTI_STRUCTURE_CORTEX_LEFT": + indices["L"] = np.asarray(model.vertex) + elif name == "CIFTI_STRUCTURE_CORTEX_RIGHT": + indices["R"] = np.asarray(model.vertex) + if "L" not in indices and "R" not in indices: + raise ValueError("CIFTI file is missing both left and right cortex structures.") + # Tolerate hemispherically split files: a missing hemisphere is an empty + # index set, which cifti_to_surface / project_fslr_to_fsaverage handle by + # leaving that hemisphere as NaN. + empty = np.array([], dtype=np.int64) + return indices.get("L", empty), indices.get("R", empty) + + +def cifti_to_surface(data, left_indices, right_indices): + """Expand CIFTI cortical grayordinates onto the full fs_LR 32k surface. + + Parameters + ---------- + data : ndarray, shape (n_grayordinates,) or (n_samples, n_grayordinates) + Cortical CIFTI data (59412 vertices for standard HCP files). Only the + cortical grayordinates are used; the left hemisphere's values come + first, followed by the right hemisphere's. + left_indices, right_indices : ndarray + Surface vertex indices per hemisphere, e.g. from + :func:`get_cifti_vertex_indices`. + + Returns + ------- + full_data : ndarray, shape (64984,) or (n_samples, 64984) + Data on the full fs_LR 32k surface, with NaN at medial-wall vertices. + """ + data = np.asarray(data) + left_indices = np.asarray(left_indices) + right_indices = np.asarray(right_indices) + n_left = len(left_indices) + n_right = len(right_indices) + right_offset = N_VERTICES_FS_LR_32K_HEM + + out_dtype = np.result_type(data.dtype, np.float32) + if data.ndim == 1: + full = np.full(N_VERTICES_FS_LR_32K, np.nan, dtype=out_dtype) + full[left_indices] = data[:n_left] + full[right_offset + right_indices] = data[n_left : n_left + n_right] + else: + full = np.full((data.shape[0], N_VERTICES_FS_LR_32K), np.nan, dtype=out_dtype) + full[:, left_indices] = data[:, :n_left] + full[:, right_offset + right_indices] = data[:, n_left : n_left + n_right] + return full + + +# --------------------------------------------------------------------------- +# Spherical barycentric resampling matrix +# --------------------------------------------------------------------------- + + +def _read_sphere(path): + """Return (points, triangles) from a sphere ``.surf.gii`` file.""" + gii = nib.load(str(path)) + pts = gii.get_arrays_from_intent("NIFTI_INTENT_POINTSET")[0].data + tris = gii.get_arrays_from_intent("NIFTI_INTENT_TRIANGLE")[0].data + return np.asarray(pts, dtype=np.float64), np.asarray(tris, dtype=np.int64) + + +def _vertex_to_triangles(triangles, n_vertices): + """List of incident triangle indices for each vertex.""" + incident = [[] for _ in range(n_vertices)] + for ti, tri in enumerate(triangles): + for v in tri: + incident[v].append(ti) + return incident + + +def _ray_triangle_bary(direction, a, b, c, tol=1e-6): + """Barycentric weights where the ray ``t*direction`` (from 0) hits triangle abc. + + Uses the Moller-Trumbore intersection. Returns ``(wa, wb, wc)`` summing to 1 + if the central projection of ``direction`` onto the plane of ``abc`` lies + inside the triangle (within ``tol``), else ``None``. Because candidate + triangles are restricted to those incident to the target's nearest source + vertices, only the correct near-side triangle is ever tested. + """ + edge1 = b - a + edge2 = c - a + pvec = np.cross(direction, edge2) + det = np.dot(edge1, pvec) + if abs(det) < 1e-12: + return None + inv_det = 1.0 / det + tvec = -a + u = np.dot(tvec, pvec) * inv_det + if u < -tol or u > 1.0 + tol: + return None + qvec = np.cross(tvec, edge1) + v = np.dot(direction, qvec) * inv_det + if v < -tol or u + v > 1.0 + tol: + return None + w = 1.0 - u - v + # Clamp tiny negative weights from edge/vertex hits, then renormalize. + weights = np.array([w, u, v], dtype=np.float64) + weights[weights < 0] = 0.0 + total = weights.sum() + if total <= 0: + return None + return tuple(weights / total) + + +def _barycentric_resample_matrix(src_sphere, tgt_sphere, k=15): + """Build a sparse barycentric resampling matrix between two spheres. + + For every target-sphere vertex, the containing triangle on the *source* + mesh is located and the target row is set to that triangle's three + barycentric weights -- exactly the interpolation + ``wb_command -metric-resample ... BARYCENTRIC`` performs. + + Parameters + ---------- + src_sphere, tgt_sphere : str or Path + Paths to the source and target sphere ``.surf.gii`` files. Both must be + in the same spherical-registration space. + k : int + Number of nearest source vertices whose incident triangles are searched + for the one containing each target vertex. + + Returns + ------- + matrix : scipy.sparse.csr_matrix, shape (n_target, n_source) + Apply with ``target_data = matrix.dot(source_data)``. Rows sum to 1. + """ + src_pts, src_tris = _read_sphere(src_sphere) + tgt_pts, _ = _read_sphere(tgt_sphere) + n_src = len(src_pts) + n_tgt = len(tgt_pts) + + # Project onto the unit sphere so the ray directions are well-conditioned. + src_unit = src_pts / np.linalg.norm(src_pts, axis=1, keepdims=True) + tgt_unit = tgt_pts / np.linalg.norm(tgt_pts, axis=1, keepdims=True) + + incident = _vertex_to_triangles(src_tris, n_src) + tree = KDTree(src_unit) + k = min(k, n_src) + _, nn = tree.query(tgt_unit, k=k) + if nn.ndim == 1: + nn = nn[:, np.newaxis] + + rows = np.empty(3 * n_tgt, dtype=np.int64) + cols = np.empty(3 * n_tgt, dtype=np.int64) + vals = np.empty(3 * n_tgt, dtype=np.float64) + n_entries = 0 + n_fallback = 0 + + for t in range(n_tgt): + direction = tgt_unit[t] + found = None + seen = set() + for v in nn[t]: + for ti in incident[v]: + if ti in seen: + continue + seen.add(ti) + ia, ib, ic = src_tris[ti] + w = _ray_triangle_bary( + direction, src_unit[ia], src_unit[ib], src_unit[ic] + ) + if w is not None: + found = (ia, ib, ic, w) + break + if found is not None: + break + + if found is None: + # No containing triangle among candidates: assign nearest vertex. + n_fallback += 1 + va = int(nn[t][0]) + rows[n_entries] = t + cols[n_entries] = va + vals[n_entries] = 1.0 + n_entries += 1 + continue + + ia, ib, ic, (wa, wb, wc) = found + rows[n_entries : n_entries + 3] = t + cols[n_entries : n_entries + 3] = (ia, ib, ic) + vals[n_entries : n_entries + 3] = (wa, wb, wc) + n_entries += 3 + + if n_fallback: + logger.warning( + "barycentric resample: %d/%d target vertices had no containing " + "source triangle within k=%d neighbours; used nearest-vertex " + "fallback.", + n_fallback, + n_tgt, + k, + ) + + matrix = coo_matrix( + (vals[:n_entries], (rows[:n_entries], cols[:n_entries])), shape=(n_tgt, n_src) + ).tocsr() + matrix.sum_duplicates() + return matrix + + +def get_fslr_to_fsaverage_matrix(hemi, target="fsaverage6", cache_dir=None, cache=True): + """Sparse fs_LR 32k -> fsaverage matrix for one hemisphere. + + Parameters + ---------- + hemi : {"L", "R"} + Hemisphere. + target : {"fsaverage5", "fsaverage6"} + Target fsaverage density. Use :func:`project_fslr_to_fsaverage` with + ``target="fsaverage"`` to reach the full-resolution surface. + cache_dir : str or Path or None + Directory for the sphere files and matrix cache. Defaults to a pycortex + user cache dir. + cache : bool + If True (default), load/save the computed matrix as an ``.npz`` file. + + Returns + ------- + matrix : scipy.sparse.csr_matrix, shape (n_target_hemi, 32492) + """ + if hemi not in ("L", "R"): + raise ValueError(f"hemi must be 'L' or 'R', got {hemi!r}") + if target not in N_VERTICES_TARGET_HEM: + raise ValueError( + f"target must be one of {sorted(N_VERTICES_TARGET_HEM)}, got {target!r}" + ) + + if cache_dir is None: + cache_dir = _default_cache_dir() + cache_dir = Path(cache_dir) + + cache_path = cache_dir / "mappers" / f"{hemi}_fs_LR_32k_to_{target}.npz" + if cache and cache_path.exists(): + logger.info("Loading cached matrix from %s", cache_path) + return load_npz(str(cache_path)) + + src_sphere = ensure_sphere_files("fs_LR_32k", hemi, cache_dir=cache_dir) + tgt_sphere = ensure_sphere_files(target, hemi, cache_dir=cache_dir) + matrix = _barycentric_resample_matrix(src_sphere, tgt_sphere) + + if cache: + cache_path.parent.mkdir(parents=True, exist_ok=True) + save_npz(str(cache_path), matrix) + logger.info("Saved matrix to %s", cache_path) + return matrix + + +# --------------------------------------------------------------------------- +# Projection +# --------------------------------------------------------------------------- + + +def _project_hemi(matrix, hemi_data, nanmean): + """Apply one hemisphere's resampling matrix with NaN-aware weighting.""" + nan_mask = np.isnan(hemi_data) + hemi_clean = np.where(nan_mask, 0.0, hemi_data) + projected = (matrix @ hemi_clean.T).T + + # Weight actually contributed by valid (non-NaN) source vertices. Matrix + # weights are non-negative (barycentric weights / unit fallback), so no abs. + valid_weight = (matrix @ (~nan_mask).astype(np.float64).T).T + all_nan = valid_weight < 1e-10 + if nanmean: + with np.errstate(invalid="ignore", divide="ignore"): + projected = np.where(all_nan, np.nan, projected / valid_weight) + else: + projected[all_nan] = np.nan + return projected + + +def project_fslr_to_fsaverage( + data, target="fsaverage", cache_dir=None, freesurfer_subjects_dir=None, nanmean=True +): + """Project fs_LR 32k data to an fsaverage surface. + + Parameters + ---------- + data : ndarray, shape (64984,) or (n_samples, 64984) + Concatenated left+right data on the full fs_LR 32k surface (medial wall + included; NaN there is fine). Use :func:`cifti_to_surface` to expand + CIFTI grayordinate data to this shape. + target : {"fsaverage", "fsaverage6", "fsaverage5"} + Destination surface. ``"fsaverage"`` (full, 327684 vertices) is reached + by projecting to fsaverage6 and then upsampling with + :func:`cortex.freesurfer.upsample_to_fsaverage`. + cache_dir : str or Path or None + Directory for sphere files and matrix caches. + freesurfer_subjects_dir : str or None + FreeSurfer ``SUBJECTS_DIR`` containing ``fsaverage``/``fsaverage6``. + Only required for ``target="fsaverage"``. If None, uses ``$SUBJECTS_DIR``. + nanmean : bool + If True (default), renormalize weights to exclude NaN sources so medial + wall / missing data do not dilute neighbouring valid values. + + Returns + ------- + projected : ndarray + Data on the target surface. Medial-wall and all-NaN-source vertices are + NaN; apply ``numpy.nan_to_num`` before handing to pycortex if needed. + """ + data = np.asarray(data) + # Preserve the input's floating precision (float32 timeseries are common and + # large); only integer inputs are promoted, to float32. + out_dtype = data.dtype if np.issubdtype(data.dtype, np.floating) else np.float32 + data = data.astype(out_dtype, copy=False) + squeeze = data.ndim == 1 + if squeeze: + data = data[np.newaxis, :] + if data.shape[1] != N_VERTICES_FS_LR_32K: + raise ValueError( + f"Expected {N_VERTICES_FS_LR_32K} fs_LR 32k vertices, got {data.shape[1]}." + ) + + matrix_target = "fsaverage6" if target == "fsaverage" else target + if matrix_target not in N_VERTICES_TARGET_HEM: + raise ValueError( + f"target must be 'fsaverage', 'fsaverage6', or 'fsaverage5', got {target!r}" + ) + + n_tgt_hemi = N_VERTICES_TARGET_HEM[matrix_target] + result = np.full((data.shape[0], 2 * n_tgt_hemi), np.nan, dtype=out_dtype) + for ih, hemi in enumerate(("L", "R")): + src_sl = slice( + ih * N_VERTICES_FS_LR_32K_HEM, (ih + 1) * N_VERTICES_FS_LR_32K_HEM + ) + tgt_sl = slice(ih * n_tgt_hemi, (ih + 1) * n_tgt_hemi) + matrix = get_fslr_to_fsaverage_matrix( + hemi, target=matrix_target, cache_dir=cache_dir + ) + result[:, tgt_sl] = _project_hemi(matrix, data[:, src_sl], nanmean) + + if target == "fsaverage": + result = upsample_to_fsaverage( + result, "fsaverage6", freesurfer_subjects_dir=freesurfer_subjects_dir + ) + + if squeeze: + result = result[0] + return result + + +def to_fsaverage( + cifti, + target="fsaverage", + cache_dir=None, + freesurfer_subjects_dir=None, + nanmean=True, +): + """Project HCP CIFTI cortical data straight to an fsaverage surface. + + Convenience wrapper chaining :func:`get_cifti_vertex_indices`, + :func:`cifti_to_surface`, and :func:`project_fslr_to_fsaverage`. + + Parameters + ---------- + cifti : str or Path or nibabel.Cifti2Image + A CIFTI-2 file (or loaded image) whose data will be projected. The + cortical grayordinates are read and expanded to the fs_LR 32k surface + before resampling. + target : {"fsaverage", "fsaverage6", "fsaverage5"} + Destination surface. + cache_dir, freesurfer_subjects_dir, nanmean + See :func:`project_fslr_to_fsaverage`. + + Returns + ------- + projected : ndarray + CIFTI data on the target fsaverage surface (NaN medial wall). + """ + if isinstance(cifti, (str, Path)): + cifti = nib.load(str(cifti)) + if not isinstance(cifti, nib.Cifti2Image): + raise TypeError("cifti must be a Cifti2Image or a path to a CIFTI-2 file.") + left_indices, right_indices = get_cifti_vertex_indices(cifti) + # CIFTI grayordinate axis is the last axis. A single-map file (e.g. a + # one-column dscalar) is squeezed to a single surface map for convenience; + # genuine multi-map/timeseries files keep their leading sample axis. + data = np.asarray(cifti.get_fdata()) + if data.ndim == 2 and data.shape[0] == 1: + data = data[0] + full = cifti_to_surface(data, left_indices, right_indices) + return project_fslr_to_fsaverage( + full, + target=target, + cache_dir=cache_dir, + freesurfer_subjects_dir=freesurfer_subjects_dir, + nanmean=nanmean, + ) diff --git a/cortex/tests/test_freesurfer.py b/cortex/tests/test_freesurfer.py index 2a16e3ad..ec8983c7 100644 --- a/cortex/tests/test_freesurfer.py +++ b/cortex/tests/test_freesurfer.py @@ -9,17 +9,13 @@ _remove_disconnected_polys, _surf2surf_nnfr_matrix, get_mri_surf2surf_matrix, + upsample_to_fsaverage, ) def test_remove_disconnected_polys_examples(): - polys = np.array([[0, 1, 2], - [0, 1, 3], - [1, 2, 4], - [5, 6, 7]]) - expected_result = np.array([[0, 1, 2], - [0, 1, 3], - [1, 2, 4]]) + polys = np.array([[0, 1, 2], [0, 1, 3], [1, 2, 4], [5, 6, 7]]) + expected_result = np.array([[0, 1, 2], [0, 1, 3], [1, 2, 4]]) result = _remove_disconnected_polys(polys) np.testing.assert_array_equal(result, expected_result) @@ -27,12 +23,12 @@ def test_remove_disconnected_polys_examples(): def test_remove_disconnected_polys_idempotence(): rng = np.random.RandomState(0) for n_polys in [10, 20, 30, 40]: - polys_0 =rng.randint(0, 100, size=3 * n_polys).reshape(-1, 3) - + polys_0 = rng.randint(0, 100, size=3 * n_polys).reshape(-1, 3) + # make sure this example filters something polys_1 = _remove_disconnected_polys(polys_0) assert len(polys_0) != len(polys_1) - + # make sure calling the function does not change anything polys_2 = _remove_disconnected_polys(polys_1) np.testing.assert_array_equal(polys_1, polys_2) @@ -42,6 +38,7 @@ def test_remove_disconnected_polys_idempotence(): # surf2surf (nnfr) matrix construction # --------------------------------------------------------------------------- + def _random_sphere(n, seed): """n points on the unit sphere (so KDTree distances behave like the real ?h.sphere.reg geometry).""" @@ -91,16 +88,20 @@ def test_surf2surf_known_averaging_example(): # Four collinear source vertices, two target vertices placed so that each # target's nearest source is an endpoint, leaving the two middle source # vertices as orphans that get folded into their nearest target. - src = np.array([[0., 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0]]) + src = np.array([[0.0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0]]) trg = np.array([[0.4, 0, 0], [2.6, 0, 0]]) m = _surf2surf_nnfr_matrix(src, trg) - expected = np.array([[0.5, 0.5, 0.0, 0.0], # mean of src 0 and 1 - [0.0, 0.0, 0.5, 0.5]]) # mean of src 2 and 3 + expected = np.array( + [ + [0.5, 0.5, 0.0, 0.0], # mean of src 0 and 1 + [0.0, 0.0, 0.5, 0.5], + ] + ) # mean of src 2 and 3 np.testing.assert_allclose(m.toarray(), expected) - data = np.array([10., 20., 30., 40.]) - np.testing.assert_allclose(m.dot(data), [15., 35.]) + data = np.array([10.0, 20.0, 30.0, 40.0]) + np.testing.assert_allclose(m.dot(data), [15.0, 35.0]) def test_get_mri_surf2surf_matrix_ignores_legacy_kwargs(monkeypatch): @@ -108,12 +109,15 @@ def test_get_mri_surf2surf_matrix_ignores_legacy_kwargs(monkeypatch): src = _random_sphere(20, seed=7) trg = _random_sphere(15, seed=8) monkeypatch.setattr( - fs, "_read_sphere_reg", - lambda subj, hemi, subjects_dir=None: src if subj == "A" else trg) + fs, + "_read_sphere_reg", + lambda subj, hemi, subjects_dir=None: src if subj == "A" else trg, + ) with pytest.warns(DeprecationWarning): - m = get_mri_surf2surf_matrix("A", "lh", target_subj="B", - n_neighbors=20, n_test_images=40) + m = get_mri_surf2surf_matrix( + "A", "lh", target_subj="B", n_neighbors=20, n_test_images=40 + ) assert m.shape == (15, 20) @@ -122,13 +126,45 @@ def test_get_mri_surf2surf_matrix_rejects_unknown_kwargs(): get_mri_surf2surf_matrix("A", "lh", target_subj="B", bogus_kwarg=1) +# --------------------------------------------------------------------------- +# upsample_to_fsaverage (bundled neighbor tables, no freesurfer needed) +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "data_space,n_src", [("fsaverage6", 81924), ("fsaverage5", 20484)] +) +def test_upsample_to_fsaverage_bundled_no_freesurfer(data_space, n_src, monkeypatch): + # fsaverage5/6 upsampling tables ship with pycortex, so this must work with + # no $SUBJECTS_DIR and no freesurfer install. + monkeypatch.delenv("SUBJECTS_DIR", raising=False) + rng = np.random.RandomState(0) + data = rng.randn(3, n_src) + out = upsample_to_fsaverage(data, data_space) + assert out.shape == (3, 327684) + # The low-resolution vertices are carried over unchanged. + np.testing.assert_array_equal(out[:, : n_src // 2], data[:, : n_src // 2]) + # A 1-D input yields a 1-D result. + assert upsample_to_fsaverage(data[0], data_space).ndim == 1 + + +def test_upsample_to_fsaverage_constant_preserved(monkeypatch): + # Nearest-neighbor fill of a constant map stays constant everywhere. + monkeypatch.delenv("SUBJECTS_DIR", raising=False) + out = upsample_to_fsaverage(np.full(81924, 2.5), "fsaverage6") + np.testing.assert_array_equal(out, np.full(327684, 2.5)) + + def _have_template(subjects_dir, name, hemi="lh"): return bool(subjects_dir) and os.path.exists( - os.path.join(subjects_dir, name, "surf", hemi + ".sphere.reg")) + os.path.join(subjects_dir, name, "surf", hemi + ".sphere.reg") + ) -@pytest.mark.skipif(shutil.which("mri_surf2surf") is None, - reason="freesurfer mri_surf2surf not available") +@pytest.mark.skipif( + shutil.which("mri_surf2surf") is None, + reason="freesurfer mri_surf2surf not available", +) def test_surf2surf_matches_freesurfer_identity(): """fsaverage -> fsaverage must be the identity, matching mri_surf2surf. @@ -140,20 +176,24 @@ def test_surf2surf_matches_freesurfer_identity(): if not _have_template(subjects_dir, "fsaverage", hemi): pytest.skip("fsaverage template with sphere.reg not found in SUBJECTS_DIR") - m = get_mri_surf2surf_matrix("fsaverage", hemi, target_subj="fsaverage", - subjects_dir=subjects_dir) + m = get_mri_surf2surf_matrix( + "fsaverage", hemi, target_subj="fsaverage", subjects_dir=subjects_dir + ) rng = np.random.RandomState(0) data = rng.randn(4, m.shape[1]).astype(np.float32) - reference = fs.mri_surf2surf(data, "fsaverage", "fsaverage", hemi, - subjects_dir=subjects_dir) + reference = fs.mri_surf2surf( + data, "fsaverage", "fsaverage", hemi, subjects_dir=subjects_dir + ) got = np.stack([m.dot(data[i]) for i in range(data.shape[0])]) - np.testing.assert_allclose(got, data, atol=1e-4) # identity - np.testing.assert_allclose(got, reference, atol=1e-4) # matches freesurfer + np.testing.assert_allclose(got, data, atol=1e-4) # identity + np.testing.assert_allclose(got, reference, atol=1e-4) # matches freesurfer -@pytest.mark.skipif(shutil.which("mri_surf2surf") is None, - reason="freesurfer mri_surf2surf not available") +@pytest.mark.skipif( + shutil.which("mri_surf2surf") is None, + reason="freesurfer mri_surf2surf not available", +) def test_surf2surf_matches_freesurfer_downsample(): """fsaverage -> fsaverage6 (icosahedral downsample) against mri_surf2surf. @@ -164,16 +204,20 @@ def test_surf2surf_matches_freesurfer_downsample(): """ subjects_dir = os.environ.get("SUBJECTS_DIR") hemi = "lh" - if not (_have_template(subjects_dir, "fsaverage", hemi) - and _have_template(subjects_dir, "fsaverage6", hemi)): + if not ( + _have_template(subjects_dir, "fsaverage", hemi) + and _have_template(subjects_dir, "fsaverage6", hemi) + ): pytest.skip("fsaverage/fsaverage6 templates not found in SUBJECTS_DIR") - m = get_mri_surf2surf_matrix("fsaverage", hemi, target_subj="fsaverage6", - subjects_dir=subjects_dir) + m = get_mri_surf2surf_matrix( + "fsaverage", hemi, target_subj="fsaverage6", subjects_dir=subjects_dir + ) rng = np.random.RandomState(0) data = rng.randn(4, m.shape[1]).astype(np.float32) - reference = fs.mri_surf2surf(data, "fsaverage", "fsaverage6", hemi, - subjects_dir=subjects_dir) + reference = fs.mri_surf2surf( + data, "fsaverage", "fsaverage6", hemi, subjects_dir=subjects_dir + ) got = np.stack([m.dot(data[i]) for i in range(data.shape[0])]) corr = np.corrcoef(reference.ravel(), got.ravel())[0, 1] diff --git a/cortex/tests/test_hcp.py b/cortex/tests/test_hcp.py new file mode 100644 index 00000000..97eb9b85 --- /dev/null +++ b/cortex/tests/test_hcp.py @@ -0,0 +1,225 @@ +import shutil +import subprocess +import tempfile +from pathlib import Path + +import numpy as np +import nibabel as nib +import pytest +from scipy.spatial import ConvexHull + +import cortex +from cortex import hcp +from cortex.hcp import ( + _barycentric_resample_matrix, + cifti_to_surface, + get_cifti_vertex_indices, +) + + +# --------------------------------------------------------------------------- +# Barycentric resampling matrix (pure geometry, no network / wb_command) +# --------------------------------------------------------------------------- + + +def _triangulated_sphere(n, seed): + """n points on the unit sphere plus a covering triangulation (convex hull).""" + rng = np.random.RandomState(seed) + pts = rng.randn(n, 3) + pts /= np.linalg.norm(pts, axis=1, keepdims=True) + tris = ConvexHull(pts).simplices + return pts, tris + + +def _write_sphere_gii(path, pts, tris): + """Write points + triangles to a ``.surf.gii`` file like the HCP spheres.""" + img = nib.gifti.GiftiImage( + darrays=[ + nib.gifti.GiftiDataArray( + np.asarray(pts, dtype=np.float32), intent="NIFTI_INTENT_POINTSET" + ), + nib.gifti.GiftiDataArray( + np.asarray(tris, dtype=np.int32), intent="NIFTI_INTENT_TRIANGLE" + ), + ] + ) + nib.save(img, str(path)) + + +def test_barycentric_identity_when_source_equals_target(tmp_path): + # Same source and target sphere: every target vertex sits on a source + # vertex, so the matrix is the identity. + pts, tris = _triangulated_sphere(60, seed=0) + src = tmp_path / "src.surf.gii" + _write_sphere_gii(src, pts, tris) + m = _barycentric_resample_matrix(src, src) + assert m.shape == (60, 60) + np.testing.assert_allclose(m.toarray(), np.eye(60), atol=1e-9) + + +def test_barycentric_rows_sum_to_one(tmp_path): + src_pts, src_tris = _triangulated_sphere(200, seed=1) + tgt_pts, tgt_tris = _triangulated_sphere(80, seed=2) + src = tmp_path / "src.surf.gii" + tgt = tmp_path / "tgt.surf.gii" + _write_sphere_gii(src, src_pts, src_tris) + _write_sphere_gii(tgt, tgt_pts, tgt_tris) + m = _barycentric_resample_matrix(src, tgt) + assert m.shape == (80, 200) + row_sums = np.asarray(m.sum(axis=1)).ravel() + np.testing.assert_allclose(row_sums, np.ones(80), atol=1e-9) + + +def test_barycentric_preserves_constants(tmp_path): + src_pts, src_tris = _triangulated_sphere(200, seed=3) + tgt_pts, tgt_tris = _triangulated_sphere(70, seed=4) + src = tmp_path / "src.surf.gii" + tgt = tmp_path / "tgt.surf.gii" + _write_sphere_gii(src, src_pts, src_tris) + _write_sphere_gii(tgt, tgt_pts, tgt_tris) + m = _barycentric_resample_matrix(src, tgt) + const = np.full(200, 3.5) + np.testing.assert_allclose(m.dot(const), np.full(70, 3.5), atol=1e-9) + + +def test_barycentric_weights_are_nonnegative(tmp_path): + src_pts, src_tris = _triangulated_sphere(150, seed=5) + tgt_pts, tgt_tris = _triangulated_sphere(60, seed=6) + src = tmp_path / "src.surf.gii" + tgt = tmp_path / "tgt.surf.gii" + _write_sphere_gii(src, src_pts, src_tris) + _write_sphere_gii(tgt, tgt_pts, tgt_tris) + m = _barycentric_resample_matrix(src, tgt) + assert m.data.min() >= 0.0 + # Barycentric interpolation uses at most three source vertices per target. + assert m.getnnz(axis=1).max() <= 3 + + +# --------------------------------------------------------------------------- +# CIFTI -> fs_LR 32k surface expansion +# --------------------------------------------------------------------------- + + +def test_cifti_to_surface_1d_places_values_and_nans_medial_wall(): + left_idx = np.array([0, 2, 5]) + right_idx = np.array([1, 3]) + data = np.array([10.0, 20.0, 30.0, 40.0, 50.0]) # 3 left then 2 right + full = cifti_to_surface(data, left_idx, right_idx) + + assert full.shape == (hcp.N_VERTICES_FS_LR_32K,) + off = hcp.N_VERTICES_FS_LR_32K_HEM + np.testing.assert_array_equal(full[left_idx], [10.0, 20.0, 30.0]) + np.testing.assert_array_equal(full[off + right_idx], [40.0, 50.0]) + # Everything else is medial wall -> NaN. + n_valid = np.sum(~np.isnan(full)) + assert n_valid == 5 + + +def test_cifti_to_surface_2d_matches_1d_per_row(): + left_idx = np.array([4, 1]) + right_idx = np.array([7]) + data = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) # 2 left, 1 right per row + full = cifti_to_surface(data, left_idx, right_idx) + assert full.shape == (2, hcp.N_VERTICES_FS_LR_32K) + off = hcp.N_VERTICES_FS_LR_32K_HEM + np.testing.assert_array_equal(full[:, left_idx], [[1.0, 2.0], [4.0, 5.0]]) + np.testing.assert_array_equal(full[:, off + right_idx], [[3.0], [6.0]]) + + +def test_get_cifti_vertex_indices_reads_brain_model_axis(): + left_idx = [0, 2, 5, 9] + right_idx = [1, 3, 8] + bm = nib.cifti2.BrainModelAxis.from_surface( + left_idx, hcp.N_VERTICES_FS_LR_32K_HEM, "CIFTI_STRUCTURE_CORTEX_LEFT" + ) + nib.cifti2.BrainModelAxis.from_surface( + right_idx, hcp.N_VERTICES_FS_LR_32K_HEM, "CIFTI_STRUCTURE_CORTEX_RIGHT" + ) + scalar = nib.cifti2.ScalarAxis(["map"]) + img = nib.cifti2.Cifti2Image( + np.zeros((1, len(left_idx) + len(right_idx))), header=(scalar, bm) + ) + + got_left, got_right = get_cifti_vertex_indices(img) + np.testing.assert_array_equal(got_left, left_idx) + np.testing.assert_array_equal(got_right, right_idx) + + +# --------------------------------------------------------------------------- +# Validation against wb_command (requires Connectome Workbench + network) +# --------------------------------------------------------------------------- + + +def _wb_metric_resample(data, src_sphere, tgt_sphere): + """Reference resample via ``wb_command -metric-resample ... BARYCENTRIC``.""" + with tempfile.TemporaryDirectory() as td: + in_p = Path(td) / "in.func.gii" + out_p = Path(td) / "out.func.gii" + img = nib.gifti.GiftiImage( + darrays=[ + nib.gifti.GiftiDataArray( + np.asarray(data, dtype=np.float32), + intent="NIFTI_INTENT_NORMAL", + datatype="NIFTI_TYPE_FLOAT32", + ) + ] + ) + nib.save(img, str(in_p)) + subprocess.run( + [ + "wb_command", + "-metric-resample", + str(in_p), + str(src_sphere), + str(tgt_sphere), + "BARYCENTRIC", + str(out_p), + ], + check=True, + capture_output=True, + text=True, + ) + return np.asarray(nib.load(str(out_p)).darrays[0].data) + + +@pytest.mark.skipif( + shutil.which("wb_command") is None, + reason="Connectome Workbench (wb_command) not installed", +) +def test_barycentric_matches_wb_command(): + # Uses the real HCP fs_LR 32k -> fsaverage5 spheres (downloaded on demand), + # so it also exercises ensure_sphere_files. Skipped if the spheres cannot + # be fetched. + cache_dir = Path(tempfile.gettempdir()) / "pycortex_hcp_test_spheres" + try: + src = hcp.ensure_sphere_files("fs_LR_32k", "L", cache_dir=cache_dir) + tgt = hcp.ensure_sphere_files("fsaverage5", "L", cache_dir=cache_dir) + except RuntimeError as exc: # network unavailable + pytest.skip(f"could not download sphere files: {exc}") + + m = _barycentric_resample_matrix(src, tgt) + rng = np.random.RandomState(0) + data = rng.randn(hcp.N_VERTICES_FS_LR_32K_HEM).astype(np.float32) + + ours = m.dot(data) + reference = _wb_metric_resample(data, src, tgt) + assert ours.shape == reference.shape + corr = np.corrcoef(ours, reference)[0, 1] + assert corr > 0.999 + # Spread of the reference is O(1); the discrepancy should be tiny. + assert np.abs(ours - reference).max() < 0.05 * reference.std() + 1e-4 + + +# --------------------------------------------------------------------------- +# Native-space rendering integration (requires the 32k_fs_LR subject) +# --------------------------------------------------------------------------- + + +@pytest.mark.skipif( + "32k_fs_LR" not in cortex.db.subjects, reason="32k_fs_LR subject not in filestore" +) +def test_native_vertex_matches_subject(): + # The subject carries the full 64984-vertex fs_LR 32k surface, so a Vertex + # built from cifti_to_surface output must line up with it. + data = np.zeros(hcp.N_VERTICES_FS_LR_32K) + v = cortex.Vertex(data, "32k_fs_LR") + assert v.data.shape[-1] == hcp.N_VERTICES_FS_LR_32K diff --git a/cortex/utils.py b/cortex/utils.py index 6ed7ab52..1613d8fc 100644 --- a/cortex/utils.py +++ b/cortex/utils.py @@ -1,5 +1,5 @@ -"""Contain utility functions -""" +"""Contain utility functions""" + import binascii import copy from importlib import import_module @@ -12,7 +12,19 @@ import urllib.request import warnings -from typing import Any, Callable, Generic, Optional, TypeVar, TYPE_CHECKING, Union, cast, overload, Literal +from typing import ( + Any, + Callable, + Generic, + Optional, + TypeVar, + TYPE_CHECKING, + Union, + cast, + overload, + Literal, +) + if sys.version_info < (3, 10): from typing_extensions import ParamSpec else: @@ -33,47 +45,67 @@ # register_cmap is deprecated in matplotlib > 3.7.0 and replaced by colormaps.register try: from matplotlib import colormaps as cm + def register_cmap(cmap): return cm.register(cmap) except ImportError: from matplotlib.cm import register_cmap -P = ParamSpec('P') -T = TypeVar('T') +P = ParamSpec("P") +T = TypeVar("T") + class DocLoader(Generic[P, T]): - def __init__(self, func, mod, package, actual_func: Optional[Callable[P, T]] = None): - self._load: Callable[[], Callable[P, T]] = lambda: getattr(import_module(mod, package), func) - self._actual_func = actual_func # stored only to resolve generic types during type checking + def __init__( + self, func, mod, package, actual_func: Optional[Callable[P, T]] = None + ): + self._load: Callable[[], Callable[P, T]] = lambda: getattr( + import_module(mod, package), func + ) + self._actual_func = ( + actual_func # stored only to resolve generic types during type checking + ) def __call__(self, *args: P.args, **kwargs: P.kwargs) -> T: return self._load()(*args, **kwargs) @overload - def __getattribute__(self, name: Literal['_load']) -> Callable[P, T]: ... + def __getattribute__(self, name: Literal["_load"]) -> Callable[P, T]: ... @overload def __getattribute__(self, name: str) -> Any: ... - def __getattribute__(self, name: Union[Literal['_load'], str]) -> Union[Any, Callable[P, T]]: + def __getattribute__( + self, name: Union[Literal["_load"], str] + ) -> Union[Any, Callable[P, T]]: if name != "_load": return getattr(self._load(), name) else: return cast(Callable[P, T], object.__getattribute__(self, name)) + if TYPE_CHECKING: from cortex.mapper import get_mapper as _get_mapper else: _get_mapper = None get_mapper = DocLoader("get_mapper", ".mapper", "cortex", actual_func=_get_mapper) + def get_roipack(*args, **kwargs): - warnings.warn('Please use db.get_overlay instead', DeprecationWarning) + warnings.warn("Please use db.get_overlay instead", DeprecationWarning) return db.get_overlay(*args, **kwargs) -def get_ctmpack(subject, types=("inflated",), method="raw", level=0, recache=False, - decimate=False, external_svg=None, - overlays_available=None): + +def get_ctmpack( + subject, + types=("inflated",), + method="raw", + level=0, + recache=False, + decimate=False, + external_svg=None, + overlays_available=None, +): """Creates ctm file for the specified input arguments. This is a cached file that specifies (1) the surfaces between which @@ -108,12 +140,10 @@ def get_ctmpack(subject, types=("inflated",), method="raw", level=0, recache=Fal ------- ctmfile : """ - lvlstr = ("%dd" if decimate else "%d")%level + lvlstr = ("%dd" if decimate else "%d") % level # Generates different cache files for each combination of disp_layers - ctmcache = "%s_[{types}]_{method}_{level}_v3.json"%subject - ctmcache = ctmcache.format(types=','.join(types), - method=method, - level=lvlstr) + ctmcache = "%s_[{types}]_{method}_{level}_v3.json" % subject + ctmcache = ctmcache.format(types=",".join(types), method=method, level=lvlstr) ctmfile = os.path.join(db.get_cache(subject), ctmcache) if os.path.exists(ctmfile) and not recache: @@ -121,20 +151,23 @@ def get_ctmpack(subject, types=("inflated",), method="raw", level=0, recache=Fal print("Generating new ctm file...") from . import brainctm - ptmap = brainctm.make_pack(ctmfile, - subject, - types=types, - method=method, - level=level, - decimate=decimate, - external_svg=external_svg, - overlays_available=overlays_available) + + ptmap = brainctm.make_pack( + ctmfile, + subject, + types=types, + method=method, + level=level, + decimate=decimate, + external_svg=external_svg, + overlays_available=overlays_available, + ) return ctmfile def get_ctmmap(subject, **kwargs): """Return a mapping from the vertices in the CTM surface to the vertices - in the freesurfer surface. + in the freesurfer surface. The mapping is a numpy array, such that `ctm2fs_left[i] = j` means that the i-th vertex in the CTM surface corresponds to the j-th vertex in the freesurfer surface. @@ -162,8 +195,9 @@ def get_ctmmap(subject, **kwargs): from scipy.spatial import cKDTree from . import brainctm + jsfile = get_ctmpack(subject, **kwargs) - ctmfile = os.path.splitext(jsfile)[0]+".ctm" + ctmfile = os.path.splitext(jsfile)[0] + ".ctm" # Load freesurfer surfaces try: @@ -209,6 +243,7 @@ def get_ctm2webgl_map(subject, **kwargs): maximum length of 65535. """ from . import brainctm + # Load CTM surfaces jsonfile = get_ctmpack(subject, **kwargs) ctmfile = os.path.splitext(jsonfile)[0] + ".ctm" @@ -290,7 +325,7 @@ def get_fs2webgl_map(subject, **kwargs): return fs2webgl_left, fs2webgl_right -def get_cortical_mask(subject, xfmname, type='nearest'): +def get_cortical_mask(subject, xfmname, type="nearest"): """Gets the cortical mask for a particular transform Parameters @@ -301,12 +336,12 @@ def get_cortical_mask(subject, xfmname, type='nearest'): Transform name type : str Mask type, one of {"cortical", "thin", "thick", "nearest", "line_nearest"}. - - 'cortical' includes voxels contained within the cortical ribbon, - between the freesurfer-estimated white matter and pial surfaces. - - 'thin' includes voxels that are < 2mm away from the fiducial surface. + - 'cortical' includes voxels contained within the cortical ribbon, + between the freesurfer-estimated white matter and pial surfaces. + - 'thin' includes voxels that are < 2mm away from the fiducial surface. - 'thick' includes voxels that are < 8mm away from the fiducial surface. - 'nearest' includes only the voxels overlapping the fiducial surface. - - 'line_nearest' includes all voxels that have any part within the cortical + - 'line_nearest' includes all voxels that have any part within the cortical ribbon. Returns @@ -316,13 +351,13 @@ def get_cortical_mask(subject, xfmname, type='nearest'): Notes ----- - "nearest" is a conservative "cortical" mask, while "line_nearest" is a liberal + "nearest" is a conservative "cortical" mask, while "line_nearest" is a liberal "cortical" mask. """ - if type == 'cortical': + if type == "cortical": ppts, polys = db.get_surf(subject, "pia", merge=True, nudge=False) wpts, polys = db.get_surf(subject, "wm", merge=True, nudge=False) - thickness = np.sqrt(((ppts - wpts)**2).sum(1)) + thickness = np.sqrt(((ppts - wpts) ** 2).sum(1)) dist, idx = get_vox_dist(subject, xfmname) cortex = np.zeros(dist.shape, dtype=bool) @@ -331,9 +366,9 @@ def get_cortical_mask(subject, xfmname, type='nearest'): mask = idx == vert cortex[mask] = dist[mask] <= thickness[vert] if i % 100 == 0: - print("%0.3f%%"%(i/float(len(verts)) * 100)) + print("%0.3f%%" % (i / float(len(verts)) * 100)) return cortex - elif type in ('thick', 'thin'): + elif type in ("thick", "thin"): dist, idx = get_vox_dist(subject, xfmname) return dist < dict(thick=8, thin=2)[type] else: @@ -379,8 +414,8 @@ def get_vox_dist(subject, xfmname, surface="fiducial", max_dist=np.inf): return dist.T, argdist.T -def get_hemi_masks(subject, xfmname, type='nearest'): - '''Returns a binary mask of the left and right hemisphere +def get_hemi_masks(subject, xfmname, type="nearest"): + """Returns a binary mask of the left and right hemisphere surface voxels for the given subject. Parameters @@ -394,11 +429,13 @@ def get_hemi_masks(subject, xfmname, type='nearest'): Returns ------- - ''' + """ return get_mapper(subject, xfmname, type=type).hemimasks -def add_roi(data, name="new_roi", open_inkscape=True, add_path=True, - overlay_file=None, **kwargs): + +def add_roi( + data, name="new_roi", open_inkscape=True, add_path=True, overlay_file=None, **kwargs +): """Add new flatmap image to the ROI file for a subject. (The subject is specified in creation of the data object) @@ -440,14 +477,16 @@ def add_roi(data, name="new_roi", open_inkscape=True, add_path=True, svg = db.get_overlay(dv.subject, overlay_file=overlay_file) fp = io.BytesIO() - quickflat.make_png(fp, dv, height=1024, with_rois=False, with_labels=False, **kwargs) + quickflat.make_png( + fp, dv, height=1024, with_rois=False, with_labels=False, **kwargs + ) fp.seek(0) - svg.rois.add_shape(name, binascii.b2a_base64(fp.read()).decode('utf-8'), add_path) + svg.rois.add_shape(name, binascii.b2a_base64(fp.read()).decode("utf-8"), add_path) if open_inkscape: - inkscape_cmd = config.get('dependency_paths', 'inkscape') - if LooseVersion(INKSCAPE_VERSION) < LooseVersion('1.0'): - cmd = [inkscape_cmd, '-f', svg.svgfile] + inkscape_cmd = config.get("dependency_paths", "inkscape") + if LooseVersion(INKSCAPE_VERSION) < LooseVersion("1.0"): + cmd = [inkscape_cmd, "-f", svg.svgfile] else: cmd = [inkscape_cmd, svg.svgfile] return sp.call(cmd) @@ -509,8 +548,8 @@ def get_roi_verts(subject, roi=None, mask=False, overlay_file=None): for name in roi: roi_idx = np.intersect1d(svg.rois.get_mask(name), goodpts) - # Now we want to include also the vertices that were removed from the flat - # surface that is, for every vertex in roi_idx we want to add the pts that are + # Now we want to include also the vertices that were removed from the flat + # surface that is, for every vertex in roi_idx we want to add the pts that are # not in goodpts but that are in pts_full # to do that, we need to find the neighboring indices from polys_full extra_idx = set() @@ -568,7 +607,7 @@ def get_roi_surf(subject, surf_type, roi, overlay_file=None): return pts[vert_idx], np.array(reindexed_polys) -def get_roi_mask(subject, xfmname, roi=None, projection='nearest'): +def get_roi_mask(subject, xfmname, roi=None, projection="nearest"): """Return a mask for the given ROI(s) Deprecated - use get_roi_masks() @@ -588,15 +627,15 @@ def get_roi_mask(subject, xfmname, roi=None, projection='nearest'): output : dict Dict of ROIs and their masks """ - warnings.warn('Deprecated! Use get_roi_masks') + warnings.warn("Deprecated! Use get_roi_masks") mapper = get_mapper(subject, xfmname, type=projection) rois = get_roi_verts(subject, roi=roi, mask=True) output = dict() for name, verts in list(rois.items()): # This is broken; unclear when/if backward mappers ever worked this way. - #left, right = mapper.backwards(vert_mask) - #output[name] = left + right + # left, right = mapper.backwards(vert_mask) + # output[name] = left + right output[name] = mapper.backwards(verts.astype(float)) # Threshold? return output @@ -643,7 +682,7 @@ def get_aseg_mask(subject, aseg_name, xfmname=None, order=1, threshold=None, **k See also get_anat(subject, type='aseg') """ - from .freesurfer import fs_aseg_dict + aseg = db.get_anat(subject, type="aseg").get_fdata().T if not isinstance(aseg_name, (list, tuple)): @@ -652,24 +691,36 @@ def get_aseg_mask(subject, aseg_name, xfmname=None, order=1, threshold=None, **k mask = np.zeros(aseg.shape) for name in aseg_name: if name in fs_aseg_dict: - tmp = aseg==fs_aseg_dict[name] + tmp = aseg == fs_aseg_dict[name] else: # Combine all masks containing `name` (e.g. all masks with 'cerebellum' in the name) keys = [k for k in fs_aseg_dict.keys() if name.lower() in k.lower()] if len(keys) == 0: - raise ValueError('Unknown aseg_name!') - tmp = np.any(np.array([aseg==fs_aseg_dict[k] for k in keys]), axis=0) + raise ValueError("Unknown aseg_name!") + tmp = np.any(np.array([aseg == fs_aseg_dict[k] for k in keys]), axis=0) mask = np.logical_or(mask, tmp) if xfmname is not None: - mask = anat2epispace(mask.astype(float), subject, xfmname, order=order, **kwargs) + mask = anat2epispace( + mask.astype(float), subject, xfmname, order=order, **kwargs + ) if threshold is not None: mask = mask > threshold return mask -def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_lr=False, - allow_overlap=False, fail_for_missing_rois=True, exclude_empty_rois=False, - threshold=None, return_dict=True, overlay_file=None): +def get_roi_masks( + subject, + xfmname, + roi_list=None, + gm_sampler="cortical", + split_lr=False, + allow_overlap=False, + fail_for_missing_rois=True, + exclude_empty_rois=False, + threshold=None, + return_dict=True, + overlay_file=None, +): """Return a dictionary of roi masks This function returns a single 3D array with a separate numerical index for each ROI, @@ -748,13 +799,17 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_ 'thin' as your `gm_sampler`. """ # Convert mapper names to pycortex sampler types - mapper_dict = {'cortical-conservative':'nearest', - 'cortical-liberal':'line_nearest'} + mapper_dict = { + "cortical-conservative": "nearest", + "cortical-liberal": "line_nearest", + } # Method use_mapper = gm_sampler in mapper_dict - use_cortex_mask = (gm_sampler in ('cortical', 'thick', 'thin')) or not isinstance(gm_sampler, str) + use_cortex_mask = (gm_sampler in ("cortical", "thick", "thin")) or not isinstance( + gm_sampler, str + ) if not (use_mapper or use_cortex_mask): - raise ValueError('Unknown gray matter sampler (gm_sampler)!') + raise ValueError("Unknown gray matter sampler (gm_sampler)!") # Initialize roi_voxels = {} pct_coverage = {} @@ -773,18 +828,28 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_ roi_verts = get_roi_verts(subject, mask=use_mapper, overlay_file=overlay_file) roi_list = list(roi_verts.keys()) else: - tmp_list = [r for r in roi_list if not r=='Cortex'] + tmp_list = [r for r in roi_list if not r == "Cortex"] try: - roi_verts = get_roi_verts(subject, roi=tmp_list, mask=use_mapper, overlay_file=overlay_file) + roi_verts = get_roi_verts( + subject, roi=tmp_list, mask=use_mapper, overlay_file=overlay_file + ) except KeyError as key: if fail_for_missing_rois: - raise KeyError("Requested ROI {} not found in overlays.svg!".format(key)) + raise KeyError( + "Requested ROI {} not found in overlays.svg!".format(key) + ) else: - roi_verts = get_roi_verts(subject, roi=None, mask=use_mapper, overlay_file=overlay_file) - missing = [r for r in roi_list if not r in roi_verts.keys()+['Cortex']] - roi_verts = dict((roi, verts) for roi, verts in roi_verts.items() if roi in roi_list) - roi_list = list(set(roi_list)-set(missing)) - print('Requested ROI(s) {} not found in overlays.svg!'.format(missing)) + roi_verts = get_roi_verts( + subject, roi=None, mask=use_mapper, overlay_file=overlay_file + ) + missing = [ + r for r in roi_list if r not in roi_verts.keys() + ["Cortex"] + ] + roi_verts = dict( + (roi, verts) for roi, verts in roi_verts.items() if roi in roi_list + ) + roi_list = list(set(roi_list) - set(missing)) + print("Requested ROI(s) {} not found in overlays.svg!".format(missing)) # Get (a) indices for nearest vertex to each voxel # and (b) distance from each voxel to nearest vertex in fiducial surface if (use_cortex_mask or split_lr) or (not return_dict): @@ -799,7 +864,7 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_ # Loop over ROIs to map vertices to volume, using mapper or cortex mask + vertex indices for roi in roi_list: if roi not in roi_verts: - if not roi=='Cortex': + if not roi == "Cortex": print("ROI {} not found...".format(roi)) continue if use_mapper: @@ -808,7 +873,9 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_ if threshold is not None: roi_voxels[roi] = roi_voxels[roi] > threshold # Check for partial / empty rois: - vert_in_scan = np.hstack([np.array((m>0).sum(1)).flatten() for m in mapper.masks]) + vert_in_scan = np.hstack( + [np.array((m > 0).sum(1)).flatten() for m in mapper.masks] + ) vert_in_scan = vert_in_scan[roi_verts[roi]] elif use_cortex_mask: vox_in_roi = np.isin(vox_idx, roi_verts[roi]) @@ -820,58 +887,69 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_ # Compute ROI coverage pct_coverage[roi] = vert_in_scan.mean() * 100 if use_mapper: - print("Found %0.2f%% of %s"%(pct_coverage[roi], roi)) + print("Found %0.2f%% of %s" % (pct_coverage[roi], roi)) # Create cortex mask all_mask = np.array(list(roi_voxels.values())).sum(0) - if 'Cortex' in roi_list: + if "Cortex" in roi_list: if use_mapper: # cortex_mask isn't defined / exactly definable if you're using a mapper - print("Cortex roi not included b/c currently not compatible with your selection for gm_sampler") - _ = roi_list.pop(roi_list.index('Cortex')) + print( + "Cortex roi not included b/c currently not compatible with your selection for gm_sampler" + ) + _ = roi_list.pop(roi_list.index("Cortex")) else: - roi_voxels['Cortex'] = (all_mask==0) & cortex_mask + roi_voxels["Cortex"] = (all_mask == 0) & cortex_mask # Optionally cull voxels assigned to > 1 ROI due to partly overlapping ROI splines # in inkscape overlays.svg file: if not allow_overlap: - print('Cutting {} overlapping voxels (should be < ~50)'.format(np.sum(all_mask > 1))) + print( + "Cutting {} overlapping voxels (should be < ~50)".format( + np.sum(all_mask > 1) + ) + ) for roi in roi_list: roi_voxels[roi][all_mask > 1] = False # Split left / right hemispheres if desired if split_lr: # Use the fiducial surface because we need to have all vertices - left_verts, _ = db.get_surf(subject, "fiducial", merge=False, nudge=True) + left_verts, _ = db.get_surf(subject, "fiducial", merge=False, nudge=True) left_mask = vox_idx < len(np.unique(left_verts[1])) right_mask = np.logical_not(left_mask) roi_voxels_lr = {} for roi in roi_list: # roi_voxels may contain float values if using a mapper, therefore we need # to manually set the voxels in the other hemisphere to False. Then we let - # numpy do the conversion False -> 0. - roi_voxels_lr[roi + '_L'] = copy.copy(roi_voxels[roi]) - roi_voxels_lr[roi + '_L'][right_mask] = False - roi_voxels_lr[roi + '_R'] = copy.copy(roi_voxels[roi]) - roi_voxels_lr[roi + '_R'][left_mask] = False + # numpy do the conversion False -> 0. + roi_voxels_lr[roi + "_L"] = copy.copy(roi_voxels[roi]) + roi_voxels_lr[roi + "_L"][right_mask] = False + roi_voxels_lr[roi + "_R"] = copy.copy(roi_voxels[roi]) + roi_voxels_lr[roi + "_R"][left_mask] = False output = roi_voxels_lr else: output = roi_voxels # Check percent coverage / optionally cull empty ROIs - for roi in set(roi_list)-set(['Cortex']): + for roi in set(roi_list) - set(["Cortex"]): if pct_coverage[roi] < 100: # if not np.any(mask) : reject ROI - if pct_coverage[roi]==0: - warnings.warn('ROI %s is entirely missing from your scan protocol!'%(roi)) + if pct_coverage[roi] == 0: + warnings.warn( + "ROI %s is entirely missing from your scan protocol!" % (roi) + ) if exclude_empty_rois: if split_lr: - _ = output.pop(roi+'_L') - _ = output.pop(roi+'_R') + _ = output.pop(roi + "_L") + _ = output.pop(roi + "_R") else: _ = output.pop(roi) else: # I think this is the only one for which this works correctly... - if gm_sampler=='cortical-conservative': - warnings.warn('ROI %s is only %0.2f%% contained in your scan protocol!'%(roi, pct_coverage[roi])) + if gm_sampler == "cortical-conservative": + warnings.warn( + "ROI %s is only %0.2f%% contained in your scan protocol!" + % (roi, pct_coverage[roi]) + ) # Support alternative outputs for backward compatibility if return_dict: @@ -885,6 +963,7 @@ def get_roi_masks(subject, xfmname, roi_list=None, gm_sampler='cortical', split_ idx_vol[left_mask] *= -1 return idx_vol, idx_labels + def get_dropout(subject: str, xfmname: str, power: float = 20): """Create a dropout Volume showing where EPI signal is very low. @@ -909,13 +988,15 @@ def get_dropout(subject: str, xfmname: str, power: float = 20): if rawdata.ndim > 3: rawdata = rawdata.mean(0) - rawdata[rawdata==0] = np.mean(rawdata[rawdata!=0]) + rawdata[rawdata == 0] = np.mean(rawdata[rawdata != 0]) normdata = (rawdata - rawdata.min()) / (rawdata.max() - rawdata.min()) normdata = (1 - normdata) ** power from .dataset import Volume + return Volume(normdata, subject, xfmname) + def make_movie(stim, outfile, fps=15, size="640x480"): """Makes an .ogv movie @@ -939,10 +1020,12 @@ def make_movie(stim, outfile, fps=15, size="640x480"): """ import shlex import subprocess as sp + cmd = "ffmpeg -r {fps} -i {infile} -b 4800k -g 30 -s {size} -vcodec libtheora {outfile}.ogv" fcmd = cmd.format(infile=stim, size=size, fps=fps, outfile=outfile) sp.call(shlex.split(fcmd)) + def vertex_to_voxel(subject): # Am I deprecated in favor of mappers??? Maybe? """ Parameters @@ -975,28 +1058,30 @@ def vertex_to_voxel(subject): # Am I deprecated in favor of mappers??? Maybe? def _set_edge_distance_graph_attribute(graph, pts, polys): - ''' + """ adds the attribute 'edge distance' to a graph - ''' + """ import networkx as nx l2_distance = lambda v1, v2: np.linalg.norm(pts[v1] - pts[v2]) - heuristic = l2_distance # A* heuristic + heuristic = l2_distance # A* heuristic - if not nx.get_edge_attributes(graph, 'distance'): # Add edge distances as an attribute to this graph if it isn't there + if not nx.get_edge_attributes( + graph, "distance" + ): # Add edge distances as an attribute to this graph if it isn't there edge_distances = dict() - for x,y,z in polys: - edge_distances[(x,y)] = heuristic(x,y) - edge_distances[(y,x)] = heuristic(y,x) - edge_distances[(y,z)] = heuristic(y,z) - edge_distances[(z,y)] = heuristic(z,y) - edge_distances[(x,z)] = heuristic(x,z) - edge_distances[(z,x)] = heuristic(z,x) - nx.set_edge_attributes(graph, edge_distances, name='distance') + for x, y, z in polys: + edge_distances[(x, y)] = heuristic(x, y) + edge_distances[(y, x)] = heuristic(y, x) + edge_distances[(y, z)] = heuristic(y, z) + edge_distances[(z, y)] = heuristic(z, y) + edge_distances[(x, z)] = heuristic(x, z) + edge_distances[(z, x)] = heuristic(z, x) + nx.set_edge_attributes(graph, edge_distances, name="distance") def get_shared_voxels(subject, xfmname, hemi="both", merge=True, use_astar=True): - '''Return voxels that are shared by multiple vertices, and for each such voxel, + """Return voxels that are shared by multiple vertices, and for each such voxel, also returns the mutually farthest pair of vertices mapping to the voxel Parameters ---------- @@ -1017,18 +1102,21 @@ def get_shared_voxels(subject, xfmname, hemi="both", merge=True, use_astar=True) vox_vert_array: np.array, array of dimensions # voxels X 3, columns being: (vox_idx, farthest_pair[0], farthest_pair[1]) - ''' + """ import networkx as nx from scipy.sparse import find as sparse_find - Lmask, Rmask = get_mapper(subject, xfmname).masks # Get masks for left and right hemisphere - if hemi == 'both': - hemispheres = ['lh', 'rh'] + + Lmask, Rmask = get_mapper( + subject, xfmname + ).masks # Get masks for left and right hemisphere + if hemi == "both": + hemispheres = ["lh", "rh"] else: hemispheres = [hemi] out = [] for hem in hemispheres: - if hem == 'lh': + if hem == "lh": mask = Lmask else: mask = Rmask @@ -1036,8 +1124,10 @@ def get_shared_voxels(subject, xfmname, hemi="both", merge=True, use_astar=True) all_voxels = mask.tolil().transpose().rows # Map from voxels to verts vert_to_vox_map = dict(zip(*(sparse_find(mask)[:2]))) # From verts to vox - pts_fid, polys_fid = db.get_surf(subject, 'fiducial', hem) # Get the fiducial surface - surf = Surface(pts_fid, polys_fid) #Get the fiducial surface + pts_fid, polys_fid = db.get_surf( + subject, "fiducial", hem + ) # Get the fiducial surface + surf = Surface(pts_fid, polys_fid) # Get the fiducial surface graph = surf.graph _set_edge_distance_graph_attribute(graph, pts_fid, polys_fid) @@ -1046,32 +1136,44 @@ def get_shared_voxels(subject, xfmname, hemi="both", merge=True, use_astar=True) heuristic = l2_distance # A* heuristic if use_astar: - shortest_path = lambda a, b: nx.astar_path(graph, a, b, heuristic=heuristic, weight='distance') # Find approximate shortest paths using A* search + shortest_path = lambda a, b: nx.astar_path( + graph, a, b, heuristic=heuristic, weight="distance" + ) # Find approximate shortest paths using A* search else: - shortest_path = surf.geodesic_path # Find shortest paths using geodesic distances + shortest_path = ( + surf.geodesic_path + ) # Find shortest paths using geodesic distances vox_vert_list = [] for vox_idx, vox in enumerate(all_voxels): if len(vox) > 1: # If the voxel maps to multiple vertices vox = np.array(vox).astype(int) - for v1 in range(vox.size-1): + for v1 in range(vox.size - 1): vert1 = vox[v1] if vert1 in vert_to_vox_map: # If the vertex is a valid vertex - for v2 in range(v1+1, vox.size): + for v2 in range(v1 + 1, vox.size): vert2 = vox[v2] - if vert2 in vert_to_vox_map: # If the vertex is a valid vertex + if ( + vert2 in vert_to_vox_map + ): # If the vertex is a valid vertex path = shortest_path(vert1, vert2) # Test whether any vertex in path goes out of the voxel - stays_in_voxel = all([(v in vert_to_vox_map) and (vert_to_vox_map[v] == vox_idx) for v in path]) + stays_in_voxel = all( + [ + (v in vert_to_vox_map) + and (vert_to_vox_map[v] == vox_idx) + for v in path + ] + ) if not stays_in_voxel: vox_vert_list.append([vox_idx, vert1, vert2]) - tmp = np.array(vox_vert_list) + tmp = np.array(vox_vert_list) # Add offset for right hem voxels - if hem=='rh': + if hem == "rh": tmp[:, 1:3] += Lmask.shape[0] out.append(tmp) - if hemi in ('lh', 'rh'): + if hemi in ("lh", "rh"): return out[0] else: if merge: @@ -1096,13 +1198,18 @@ def load_sparse_array(fname, varname): conventions, so cannot be used to load arbitrary sparse arrays. """ import scipy.sparse + with h5py.File(fname) as hf: - data = (hf['%s_data'%varname], hf['%s_indices'%varname], hf['%s_indptr'%varname]) - sparsemat = scipy.sparse.csr_matrix(data, shape=hf['%s_shape'%varname]) + data = ( + hf["%s_data" % varname], + hf["%s_indices" % varname], + hf["%s_indptr" % varname], + ) + sparsemat = scipy.sparse.csr_matrix(data, shape=hf["%s_shape" % varname]) return sparsemat -def save_sparse_array(fname, data, varname, mode='a'): +def save_sparse_array(fname, data, varname, mode="a"): """Save a numpy sparse array to an hdf file Results in relatively smaller file size than numpy.savez @@ -1119,19 +1226,20 @@ def save_sparse_array(fname, data, varname, mode='a'): write / append mode set, one of ['w','a'] (passed to h5py.File()) """ import scipy.sparse + if not isinstance(data, scipy.sparse.csr.csr_matrix): data_ = scipy.sparse.csr_matrix(data) else: data_ = data with h5py.File(fname, mode=mode) as hf: # Save indices - hf.create_dataset(varname + '_indices', data=data_.indices, compression='gzip') + hf.create_dataset(varname + "_indices", data=data_.indices, compression="gzip") # Save data - hf.create_dataset(varname + '_data', data=data_.data, compression='gzip') + hf.create_dataset(varname + "_data", data=data_.data, compression="gzip") # Save indptr - hf.create_dataset(varname + '_indptr', data=data_.indptr, compression='gzip') + hf.create_dataset(varname + "_indptr", data=data_.indptr, compression="gzip") # Save shape - hf.create_dataset(varname + '_shape', data=data_.shape, compression='gzip') + hf.create_dataset(varname + "_shape", data=data_.shape, compression="gzip") def get_cmap(name): @@ -1149,10 +1257,11 @@ def get_cmap(name): """ import matplotlib.pyplot as plt from matplotlib import colors + # unknown colormap, test whether it's in pycortex colormaps - cmapdir = config.get('webgl', 'colormaps') + cmapdir = config.get("webgl", "colormaps") colormaps = os.listdir(cmapdir) - colormaps = sorted([c for c in colormaps if '.png' in c]) + colormaps = sorted([c for c in colormaps if ".png" in c]) colormaps = dict((c[:-4], os.path.join(cmapdir, c)) for c in colormaps) if name in colormaps: I = plt.imread(colormaps[name]) @@ -1165,15 +1274,16 @@ def get_cmap(name): try: cmap = plt.cm.get_cmap(name) except: - raise Exception('Unkown color map!') + raise Exception("Unkown color map!") return cmap + def add_cmap(cmap, name, cmapdir=None): """Add a colormap to pycortex. This stores a matplotlib colormap in the pycortex filestore, such that it can - be used in the webgl viewer in pycortex. See - https://matplotlib.org/stable/users/explain/colors/colormap-manipulation.html + be used in the webgl viewer in pycortex. See + https://matplotlib.org/stable/users/explain/colors/colormap-manipulation.html for more information about how to generate colormaps in matplotlib. Parameters @@ -1182,8 +1292,8 @@ def add_cmap(cmap, name, cmapdir=None): Color map to be saved name : str Name for colormap, e.g. 'jet', 'blue_to_yellow', etc. The name will be used - to generate a filename for the colormap stored in the pycortex store, - so avoid illegal characters for a filename. This name will also be used to + to generate a filename for the colormap stored in the pycortex store, + so avoid illegal characters for a filename. This name will also be used to specify this colormap in future calls to `cortex.quickflat.make_figure()` or `cortex.webgl.show()`. """ @@ -1199,8 +1309,9 @@ def add_cmap(cmap, name, cmapdir=None): plt.imsave(os.path.join(cmapdir, name), cmap_im, format="png") -def download_subject(subject_id='fsaverage', url=None, pycortex_store=None, - download_again=False): +def download_subject( + subject_id="fsaverage", url=None, pycortex_store=None, download_again=False +): """Download subjects to pycortex store Parameters @@ -1224,15 +1335,20 @@ def download_subject(subject_id='fsaverage', url=None, pycortex_store=None, warnings.warn( "{} is already present in the database. " "Set download_again to True if you wish to download " - "the subject again.".format(subject_id)) + "the subject again.".format(subject_id) + ) return # Map codes to URLs; more coming eventually id_to_url = dict( - fsaverage='https://ndownloader.figshare.com/files/17827577?private_link=4871247dce31e188e758', + fsaverage="https://ndownloader.figshare.com/files/17827577?private_link=4871247dce31e188e758", + # HCP fs_LR 32k pycortex subject (derived from HCP S1200 group-average + # Open Access surfaces; see cortex.hcp). The tarball contains a + # top-level ``32k_fs_LR/`` directory. + **{"32k_fs_LR": "https://ndownloader.figshare.com/files/66442049"}, ) if url is None: if subject_id not in id_to_url: - raise ValueError('Unknown subject_id!') + raise ValueError("Unknown subject_id!") url = id_to_url[subject_id] # Setup pycortex store location if pycortex_store is None: @@ -1242,12 +1358,11 @@ def download_subject(subject_id='fsaverage', url=None, pycortex_store=None, # Download to temp dir print("Downloading from: {}".format(url)) with tempfile.TemporaryDirectory() as tmp_dir: - print('Downloading subject {} to {}'.format(subject_id, tmp_dir)) + print("Downloading subject {} to {}".format(subject_id, tmp_dir)) fnout, _ = urllib.request.urlretrieve( - url, - os.path.join(tmp_dir, f"{subject_id}.tar.gz") + url, os.path.join(tmp_dir, f"{subject_id}.tar.gz") ) - print(f'Done downloading to {fnout}') + print(f"Done downloading to {fnout}") # Un-tar to pycortex store with tarfile.open(fnout, "r:gz") as tar: print("Extracting subject {} to {}".format(subject_id, pycortex_store)) @@ -1259,51 +1374,56 @@ def download_subject(subject_id='fsaverage', url=None, pycortex_store=None, def rotate_flatmap(surf_id, theta, plot=False): """Rotate flatmap to be less V-shaped - + Parameters ---------- surf_id : str pycortex surface identifier theta : scalar - angle in degrees to rotate flatmaps (rotation is clockwise + angle in degrees to rotate flatmaps (rotation is clockwise for right hemisphere and counter-clockwise for left) plot : bool Whether to make a coarse plot to visualize the changes """ # Lazy load of matplotlib import matplotlib.pyplot as plt - paths = db.get_paths(surf_id)['surfs']['flat'] + + paths = db.get_paths(surf_id)["surfs"]["flat"] theta = np.radians(theta) if plot: fig, axs = plt.subplots(2, 2) - for j, hem in enumerate(('lh','rh')): + for j, hem in enumerate(("lh", "rh")): this_file = paths[hem] pts, polys = formats.read_gii(this_file) # Rotate clockwise (- rotation) for RH, counter-clockwise (+ rotation) for LH - if hem == 'rh': - rtheta = - theta + if hem == "rh": + rtheta = -theta else: rtheta = copy.copy(theta) - rotation_mat = np.array([[np.cos(rtheta), -np.sin(rtheta)], [np.sin(rtheta), np.cos(rtheta)]]) + rotation_mat = np.array( + [[np.cos(rtheta), -np.sin(rtheta)], [np.sin(rtheta), np.cos(rtheta)]] + ) rotated = rotation_mat.dot(pts[:, :2].T).T pts_new = pts.copy() pts_new[:, :2] = rotated new_file, bkup_num = copy.copy(this_file), 0 while os.path.exists(new_file): - new_file = this_file.replace('.gii', '_rotbkup%02d.gii'%bkup_num) + new_file = this_file.replace(".gii", "_rotbkup%02d.gii" % bkup_num) bkup_num += 1 - print('Backing up file at %s...' % new_file) + print("Backing up file at %s..." % new_file) shutil.copy(this_file, new_file) formats.write_gii(this_file, pts_new, polys) - print('Overwriting %s...' % this_file) + print("Overwriting %s..." % this_file) if plot: - axs[0,j].plot(*pts[::100, :2].T, marker='r.') - axs[0,j].axis('equal') - axs[1,j].plot(*pts_new[::100, :2].T, marker='b.') - axs[1,j].axis('equal') + axs[0, j].plot(*pts[::100, :2].T, marker="r.") + axs[0, j].axis("equal") + axs[1, j].plot(*pts_new[::100, :2].T, marker="b.") + axs[1, j].axis("equal") # Remove and back up overlays file - overlay_file = db.get_paths(surf_id)['overlays'] - shutil.copy(overlay_file, overlay_file.replace('.svg', '_rotbkup%02d.svg'%bkup_num)) + overlay_file = db.get_paths(surf_id)["overlays"] + shutil.copy( + overlay_file, overlay_file.replace(".svg", "_rotbkup%02d.svg" % bkup_num) + ) os.unlink(overlay_file) # Regenerate file svg = db.get_overlay(surf_id) diff --git a/examples/hcp/README.txt b/examples/hcp/README.txt new file mode 100644 index 00000000..b6f25b1f --- /dev/null +++ b/examples/hcp/README.txt @@ -0,0 +1,5 @@ +Examples with HCP fs_LR data +------------------------------ + +Examples showing how to use PyCortex to visualize HCP fs_LR 32k data, both on +the HCP template and resampled to fsaverage. diff --git a/examples/hcp/plot_hcp_fs_lr.py b/examples/hcp/plot_hcp_fs_lr.py new file mode 100644 index 00000000..f24d9cae --- /dev/null +++ b/examples/hcp/plot_hcp_fs_lr.py @@ -0,0 +1,77 @@ +""" +================================================================== +Visualize HCP fs_LR data on the HCP template and on fsaverage +================================================================== + +This example shows how to: + +a) visualize data defined on the HCP fs_LR 32k surface directly on the HCP + template (the ``32k_fs_LR`` pycortex subject), and + +b) resample that data to fsaverage6 with :mod:`cortex.hcp` and visualize it as + a flatmap. + +Both steps produce a flatmap. + +Notes +----- +* The first time you run this, the ``32k_fs_LR`` and ``fsaverage`` subjects are + downloaded into your pycortex filestore. +* The resampling matrix is built from the HCP standard-mesh spheres, which are + downloaded on demand (no Connectome Workbench needed). +* fsaverage6 is not itself a renderable pycortex subject, so the fsaverage6 + result is upsampled to the full ``fsaverage`` surface for display using + :func:`cortex.freesurfer.upsample_to_fsaverage`. The fsaverage5/6 upsampling + tables ship with pycortex, so no FreeSurfer installation is needed. +""" + +import matplotlib.pyplot as plt +import numpy as np + +import cortex +import cortex.hcp + +# --------------------------------------------------------------------------- +# a) Visualize data on the HCP fs_LR 32k template +# --------------------------------------------------------------------------- + +hcp_subject = "32k_fs_LR" + +# Make sure the HCP template is in the pycortex filestore, downloading if not. +if hcp_subject not in cortex.db.subjects: + cortex.hcp.download_fs_lr() + +# Create some smooth demo data on the full fs_LR 32k surface (64984 vertices, +# both hemispheres, medial wall included). Here we use the anterior-posterior +# (y) coordinate of the inflated surface as a smooth gradient. For real HCP +# data stored as CIFTI grayordinates, expand it to the full surface first with +# ``cortex.hcp.cifti_to_surface`` (or go straight to fsaverage with +# ``cortex.hcp.to_fsaverage``). +pts, _ = cortex.db.get_surf(hcp_subject, "inflated", merge=True) +data_fslr = pts[:, 1].astype(float) +vmin, vmax = np.percentile(data_fslr, [1, 99]) + +# Visualize on the HCP template just like any other vertex dataset. +vtx_hcp = cortex.Vertex(data_fslr, hcp_subject, vmin=vmin, vmax=vmax, cmap="turbo") +cortex.quickshow(vtx_hcp, with_curvature=True, with_rois=False, with_labels=False) +plt.gcf().suptitle("HCP fs_LR 32k (native)") + +# --------------------------------------------------------------------------- +# b) Resample the same data to fsaverage6 and visualize it +# --------------------------------------------------------------------------- + +# Project fs_LR 32k -> fsaverage6 (81924 vertices, both hemispheres). +data_fs6 = cortex.hcp.project_fslr_to_fsaverage(data_fslr, target="fsaverage6") +print("fsaverage6 data shape:", data_fs6.shape) + +# fsaverage6 is not a renderable pycortex subject, so upsample the fsaverage6 +# result to the full fsaverage surface for the flatmap. +if "fsaverage" not in cortex.db.subjects: + cortex.download_subject("fsaverage") +data_fs = cortex.freesurfer.upsample_to_fsaverage(np.nan_to_num(data_fs6), "fsaverage6") + +vtx_fs = cortex.Vertex(data_fs, "fsaverage", vmin=vmin, vmax=vmax, cmap="turbo") +cortex.quickshow(vtx_fs, with_curvature=True, with_rois=False, with_labels=False) +plt.gcf().suptitle("Resampled to fsaverage6 (shown on fsaverage)") + +plt.show() diff --git a/pyproject.toml b/pyproject.toml index 6acd93aa..29a3c476 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,9 @@ test = [ "pytest-cov", "pytest-timeout", ] +# Note: the cortex.hcp barycentric-resampling validation test additionally +# needs the external Connectome Workbench ``wb_command`` binary on PATH (not a +# Python package). That test is skipped automatically when it is absent. [project.optional-dependencies] headless = [ diff --git a/setup.py b/setup.py index 86823c57..3737f18e 100644 --- a/setup.py +++ b/setup.py @@ -17,17 +17,17 @@ def set_default_filestore(prefix, optfile): config.read(optfile) config.set("basic", "filestore", os.path.join(prefix, "db")) config.set("webgl", "colormaps", os.path.join(prefix, "colormaps")) - with open(optfile, 'w') as fp: + with open(optfile, "w") as fp: config.write(fp) class my_install(install): def run(self): install.run(self) - optfile = [f for f in self.get_outputs() if 'defaults.cfg' in f] + optfile = [f for f in self.get_outputs() if "defaults.cfg" in f] prefix = os.path.join(self.install_base, "share", "pycortex") set_default_filestore(prefix, optfile[0]) - self.copy_tree('filestore', prefix) + self.copy_tree("filestore", prefix) for root, folders, files in os.walk(prefix): for folder in folders: os.chmod(os.path.join(root, folder), 511) @@ -42,110 +42,109 @@ def run(self): # ] data_files = [ # [10:] to remove "filestore/" - ('share/pycortex/' + os.path.dirname(file)[10:], [file]) - for file in glob('filestore/**/*', recursive=True) + ("share/pycortex/" + os.path.dirname(file)[10:], [file]) + for file in glob("filestore/**/*", recursive=True) if os.path.isfile(file) ] +ctm = Extension( + "cortex.openctm", + [ + "cortex/openctm.pyx", + "OpenCTM-1.0.3/lib/openctm.c", + "OpenCTM-1.0.3/lib/stream.c", + "OpenCTM-1.0.3/lib/compressRAW.c", + "OpenCTM-1.0.3/lib/compressMG1.c", + "OpenCTM-1.0.3/lib/compressMG2.c", + "OpenCTM-1.0.3/lib/liblzma/Alloc.c", + "OpenCTM-1.0.3/lib/liblzma/LzFind.c", + "OpenCTM-1.0.3/lib/liblzma/LzmaDec.c", + "OpenCTM-1.0.3/lib/liblzma/LzmaEnc.c", + "OpenCTM-1.0.3/lib/liblzma/LzmaLib.c", + ], + libraries=["m"], + include_dirs=[ + "OpenCTM-1.0.3/lib/", + "OpenCTM-1.0.3/lib/liblzma/", + numpy.get_include(), + ], + define_macros=[ + ("LZMA_PREFIX_CTM", None), + ("OPENCTM_BUILD", None), + ("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION"), + # ('__DEBUG_', None), + ], +) +formats = Extension( + "cortex.formats", ["cortex/formats.pyx"], include_dirs=[numpy.get_include()] +) - - -ctm = Extension('cortex.openctm', [ - 'cortex/openctm.pyx', - 'OpenCTM-1.0.3/lib/openctm.c', - 'OpenCTM-1.0.3/lib/stream.c', - 'OpenCTM-1.0.3/lib/compressRAW.c', - 'OpenCTM-1.0.3/lib/compressMG1.c', - 'OpenCTM-1.0.3/lib/compressMG2.c', - 'OpenCTM-1.0.3/lib/liblzma/Alloc.c', - 'OpenCTM-1.0.3/lib/liblzma/LzFind.c', - 'OpenCTM-1.0.3/lib/liblzma/LzmaDec.c', - 'OpenCTM-1.0.3/lib/liblzma/LzmaEnc.c', - 'OpenCTM-1.0.3/lib/liblzma/LzmaLib.c',], - libraries=['m'], include_dirs=[ - 'OpenCTM-1.0.3/lib/', - 'OpenCTM-1.0.3/lib/liblzma/', numpy.get_include()], - define_macros=[ - ('LZMA_PREFIX_CTM', None), - ('OPENCTM_BUILD', None), - ('NPY_NO_DEPRECATED_API', 'NPY_1_7_API_VERSION'), - #('__DEBUG_', None), - ] - ) -formats = Extension('cortex.formats', ['cortex/formats.pyx'], - include_dirs=[numpy.get_include()]) - -DISTNAME = 'pycortex' +DISTNAME = "pycortex" # VERSION is now automatically derived from git tags via setuptools-scm -DESCRIPTION = 'Python Cortical mapping software for fMRI data' -with open('README.md') as f: +DESCRIPTION = "Python Cortical mapping software for fMRI data" +with open("README.md") as f: LONG_DESCRIPTION = f.read() -AUTHOR = 'James Gao' -AUTHOR_EMAIL = 'james@jamesgao.com' -LICENSE = '2-clause BSD license' -URL = 'http://gallantlab.github.io/pycortex' +AUTHOR = "James Gao" +AUTHOR_EMAIL = "james@jamesgao.com" +LICENSE = "2-clause BSD license" +URL = "http://gallantlab.github.io/pycortex" DOWNLOAD_URL = URL -with open('requirements.txt') as f: +with open("requirements.txt") as f: INSTALL_REQUIRES = [ - line.strip() - for line in f - if line.strip() and not line.lstrip().startswith('#') + line.strip() for line in f if line.strip() and not line.lstrip().startswith("#") ] -setup(name=DISTNAME, - description=DESCRIPTION, - long_description=LONG_DESCRIPTION, - long_description_content_type='text/markdown', - author=AUTHOR, - author_email=AUTHOR_EMAIL, - license=LICENSE, - url=URL, - download_url=DOWNLOAD_URL, - packages=[ - 'cortex', - 'cortex.webgl', - 'cortex.mapper', - 'cortex.dataset', - 'cortex.blender', - 'cortex.tests', - 'cortex.quickflat', - 'cortex.polyutils', - 'cortex.export' - ], - data_files=data_files, - ext_modules=cythonize([ctm, formats]), - package_data={ - 'cortex': [ - 'svgbase.xml', - 'defaults.cfg', - 'bbr.sch' - ], - 'cortex.webgl': [ - '*.html', - 'favicon.ico', - 'resources/js/*.js', - 'resources/js/ctm/*.js', - 'resources/css/*.css', - 'resources/css/images/*', - 'resources/css/ui-lightness/*.css', - 'resources/css/ui-lightness/images/*', - 'resources/images/*' - ] - }, - setup_requires=['Cython', 'numpy'], - install_requires=INSTALL_REQUIRES, - # Don't use `extras_require` here. Put them in pyproject.toml . - cmdclass=dict(install=my_install), - include_package_data=True, - classifiers=[ - 'Development Status :: 6 - Mature', - 'Intended Audience :: Science/Research', - 'License :: OSI Approved :: BSD License', - 'Operating System :: OS Independent', - 'Programming Language :: Python', - 'Programming Language :: Python :: Implementation :: CPython', - 'Topic :: Scientific/Engineering :: Visualization' - ] +setup( + name=DISTNAME, + description=DESCRIPTION, + long_description=LONG_DESCRIPTION, + long_description_content_type="text/markdown", + author=AUTHOR, + author_email=AUTHOR_EMAIL, + license=LICENSE, + url=URL, + download_url=DOWNLOAD_URL, + packages=[ + "cortex", + "cortex.webgl", + "cortex.mapper", + "cortex.dataset", + "cortex.blender", + "cortex.tests", + "cortex.quickflat", + "cortex.polyutils", + "cortex.export", + ], + data_files=data_files, + ext_modules=cythonize([ctm, formats]), + package_data={ + "cortex": ["svgbase.xml", "defaults.cfg", "bbr.sch", "data/*.npz"], + "cortex.webgl": [ + "*.html", + "favicon.ico", + "resources/js/*.js", + "resources/js/ctm/*.js", + "resources/css/*.css", + "resources/css/images/*", + "resources/css/ui-lightness/*.css", + "resources/css/ui-lightness/images/*", + "resources/images/*", + ], + }, + setup_requires=["Cython", "numpy"], + install_requires=INSTALL_REQUIRES, + # Don't use `extras_require` here. Put them in pyproject.toml . + cmdclass=dict(install=my_install), + include_package_data=True, + classifiers=[ + "Development Status :: 6 - Mature", + "Intended Audience :: Science/Research", + "License :: OSI Approved :: BSD License", + "Operating System :: OS Independent", + "Programming Language :: Python", + "Programming Language :: Python :: Implementation :: CPython", + "Topic :: Scientific/Engineering :: Visualization", + ], )