From 3ad1939c2185eae4d99d07e7e7b5a8ef358e9098 Mon Sep 17 00:00:00 2001 From: MShinkle <49564869+MShinkle@users.noreply.github.com> Date: Thu, 26 Sep 2019 21:14:07 -0700 Subject: [PATCH] Fsaverage transform (#334) Related to fsaverage transform: * ENH add mri_surf2surf matrix * ENH: allow setting custom subjects_dir in mri_surf2surf * ENH: informative outputs for calls to freesurfer functions * ENH: added caching of surf2surf transforms in pycortex filestore (in surf2surf directory in each subject) * ENH: added Vertex.map() function to map vertex object from one subject (using cached transforms) to another subject via freesurfer's sphere registration. Not related to fsaverage, misc fixes included because this PR went on too long: * FIX: misc python3 fixes, incl removal of fstrings and whitespace * FIX: dict comprehension fix in anonymization code. * ENH: update to freesurfer alignment to allow optional reference image input from intial alignment * ENH: added output text to fs_manual() to clarify how to save transform, option for which surface to use for fs aligner. * FIX: killed string/bytes bug w/ addition of data to cutting in blender * ENH: Update to cut_surface to check whether fiducial vertex count matches inflated vertex count; avoids wasted work segmenting surfaces that won't flatten * FIX: potential config-stored cache was not correctly handled in freesurfer re-import of flatmaps. * minor changes (cleanup, simplification) to fs_align and xfm.py, minor changes to Volume.map() function (preserving vmin, vmax, cmap) --- cortex/align.py | 100 ++++++++++++----- cortex/blender/__init__.py | 5 + cortex/database.py | 67 ++++++++++- cortex/dataset/braindata.py | 6 + cortex/dataset/views.py | 35 ++++++ cortex/freesurfer.py | 218 +++++++++++++++++++++++++++++++++++- cortex/segment.py | 9 ++ cortex/utils.py | 56 +++++++++ cortex/xfm.py | 52 ++++----- 9 files changed, 482 insertions(+), 66 deletions(-) diff --git a/cortex/align.py b/cortex/align.py index fe8202d98..898397122 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -82,7 +82,8 @@ def view_callback(aligner): return m -def fs_manual(subject, xfmname, output_name="register.lta", wm_color="blue", pial_color="red", noclean=False): +def fs_manual(subject, xfmname, output_name="register.lta", wm_color="yellow", + pial_color="blue", wm_surface='white', noclean=False, reference=None, inspect_only=False): """Open Freesurfer FreeView GUI for manually aligning/adjusting a functional volume to the cortical surface for `subject`. This creates a new transform called `xfmname`. The name of a nibabel-readable file (e.g. NIfTI) should be @@ -90,6 +91,16 @@ def fs_manual(subject, xfmname, output_name="register.lta", wm_color="blue", pia IMPORTANT: This function assumes that the resulting .lta file is saved as: "{default folder chosen by FreeView (should be /tmp/fsalign_xxx)}/{output_name}". + + NOTE: Half-fixed some potential bugs in here, related to assumptions about how + results from mri_info calls would be formatted. IFF .dat files are written + based on nii files that have been stripped of their headers, then there will + be an extra line at the top stating that the coordinates are assumed to be in mm. + Without this line, the code here fails. Seems brittle, ripe for future bugs. + + ALSO: all the freesurfer environment stuff shouldn't be necessary, except that + I don't know what vox2ras-tkr is doing. + Parameters ---------- @@ -109,6 +120,15 @@ def fs_manual(subject, xfmname, output_name="register.lta", wm_color="blue", pia If True, intermediate files will not be removed from /tmp/fsalign_xxx (this is useful for debugging things), and the returned value will be the name of the temp directory. Default False. + reference : str + name of reference (generally, functional) volume. Only provide this if + you are working from scratch (if no transform exists already), else + it will throw an error. + inspect_only : boolean | False + Whether to open transform to view only (if True, nothing is saved + when freeview is closed) + wm_surface : string + name for white matter surface to use. 'white' or 'smoothwm' Returns ------- @@ -130,48 +150,68 @@ def fs_manual(subject, xfmname, output_name="register.lta", wm_color="blue", pia # if masks have been cached, quit! user must remove them by hand from glob import glob - if len(glob(db.get_paths(subject)['masks'].format(xfmname=xfmname, type='*'))): + masks_exist = len(glob(db.get_paths(subject)['masks'].format(xfmname=xfmname, type='*'))) + if masks_exist and not inspect_only: print('Refusing to overwrite existing transform %s because there are cached masks. Delete the masks manually if you want to modify the transform.' % xfmname) raise ValueError('Exiting...') + if reference is not None: + raise ValueError('Refusing to overwrite extant reference for transform') except IOError: - print("Transform does not exist!") - - # Load transform-relevant things - reference = sub_xfm.reference.get_filename() - xfm_dir = os.path.dirname(reference) - _ = sub_xfm.to_freesurfer(os.path.join(cache, "register.dat"), subject) # Transform in freesurfer .dat format - - # Command for FreeView and run - cmd = ("freeview -v $SUBJECTS_DIR/{sub}/mri/orig.mgz " - "{ref}:reg={reg} " - "-f $SUBJECTS_DIR/{sub}/surf/lh.white:edgecolor={wmc} $SUBJECTS_DIR/{sub}/surf/rh.white:edgecolor={wmc} " - "$SUBJECTS_DIR/{sub}/surf/lh.pial:edgecolor={pialc} $SUBJECTS_DIR/{sub}/surf/rh.pial:edgecolor={pialc}") - cmd = cmd.format(sub=subject, ref=reference, reg=os.path.join(cache, "register.dat"), - wmc=wm_color, pialc=pial_color) + if reference is None: + print("Transform does not exist!") + + if reference is None: + # Load load extant transform-relevant things + reference = sub_xfm.reference.get_filename() + _ = sub_xfm.to_freesurfer(os.path.join(cache, "register.dat"), subject) # Transform in freesurfer .dat format + # Command for FreeView and run + cmd = ("freeview -v $SUBJECTS_DIR/{sub}/mri/orig.mgz " + "{ref}:reg={reg} " + "-f $SUBJECTS_DIR/{sub}/surf/lh.{wms}:edgecolor={wmc} $SUBJECTS_DIR/{sub}/surf/rh.{wms}:edgecolor={wmc} " + "$SUBJECTS_DIR/{sub}/surf/lh.pial:edgecolor={pialc} $SUBJECTS_DIR/{sub}/surf/rh.pial:edgecolor={pialc}") + cmd = cmd.format(sub=subject, ref=reference, reg=os.path.join(cache, "register.dat"), + wmc=wm_color, pialc=pial_color, wms=wm_surface) + print('=== Calling (NO REFERENCE PROVIDED): ===') + print(cmd) + else: + # Command for FreeView and run + cmd = ("freeview -v $SUBJECTS_DIR/{sub}/mri/orig.mgz " + "{ref} " + "-f $SUBJECTS_DIR/{sub}/surf/lh.{wms}:edgecolor={wmc} $SUBJECTS_DIR/{sub}/surf/rh.{wms}:edgecolor={wmc} " + "$SUBJECTS_DIR/{sub}/surf/lh.pial:edgecolor={pialc} $SUBJECTS_DIR/{sub}/surf/rh.pial:edgecolor={pialc}") + cmd = cmd.format(sub=subject, ref=reference, + wmc=wm_color, pialc=pial_color, + wms=wm_surface) + print('=== Calling: ===') + print(cmd) + + if not inspect_only: + sfile = os.path.join(cache, output_name) + print('\nREGISTRATION MUST BE SAVED AS:\n\n{}'.format(sfile)) # Run and save transform when user is done editing if sp.call(cmd, shell=True) != 0: raise IOError("Problem with FreeView!") else: - # Convert transform into .dat format - reg_dat = os.path.join(cache, os.path.splitext(output_name)[0] + ".dat") - cmd = "lta_convert --inlta {inlta} --outreg {regdat}" - cmd = cmd.format(inlta=os.path.join(cache, output_name), regdat=reg_dat) - if sp.call(cmd, shell=True) != 0: - raise IOError("Error converting lta into dat!") - - # Save transform to pycortex - xfm = Transform.from_freesurfer(reg_dat, reference, subject) - db.save_xfm(subject, xfmname, xfm.xfm, xfmtype='coord', reference=reference) - - print("saved xfm:", subject, xfmname) - + if not inspect_only: + # Convert transform into .dat format + # Unclear why we're not just saving in .dat format above...? + reg_dat = os.path.join(cache, os.path.splitext(output_name)[0] + ".dat") + cmd = "lta_convert --inlta {inlta} --outreg {regdat}" + cmd = cmd.format(inlta=os.path.join(cache, output_name), regdat=reg_dat) + if sp.call(cmd, shell=True) != 0: + raise IOError("Error converting lta into dat!") + # Save transform to pycortex + xfm = Transform.from_freesurfer(reg_dat, reference, subject) + db.save_xfm(subject, xfmname, xfm.xfm, xfmtype='coord', reference=reference) + print("saved xfm") + except Exception as e: + raise(e) finally: if not noclean: shutil.rmtree(cache) else: retval = cache - return retval diff --git a/cortex/blender/__init__.py b/cortex/blender/__init__.py index 40926dc60..48a24c658 100644 --- a/cortex/blender/__init__.py +++ b/cortex/blender/__init__.py @@ -1,4 +1,5 @@ import os +import six import shlex import xdrlib import tempfile @@ -76,6 +77,10 @@ def add_cutdata(fname, braindata, name="retinotopy", projection="nearest", mesh= rcolor = cmap((right - vmin) / (vmax - vmin))[:,:3] p = xdrlib.Packer() + if six.PY3: + mesh = mesh.encode() + name = name.encode() + p.pack_string(mesh) p.pack_string(name) p.pack_array(lcolor.ravel(), p.pack_double) diff --git a/cortex/database.py b/cortex/database.py index b613b7cb3..2b87f21f0 100644 --- a/cortex/database.py +++ b/cortex/database.py @@ -168,7 +168,8 @@ def __getattr__(self, attr): def __dir__(self): return ["save_xfm","get_xfm", "get_surf", "get_anat", "get_surfinfo", "subjects", # "get_paths", # Add? - "get_mask", "get_overlay","get_cache", "get_view", "save_view", "get_mnixfm"] + list(self.subjects.keys()) + "get_mask", "get_overlay","get_cache", "get_view", "save_view", "get_mnixfm", + 'get_mri_surf2surf_matrix'] + list(self.subjects.keys()) @property def subjects(self): @@ -269,6 +270,69 @@ def get_surfinfo(self, subject, type="curvature", recache=False, **kwargs): npz.close() return Vertex(verts, subject) return npz + + def get_mri_surf2surf_matrix(self, subject, surface_type, hemi='both', + fs_subj=None, target_subj='fsaverage', + **kwargs): + """Get matrix generated by surf2surf to map one subject's surface to another's + + Parameters + ---------- + subject : pycortex + pycortex subject ID + surface_type : str + type of surface; one of ['fiducial', 'inflated', 'white', + 'pial'] ... more? + hemi : str or list + one of: 'both', 'lh', 'rh' + fs_subj : str + string ID for freesurfer subject (if different from + pycortex subject ID; None assumes they are the same) + target_subj : str + string ID for freesurfer subject to which to map surface + from `subject` (or `fs_subj`) + + Other Parameters + ---------------- + subjects_dir : str + Custom path to freesurfer subjects dir can be specified in the + keyword args + n_neighbors : int + ... + random_state : scalar + ... + n_test_images : int + ... + coef_threshold : scalar or None + ... + renormalize : bool + ... + + """ + from .freesurfer import get_mri_surf2surf_matrix as mri_s2s + from .utils import load_sparse_array, save_sparse_array + if fs_subj is None: + fs_subj = subject + fpath = self.get_paths(subject)['surf2surf'].format(source=fs_subj, target=target_subj) + # Backward compatiblity + fdir, _ = os.path.split(fpath) + if not os.path.exists(fdir): + print("Creating surf2surf directory for subject %s"%(subject)) + os.makedirs(fdir) + if hemi == 'both': + hemis = ['lh', 'rh'] + else: + hemis = [hemi] + if os.path.exists(fpath): + mats = [load_sparse_array(fpath, h) for h in hemis] + else: + mats = [] + for h in hemis: + tmp = mri_s2s(fs_subj, h, surface_type, + target_subj=target_subj, **kwargs) + mats.append(tmp) + save_sparse_array(fpath, tmp, h, mode='a') + return mats def get_overlay(self, subject, overlay_file=None, **kwargs): from . import svgoverlay @@ -596,6 +660,7 @@ def get_paths(self, subject): rois=os.path.join(self.filestore, subject, "rois.svg").format(subj=subject), overlays=os.path.join(self.filestore, subject, "overlays.svg").format(subj=subject), views=sorted([os.path.splitext(f)[0] for f in views]), + surf2surf=os.path.join(self.filestore, subject, "surf2surf", "{source}_to_{target}", "matrices.hdf"), ) return filenames diff --git a/cortex/dataset/braindata.py b/cortex/dataset/braindata.py index 2ebf3a433..45d29b358 100644 --- a/cortex/dataset/braindata.py +++ b/cortex/dataset/braindata.py @@ -256,6 +256,12 @@ def map(self, projection="nearest"): from cortex import utils mapper = utils.get_mapper(self.subject, self.xfmname, projection) data = mapper(self) + # Note: this is OK, because VolumeRGB and Volume2D objects (which + # have different requirements for vmin, vmax, cmap) do not inherit + # from VolumeData, and so do not have this method. + data.vmin = self.vmin + data.vmax = self.vmax + data.cmap = self.cmap return data def __repr__(self): diff --git a/cortex/dataset/views.py b/cortex/dataset/views.py index f18032263..6df7c3d30 100644 --- a/cortex/dataset/views.py +++ b/cortex/dataset/views.py @@ -333,7 +333,42 @@ def raw(self): description=self.description, state=self.state, **self.attrs) + def map(self, target_subj, surface_type='fiducial', + hemi='both', fs_subj=None, **kwargs): + """Map this data from this surface to another surface + + Calls `cortex.freesurfer.vertex_to_vertex()` with this + vertex object as the first argument. + + NOTE: Requires either previous computation of mapping matrices + (with `cortex.db.get_mri_surf2surf_matrix`) or active + freesurfer environment. + Parameters + ---------- + target_subj : str + freesurfer subject to which to map + + Other Parameters + ---------------- + kwargs map to `cortex.freesurfer.vertex_to_vertex()` + """ + # Input check + if hemi not in ['lh', 'rh', 'both']: + raise ValueError("`hemi` kwarg must be 'lh', 'rh', or 'both'") + mats = db.get_mri_surf2surf_matrix(self.subject, surface_type, + hemi='both', target_subj=target_subj, fs_subj=fs_subj, + **kwargs) + new_data = [mats[0].dot(self.left), mats[1].dot(self.right)] + if hemi == 'both': + new_data = np.hstack(new_data) + elif hemi == 'lh': + new_data = np.hstack([new_data[0], np.nan * np.zeros(new_data[1].shape)]) + elif hemi == 'rh': + new_data = np.hstack([np.nan * np.zeros(new_data[0].shape), new_data[1]]) + vx = Vertex(new_data, target_subj, vmin=self.vmin, vmax=self.vmax, cmap=self.cmap) + return vx + def u(s, encoding='utf8'): try: return s.decode(encoding) diff --git a/cortex/freesurfer.py b/cortex/freesurfer.py index a43fea4d6..8086b62c0 100644 --- a/cortex/freesurfer.py +++ b/cortex/freesurfer.py @@ -1,5 +1,6 @@ """Contains functions for interfacing with freesurfer """ +from __future__ import print_function import os import copy import shutil @@ -11,6 +12,13 @@ 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 @@ -248,11 +256,11 @@ def import_flat(subject, patch, hemis=['lh', 'rh'], sname=None, cache = os.path.join(database.default_filestore, sname, "cache") shutil.rmtree(cache) os.makedirs(cache) - # clear config-specified cache - from .options import config - config_cache = os.path.expanduser(os.path.join(config.get('basic', 'cache'), sname, 'cache')) - shutil.rmtree(config_cache) - os.makedirs(config_cache) + # clear config-specified cache, if different + config_cache = database.db.get_cache(sname) + if config_cache != cache: + shutil.rmtree(config_cache) + os.makedirs(config_cache) def _remove_disconnected_polys(polys): @@ -530,6 +538,206 @@ def get_label(subject, label, fs_subject=None, fs_dir=None, src_subject='fsavera return idx, values +def _mri_surf2surf_command(src_subj, trg_subj, input_file, output_file, hemi): + # mri_surf2surf --srcsubject --srcsurfval + # --trgsubject --trgsurfval --hemi + + 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") + 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) + 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(tf_out.name) + output_data = np.array([da.data 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 + + def get_curv(subject, hemi, type='wm', freesurfer_subject_dir=None): """Load freesurfer curv file """ diff --git a/cortex/segment.py b/cortex/segment.py index 2fa08a1fa..ba3b18539 100644 --- a/cortex/segment.py +++ b/cortex/segment.py @@ -143,6 +143,15 @@ def cut_surface(cx_subject, hemi, name='flatten', fs_subject=None, data=None, fs_subject = cx_subject opts = "[hemi=%s,name=%s]"%(hemi, name) fname = db.get_paths(cx_subject)['anats'].format(type='cutsurf', opts=opts, ext='blend') + # Double-check that fiducial and inflated vertex counts match + # (these may not match if a subject is initially imported from freesurfer to pycortex, + # and then edited further for a better segmentation and not re-imported) + ipt, ipoly, inrm = freesurfer.get_surf(fs_subject, hemi, 'inflated') + fpt, fpoly, fnrm = freesurfer.get_surf(fs_subject, hemi, 'fiducial') + if ipt.shape[0] != fpt.shape[0]: + raise ValueError("Please re-import subject - fiducial and inflated vertex counts don't match!") + else: + print('Vert check ok!') if not os.path.exists(fname): blender.fs_cut(fname, fs_subject, hemi, freesurfer_subject_dir) # Add localizer data to facilitate cutting diff --git a/cortex/utils.py b/cortex/utils.py index 5514cf489..918ec38bd 100644 --- a/cortex/utils.py +++ b/cortex/utils.py @@ -2,6 +2,7 @@ """ import io import os +import h5py import copy import binascii import warnings @@ -862,6 +863,61 @@ def get_shared_voxels(subject, xfmname, hemi="both", merge=True, use_astar=True) else: return tuple(out) + +def load_sparse_array(fname, varname): + """Load a numpy sparse array from an hdf file + + Parameters + ---------- + fname: string + file name containing array to be loaded + varname: string + name of variable to be loaded + + Notes + ----- + This function relies on variables being stored with specific naming + 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]) + return sparsemat + + +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 + + Parameters + ---------- + fname : string + file name to save + data : sparse array + data to save + varname : string + name of variable to save + mode : string + 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') + # Save data + hf.create_dataset(varname + '_data', data=data_.data, compression='gzip') + # Save indptr + hf.create_dataset(varname + '_indptr', data=data_.indptr, compression='gzip') + # Save shape + hf.create_dataset(varname + '_shape', data=data_.shape, compression='gzip') + + def get_cmap(name): """Gets a colormaps diff --git a/cortex/xfm.py b/cortex/xfm.py index 380e2d299..c829c39a4 100644 --- a/cortex/xfm.py +++ b/cortex/xfm.py @@ -300,7 +300,7 @@ def to_freesurfer(self, fs_register, subject, freesurfer_subject_dir=None): """Converts a pycortex transform to a FreeSurfer transform. Converts a transform stored in pycortex xfm object to the FreeSurfer format - (i.e., register.dat format) + (i.e., register.dat format: https://surfer.nmr.mgh.harvard.edu/fswiki/RegisterDat) Parameters ---------- @@ -310,7 +310,7 @@ def to_freesurfer(self, fs_register, subject, freesurfer_subject_dir=None): FreeSurfer subject name from which the pycortex subject was imported freesurfer_subject_dir : str | None - Directory of FreeSurfer subjects. Defaults to the value for + Directory of FreeSurfer subjects. If None, defaults to the value for the environment variable 'SUBJECTS_DIR' (which should be set by freesurfer) @@ -319,32 +319,18 @@ def to_freesurfer(self, fs_register, subject, freesurfer_subject_dir=None): import subprocess import nibabel import numpy.linalg as npl + from .database import db inv = npl.inv - # Set FreeSurfer subject directory - if freesurfer_subject_dir is None: - freesurfer_subject_dir = os.environ['SUBJECTS_DIR'] - # Set path to the anatomical volume for the FreeSurfer subject - anat_mgz = os.path.join(freesurfer_subject_dir, subject, 'mri', 'orig.mgz') - - # Write out the functional volume as Nifti file - tmp = tempfile.mkstemp(suffix='.nii') - func_nii = tmp[1] - nibabel.save(self.reference, func_nii) + anat = db.get_anat(subject, type='raw') # Read vox2ras transform for the anatomical volume - try: - cmd = ('mri_info', '--vox2ras', anat_mgz) - L = decode(subprocess.check_output(cmd)).splitlines() - anat_vox2ras = np.array([[np.float(s) for s in ll.split() if s] for ll in L]) - except OSError: - print ("Error occured while executing:\n{}".format(' '.join(cmd))) - raise + anat_vox2ras = anat.affine # Read tkrvox2ras transform for the anatomical volume try: - cmd = ('mri_info', '--vox2ras-tkr', anat_mgz) + cmd = ('mri_info', '--vox2ras-tkr', anat.get_filename()) L = decode(subprocess.check_output(cmd)).splitlines() anat_tkrvox2ras = np.array([[np.float(s) for s in ll.split() if s] for ll in L]) except OSError: @@ -353,25 +339,31 @@ def to_freesurfer(self, fs_register, subject, freesurfer_subject_dir=None): # Read tkvox2ras transform for the functional volume try: - cmd = ('mri_info', '--vox2ras-tkr', func_nii) + cmd = ('mri_info', '--vox2ras-tkr', self.reference.get_filename()) + # This next line is potentially problematic. The [1:] index + # skips a first line that is present only if an error occurs + # or some info is missing when the transform is created - i.e. + # not in all cases, just in the case that the transform is + # created exactly as it is now. Better would be to check the first + # line e.g. with a regular expression. L = decode(subprocess.check_output(cmd)).splitlines()[1:] func_tkrvox2ras = np.array([[np.float(s) for s in ll.split() if s] for ll in L]) except OSError: print ("Error occured while executing:\n{}".format(' '.join(cmd))) raise + # Debugging code + #print('Shape of func_tkrvox2ras is: [SHOULD BE (4,4) !!]') + #print(func_tkrvox2ras.shape) # Read voxel resolution of the functional volume - try: - cmd = ('mri_info', '--res', func_nii) - ll = decode(subprocess.check_output(cmd)).split("\n")[1] - func_voxres = np.array([np.float(s) for s in ll.split() if s]) - except OSError: - print ("Error occured while executing:\n{}".format(' '.join(cmd))) - raise + func_voxres = self.reference.header.get_zooms() # Calculate FreeSurfer transform fs_anat2func = np.dot(func_tkrvox2ras, np.dot(self.xfm, np.dot(anat_vox2ras, inv(anat_tkrvox2ras)))) + # Debugging code + #if not fs_anat2func.shape == (4, 4): + # raise Exception("bad assumptions led to bad transformation matrix.") # Write out to `fs_register` in register.dat format with open(fs_register, 'w') as fid: fid.write('{}\n'.format(subject)) @@ -380,8 +372,8 @@ def to_freesurfer(self, fs_register, subject, freesurfer_subject_dir=None): fid.write('0.150000\n') for row in fs_anat2func: fid.write(' '.join(['{:.15e}'.format(x) for x in row]) + '\n') - - os.remove(func_nii) + print('Wrote:') + subprocess.call(('cat', fs_register)) return fs_anat2func