diff --git a/cortex/align.py b/cortex/align.py index 2a49a2435..6307ac39f 100644 --- a/cortex/align.py +++ b/cortex/align.py @@ -200,3 +200,76 @@ def autotweak(subject, xfmname): print('Saved transform as (%s, %s)'%(subject, xfmname+'_auto')) finally: shutil.rmtree(cache) + + +def automatic_bbregister(subject, xfmname, reference, noclean=False): + """Create an automatic alignment using the FLIRT boundary-based alignment (BBR) from FSL. + + If `noclean`, intermediate files will not be removed from /tmp. The `reference` image and resulting + transform called `xfmname` will be automatically stored in the database. + + It's good practice to open up this transform afterward in the manual aligner and check how it worked. + Do that using the following (with the same `subject` and `xfmname` used here, no need for `reference`): + > align.manual(subject, xfmname) + + Parameters + ---------- + subject : str + Subject identifier. + xfmname : str + String identifying the transform to be created. + reference : str + Path to a nibabel-readable image that will be used as the reference for this transform. + Usually, this is a single (3D) functional data volume. + noclean : bool, optional + If True intermediate files will not be removed from /tmp (this is useful for debugging things), + and the returned value will be the name of the temp directory. Default False. + + Returns + ------- + Nothing unless `noclean` is True. + """ + import shlex + import shutil + import tempfile + import subprocess as sp + + from .database import db + from .xfm import Transform + from .options import config + + fsl_prefix = config.get("basic", "fsl_prefix") + schfile = os.path.join(os.path.split(os.path.abspath(__file__))[0], "bbr.sch") + + retval = None + try: + cache = tempfile.mkdtemp() + absreference = os.path.abspath(reference) + raw = db.get_anat(subject, type='raw').get_filename() + bet = db.get_anat(subject, type='brainmask').get_filename() + wmseg = db.get_anat(subject, type='whitematter').get_filename() + # Compute anatomical-to-epi transform + + print('Running bbregister') + cmd = 'bbregister --s {subject} --mov {absref} --bold --init-fsl --reg {cache}/register.dat --fslmat {cache}/out.mat' + cmd = cmd.format(cache=cache, absref=absreference, subject=subject) + + if sp.call(cmd, shell=True) != 0: + raise IOError('Error calling bbregister') + + x = np.loadtxt(os.path.join(cache, "out.mat")) + # Pass transform as FROM epi TO anat; transform will be inverted + # back to anat-to-epi, standard direction for pycortex internal + # storage by from_fsl + xfm = Transform.from_fsl(x,absreference,raw) + # Save as pycortex 'coord' transform + xfm.save(subject,xfmname,'coord') + print('Success') + + finally: + if not noclean: + shutil.rmtree(cache) + else: + retval = cache + + return retval \ No newline at end of file diff --git a/cortex/polyutils.py b/cortex/polyutils.py index 4832ced3c..788038682 100644 --- a/cortex/polyutils.py +++ b/cortex/polyutils.py @@ -43,6 +43,55 @@ def __init__(self, pts, polys): self._rlfac_solvers = dict() self._nLC_solvers = dict() + @property + def poly_graph(self): + """NetworkX undirected graph representing polygons of this Surface. + """ + import networkx as nx + from collections import defaultdict + edges = defaultdict(list) + for ii,(a,b,c) in enumerate(self.polys): + edges[frozenset([a,b])].append(ii) + edges[frozenset([a,c])].append(ii) + edges[frozenset([b,c])].append(ii) + + #nedges = len(edges) + #ii,jj = np.vstack(edges.values()).T + #polymat = sparse.coo_matrix((np.ones((nedges,)), (ii, jj)), shape=[len(self.polys)]*2) + polygraph = nx.Graph() + polygraph.add_edges_from(((p[0], p[1], dict(verts=k)) for k,p in edges.iteritems())) + return polygraph + + def get_boundary(self, vertices, remove_danglers=False): + """Return interior and exterior boundary sets for `vertices`. + If `remove_danglers` is True vertices in the internal boundary with + only one neighbor in the internal boundary will be moved to the external + boundary. + """ + if not len(vertices): + return [], [] + + import networkx as nx + + # Use networkx to get external boundary + external_boundary = set(nx.node_boundary(self.graph, vertices)) + + # Find adjacent vertices to get inner boundary + internal_boundary = set.union(*[set(self.graph[v].keys()) + for v in external_boundary]).intersection(set(vertices)) + + if remove_danglers: + ingraph = self.graph.subgraph(internal_boundary) + danglers = [n for n,d in ingraph.degree().items() if d==1] + while danglers: + internal_boundary -= set(danglers) + external_boundary |= set(danglers) + + ingraph = self.graph.subgraph(internal_boundary) + danglers = [n for n,d in ingraph.degree().items() if d<2] + + return list(internal_boundary), list(external_boundary) + @property @_memo def ppts(self): diff --git a/cortex/rois.py b/cortex/rois.py index 6a762fb25..2bc55f29b 100644 --- a/cortex/rois.py +++ b/cortex/rois.py @@ -7,11 +7,12 @@ import networkx as nx from db import surfs + +from dataset import db from svgroi import get_roipack, _make_layer, _find_layer, parser from lxml import etree -from dataset import VertexData -from polyutils import Surface, boundary_edges -from utils import get_curvature, add_roi +from polyutils import Surface + import quickflat class ROIpack(object): @@ -29,9 +30,12 @@ def load_roifile(self): print("ROI file %s doesn't exist.." % self.roifile) return - # Create basic VertexData to avoid expensive initialization.. - empty = VertexData(None, self.subject) - + #Used to create basic VertexData to avoid expensive initialization.. + #empty = VertexData(None, self.subject) + + #hack just gets curvature VertexData to prevent error from old way + empty = db.get_surfinfo(self.subject) + # Load ROIs from file if self.roifile.endswith("npz"): roidata = np.load(self.roifile) @@ -79,7 +83,7 @@ def to_svg(self, open_inkscape=False, filename=None): # Add default layers # Add curvature from matplotlib import cm - curv = VertexData(np.hstack(get_curvature(self.subject)), self.subject) + curv = db.get_surfinfo(self.subject) fp = cStringIO.StringIO() curvim = quickflat.make_png(fp, curv, height=1024, with_rois=False, with_labels=False, with_colorbar=False, cmap=cm.gray,recache=True)