merge

.gitignore

@@ -34,3 +34,9 @@ nosetests.xml
 .mr.developer.cfg
 .project
 .pydevproject
+
+# OS X bullshit
+.DS_store
+
+# vim temp files
+*~

cortex/blender.py

@@ -0,0 +1,127 @@
+import os
+import struct
+import numpy as np
+from matplotlib import cm, colors
+
+import bpy.ops
+from bpy import context as C
+from bpy import data as D
+
+def make_object(pts, polys, name="mesh"):
+    mesh = D.meshes.new(name)
+    mesh.from_pydata(pts.tolist(), [], polys.tolist())
+    obj = D.objects.new(name, mesh)
+    C.scene.objects.link(obj)
+    return obj, mesh
+
+def add_vcolor(color, mesh=None, name='color', cmap=cm.RdBu, vmin=None, vmax=None):
+    if mesh is None:
+        mesh = C.scene.objects.active.data
+    elif isinstance(mesh, str):
+        mesh = D.meshes[mesh]
+
+    bpy.ops.object.mode_set(mode='OBJECT')
+    if color.ndim == 1:
+        if vmin is None:
+            vmin = color.min()
+        if vmax is None:
+            vmax = color.max()
+        color = cmap((color - vmin) / (vmax - vmin))[:,:3]
+
+    loopidx = np.zeros((len(mesh.loops),), dtype=np.uint32)
+    mesh.loops.foreach_get('vertex_index', loopidx)
+
+    vcolor = mesh.vertex_colors.new(name)
+    for i, j in enumerate(loopidx):
+        vcolor.data[i].color = list(color[j])
+    return vcolor
+
+def add_shapekey(shape, name=None):
+    bpy.ops.object.shape_key_add()
+    key = D.shape_keys[-1].key_blocks[-1]
+    if name is not None:
+        key.name = name
+
+    for i in range(len(key.data)):
+        key.data[i].co = shape[i]
+    return key
+
+def cut_data(volumedata, name="retinotopy", projection="nearest", cmap=cm.RdBu, vmin=None, vmax=None, mesh="hemi"):
+    if isinstance(mesh, str):
+        mesh = D.meshes[mesh]
+
+    mapped = volumedata.map(projection)
+    if mapped.llen == len(mesh.vertices):
+        print("left hemisphere")
+        vcolor = mapped.left
+    else:
+        print ("right hemisphere")
+        vcolor = mapped.right
+
+    return add_vcolor(vcolor, mesh=mesh, name=name, cmap=cmap, vmin=vmin, vmax=vmax)
+
+def fs_cut(subject, hemi):
+    from .freesurfer import get_surf
+    wpts, polys, curv = get_surf(subject, hemi, 'smoothwm')
+    ipts, _, _ = get_surf(subject, hemi, 'inflated')
+
+    obj, mesh = make_object(wpts, polys, name='hemi')
+    obj.scale = .1, .1, .1
+    C.scene.objects.active = obj
+    bpy.ops.object.shape_key_add()
+    add_vcolor(curv, mesh, vmin=-.6, vmax=.6, name='curvature')
+    add_shapekey(ipts, name='inflated')
+    obj.use_shape_key_edit_mode = True
+    return mesh
+
+def write_patch(subject, hemi, name, mesh='hemi'):
+    from .freesurfer import get_paths, write_patch
+    if isinstance(mesh, str):
+        mesh = D.meshes[mesh]
+
+    bpy.ops.object.mode_set(mode='OBJECT')
+    bpy.ops.object.mode_set(mode='EDIT')
+    bpy.ops.mesh.select_all(action='DESELECT')
+    C.tool_settings.mesh_select_mode = False, True, False
+    bpy.ops.mesh.select_non_manifold()
+    bpy.ops.object.mode_set(mode='OBJECT')
+
+    mwall_edge = set()
+    for edge in mesh.edges:
+        if edge.select:
+            mwall_edge.add(edge.vertices[0])
+            mwall_edge.add(edge.vertices[1])
+
+    bpy.ops.object.mode_set(mode='EDIT') 
+    C.tool_settings.mesh_select_mode = True, False, False
+    bpy.ops.mesh.select_all(action='DESELECT')
+    bpy.ops.object.mode_set(mode='OBJECT')
+    seam = set()
+    for edge in mesh.edges:
+        if edge.use_seam:
+            seam.add(edge.vertices[0])
+            seam.add(edge.vertices[1])
+            edge.select = True
+
+    bpy.ops.object.mode_set(mode='EDIT') 
+    bpy.ops.mesh.select_more()
+    bpy.ops.object.mode_set(mode='OBJECT')
+    smore = set()
+    for i, vert in enumerate(mesh.vertices):
+        if vert.select:
+            smore.add(i)
+
+    bpy.ops.object.mode_set(mode='EDIT') 
+    bpy.ops.mesh.select_all(action='DESELECT')
+    bpy.ops.object.mode_set(mode='OBJECT')
+
+    fverts = set()
+    for face in mesh.polygons:
+        fverts.add(face.vertices[0])
+        fverts.add(face.vertices[1])
+        fverts.add(face.vertices[2])
+
+    edges = mwall_edge | (smore - seam)
+    verts = fverts - seam
+    pts = [(v, D.shape_keys['Key'].key_blocks['inflated'].data[v].co) for v in verts]
+    write_patch(get_paths(subject, hemi).format(name=name), pts, edges)

