Source code for cortex.freesurfer

"""Contains functions for interfacing with freesurfer
from __future__ import print_function
import os
import copy
import shutil
import struct
import tempfile
import warnings
import shlex
import subprocess as sp
from builtins import input

import numpy as np
import nibabel
from nibabel import gifti
from tempfile import NamedTemporaryFile
from scipy.spatial.kdtree import KDTree
from scipy.linalg import lstsq
from scipy.sparse import coo_matrix

from . import database
from . import anat

[docs]def get_paths(fs_subject, hemi, type="patch", freesurfer_subject_dir=None): """Retrive paths for all surfaces for a subject processed by freesurfer Parameters ---------- subject : string Subject ID for freesurfer subject hemi : string ['lh'|'rh'] Left ('lh') or right ('rh') hemisphere type : string ['patch'|'surf'|'curv'] Which type of files to return freesurfer_subject_dir : string | None Directory of freesurfer subjects. Defaults to the value for the environment variable 'SUBJECTS_DIR' (which should be set by freesurfer) """ if freesurfer_subject_dir is None: 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") elif type == "surf": return os.path.join(base, "surf", hemi+".{name}") elif type == "curv": return os.path.join(base, "surf", hemi+".curv{name}") elif type == "slim": return os.path.join(base, "surf", hemi+".{name}_slim.obj")
[docs]def autorecon(fs_subject, type="all", parallel=True, n_cores=None): """Run Freesurfer's autorecon-all command for a given freesurfer subject Parameters ---------- fs_subject : string Freesurfer subject ID (should be a folder in your freesurfer $SUBJECTS_DIR) type : string Which steps of autorecon-all to perform. {'all', '1','2','3','cp','wm', 'pia'} """ 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" } if str(type) in times: 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'): # 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) print("Calling:\n{cmd}".format(cmd=cmd)) sp.check_call(shlex.split(cmd))
[docs]def flatten(fs_subject, hemi, patch, freesurfer_subject_dir=None, save_every=None): """Perform flattening of a brain using freesurfer Parameters ---------- fs_subject : str Freesurfer subject ID hemi : str ['lh' | 'rh'] hemisphere to flatten patch : str name for freesurfer patch (used as `name` argument to format output of `get_paths()`) freesurfer_subject_dir : str Freesurfer subjects directory location. None defaults to $SUBJECTS_DIR save_every: int If not None, this saves a version of the mesh every `save_every` iterations of the flattening process. Useful for determining why a flattening fails. Returns ------- Notes ----- To look into: link below shows how to give continuous output for a subprocess. There maybe indications that a flattening is going badly that we could detect in the stdout; perhaps even continuously update a visualization of the generated files using segment.show_surface() with the outputs (triggered to update once stdout shows that a flattening iteration has completed) """ 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 = '' 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) print("Calling: ") print(cmd) sp.check_call(shlex.split(cmd)) return True else: print("Not going to flatten...") return False
[docs]def import_subj(fs_subject, cx_subject=None, freesurfer_subject_dir=None, whitematter_surf='smoothwm'): """Imports a subject from freesurfer 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. Parameters ---------- fs_subject : string Freesurfer subject name cx_subject : string, optional Pycortex subject name (These variable names should be changed). By default uses the same name as the freesurfer subject. Best to stick to that convention, if possible (your life will go more smoothly.) This optional kwarg is for edge cases. freesurfer_subject_dir : string, optional Freesurfer subject directory to pull data from. By default uses the directory given by the environment variable $SUBJECTS_DIR. whitematter_surf : string, optional 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'. """ if cx_subject is None: cx_subject = fs_subject # Create and/or replace extant subject. Throws a warning that this will happen. database.db.make_subj(cx_subject) import nibabel surfs = os.path.join(database.default_filestore, cx_subject, "surfaces", "{name}_{hemi}.gii") anats = os.path.join(database.default_filestore, cx_subject, "anatomicals", "{name}.nii.gz") surfinfo = os.path.join(database.default_filestore, cx_subject, "surface-info", "{name}.npz") if freesurfer_subject_dir is None: freesurfer_subject_dir = os.environ['SUBJECTS_DIR'] fspath = os.path.join(freesurfer_subject_dir, fs_subject, 'mri') curvs = os.path.join(freesurfer_subject_dir, fs_subject, 'surf', '{hemi}.{name}') #import anatomicals for fsname, name in dict(T1="raw", aseg="aseg", wm="raw_wm").items(): path = os.path.join(fspath, "{fsname}.mgz").format(fsname=fsname) out = anats.format(subj=cx_subject, name=name) cmd = "mri_convert {path} {out}".format(path=path, out=out) sp.check_output(shlex.split(cmd)) # (Re-)Make the fiducial files # NOTE: these are IN THE FREESURFER $SUBJECTS_DIR !! which can cause confusion. make_fiducial(fs_subject, freesurfer_subject_dir=freesurfer_subject_dir) # Freesurfer uses FOV/2 for center, let's set the surfaces to use the # magnet isocenter trans = nibabel.load(out).get_affine()[:3, -1] surfmove = trans - np.sign(trans) * [128, 128, 128] from . import formats for fsname, name in [(whitematter_surf, "wm"), ('pial', "pia"), ('inflated', "inflated")]: for hemi in ("lh", "rh"): pts, polys, _ = get_surf(fs_subject, hemi, fsname, freesurfer_subject_dir=freesurfer_subject_dir) fname = str(surfs.format(subj=cx_subject, name=name, hemi=hemi)) formats.write_gii(fname, pts=pts + surfmove, polys=polys) for curv, info in dict(sulc="sulcaldepth", thickness="thickness", curv="curvature").items(): lh, rh = [parse_curv(curvs.format(hemi=hemi, name=curv)) for hemi in ['lh', 'rh']] np.savez(surfinfo.format(subj=cx_subject, name=info), left=-lh, right=-rh) database.db = database.Database()
[docs]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). Parameters ---------- fs_subject : str Freesurfer subject name patch : str Name of flat.patch.3d file; e.g., "flattenv01" hemis : list List of hemispheres to import. Defaults to both hemispheres. cx_subject : str Pycortex subject name freesurfer_subject_dir : str directory for freesurfer subjects. None defaults to evironment variable $SUBJECTS_DIR clean : bool If True, the flat surface is cleaned to remove the disconnected polys. Returns ------- """ 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', '']: 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") from . import formats for hemi in hemis: if flat_type == 'freesurfer': pts, polys, _ = get_surf(fs_subject, hemi, "patch", patch+".flat", freesurfer_subject_dir=freesurfer_subject_dir) # Reorder axes: X, Y, Z instead of Y, X, Z flat = pts[:, [1, 0, 2]] # Flip Y axis upside down flat[:, 1] = -flat[:, 1] 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) if clean: polys = _remove_disconnected_polys(polys) flat = _move_disconnect_points_to_zero(flat, polys) fname = surfs.format(hemi=hemi) 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'] os.unlink(overlays_file)
# 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. """ n_points = np.max(polys) + 1 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] ]) data = np.ones(len(col), 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) non_trivial_components = unique_labels[np.where(counts > 1)[0]] main_component = unique_labels[np.argmax(counts)] extra_components = non_trivial_components[non_trivial_components != main_component] # filter all components not in the largest component disconnected_pts = np.where(np.isin(labels, extra_components))[0] disconnected_polys_mask = np.isin(polys[:, 0], disconnected_pts) return polys[~disconnected_polys_mask] 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). """ mask = np.zeros(len(pts), dtype=bool) mask[np.unique(polys)] = True pts[~mask] = 0 return pts
[docs]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") write_surf(fname, (spts + ppts) / 2, polys)
[docs]def parse_surf(filename): """ """ with open(filename, 'rb') as fp: #skip magic comment = fp.readline() fp.readline() print(comment) verts, faces = struct.unpack('>2I', pts = np.fromstring(*3*verts), dtype='f4').byteswap() polys = np.fromstring(*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))) fp.write(pts.astype(np.float32).byteswap().tostring()) fp.write(polys.astype(np.uint32).byteswap().tostring()) fp.write(b'\n') def write_patch(filename, pts, edges=None): """Writes a patch file that is readable by freesurfer. Note this function is duplicated here and in blendlib. This function writes freesurfer format, so seems natural to place here, but it also needs to be called from blender, and the blendlib functions are the only ones currently that can easily be called in a running blender session. Parameters ---------- filename : name for patch to write. Should be of the form <subject>.flatten.3d pts : array-like points in the mesh edges : array-like edges in the mesh. """ if edges is None: edges = set() 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)) else: fp.write(struct.pack('>i3f', i+1, *pt))
[docs]def parse_curv(filename): """ """ with open(filename, 'rb') as fp: return np.fromstring(, dtype='>f4').byteswap().newbyteorder()
[docs]def parse_patch(filename): """ """ with open(filename, 'rb') as fp: header, = struct.unpack('>i', nverts, = struct.unpack('>i', data = np.fromstring(, dtype=[('vert', '>i4'), ('x', '>f4'), ('y', '>f4'), ('z', '>f4')]) assert len(data) == nverts return data
[docs]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') else: 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) if flatten_step is not None: patch_file += '%04d'%flatten_step patch = parse_patch(patch_file) verts = patch[patch['vert'] > 0]['vert'] - 1 edges = -patch[patch['vert'] < 0]['vert'] - 1 idx = np.zeros((len(pts),), dtype=bool) idx[verts] = True idx[edges] = True valid = idx[polys.ravel()].reshape(-1, 3).all(1) polys = polys[valid] idx = np.zeros((len(pts),)) idx[verts] = 1 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] return pts, polys, idx 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'): """subject is a freesurfer subject""" if fs_dir is None: 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)) if not os.path.exists(srclabel): raise ValueError("Label {} doesn't exist!".format(srclabel)) 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) 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) stdout, stderr = proc.communicate() if stderr not in ('', b''): raise Exception("Error in freesurfer function call:\n{}".format(stderr)) print("Labels transferred") def _parse_labels(label_files, cx_subject): """Extract values from freesurfer label file(s) and map to vertices Parameters ---------- label_files : str or list full paths to label file or files to load cx_subject : str pycortex subject ID """ if not isinstance(label_files, (list, tuple)): label_files = [label_files] verts = [] values = [] 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:]] 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]) 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): """Get data from a label file for fsaverage subject Parameters ---------- cx_subject : str A pycortex subject ID label : str Label name fs_subject : str Freesurfer subject ID, if different from pycortex subject ID src_subject : str Freesurfer subject ID from which to transfer the label. fs_dir : str Freesurfer subject directory; None defaults to OS environment variable hemisphere : list | tuple """ if fs_dir is None: fs_dir = os.environ['SUBJECTS_DIR'] else: 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']: # If label file doesn't exist, try to move it there 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) verts, values = _parse_labels(label_files, cx_subject) idx = verts.astype( return idx, values def _mri_surf2surf_command(src_subj, trg_subj, input_file, output_file, hemi): # mri_surf2surf --srcsubject <source subject name> --srcsurfval # <sourcefile> --trgsubject <target suhject name> --trgsurfval <target # file> --hemi <hemifield> cmd = ["mri_surf2surf", "--srcsubject", src_subj, "--sval", input_file, "--trgsubject", trg_subj, "--tval", output_file, "--hemi", hemi, ] return cmd 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. """ data_arrays = [gifti.GiftiDataArray(d) for d in data] gifti_image = gifti.GiftiImage(darrays=data_arrays) tf_in = NamedTemporaryFile(suffix=".gii"), tf_out = NamedTemporaryFile(suffix='.gii') cmd = _mri_surf2surf_command(source_subj, target_subj,,, hemi) if subjects_dir is not None: env = os.environ.copy() env['SUBJECTS_DIR'] = subjects_dir else: env = None 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)) tf_in.close() output_img = nibabel.load( output_data = np.array([ for da in output_img.darrays]) tf_out.close() return output_data def get_mri_surf2surf_matrix(source_subj, hemi, surface_type, target_subj='fsaverage', subjects_dir=None, n_neighbors=20, random_state=0, n_test_images=40, coef_threshold=None, renormalize=True): """Creates a matrix implementing freesurfer mri_surf2surf command. A surface-to-surface transform is a linear transform between vertex spaces. Such a transform must be highly localized in the sense that a vertex in the target surface only draws its values from very few source vertices. This function exploits the localization to create an inverse problem for each vertex. The source neighborhoods for each target vertex are found by using mri_surf2surf to transform the three coordinate maps from the source surface to the target surface, yielding three coordinate values for each target vertex, for which we find the nearest neighbors in the source space. A small number of test images is transformed from source surface to target surface. For each target vertex in the transformed test images, a regression is performed using only the corresponding source image neighborhood, yielding the entries for a sparse matrix encoding the transform. Parameters ========== source_subj: str Freesurfer name of source subject hemi: str in ("lh", "rh") Indicator for hemisphere surface_type: str in ("white", "pial", ...) Indicator for surface layer target_subj: str, default "fsaverage" Freesurfer name of target subject subjects_dir: str, default os.environ["SUBJECTS_DIR"] The freesurfer subjects directory n_neighbors: int, default 20 The size of the neighborhood to take into account when estimating the source support of a vertex random_state: int, default 0 Random number generator or seed for generating test images n_test_images: int, default 40 Number of test images transformed to compute inverse problem. This should be greater than n_neighbors or equal. coef_treshold: float, default 1 / (10 * n_neighbors) Value under which to set a weight to zero in the inverse problem. renormalize: boolean, default True Determines whether the rows of the output matrix should add to 1, implementing what is sensible: a weighted averaging Notes ===== It turns out that freesurfer seems to do the following: For each target vertex, find, on the sphere, the nearest source vertices, and average their values. Try to be as one-to-one as possible. """ source_verts, _, _ = get_surf(source_subj, hemi, surface_type, freesurfer_subject_dir=subjects_dir) transformed_coords = mri_surf2surf(source_verts.T, source_subj, target_subj, hemi, subjects_dir=subjects_dir) kdt = KDTree(source_verts) print("Getting nearest neighbors") distances, indices = kdt.query(transformed_coords.T, k=n_neighbors) print("Done") rng = (np.random.RandomState(random_state) if isinstance(random_state, int) else random_state) test_images = rng.randn(n_test_images, len(source_verts)) transformed_test_images = mri_surf2surf(test_images, source_subj, target_subj, hemi, subjects_dir=subjects_dir) # Solve linear problems to get coefficients all_coefs = [] residuals = [] print("Computing coefficients") i = 0 for target_activation, source_inds in zip( transformed_test_images.T, indices): i += 1 print("{i}".format(i=i), end="\r") source_values = test_images[:, source_inds] r = lstsq(source_values, target_activation, overwrite_a=True, overwrite_b=True) all_coefs.append(r[0]) residuals.append(r[1]) print("Done") all_coefs = np.array(all_coefs) if coef_threshold is None: # we know now that coefs are doing averages coef_threshold = (1 / 10. / n_neighbors ) all_coefs[np.abs(all_coefs) < coef_threshold] = 0 if renormalize: all_coefs /= np.abs(all_coefs).sum(axis=1)[:, np.newaxis] + 1e-10 # there seem to be like 7 vertices that don't constitute an average over # 20 vertices or less, but all the others are such an average. # Let's make a matrix that does the transform: col_indices = indices.ravel() row_indices = (np.arange(indices.shape[0])[:, np.newaxis] * np.ones(indices.shape[1], dtype='int')).ravel() data = all_coefs.ravel() shape = (transformed_coords.shape[1], source_verts.shape[0]) matrix = coo_matrix((data, (row_indices, col_indices)), shape=shape) return matrix
[docs]def get_curv(fs_subject, hemi, type='wm', freesurfer_subject_dir=None): """Load freesurfer curv file for a freesurfer subject Parameters ---------- fs_subject : str freesurfer subject identifier hemi : str 'lh' or 'rh' for left or right hemisphere, respectively 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 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='') else: curv_file = get_paths(fs_subject, hemi, 'curv', freesurfer_subject_dir=freesurfer_subject_dir).format(name='.'+type) return parse_curv(curv_file)
[docs]def show_surf(subject, hemi, type, patch=None, curv=True, freesurfer_subject_dir=None): """Show a surface from a Freesurfer subject directory Parameters ---------- subject : str Freesurfer subject name hemi : str ['lh' | 'rh'] Left or right hemisphere type : patch : curv : bool 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!')) from mayavi import mlab from tvtk.api import tvtk 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) 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) 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'] fig.scene.background = (0,0,0) fig.scene.interactor.interactor_style = tvtk.InteractorStyleTerrain() path = os.path.join(os.environ['SUBJECTS_DIR'], subject) def picker_callback(picker): if in npts = np.append(, [pts[picker.point_id]], axis=0) = 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)) picker = fig.on_mouse_pick(picker_callback) picker.tolerance = 0.01 return fig, surf
[docs]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 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('node [shape=point,label=""];\n') for a, b in graph.edges_iter(): 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("}")
[docs]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: fp.readline() fp.readline() fp.readline() fp.readline() 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(' ') return data
[docs]def write_decimated(path, pts, polys): """ """ from .polyutils import decimate, boundary_edges dpts, dpolys = decimate(pts, polys) 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))) fp.write(data.tostring())
[docs]class SpringLayout(object): """ """
[docs] def __init__(self, pts, polys, dpts=None, pins=None, stepsize=1, neighborhood=0): self.pts = pts self.polys = polys self.stepsize = stepsize pinmask = np.zeros((len(pts),), dtype=bool) if isinstance(pins, (list, set, np.ndarray)): pinmask[pins] = True self.pins = pinmask self.neighbors = [set() for _ in range(len(pts))] for i, j, k in polys: self.neighbors[i].add(j) self.neighbors[i].add(k) self.neighbors[j].add(i) self.neighbors[j].add(k) self.neighbors[k].add(i) self.neighbors[k].add(j) for _ in range(neighborhood): _neighbors = copy.deepcopy(self.neighbors) for v, neighbors in enumerate(self.neighbors): for n in neighbors: _neighbors[v] |= self.neighbors[n] self.neighbors = _neighbors for i in range(len(self.neighbors)): self.neighbors[i] = list(self.neighbors[i] - set([i])) if dpts is None: dpts = 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._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._num[i] = len(n) self._dists = np.hstack(self._dists) self._idx = np.hstack(self._idx).astype(np.uint) self._neigh = np.hstack(self.neighbors).astype(np.uint) self.figure = None
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] fvec = force[:, np.newaxis] * svec self._move[self._mask] = self.stepsize * fvec return self._move.sum(1) / self._num[:, np.newaxis] def _estatic(self, idx): dist, neighbors = self.kdt.query(self.pts[idx], k=20) valid = dist > 0 mag = self.stepsize * (1 / dist) diff = self.pts[neighbors] - self.pts[idx] return (mag[valid] * diff[valid].T).T.mean(0) def step(self): move = self._spring()[~self.pins] 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) def run(self, n=1000): for _ in range(n): self.step() print(_) 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.step() self.figure.mlab_source.set(x=self.pts[:,0], y=self.pts[:,1], z=self.pts[:,2])
[docs]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]) 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)
# 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} if __name__ == "__main__": import sys show_surf(sys.argv[1], sys.argv[2], sys.argv[3])