Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

bumpy flatmap fixes #310

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

bumpy flatmap fixes #310

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2bc274b
fixed left hem smoothing, tweaked lights
alexhuth Jan 25, 2019
30c03a5
resurrected the ROIpack
alexhuth Aug 8, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
89 changes: 70 additions & 19 deletions cortex/rois.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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"):
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)

17 changes: 9 additions & 8 deletions cortex/webgl/resources/js/axes3d.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = [];

Expand Down
18 changes: 11 additions & 7 deletions cortex/webgl/resources/js/mriview_surface.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);
Expand All @@ -179,18 +187,14 @@ 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;});
// flatareas = module.iterativelySmoothVertexData(hemi.attributes.position, hemi.attributes.index, hemi.offsets, flatareas, smoothfactor, smoothiter);
// 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;
Expand Down