cortex/blender/__init__.py

@@ -71,6 +71,39 @@ def add_cutdata(fname, dataview, name="retinotopy", projection="nearest", mesh="
 
     return 
 
+
+def gii_cut(fname, subject, hemi):
+    '''
+    Add gifti surface to blender
+    '''
+    from ..database import db
+    hemis = dict(lh='left',
+                 rh='right')
+
+    wpts, polys = db.get_surf(subject, 'wm', hemi)
+    ipts, _ = db.get_surf(subject, 'very_inflated', hemi)
+    curvature = db.getSurfInfo(subject, 'curvature')
+    rcurv = curvature.__getattribute__(hemis[hemi])
+
+    p = xdrlib.Packer()
+    p.pack_array(wpts.ravel(), p.pack_double)
+    p.pack_array(ipts.ravel(), p.pack_double)
+    p.pack_array(polys.ravel(), p.pack_uint)
+    p.pack_array(rcurv.ravel(), p.pack_double)
+    with tempfile.NamedTemporaryFile() as tf:
+        tf.write(p.get_buffer())
+        tf.flush()
+        code = """with open('{tfname}', 'rb') as fp:
+            u = xdrlib.Unpacker(fp.read())
+            wpts = u.unpack_array(u.unpack_double)
+            ipts = u.unpack_array(u.unpack_double)
+            polys = u.unpack_array(u.unpack_uint)
+            curv = u.unpack_array(u.unpack_double)
+            blendlib.init_subject(wpts, ipts, polys, curv)
+        """.format(tfname=tf.name)
+        _call_blender(fname, code)
+
+
 def fs_cut(fname, subject, hemi):
     wpts, polys, curv = freesurfer.get_surf(subject, hemi, 'smoothwm')
     ipts, _, _ = freesurfer.get_surf(subject, hemi, 'inflated')
@@ -109,4 +142,5 @@ def write_patch(bname, pname, mesh="hemi"):
             mesh = u.unpack_string().decode('utf-8')
             blendlib.save_patch(pname, mesh)
         """.format(tfname=tf.name)
-        _call_blender(bname, code)
+        _call_blender(bname, code)
+

cortex/brainctm.py

@@ -89,13 +89,23 @@ def addCurvature(self, **kwargs):
             self.left.aux[:,1] = npz.left
             self.right.aux[:,1] = npz.right
 
-    def save(self, path, method='mg2', **kwargs):
+    def save(self, path, method='mg2', disp_layers=['rois'], extra_disp=None, **kwargs):
+        """Save CTM file for static html display. 
+
+        Parameters
+        ----------
+        path : string filepath
+        method : idkwtf
+        disp_layers : tuple|list
+            list of strings; names of layers from rois.svg to show
+        extra_disp : tuple
+            (svgfile,[layers]) - tuple of (external display .svg filename, [list_of,layers_in_file,to_display])
+        """
         ctmname = path+".ctm"
         svgname = path+".svg"
         jsname = path+".json"
-        #ptmapname = path+".npz"
 
-        ##### Save CTM concatenation
+        # Save CTM concatenation
         (lpts, _, _), lbin = self.left.save(method=method, **kwargs)
         (rpts, _, _), rbin = self.right.save(method=method, **kwargs)
 
@@ -105,50 +115,62 @@ def save(self, path, method='mg2', **kwargs):
             offsets.append(fp.tell())
             fp.write(rbin)
 
-        ##### Save the JSON descriptor
-        jsdict = dict(rois=os.path.split(svgname)[1], data=os.path.split(ctmname)[1], names=self.types, 
-            materials=[], offsets=offsets)
+        # Save the JSON descriptor | Need to add to this for extra_disp?
+        jsdict = dict(rois=os.path.split(svgname)[1],
+                      data=os.path.split(ctmname)[1],
+                      names=self.types, 
+                      materials=[],
+                      offsets=offsets)
         if self.flatlims is not None:
             jsdict['flatlims'] = self.flatlims
         json.dump(jsdict, open(jsname, 'w'))
 
-        ##### Compute and save the index map
+        # Compute and save the index map
         if method != 'raw':
             ptmap, inverse = [], []
             for hemi, pts in zip([self.left, self.right], [lpts, rpts]):
                 kdt = cKDTree(hemi.pts)
                 diff, idx = kdt.query(pts)
                 ptmap.append(idx)
                 inverse.append(idx.argsort())
-
-            # np.savez(ptmapname, left=ptmap[0], right=ptmap[1])
         else:
             ptmap = inverse = np.arange(len(self.left.ctm)), np.arange(len(self.right.ctm))
 
-        ##### Save the SVG with remapped indices
+        # Save the SVG with remapped indices (map 2D flatmap locations to vertices)
         if self.left.flat is not None:
+            # add sulci & display layers
             flatpts = np.vstack([self.left.flat, self.right.flat])
-            roipack = db.get_overlay(self.subject, pts=flatpts)
-            layer = roipack.setup_labels()
+            roipack = db.get_overlay(self.subject, pts=flatpts, otype=disp_layers)
+            # optionally add extra display layers
+            if not extra_disp is None:
+                esvgfile,elayerlist = extra_disp
+                eroipack = db.get_overlay(self.subject,pts=flatpts, otype='external',
+                    svgfile=esvgfile,layer=elayerlist)
+                roipack = roipack + eroipack
+            layers = roipack.setup_labels()
+            if not isinstance(layers, (tuple, list)):
+                layers = (layers,)
+
+            # assign coordinates in left hemisphere negative values
             with open(svgname, "w") as fp:
-                for element in layer.findall(".//{http://www.w3.org/2000/svg}text"):
-                    idx = int(element.attrib["data-ptidx"])
-                    if idx < len(inverse[0]):
-                        idx = inverse[0][idx]
-                    else:
-                        idx -= len(inverse[0])
-                        idx = inverse[1][idx] + len(inverse[0])
-                    element.attrib["data-ptidx"] = str(idx)
+                for layer in layers:
+                    for element in layer.findall(".//{http://www.w3.org/2000/svg}text"):
+                        idx = int(element.attrib["data-ptidx"])
+                        if idx < len(inverse[0]):
+                            idx = inverse[0][idx]
+                        else:
+                            idx -= len(inverse[0])
+                            idx = inverse[1][idx] + len(inverse[0])
+                        element.attrib["data-ptidx"] = str(idx)
                 fp.write(roipack.toxml())
-
         return ptmap
 
 class Hemi(object):
     def __init__(self, pts, polys, norms=None):
         self.tf = tempfile.NamedTemporaryFile()
         self.ctm = CTMfile(self.tf.name, "w")
 
-        self.ctm.setMesh(pts.astype(np.float32), polys, norms=norms)
+        self.ctm.setMesh(pts.astype(np.float32), polys.astype(np.uint32), norms=norms)
 
         self.pts = pts
         self.polys = polys
@@ -216,7 +238,10 @@ def setFlat(self, pts):
     def addSurf(self, pts, **kwargs):
         super(DecimatedHemi, self).addSurf(pts[self.mask], **kwargs)
 
-def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, decimate=False):
+def make_pack(outfile, subj, types=("inflated",), method='raw', level=0,
+              decimate=False, disp_layers=['rois'],extra_disp=None):
+    """Generates a cached CTM file"""
+
     ctm = BrainCTM(subj, decimate=decimate)
     ctm.addCurvature()
     for name in types:
@@ -225,7 +250,11 @@ def make_pack(outfile, subj, types=("inflated",), method='raw', level=0, decimat
     if not os.path.exists(os.path.split(outfile)[0]):
         os.makedirs(os.path.split(outfile)[0])
 
-    return ctm.save(os.path.splitext(outfile)[0], method=method, level=level)
+    return ctm.save(os.path.splitext(outfile)[0],
+                    method=method,
+                    level=level,
+                    disp_layers=disp_layers,
+                    extra_disp=extra_disp)
 
 def read_pack(ctmfile):
     fname = os.path.splitext(ctmfile)[0]
@@ -247,4 +276,4 @@ def read_pack(ctmfile):
             pts, polys, norms = ctm.getMesh()
             meshes.append((pts, polys))
 
-    return meshes
+    return meshes