From 2bc274bbadff4567fe759da6a4017f9498140bb1 Mon Sep 17 00:00:00 2001 From: alexhuth Date: Fri, 25 Jan 2019 07:41:48 -0800 Subject: [PATCH 1/2] fixed left hem smoothing, tweaked lights --- cortex/webgl/resources/js/axes3d.js | 17 +++++++++-------- cortex/webgl/resources/js/mriview_surface.js | 18 +++++++++++------- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/cortex/webgl/resources/js/axes3d.js b/cortex/webgl/resources/js/axes3d.js index 89ab251d9..a1b61eba2 100644 --- a/cortex/webgl/resources/js/axes3d.js +++ b/cortex/webgl/resources/js/axes3d.js @@ -22,16 +22,17 @@ var jsplot = (function (module) { //These lights approximately match what's done by vtk this.lights = [ - new THREE.DirectionalLight(0xffffff, .47), - new THREE.DirectionalLight(0xffffff, .29), - new THREE.DirectionalLight(0xffffff, .24) + //new THREE.DirectionalLight(0xffffff, 1.0), + new THREE.DirectionalLight(0xffffff, 1.0), + // new THREE.DirectionalLight(0xffffff, .29), + // new THREE.DirectionalLight(0xffffff, .24) ]; - this.lights[0].position.set( 1, -1, -1 ).normalize(); - this.lights[1].position.set( -1, -.25, .75 ).normalize(); - this.lights[2].position.set( 1, -.25, .75 ).normalize(); + this.lights[0].position.set( -100, 200, 10 ); + // this.lights[1].position.set( -1, -.25, .75 ).normalize(); + // this.lights[2].position.set( 1, -.25, .75 ).normalize(); this.camera.add( this.lights[0] ); - this.camera.add( this.lights[1] ); - this.camera.add( this.lights[2] ); + // this.camera.add( this.lights[1] ); + // this.camera.add( this.lights[2] ); this.views = []; diff --git a/cortex/webgl/resources/js/mriview_surface.js b/cortex/webgl/resources/js/mriview_surface.js index de14d8c2f..2987fb9fa 100644 --- a/cortex/webgl/resources/js/mriview_surface.js +++ b/cortex/webgl/resources/js/mriview_surface.js @@ -131,6 +131,14 @@ var mriview = (function(module) { left :{positions:[], normals:[]}, right:{positions:[], normals:[]}, }; + + + var areasmoothfactor = 0.1; + var areasmoothiter = 5; + + var distsmoothfactor = 0.1; + var distsmoothiter = 20; + for (var name in names) { var hemi = geometries[names[name]]; posdata[name].map = hemi.indexMap; @@ -155,11 +163,11 @@ var mriview = (function(module) { } //Setup flatmap mix var wmareas = module.computeAreas(hemi.attributes.wm, hemi.attributes.index, hemi.offsets); - wmareas = module.iterativelySmoothVertexData(hemi.attributes.wm, hemi.attributes.index, hemi.offsets, wmareas, smoothfactor, smoothiter); + wmareas = module.iterativelySmoothVertexData(hemi.attributes.wm, hemi.attributes.index, hemi.offsets, wmareas, areasmoothfactor, areasmoothiter); hemi.wmareas = wmareas; var pialareas = module.computeAreas(hemi.attributes.position, hemi.attributes.index, hemi.offsets); - pialareas = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, pialareas, smoothfactor, smoothiter); + pialareas = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, pialareas, areasmoothfactor, areasmoothiter); hemi.pialareas = pialareas; var pialarea_attr = new THREE.BufferAttribute(pialareas, 1); @@ -179,10 +187,6 @@ var mriview = (function(module) { posdata[name].positions.push(hemi.attributes['mixSurfs'+json.names.length]); posdata[name].normals.push(hemi.attributes['mixNorms'+json.names.length]); - - var smoothfactor = 0.1; - var smoothiter = 50; - // var flatareas = module.computeAreas(hemi.attributes.mixSurfs1, hemi.culled.index, hemi.culled.offsets); // var flatareascale = flatscale ** 2; // flatareas = flatareas.map(function (a) { return a / flatareascale;}); @@ -190,7 +194,7 @@ var mriview = (function(module) { // hemi.flatareas = flatareas; var dists = module.computeDist(hemi.attributes.position, hemi.attributes.wm); - dists = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, dists, smoothfactor, smoothiter); + dists = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, dists, distsmoothfactor, distsmoothiter); var vertexvolumes = module.computeVertexPrismVolume(wmareas, pialareas, dists); hemi.vertexvolumes = vertexvolumes; From 30c03a52f49b1bc2a952e9441931e38243d69ddc Mon Sep 17 00:00:00 2001 From: alexhuth Date: Thu, 8 Aug 2019 12:10:20 -0700 Subject: [PATCH 2/2] resurrected the ROIpack --- cortex/rois.py | 89 +++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 70 insertions(+), 19 deletions(-) diff --git a/cortex/rois.py b/cortex/rois.py index 333df71b6..30a8ed344 100644 --- a/cortex/rois.py +++ b/cortex/rois.py @@ -7,12 +7,12 @@ import networkx as nx -from db import surfs -from svgroi import get_roipack, _make_layer, _find_layer, parser +from cortex import db +from svgoverlay import get_overlay, _make_layer, _find_layer, parser from lxml import etree -from dataset import VertexData +from dataset import Vertex from polyutils import Surface, boundary_edges -from utils import get_curvature, add_roi +from utils import add_roi import quickflat class ROIpack(object): @@ -30,8 +30,8 @@ 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) + # Create basic Vertex to avoid expensive initialization.. + empty = Vertex(None, self.subject) # Load ROIs from file if self.roifile.endswith("npz"): @@ -43,7 +43,7 @@ def load_roifile(self): elif self.roifile.endswith("svg"): pts, polys = surfs.getSurf(self.subject, "flat", merge=True, nudge=True) npts = len(pts) - svgroipack = get_roipack(self.roifile, pts, polys) + svgroipack = get_overlay(self.subject, self.roifile, pts, polys) for name in svgroipack.names: roimask = np.zeros((npts,)) roimask[svgroipack.get_roi(name)] = 1 @@ -68,24 +68,25 @@ def to_svg(self, open_inkscape=False, filename=None): if filename is None: filename = tempfile.mktemp(suffix=".svg", prefix=self.subject+"-rois-") - mpts, mpolys = surfs.getSurf(self.subject, "flat", merge=True, nudge=True) + mpts, mpolys = db.get_surf(self.subject, "flat", merge=True, nudge=True) svgmpts = mpts[:,:2].copy() svgmpts -= svgmpts.min(0) svgmpts *= 1024 / svgmpts.max(0)[1] svgmpts[:,1] = 1024 - svgmpts[:,1] npts = len(mpts) - svgroipack = get_roipack(filename, mpts, mpolys) + svgroipack = get_overlay(self.subject, filename, mpts, mpolys) # Add default layers # Add curvature from matplotlib import cm - curv = VertexData(np.hstack(get_curvature(self.subject)), self.subject) - fp = io.BytesIO() - curvim = quickflat.make_png(fp, curv, height=1024, with_rois=False, with_labels=False, - with_colorbar=False, cmap=cm.gray,recache=True) - fp.seek(0) - svgroipack.add_roi("curvature", binascii.b2a_base64(fp.read()), add_path=False) + #curv = Vertex(np.hstack(get_curvature(self.subject)), self.subject) + #curv = db.get_surfinfo(self.subject, 'curvature') + #fp = io.BytesIO() + #curvim = quickflat.make_png(fp, curv, height=1024, with_rois=False, with_labels=False, + #with_colorbar=False, cmap=cm.gray,recache=True) + #fp.seek(0) + #svgroipack.add_roi("curvature", binascii.b2a_base64(fp.read()), add_path=False) # Add thickness @@ -94,12 +95,12 @@ def to_svg(self, open_inkscape=False, filename=None): svg = etree.parse(svgroipack.svgfile, parser=parser) # Find boundary vertices for each ROI - lsurf, rsurf = [Surface(*pp) for pp in surfs.getSurf(self.subject, "fiducial")] - flsurf, frsurf = [Surface(*pp) for pp in surfs.getSurf(self.subject, "flat")] + lsurf, rsurf = [Surface(*pp) for pp in db.get_surf(self.subject, "fiducial")] + flsurf, frsurf = [Surface(*pp) for pp in db.get_surf(self.subject, "flat")] valids = [set(np.unique(flsurf.polys)), set(np.unique(frsurf.polys))] # Construct polygon adjacency graph for each surface - polygraphs = [lsurf.poly_graph, rsurf.poly_graph] + polygraphs = [poly_graph(lsurf), poly_graph(rsurf)] for roi in self.rois.keys(): print("Adding %s.." % roi) masks = self.rois[roi].left, self.rois[roi].right @@ -111,7 +112,7 @@ def to_svg(self, open_inkscape=False, filename=None): continue # Find bounds - inbound, exbound = surf.get_boundary(np.nonzero(mask)[0]) + inbound, exbound = get_boundary(surf, np.nonzero(mask)[0]) # Find polys allbpolys = np.unique(surf.connected[inbound+exbound].indices) @@ -161,3 +162,53 @@ def to_svg(self, open_inkscape=False, filename=None): with open(svgroipack.svgfile, "w") as xml: xml.write(etree.tostring(svg, pretty_print=True)) + + +def poly_graph(surf): + """NetworkX undirected graph representing polygons of a Surface. + """ + import networkx as nx + from collections import defaultdict + edges = defaultdict(list) + for ii,(a,b,c) in enumerate(surf.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(surf, 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(surf.graph, vertices)) + + # Find adjacent vertices to get inner boundary + internal_boundary = set.union(*[set(surf.graph[v].keys()) + for v in external_boundary]).intersection(set(vertices)) + + if remove_danglers: + ingraph = surf.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 = surf.graph.subgraph(internal_boundary) + danglers = [n for n,d in ingraph.degree().items() if d<2] + + return list(internal_boundary), list(external_boundary) +