diff --git a/cortex/__init__.py b/cortex/__init__.py index 0c0bc11c..a1a61067 100644 --- a/cortex/__init__.py +++ b/cortex/__init__.py @@ -29,6 +29,11 @@ except ImportError: pass +try: + from cortex import hcp +except ImportError: + pass + # Create deprecated interface for database class dep(object): def __getattr__(self, name): 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..6a9c4100 100644 --- a/cortex/freesurfer.py +++ b/cortex/freesurfer.py @@ -3,6 +3,7 @@ from __future__ import print_function import copy +import functools import os import shlex import shutil @@ -19,7 +20,7 @@ from scipy.sparse import coo_matrix from scipy.spatial import KDTree -from . import anat, database +from . import anat, appdirs, database def get_paths(fs_subject, hemi, type="patch", freesurfer_subject_dir=None): @@ -1114,6 +1115,77 @@ def stretch_mwall(pts, polys, mwall): return SpringLayout(pts, polys, inflated, pins=mwall) +def _n_vertices_ico(icoorder): + return 4 ** icoorder * 10 + 2 + + +# Precomputed upsampling neighbor arrays bundled with pycortex (see +# ``_get_upsample_neighbors``). ``NpzFile`` keyed by ``"{hemi}_ico{order}"``. +_UPSAMPLE_NEIGHBORS_FILE = os.path.join( + os.path.dirname(__file__), "data", "upsample_fsaverage_neighbors.npz") + + +@functools.lru_cache(maxsize=1) +def _load_upsample_neighbors_bundle(): + """Load (once) the neighbor arrays bundled with pycortex, or {} if absent.""" + if os.path.exists(_UPSAMPLE_NEIGHBORS_FILE): + return np.load(_UPSAMPLE_NEIGHBORS_FILE) + return {} + + +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. + """ + key = f"{hemi}_ico{ico_order}" + + # 1. Bundled with the package. + bundle = _load_upsample_neighbors_bundle() + bundle_keys = getattr(bundle, "files", list(bundle)) + if key in bundle_keys: + return np.asarray(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 ): @@ -1129,8 +1201,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 +1214,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,28 +1236,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 + 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) ) diff --git a/cortex/hcp.py b/cortex/hcp.py new file mode 100644 index 00000000..ee46ee0e --- /dev/null +++ b/cortex/hcp.py @@ -0,0 +1,633 @@ +"""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 functions: + +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. :func:`to_fsaverage` -- resample fs_LR 32k data to ``fsaverage`` (the + lower-level :func:`project_fslr_to_fsaverage` does the surface-to-surface + step). 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", + }, +} + + +_HEMI_ALIASES = {"L": "L", "lh": "L", "R": "R", "rh": "R"} + + +def _normalize_hemi(hemi): + """Normalize a hemisphere name to ``"L"``/``"R"`` (accepts ``lh``/``rh``).""" + try: + return _HEMI_ALIASES[hemi] + except (KeyError, TypeError): + raise ValueError( + f"hemi must be 'L'/'lh' or 'R'/'rh', got {hemi!r}" + ) from None + + +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)}" + ) + hemi = _normalize_hemi(hemi) + + 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 + + n_grayordinates = n_left + n_right + if data.shape[-1] != n_grayordinates: + raise ValueError( + f"data has {data.shape[-1]} grayordinates on its last axis, but the " + f"indices specify {n_grayordinates} ({n_left} left + {n_right} right)." + ) + + # Ellipsis indexing supports 1-D, 2-D, and higher-dimensional inputs + # (e.g. subjects x time x grayordinates). + out_dtype = np.result_type(data.dtype, np.float32) + full = np.full(data.shape[:-1] + (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) + """ + hemi = _normalize_hemi(hemi) + 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) + if data.shape[-1] != N_VERTICES_FS_LR_32K: + raise ValueError( + f"Expected {N_VERTICES_FS_LR_32K} fs_LR 32k vertices on the last " + f"axis, got {data.shape[-1]}." + ) + # Preserve the input's floating precision (float32 timeseries are common and + # large); only integer inputs are promoted, to float32. Flatten any leading + # dimensions so 1-D, 2-D, and N-D inputs are all handled uniformly. + out_dtype = data.dtype if np.issubdtype(data.dtype, np.floating) else np.float32 + leading_shape = data.shape[:-1] + data_2d = data.reshape(-1, N_VERTICES_FS_LR_32K).astype(out_dtype, copy=False) + + 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_2d.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_2d[:, src_sl], nanmean) + + if target == "fsaverage": + result = upsample_to_fsaverage( + result, "fsaverage6", freesurfer_subjects_dir=freesurfer_subjects_dir + ) + + return result.reshape(leading_shape + (result.shape[-1],)) + + +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..19d7283f 100644 --- a/cortex/tests/test_freesurfer.py +++ b/cortex/tests/test_freesurfer.py @@ -9,6 +9,7 @@ _remove_disconnected_polys, _surf2surf_nnfr_matrix, get_mri_surf2surf_matrix, + upsample_to_fsaverage, ) @@ -122,6 +123,34 @@ 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")) diff --git a/cortex/tests/test_hcp.py b/cortex/tests/test_hcp.py new file mode 100644 index 00000000..a031df23 --- /dev/null +++ b/cortex/tests/test_hcp.py @@ -0,0 +1,249 @@ +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, + _normalize_hemi, + cifti_to_surface, + get_cifti_vertex_indices, +) + + +def test_normalize_hemi_accepts_aliases(): + assert _normalize_hemi("L") == "L" and _normalize_hemi("lh") == "L" + assert _normalize_hemi("R") == "R" and _normalize_hemi("rh") == "R" + with pytest.raises(ValueError): + _normalize_hemi("left") + + +def test_cifti_to_surface_ndim_and_validation(): + left_idx = [0, 2, 5] + right_idx = [1, 3] + # 3-D input (subjects x time x grayordinates), 5 grayordinates on last axis. + data = np.arange(2 * 4 * 5, dtype=float).reshape(2, 4, 5) + full = cifti_to_surface(data, left_idx, right_idx) + assert full.shape == (2, 4, hcp.N_VERTICES_FS_LR_32K) + np.testing.assert_array_equal(full[..., left_idx], data[..., :3]) + np.testing.assert_array_equal( + full[..., hcp.N_VERTICES_FS_LR_32K_HEM + np.array(right_idx)], data[..., 3:5] + ) + # Wrong number of grayordinates on the last axis is rejected. + with pytest.raises(ValueError): + cifti_to_surface(np.zeros((2, 7)), left_idx, right_idx) + + +# --------------------------------------------------------------------------- +# 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..4badb51b 100644 --- a/cortex/utils.py +++ b/cortex/utils.py @@ -1165,7 +1165,7 @@ def get_cmap(name): try: cmap = plt.cm.get_cmap(name) except: - raise Exception('Unkown color map!') + raise Exception('Unknown color map!') return cmap def add_cmap(cmap, name, cmapdir=None): @@ -1229,6 +1229,10 @@ def download_subject(subject_id='fsaverage', url=None, pycortex_store=None, # Map codes to URLs; more coming eventually id_to_url = dict( 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: 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/setup.py b/setup.py index 86823c57..c9d36fc4 100644 --- a/setup.py +++ b/setup.py @@ -120,7 +120,8 @@ def run(self): 'cortex': [ 'svgbase.xml', 'defaults.cfg', - 'bbr.sch' + 'bbr.sch', + 'data/*.npz' ], 'cortex.webgl': [ '*.html',