"""
Contains functions for computing various surface properties. Mostly wrappers
for functions in `polyutils.Surface` and `polyutils.Distortion`.
"""
import os
import shlex
import shutil
import tempfile
import subprocess as sp
import numpy as np
from . import utils
from . import polyutils
from .database import db
from .xfm import Transform
[docs]
def curvature(outfile, subject, smooth=20):
"""
Compute smoothed mean curvature of the fiducial surface for the given
subject and save it to `outfile`.
Parameters
----------
outfile : str
Path where the curvature map will be saved as an npz file.
subject : str
Subject in the pycortex database for whom curvature will be computed.
smooth : float, optional
Amount of smoothing to apply to the curvature map. Default 20.
"""
curvs = []
for pts, polys in db.get_surf(subject, "fiducial"):
surf = polyutils.Surface(pts, polys)
curv = surf.smooth(surf.mean_curvature(), smooth)
curvs.append(curv)
np.savez(outfile, left=curvs[0], right=curvs[1])
[docs]
def distortion(outfile, subject, dist_type='areal', smooth=20):
"""
Compute distortion of flatmap relative to fiducial surface and save it
at `outfile`. Several different types of distortion are available:
'areal': computes the areal distortion for each triangle in the flatmap, defined as the
log ratio of the area in the fiducial mesh to the area in the flat mesh. Returns
a per-vertex value that is the average of the neighboring triangles.
See: http://brainvis.wustl.edu/wiki/index.php/Caret:Operations/Morphing
'metric': computes the linear distortion for each vertex in the flatmap, defined as
the mean squared difference between distances in the fiducial map and distances in
the flatmap, for each pair of neighboring vertices. See Fishl, Sereno, and Dale, 1999.
Parameters
----------
outfile : str
Path where the distortion map will be saved as an npz file.
subject : str
Subject in the pycortex database for whom distortion will be computed.
dist_type : ['areal', 'metric'], optional
Type of distortion to compute. Default 'areal'.
smooth : float, optional
Amount of smoothing to apply to the distortion map before returning.
Default 20.
"""
distortions = []
for hem in ["lh", "rh"]:
fidvert, fidtri = db.get_surf(subject, "fiducial", hem)
flatvert, flattri = db.get_surf(subject, "flat", hem)
surf = polyutils.Surface(fidvert, fidtri)
dist = getattr(polyutils.Distortion(flatvert, fidvert, flattri), dist_type)
smdist = surf.smooth(dist, smooth)
distortions.append(smdist)
np.savez(outfile, left=distortions[0], right=distortions[1])
[docs]
def thickness(outfile, subject):
"""
Compute cortical thickness as the distance between corresponding pial and
white matter vertices for the given subject. Note that this is slightly
different than the method used by Freesurfer, and will yield ever-so-slightly
different results.
Parameters
----------
outfile : str
Path where the thickness map will be saved.
subject : str
Subject in the pycortex database for whom cortical thickness will be
computed.
"""
pl, pr = db.get_surf(subject, "pia")
wl, wr = db.get_surf(subject, "wm")
left = np.sqrt(((pl[0] - wl[0])**2).sum(1))
right = np.sqrt(((pr[0] - wr[0])**2).sum(1))
np.savez(outfile, left=left, right=right)
[docs]
def tissots_indicatrix(outfile, sub, radius=10, spacing=50):
"""
Compute a Tissot's indicatrix for the given subject and save the result to
a file. This involves randomly filling in discs of fixed geodesic radius
on the fiducial surface.
See https://en.wikipedia.org/wiki/Tissot's_indicatrix for more info.
Parameters
----------
outfile : str
Path where the indicatrix map will be saved.
sub : str
Subject in the pycortex database for whom the indicatrix will be
computed.
radius : float, optional
The geodesic radius of each disc in mm. Default 10.
spacing : float, optional
The minimum distance between disc centers in mm. Default 50.
"""
tissots = []
allcenters = []
for hem in ["lh", "rh"]:
fidpts, fidpolys = db.get_surf(sub, "fiducial", hem)
#G = make_surface_graph(fidtri)
surf = polyutils.Surface(fidpts, fidpolys)
nvert = fidpts.shape[0]
tissot_array = np.zeros((nvert,))
centers = [np.random.randint(nvert)]
cdists = [surf.geodesic_distance(centers)]
while True:
## Find possible vertices
mcdist = np.vstack(cdists).min(0)
possverts = np.nonzero(mcdist > spacing)[0]
#possverts = np.nonzero(surf.geodesic_distance(centers) > spacing)[0]
if not len(possverts):
break
## Pick random vertex
centervert = possverts[np.random.randint(len(possverts))]
centers.append(centervert)
print("Adding vertex %d.." % centervert)
dists = surf.geodesic_distance([centervert])
cdists.append(dists)
## Find appropriate set of vertices
selverts = dists < radius
tissot_array[selverts] = 1
tissots.append(tissot_array)
allcenters.append(np.array(centers))
# make an array of objects to allow different lengths for each hemisphere
allcenters = np.array(allcenters, dtype="object")
np.savez(outfile, left=tissots[0], right=tissots[1], centers=allcenters)
[docs]
def flat_border(outfile, subject):
"""
Compute the boundary of the flatmap as a set of line segments, split into
segments that lie along the medial wall and segments that don't, and
save the result to `outfile`.
Parameters
----------
outfile : str
Path where the border map will be saved as an npz file, with `lines`
(a list of boundary line segments, each an array of 2D points) and
`ismwalls` (a boolean flag for each segment indicating whether it
lies along the medial wall) arrays.
subject : str
Subject in the pycortex database for whom the flatmap border will be
computed.
"""
flatpts, flatpolys = db.get_surf(subject, "flat", merge=True, nudge=True)
flatpolyset = set([tuple(x) for x in flatpolys])
fidpts, fidpolys = db.get_surf(subject, "fiducial", merge=True, nudge=True)
fidpolyset = set([tuple(x) for x in fidpolys])
fidonlypolys = fidpolyset - flatpolyset
fidonlypolyverts = np.unique(np.array(list(fidonlypolys)).ravel())
fidonlyverts = np.setdiff1d(fidpolys.ravel(), flatpolys.ravel())
import networkx as nx
def iter_surfedges(tris):
for a,b,c in tris:
yield a,b
yield b,c
yield a,c
def make_surface_graph(tris):
graph = nx.Graph()
graph.add_edges_from(iter_surfedges(tris))
return graph
bounds = [p for p in polyutils.trace_poly(polyutils.boundary_edges(flatpolys))]
allbounds = np.hstack(bounds)
g = make_surface_graph(fidonlypolys)
fog = g.subgraph(fidonlyverts)
badverts = np.array([v for v,d in fog.degree().items() if d<2])
g.remove_nodes_from(badverts)
fog.remove_nodes_from(badverts)
mwallset = set.union(*(set(g[v]) for v in fog.nodes())) & set(allbounds)
#cutset = (set(g.nodes()) - mwallset) & set(allbounds)
mwallbounds = [np.isin(b, mwallset) for b in bounds]
changes = [np.nonzero(np.diff(b.astype(float))!=0)[0]+1 for b in mwallbounds]
#splitbounds = [np.split(b, c) for b,c in zip(bounds, changes)]
splitbounds = []
for b,c in zip(bounds, changes):
sb = []
rb = [b[-1]] + b
rc = [1] + (c + 1).tolist() + [len(b)]
for ii in range(len(rc)-1):
sb.append(rb[rc[ii]-1 : rc[ii+1]])
splitbounds.append(sb)
ismwall = [[s.mean()>0.5 for s in np.split(mwb, c)] for mwb,c in zip(mwallbounds, changes)]
# FIXME: `height` is undefined -- this raises NameError as written.
# It was presumably meant to be a parameter (see the commented-out
# Image.new(...) call below, which used it as the output image height
# in pixels). Needs a reviewer decision on the right fix/default.
aspect = (height / (flatpts.max(0) - flatpts.min(0))[1])
lpts = (flatpts - flatpts.min(0)) * aspect
rpts = (flatpts - flatpts.min(0)) * aspect
#im = Image.new('RGBA', (int(aspect * (flatpts.max(0) - flatpts.min(0))[0]), height))
#draw = ImageDraw.Draw(im)
ismwalls = []
lines = []
for bnds, mw, pts in zip(splitbounds, ismwall, [lpts, rpts]):
for pbnd, pmw in zip(bnds, mw):
#color = {True:(0,0,255,255), False:(255,0,0,255)}[pmw]
#draw.line(pts[pbnd,:2].ravel().tolist(), fill=color, width=2)
ismwalls.append(pmw)
lines.append(pts[pbnd,:2])
np.savez(outfile, lines=lines, ismwalls=ismwalls)
[docs]
def bumpy_flatmap(outfile, subject, **kwargs):
"""Give the flatmap the relief of the cortical slab, and cache it.
Hemispheres with no flat surface get an array of zeros rather than an error.
Parameters
----------
outfile : str
Path where the offsets will be saved as an npz file.
subject : str
Subject in the pycortex database whose flatmap gets the relief.
**kwargs
Passed to `cortex.polyutils.FlatSlab`.
Notes
-----
Stored under ``bump_left`` and ``bump_right`` rather than ``left`` and
``right`` on purpose: `get_surfinfo` turns a file with the latter into a
`Vertex` by concatenating them, which assumes one value per vertex and
would mangle these three-component offsets.
"""
offsets = []
for hemi in ("lh", "rh"):
wm, _ = db.get_surf(subject, "wm", hemi)
pia, _ = db.get_surf(subject, "pia", hemi)
try:
flat, flatpolys = db.get_surf(subject, "flat", hemi)
except IOError:
offsets.append(np.zeros_like(wm))
continue
offsets.append(polyutils.FlatSlab(flat, wm, pia, flatpolys,
**kwargs).relaxed)
# Compressed, and single precision: these are millimetre offsets that end
# up in a float32 vertex attribute anyway, so the extra digits are only
# taking up room in the filestore.
np.savez_compressed(outfile,
bump_left=offsets[0].astype(np.float32),
bump_right=offsets[1].astype(np.float32))
[docs]
def equivolume_areas(outfile, subject, smooth=1.0):
"""
Compute smoothed vertex areas on the white matter and pial surfaces.
These are what the webgl viewer's equivolume depth sampling needs in order to
turn a requested volume fraction through the cortical sheet into a position
between the two surfaces. They used to be recomputed in javascript on every
viewer load with a uniform "umbrella" smoothing; computing them here instead
means they are cached, and lets them be smoothed with the cotangent-weighted
operator in `cortex.polyutils.Surface.smooth`, which respects the varying
size of the triangles rather than treating every neighbour equally.
Parameters
----------
outfile : str
Path where the areas will be saved as an npz file.
subject : str
Subject in the pycortex database for whom the areas will be computed.
smooth : float, optional
Amount of smoothing to apply. Default 1.0. Pass 0 for the raw
barycentric vertex areas.
Notes
-----
Stored under ``wm_left`` / ``wm_right`` / ``pia_left`` / ``pia_right``; see
the note in `bumpy_flatmap` for why ``left`` and ``right`` are avoided.
Changing the smoothing operator does move the depths the viewer samples at,
slightly: against the five umbrella iterations the javascript used, the
depth equivolume sampling picks for a requested fraction of 0.5 shifts by
about 0.01 of the cortical thickness at the median on S1. The tail is larger,
but it is concentrated where the two areas are nearly equal and the depth is
poorly determined either way.
"""
areas = dict()
for hemi, side in zip(["lh", "rh"], ["left", "right"]):
for name in ["wm", "pia"]:
pts, polys = db.get_surf(subject, name, hemi)
surf = polyutils.Surface(pts, polys)
# The lumped mass matrix of the Laplace-Beltrami operator is exactly
# the barycentric vertex area, i.e. a third of each incident face.
_, vertex_area, _, _ = surf.laplace_operator
areas["%s_%s" % (name, side)] = surf.smooth(vertex_area, smooth)
np.savez(outfile, **areas)