Skip to content
Permalink
Browse files

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)
  • Loading branch information...
MShinkle authored and marklescroart committed Sep 27, 2019
1 parent 4a6e4ae commit 3ad1939c2185eae4d99d07e7e7b5a8ef358e9098
Showing with 482 additions and 66 deletions.
  1. +70 −30 cortex/align.py
  2. +5 −0 cortex/blender/__init__.py
  3. +66 −1 cortex/database.py
  4. +6 −0 cortex/dataset/braindata.py
  5. +35 −0 cortex/dataset/views.py
  6. +213 −5 cortex/freesurfer.py
  7. +9 −0 cortex/segment.py
  8. +56 −0 cortex/utils.py
  9. +22 −30 cortex/xfm.py
@@ -82,14 +82,25 @@ 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
supplied as `reference`. This image will be copied into the database.
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


@@ -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)
@@ -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
@@ -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):
@@ -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)

0 comments on commit 3ad1939

Please sign in to comment.
You can’t perform that action at this time.