####################################################################################################
# neuropythy/graphics/core.py
# Core implementation of the neuropythy graphics library for plotting cortical surfaces
import numpy as np
import numpy.linalg as npla
import scipy as sp
import scipy.sparse as sps
import scipy.spatial as spspace
import neuropythy.geometry as geo
import neuropythy.mri as mri
import neuropythy.io as nyio
import os, sys, six, itertools, atexit, shutil, tempfile, warnings, pimms
from ..util import (times, zdivide, plus, minus, to_hemi_str, nanlog)
from ..vision import (visual_area_names, visual_area_numbers)
# 2D Graphics ######################################################################################
# Below are various graphics-related functions; some are wrapped in a try block in case matplotlib
# is not installed.
def _matplotlib_load_error(*args, **kwargs):
raise RuntimeError('load failure: the requested object could not be loaded, probably because ' +
'you do not have matplotlib installed correctly')
cmap_curvature = _matplotlib_load_error
cmap_polar_angle_sym = _matplotlib_load_error
cmap_polar_angle_lh = _matplotlib_load_error
cmap_polar_angle_rh = _matplotlib_load_error
cmap_polar_angle = _matplotlib_load_error
cmap_theta_sym = _matplotlib_load_error
cmap_theta_lh = _matplotlib_load_error
cmap_theta_rh = _matplotlib_load_error
cmap_theta = _matplotlib_load_error
cmap_eccentricity = _matplotlib_load_error
cmap_log_eccentricity = _matplotlib_load_error
cmap_radius = _matplotlib_load_error
cmap_log_radius = _matplotlib_load_error
cmap_cmag = _matplotlib_load_error
cmap_log_cmag = _matplotlib_load_error
cmap_eccenflat = _matplotlib_load_error
label_cmap = _matplotlib_load_error
try:
import matplotlib, matplotlib.pyplot, matplotlib.tri, matplotlib.colors
# we use this below
blend_cmap = matplotlib.colors.LinearSegmentedColormap.from_list
cmap_curvature = matplotlib.colors.LinearSegmentedColormap(
'curv',
{name: ((0.0, 0.0, 0.5), (0.5, 0.5, 0.2), (1.0, 0.2, 0.0))
for name in ['red', 'green', 'blue']})
cmap_curvature.__doc__ = '''
cmap_curvature is a colormap for plotting the curvature of a vertex; note that this colormap
assumes FreeSurfer's standard way of storing curvature where negative values indicate gyri and
positive values indicate sulci.
Values passed to cmap_curvature should be scaled such that (-1,1) -> (0,1).
'''
cmap_polar_angle_sym = blend_cmap(
'polar_angle_sym',
[(0.5,0,0), (1,1,0), (0,0.5,0), (0,1,1), (0,0,0.5), (0,1,1), (0,0.5,0), (1,1,0), (0.5,0,0)])
cmap_polar_angle_sym.__doc__ = '''
cmap_polar_angle_sym is a colormap for plotting the pRF polar angle of a vertex.
Values passed to cmap_polar_angle_sym should be scaled such that (-180,180 deg) -> (0,1).
cmap_polar_angle_sym is a circular colormap that is left-right symmetric with green representing
the horizontal meridian, blue representing the upper vertical meridian, and red representing the
lower vertical meridian.
'''
cmap_polar_angle_lh = blend_cmap(
'polar_angle_lh',
[(0.5,0,0), (0.75,0,0.5), (1,0,1), (0.5,0,0.75), (0,0,0.5), (0,1,1), (0,0.5,0), (1,1,0),
(0.5,0,0)])
cmap_polar_angle_lh.__doc__ = '''
cmap_polar_angle_lh is a colormap for plotting the pRF polar angle of a vertex.
Values passed to cmap_polar_angle_lh should be scaled such that (-180,180 deg) -> (0,1).
cmap_polar_angle_lh is a circular colormap that emphasizes colors in the right visual field; the
left visual field appears mostly magenta.
'''
cmap_polar_angle_rh = blend_cmap(
'polar_angle_rh',
[(0.5,0,0), (1,1,0), (0,0.5,0), (0,1,1), (0,0,0.5), (0.5,0,0.75), (1,0,1), (0.75,0,0.5),
(0.5,0,0)])
cmap_polar_angle_rh.__doc__ = '''
cmap_polar_angle_rh is a colormap for plotting the pRF polar angle of a vertex.
Values passed to cmap_polar_angle_rh should be scaled such that (-180,180 deg) -> (0,1).
cmap_polar_angle_rh is a circular colormap that emphasizes colors in the left visual field; the
right visual field appears mostly magenta.
'''
cmap_polar_angle = blend_cmap(
'polar_angle',
[(0.5,0,0), (1,0,1), (0,0,0.5), (0,1,1), (0,0.5,0), (1,1,0), (0.5,0,0)])
cmap_polar_angle.__doc__ = '''
cmap_polar_angle is a colormap for plotting the pRF polar angle of a vertex.
Values passed to cmap_polar_angle should be scaled such that (-180,180 deg) -> (0,1).
cmap_polar_angle is a 6-pronged circular colormap; note that it does not have dark or bright
values at the horizontal meridia as cmap_polar_angle_sym, cmap_polar_angle_lh, and
cmap_polar_angle_rh do.
'''
cmap_theta_sym = blend_cmap(
'theta_sym',
[(0.5,0,0), (1,1,0), (0,0.5,0), (0,1,1), (0,0,0.5), (0,1,1), (0,0.5,0), (1,1,0), (0.5,0,0)])
cmap_theta_sym.__doc__ = '''
cmap_theta_sym is a colormap for plotting the pRF theta of a vertex.
Values passed to cmap_theta_sym should be scaled such that (-pi,pi rad) -> (0,1). Note that 0 is
the right right horizontal meridian and positive is the counter-clockwise direction.
cmap_theta_sym is a circular colormap that is left-right symmetric with green representing
the horizontal meridian, blue representing the upper vertical meridian, and red representing the
lower vertical meridian.
'''
cmap_theta_lh = blend_cmap(
'theta_lh',
[(0.5,0,0), (1,1,0), (0,0.5,0), (0,1,1), (0,0,0.5),
(0.5,0,0.75), (1,0,1), (0.75,0,0.5), (0.5,0,0)])
cmap_theta_lh.__doc__ = '''
cmap_theta_lh is a colormap for plotting the pRF theta of a vertex.
Values passed to cmap_theta_lh should be scaled such that (-pi,pi rad) -> (0,1). Note that 0 is
the right right horizontal meridian and positive is the counter-clockwise direction.
cmap_theta_lh is a circular colormap that emphasizes colors in the right visual field; the
left visual field appears mostly magenta.
'''
cmap_theta_rh = blend_cmap(
'theta_rh',
[(0.5,0,0), (0.75,0,0.5), (0.5,0,0.75), (1,0,1),
(0,0,0.5), (0,1,1), (0,0.5,0), (1,1,0), (0.5,0,0)])
cmap_theta_rh.__doc__ = '''
cmap_theta_rh is a colormap for plotting the pRF theta of a vertex.
Values passed to cmap_theta_rh should be scaled such that (-pi,pi rad) -> (0,1). Note that 0 is
the right right horizontal meridian and positive is the counter-clockwise direction.
cmap_theta_rh is a circular colormap that emphasizes colors in the left visual field; the
right visual field appears mostly magenta.
'''
cmap_theta = blend_cmap(
'theta',
[(0, cmap_polar_angle(0.25)),
( 1.0/12, ( 1, 0, 1)),
( 3.0/12, (0.5, 0, 0)),
( 5.0/12, ( 1, 1, 0)),
( 7.0/12, ( 0,0.5, 0)),
( 9.0/12, ( 0, 1, 1)),
(11.0/12, ( 0, 0,0.5)),
(1, cmap_polar_angle(0.25))])
cmap_theta.__doc__ = '''
cmap_theta is a colormap for plotting the pRF theta of a vertex.
Values passed to cmap_theta should be scaled such that (-pi,pi rad) -> (0,1). Note that 0 is the
right right horizontal meridian and positive is the counter-clockwise direction.
cmap_theta is a 6-pronged circular colormap; note that it does not have dark or bright
values at the horizontal meridia as cmap_theta_sym, cmap_theta_lh, and
cmap_theta_rh do.
'''
cmap_eccentricity = blend_cmap(
'eccentricity',
[(0, ( 0, 0, 0)),
(1.25/90, ( 0, 0,0.5)),
(2.5/90, ( 1, 0, 1)),
(5.0/90, (0.5, 0, 0)),
(10.0/90, ( 1, 1, 0)),
(20.0/90, ( 0,0.5, 0)),
(40.0/90, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_eccentricity.__doc__ = '''
cmap_eccentricity is a colormap for plotting the pRF eccentricity of a vertex.
Values passed to cmap_eccentricity should be scaled such that (0,90 deg) -> (0,1).
'''
cmap_log_eccentricity = blend_cmap(
'log_eccentricity',
[(0, ( 0, 0, 0)),
(1.0/7, ( 0, 0,0.5)),
(2.0/7, ( 1, 0, 1)),
(3.0/7, (0.5, 0, 0)),
(4.0/7, ( 1, 1, 0)),
(5.0/7, ( 0,0.5, 0)),
(6.0/7, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_log_eccentricity.__doc__ = '''
cmap_log_eccentricity is a colormap for plotting the log of eccentricity.
Values passed to cmap_log_cmag should be scaled however desired, but note that the colormap
itself runs linearly from 0 to 1, so eccentricity data should be log-transformed
then scaled before being passed.
'''
cmap_radius = blend_cmap(
'radius',
[(0, ( 0, 0, 0)),
(1.25/60, ( 0, 0,0.5)),
(1.25/30, ( 1, 0, 1)),
( 2.5/30, (0.5, 0, 0)),
( 5.0/30, ( 1, 1, 0)),
(10.0/30, ( 0,0.5, 0)),
(20.0/30, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_radius.__doc__ = '''
cmap_radius is a colormap for plotting the pRF radius (sigma) of a vertex.
Values passed to cmap_radius should be scaled such that (0,30 deg) -> (0,1).
'''
cmap_log_radius = blend_cmap(
'log_radius',
[(0, ( 0, 0, 0)),
(1.0/7, ( 0, 0,0.5)),
(2.0/7, ( 1, 0, 1)),
(3.0/7, (0.5, 0, 0)),
(4.0/7, ( 1, 1, 0)),
(5.0/7, ( 0,0.5, 0)),
(6.0/7, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_log_radius.__doc__ = '''
cmap_log_radius is a colormap for plotting the log of radius.
Values passed to cmap_log_cmag should be scaled however desired, but note that the colormap
itself runs linearly from 0 to 1, so radius data should be log-transformed then scaled before
being passed.
'''
cmap_cmag = blend_cmap(
'cmag',
[(0, ( 0, 0, 0)),
(0.25/256, ( 0, 0, 0)),
(1.23/256, ( 0, 0,0.5)),
(2.97/256, ( 1, 0, 1)),
(7.25/256, (0.5, 0, 0)),
(17.69/256, ( 1, 1, 0)),
(43.13/256, ( 0,0.5, 0)),
(105.0/256, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_cmag.__doc__ = '''
cmap_cmag is a colormap for plotting cortical magnification.
It is generally advised to use the cmap_log_cmag colormap, as it will give better results for
the logarithmically-distributed corotical magnification colormap.
'''
cmap_log_cmag = blend_cmap(
'log_cmag',
[(0, ( 0, 0, 0)),
(1.0/7, ( 0, 0,0.5)),
(2.0/7, ( 1, 0, 1)),
(3.0/7, (0.5, 0, 0)),
(4.0/7, ( 1, 1, 0)),
(5.0/7, ( 0,0.5, 0)),
(6.0/7, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_log_cmag.__doc__ = '''
cmap_log_cmag is a colormap for plotting the log of cortical magnification.
Values passed to cmap_log_cmag should be scaled however desired, but note that the colormap
itself runs linearly from 0 to 1, so cortical magnification data should be log-transformed
then scaled before being passed.
'''
cmap_eccenflat = blend_cmap(
'eccenflat',
[(0, ( 0, 0, 0)),
(1.0/7, ( 0, 0,0.5)),
(2.0/7, ( 1, 0, 1)),
(3.0/7, (0.5, 0, 0)),
(4.0/7, ( 1, 1, 0)),
(5.0/7, ( 0,0.5, 0)),
(6.0/7, ( 0, 1, 1)),
(1, ( 1, 1, 1))])
cmap_eccenflat.__doc__ = '''
cmap_eccenflat is a colormap for plotting the log of eccentricity; in fact, it is identical
to the cmap_log_eccentricity colorscale, but when used with neuropythy it is assumed to be
linear instead of logarithmic.
'''
# A few other handy colormaps:
cmap_temperature_dark = blend_cmap('temperature_dark',
[(0,1,1), (0,0,1), (0,0,0), (1,0,0), (1,1,0)])
cmap_temperature = blend_cmap('temperature',
[(0,0,1), (0,1,1), (1,1,1), (1,1,0), (1,0,0)])
cmap_lightsaber_dark = blend_cmap('lightsaber_dark',
[(0,1,1), (0,1,0), (0,0,0), (1,0,0), (1,1,0)])
cmap_lightsaber = blend_cmap('lightsaber',
[(0,1,0), (0,1,1), (1,1,1), (1,1,0), (1,0,0)])
cmap_electricity_dark = blend_cmap('electricity_dark',
[(0,1,1), (0,0,1), (0,0,0), (0,1,0), (1,1,0)])
cmap_electricity = blend_cmap('electricity',
[(0,0,1), (0,1,1), (1,1,1), (1,1,0), (0,1,0)])
cmap_reddish = blend_cmap('reddish',
[(1,1,1), (1,1,0), (1,0,0), (0.5, 0, 0.25)])
cmap_reddish_dark = blend_cmap('reddish_dark',
[(0,0,0), (1,0,0), (1,1,0), (1, 1, 0.25)])
cmap_bluish = blend_cmap('bluish',
[(1,1,1), (0,1,1), (0,0,1), (0.25, 0, 0.5)])
cmap_bluish_dark = blend_cmap('bluish_dark',
[(0,0,0), (0,0,1), (0,1,1), (0.25, 1, 1)])
cmap_greenish = blend_cmap('greenish',
[(1,1,1), (0,1,1), (0,1,0), (0.25, 0.5, 0)])
cmap_greenish_dark = blend_cmap('greenish_dark',
[(0,0,0), (0,1,0), (1,1,0), (1, 1, 0.5)])
colormaps = {
'curvature': (cmap_curvature, (-1,1)),
'polar_angle': (cmap_polar_angle, (-180,180), 'deg'),
'polar_angle_sym': (cmap_polar_angle_sym, (-180,180), 'deg'),
'polar_angle_lh': (cmap_polar_angle_lh, (-180,180), 'deg'),
'polar_angle_rh': (cmap_polar_angle_rh, (-180,180), 'deg'),
'theta': (cmap_theta, (-np.pi,np.pi), 'rad'),
'theta_sym': (cmap_theta_sym, (-np.pi,np.pi), 'rad'),
'theta_lh': (cmap_theta_lh, (-np.pi,np.pi), 'rad'),
'theta_rh': (cmap_theta_rh, (-np.pi,np.pi), 'rad'),
'eccentricity': (cmap_eccentricity, (0,90), 'deg'),
'log_eccentricity': (cmap_log_eccentricity, (np.log(0.75), np.log(90.75)), 'deg'),
'radius': (cmap_radius, (0, 30), 'deg'),
'log_radius': (cmap_log_radius, (np.log(0.25), np.log(30.25)), 'deg'),
'cmag2': (cmap_cmag, (0.25, 512.25), 'mm**2/deg**2'),
'log_cmag2': (cmap_log_cmag, (np.log(0.25), np.log(512.25)), 'mm**2/deg**2'),
'cmag': (cmap_cmag, (0.5, 256.5), 'mm/deg'),
'log_cmag': (cmap_log_cmag, (np.log(0.5), np.log(256.5)), 'mm/deg'),
'eccenflat': (cmap_eccenflat, (0, 1)),
# the handy but non-neuroscience-based ones:
'temperature': (cmap_temperature, (-1,1)),
'temperature_dark': (cmap_temperature_dark, (-1,1)),
'lightsaber': (cmap_lightsaber, (-1,1)),
'lightsaber_dark': (cmap_lightsaber_dark, (-1,1)),
'electricity': (cmap_electricity, (-1,1)),
'electricity_dark': (cmap_electricity_dark, (-1,1)),
'reddish': (cmap_reddish, (0,1)),
'reddish_dark': (cmap_reddish_dark, (0,1)),
'greenish': (cmap_greenish, (0,1)),
'greenish_dark': (cmap_greenish_dark, (0,1)),
'bluish': (cmap_bluish, (0,1)),
'bluish_dark': (cmap_bluish_dark, (0,1))}
for (k,cmdat) in six.iteritems(colormaps): matplotlib.cm.register_cmap(k, cmdat[0])
def _diff_order(n):
u0 = np.arange(n)
d = int(np.ceil(np.sqrt(n)))
mtx = np.reshape(np.pad(u0, [(0,d*d-n)], 'constant', constant_values=[(0,-1)]), (d,d))
h = int((d+1)/2)
u = np.vstack([mtx[::2], mtx[1::2]]).T.flatten()
return u[u >= 0]
def label_cmap(colors, cmap=None, name=Ellipsis):
'''
label_cmap(n) yields a colormap with n color-steps that should be optimized such that each
colors is relatively different from its neighbors on the color spectrum. This is generally
well-suited for discrete catgory/label colormaps.
Note that this function uses a heuristic and is not guaranteed to be optimal in any way for any
value of n--but it generally works well enough for most common purposes.
The following optional arguments may be given:
* cmap (default: None) specifies a colormap to use as a base. If this is None, then a varianct
of 'hsv' is used.
* name (default: None) specifies a name (string) that will be used for the colormap name; if
Ellipsis is given, then ('label%d' % colors) is used.
'''
if not pimms.is_int(colors):
(lbls,ris) = np.unique(colors, return_inverse=True)
if lbls[0] == 0:
lbls = lbls[1:]
cm = label_cmap(len(lbls), cmap=cmap)
mx = np.max(lbls)
clrs = cm(np.linspace(0, 1, len(lbls)))
return blend_cmap(name, [(0, [0,0,0,0])] + list(zip(lbls/mx, clrs)))
if name is Ellipsis: name = 'label%d' % colors
if pimms.is_str(cmap): cmap = getattr(mpl.cm, cmap)
# get a diff-ordering
u = _diff_order(colors)
cm = matplotlib.cm.hsv if cmap is None else cmap
clrs = cm(u / float(colors))
if cmap is None and len(clrs) > 9:
# we use the hsv-like map: equivalent to using hsv then modifying it with a saturation
# and value label-cmap; this shouldn't go too low to prevent colors getting washed out
d = int(np.ceil(np.sqrt(colors)))
uu = _diff_order(d) / float(d-1)
uu = (uu + 1) / 2
ii = np.where(uu == 1)[0][0]
uu = np.roll(uu, -ii, 0)
uu = np.asarray([[x]*d for x in uu]).flatten()[:colors]
# make sure the highest value comes first
clrs[:,:3] *= np.reshape(uu, (-1,1))
cm = blend_cmap(name, list(zip(np.linspace(0,1,colors), clrs)))
return cm
except Exception: pass
def scale_for_cmap(cmap, x, vmin=Ellipsis, vmax=Ellipsis, unit=Ellipsis):
'''
scale_for_cmap(cmap, x) yields the values in x rescaled to be appropriate for the given
colormap cmap. The cmap must be the name of a colormap or a colormap object.
For a given cmap argument, if the object is a colormap itself, it is treated as cmap.name.
If the cmap names a colormap known to neuropythy, neuropythy will rescale the values in x
according to a heuristic.
'''
import matplotlib as mpl
if isinstance(cmap, mpl.colors.Colormap): cmap = cmap.name
(name, cm) = (None, None)
if cmap not in colormaps:
for (k,v) in six.iteritems(colormaps):
if cmap in k:
(name, cm) = (k, v)
break
else: (name, cm) = (cmap, colormaps[cmap])
if cm is not None:
cm = cm if len(cm) == 3 else (cm + (None,))
(cm, (mn,mx), uu) = cm
if vmin is Ellipsis: vmin = mn
if vmax is Ellipsis: vmax = mx
if unit is Ellipsis: unit = uu
if vmin is Ellipsis: vmin = None
if vmax is Ellipsis: vmax = None
if unit is Ellipsis: unit = None
x = pimms.mag(x) if unit is None else pimms.mag(x, unit)
if name is not None and name.startswith('log_'):
emn = np.exp(vmin)
x = np.log(x + emn)
vmin = np.nanmin(x) if vmin is None else vmin
vmax = np.nanmax(x) if vmax is None else vmax
return zdivide(x - vmin, vmax - vmin, null=np.nan)
def visual_field_legend(cmap, on=Ellipsis, max_eccentricity=12, transform=Ellipsis, pixels=288,
background=None, boundary_pixels=0):
'''
visual_field_map('polar_angle') yields an image array of a legend of the polar angle colormap,
plotted in the visual-field.
visual_field_map('eccentricity') yields an image array of a legend of the eccentricity colormap,
plotted in the visual-field.
visual_field_map(colormap, property) is used to color an arbitrary colormap on the given
property. The prorperty may be 'polar_angle', 'eccentricity', 'theta', 'rho', 'x', or 'y'. For
explanations of these properties, see below.
Note: the call visual_field_legend('log_eccentricity') will plot an eccentricity legend using
the log_eccentricity colormap; this is not a plot of the log-eccentricity, but rather is usually
just a smoother/nicer version of the 'eccentricity' plot. This automatic scaling is performed
only if the second parameter is not provided or is given as Ellipsis.
The following optional arguments are accepted:
* max_eccentricity (default: 12) specifies the maximum eccentricity of the plot; for polar
angle plots this is not important, but for eccentricity, this will limit the range of the
plot.
* transform (default: Ellipsis) specifies a transform function f(x) that is applied to the
given property before being passed to the colormap. Note that this function should accept
and return a vector of values between 0 and 1 (for use by the colormap). A value of None
indicates that no transformation should be used. The default value of Ellipsis indicates
that the transformation should be detected automatically based on the plotted property; see
below for default scaling.
* pixels (default: 288) specifies the number of pixels in the width/height of the image.
* background (default: None) specifies the background on which to plot the image. The default
value of None is equivalent to (1,1,1,0) or a transparent background. Anything that, when
passed to to_rgba() yields a single color-vector, can be passed to background.
* boundary_pixels (default: 0) specifies the number of pixels around the edge of the image to
reserve for the boundary. This may be specified as (left, right, bottom, top) or as
(sides, top-bottom). In all cases, the boundary_pixels value does not change the width of
the image but rather makes the plot smaller.
'''
# Start by parsing arguments:
if max_eccentricity is None: max_eccentricity = 90
if pimms.is_str(cmap):
cmap = cmap.replace('-', '_').replace(' ', '_').lower()
if on is Ellipsis:
if 'polar_angle' in cmap: on = 'polar_angle'
elif 'theta' in cmap: on = 'theta'
elif 'log_eccentricity' in cmap:
on = 'eccentricity'
(mn,mx) = np.log([0.75, 0.75 + max_eccentricity])
if transform is Ellipsis:
transform = lambda x: (np.log(x + 0.75) - mn) / (mx - mn)
elif 'eccentricity' in cmap: on = 'eccentricity'
if cmap in colormaps:
(cm, (mn,mx)) = colormaps[cmap][:2]
if transform is Ellipsis:
if cmap.startswith('log'):
emn = np.exp(mn)
transform = lambda x: (np.log(x+emn) - mn) / (mx - mn)
else:
transform = lambda x: (x - mn) / (mx - mn)
cmap = cm
else:
cmap = getattr(matplotlib.cm, cmap)
if transform is Ellipsis: transform = None
if background is None: background = (1,1,1,0)
else: background = to_rgba(background)
# go ahead and make the image and the image-portion we are using
whole_im = np.zeros((pixels, pixels, 4))
whole_im[:,:,:] = np.reshape(background, (1,1,4))
if boundary_pixels is None or boundary_pixels <= 0: im = whole_im
else: im = whole_im[boundary_pixels:-boundary_pixels, boundary_pixels:-boundary_pixels, :]
# note that this is visual field pixels:
mid = im.shape[0] / 2
(x_im, y_im) = np.meshgrid(np.arange(im.shape[0]) - mid, mid - np.arange(im.shape[1]))
r_im = np.sqrt(x_im**2 + y_im**2)
ii = np.where(r_im <= mid)
x = x_im[ii] / mid * max_eccentricity
y = y_im[ii] / mid * max_eccentricity
# what property are we plotting on?
if on is Ellipsis: raise ValueError('could not deduce property on which to operate')
elif not pimms.is_str(on): raise ValueError('plot-property must be a string')
on = on.replace('-', '_').replace(' ', '_').lower()
if on == 'polar_angle': p = np.mod(90 - 180/np.pi*np.arctan2(y, x) + 180, 360) - 180
elif on == 'theta': p = np.mod(np.arctan2(y, x) + np.pi, 2*np.pi) - np.pi
elif on == 'eccentricity': p = np.sqrt(x**2 + y**2)
elif on == 'rho': p = np.sqrt(x**2 + y**2) * np.pi/180
elif on == 'x': p = x
elif on == 'y': p = y
else: raise ValueError('unrecognized plot parameter: %s' % (on,))
# okay, transform p:
p = np.asarray(p if transform is None else transform(p))
# if there are nans, let's get rid of them
p[np.isnan(p)] = -np.inf
# and run it through the cmap...
clrs = cmap(p)
# and set the appropriate pixels...
im[ii] = clrs
# check for undersize requests
if boundary_pixels < 0:
bp = -boundary_pixels
whole_im = whole_im[bp:-bp, bp:-bp, :]
# and return the whole thing!
return whole_im
def to_rgba(val):
'''
to_rgba(val) is identical to matplotlib.colors.to_rgba(val) except that it operates over lists
as well as individual elements to yield matrices of rgba values. In addition, it always yields
numpy vectors or matrices.
'''
if pimms.is_npmatrix(val) and val.shape[1] == 4: return val
try: return np.asarray(matplotlib.colors.to_rgba(val))
except Exception: return np.asarray([matplotlib.colors.to_rgba(u) for u in val])
def color_overlap(color1, *args):
'''
color_overlap(color1, color2...) yields the rgba value associated with overlaying color2 on top
of color1 followed by any additional colors (overlaid left to right). This respects alpha
values when calculating the results.
Note that colors may be lists of colors, in which case a matrix of RGBA values is yielded.
'''
args = list(args)
args.insert(0, color1)
rgba = np.asarray([0.5,0.5,0.5,0])
for c in args:
c = to_rgba(c)
a = c[...,3]
a0 = rgba[...,3]
if np.isclose(a0, 0).all(): rgba = np.ones(rgba.shape) * c
elif np.isclose(a, 0).all(): continue
else: rgba = times(a, c) + times(1-a, rgba)
return rgba
_vertex_angle_empirical_prefixes = ['prf_', 'measured_', 'empiirical_']
_vertex_angle_model_prefixes = ['model_', 'predicted_', 'inferred_', 'template_', 'atlas_',
'benson14_', 'benson17_']
_vertex_angle_prefixes = ([''] + _vertex_angle_model_prefixes + _vertex_angle_model_prefixes)
def vertex_curvature_color(m):
return [0.2,0.2,0.2,1.0] if m['curvature'] > -0.025 else [0.7,0.7,0.7,1.0]
def vertex_prop(m, property=None):
if property is None or property is Ellipsis: return None
elif pimms.is_vector(property): return property
elif pimms.is_number(property): return property
elif not pimms.is_str(property): raise ValueError('argument is not property-like: %s'%property)
elif property in m: return m[property]
else: return None
def vertex_weight(m, property=None):
p = vertex_prop(m, property)
if p is None:
return next((m[k]
for name in ['weight', 'variance_explained']
for pref in ([''] + _vertex_angle_empirical_prefixes)
for k in [pref + name]
if k in m),
1.0)
return p
def vertex_angle(m, property=None):
p = vertex_prop(m, property)
if p is None:
ang0 = next((m[kk] for k in _vertex_angle_prefixes for kk in [k+'polar_angle'] if kk in m),
None)
if ang0 is not None: return ang0
ang0 = next((m[kk]
for name in ['angle', 'theta']
for k in _vertex_angle_prefixes for kk in [k + name]
if kk in m),
None)
if ang0 is not None:
return np.mod(90.0 - 180.0/np.pi*ang0 + 180, 360) - 180
return None
return p
def vertex_eccen(m, property=None):
p = vertex_prop(m, property)
if p is None:
ecc0 = next((m[k]
for kk in _vertex_angle_prefixes
for k in [kk + 'eccentricity']
if k in m),
None)
if ecc0 is not None: return ecc0
ecc0 = next((m[k]
for kk in _vertex_angle_prefixes
for k in [kk + 'rho']
if k in m),
None)
if ecc0 is not None:
return 180.0/np.pi*ecc0
return None
return p
def vertex_sigma(m, property=None):
p = vertex_prop(m, property)
if p is not None: return p
return next((m[k]
for kk in _vertex_angle_prefixes
for nm in ['sigma', 'radius', 'size']
for k in [kk + nm]
if k in m),
None)
def vertex_varea(m, property=None):
p = vertex_prop(m, property)
if p is not None: return p
return next((m[k]
for kk in _vertex_angle_prefixes
for nm in ['visual_area', 'varea', 'label']
for k in [kk + nm]
if k in m),
None)
def vertex_angle_color(m, weight_min=0.2, weighted=True, hemi=None, property=Ellipsis,
weight=Ellipsis, null_color='curvature'):
if m is Ellipsis:
kw0 = {'weight_min':weight_min, 'weighted':weighted, 'hemi':hemi,
'property':property, 'weight':weight, 'null_color':null_color}
def _pass_fn(x, **kwargs):
return vertex_angle_color(x, **{k:(kwargs[k] if k in kwargs else kw0[k])
for k in list(set(kwargs.keys() + kw0.keys()))})
return _pass_fn
if isinstance(null_color, basestring):
null_color = null_color.lower()
if null_color == 'curvature' or null_color == 'curv':
nullColor = np.asarray(vertex_curvature_color(m))
else:
raise ValueError('bad null color: %s' % null_color)
else:
nullColor = np.asarray(null_color)
ang = vertex_angle(m, property)
if ang is None: return nullColor
if isinstance(hemi, basestring):
hemi = hemi.lower()
if hemi == 'lh' or hemi == 'left':
ang = ang
elif hemi == 'rh' or hemi == 'right':
ang = -ang
elif hemi == 'abs':
ang = np.abs(ang)
else: raise ValueError('bad hemi argument: %s' % hemi)
w = vertex_weight(m, property=weight) if weighted else None
if weighted and (not pimms.is_number(w) or w < weight_min):
return nullColor
angColor = np.asarray(cmap_polar_angle_sym((ang + 180.0) / 360.0))
if weighted:
return angColor*w + nullColor*(1-w)
else:
return angColor
def vertex_eccen_color(m, weight_min=0.1, weighted=True, hemi=None,
property=Ellipsis, null_color='curvature', weight=Ellipsis):
if m is Ellipsis:
kw0 = {'weight_min':weight_min, 'weighted':weighted, 'hemi':hemi,
'property':property, 'weight':weight, 'null_color':null_color}
def _pass_fn(x, **kwargs):
return vertex_eccen_color(x, **{k:(kwargs[k] if k in kwargs else kw0[k])
for k in list(set(kwargs.keys() + kw0.keys()))})
return _pass_fn
if isinstance(null_color, basestring):
null_color = null_color.lower()
if null_color == 'curvature' or null_color == 'curv':
nullColor = np.asarray(vertex_curvature_color(m))
else:
raise ValueError('bad null color: %s' % null_color)
else:
nullColor = np.asarray(null_color)
ecc = vertex_eccen(m, property)
if ecc is None: return nullColor
w = vertex_weight(m, property=weight) if weighted else None
if weighted and (not pimms.is_number(w) or w < weight_min):
return nullColor
eccColor = np.asarray(cmap_eccentricity((ecc if 0 < ecc < 90 else 0 if ecc < 0 else 90)/90.0))
if weighted:
return eccColor*w + nullColor*(1-w)
else:
return eccColor
def vertex_sigma_color(m, weight_min=0.1, weighted=True, hemi=None,
property=Ellipsis, null_color='curvature', weight=Ellipsis):
if m is Ellipsis:
kw0 = {'weight_min':weight_min, 'weighted':weighted, 'hemi':hemi,
'property':property, 'weight':weight, 'null_color':null_color}
def _pass_fn(x, **kwargs):
return vertex_sigma_color(x, **{k:(kwargs[k] if k in kwargs else kw0[k])
for k in list(set(kwargs.keys() + kw0.keys()))})
return _pass_fn
if isinstance(null_color, basestring):
null_color = null_color.lower()
if null_color == 'curvature' or null_color == 'curv':
nullColor = np.asarray(vertex_curvature_color(m))
else:
raise ValueError('bad null color: %s' % null_color)
else:
nullColor = np.asarray(null_color)
sig = vertex_sigma(m, property)
if sig is None: return nullColor
w = vertex_weight(m, property=weight) if weighted else None
if weighted and (not pimms.is_number(w) or w < weight_min):
return nullColor
sigColor = np.asarray(cmap_radius(sig/30.0))
if weighted:
return sigColor*w + nullColor*(1-w)
else:
return sigColor
_varea_colors = {'V1': (1,0,0), 'V2': (0,1,0), 'V3': (0,0,1),
'hV4':(0,1,1), 'VO1':(0,0.5,1), 'VO2':(0,1,0.5),
'LO1':(1,0,1), 'LO2':(0.5,0,1), 'TO1':(1,0,0.5), 'TO2':(0.5,0,0.5),
'V3a':(1,1,0), 'V3a':(0.5,1,0)}
def vertex_varea_color(m, property=Ellipsis, null_color='curvature', hemi=None,
weight=Ellipsis, weight_min=0.1, weighted=False):
if m is Ellipsis:
kw0 = {'weight_min':weight_min, 'weighted':weighted, 'hemi':hemi,
'property':property, 'weight':weight, 'null_color':null_color}
def _pass_fn(x, **kwargs):
return vertex_varea_color(x, **{k:(kwargs[k] if k in kwargs else kw0[k])
for k in list(set(kwargs.keys() + kw0.keys()))})
return _pass_fn
if isinstance(null_color, basestring):
null_color = null_color.lower()
if null_color == 'curvature' or null_color == 'curv':
nullColor = np.asarray(vertex_curvature_color(m))
else:
raise ValueError('bad null color: %s' % null_color)
else:
nullColor = np.asarray(null_color)
lbl = vertex_varea(m, property)
if lbl is None or lbl == 0: return nullColor
w = vertex_weight(m, property=weight) if weighted else None
if weighted and (not pimms.is_number(w) or w < weight_min):
return nullColor
if not pimms.is_str(lbl):
lbl = None if lbl is None or lbl < 0 or lbl > len(visual_area_names) else \
visual_area_names[lbl]
lblColor = np.asarray(_varea_colors.get(lbl, None))
if lblColor is None: return nullColor
if weighted: return lblColor*w + nullColor*(1-w)
else: return lblColor
def curvature_colors(m):
'''
curvature_colors(m) yields an array of curvature colors for the vertices of the given
property-bearing object m.
'''
return np.asarray(m.map(vertex_curvature_color))
def retino_colors(vcolorfn, *args, **kwargs):
'See eccen_colors, angle_colors, sigma_colors, and varea_colors.'
if len(args) == 0:
def _retino_color_pass(*args, **new_kwargs):
return retino_colors(vcolorfn, *args,
**{k:(new_kwargs[k] if k in new_kwargs else kwargs[k])
for k in set(kwargs.keys() + new_kwargs.keys())})
return _retino_color_pass
elif len(args) > 1:
raise ValueError('retinotopy color functions accepts at most one argument')
m = args[0]
# we need to handle the arguments
if isinstance(m, (geo.VertexSet, pimms.ITable)):
tbl = m.properties if isinstance(m, geo.VertexSet) else m
n = tbl.row_count
# if the weight or property arguments are lists, we need to thread these along
if 'property' in kwargs:
props = kwargs['property']
del kwargs['property']
if not (pimms.is_vector(props) or pimms.is_matrix(props)):
props = [props for _ in range(n)]
else: props = None
if 'weight' in kwargs:
ws = kwargs['weight']
del kwargs['weight']
if not pimms.is_vector(ws) and not pimms.is_matrix(ws): ws = [ws for _ in range(n)]
else: ws = None
vcolorfn0 = vcolorfn(Ellipsis, **kwargs) if len(kwargs) > 0 else vcolorfn
if props is None and ws is None: vcfn = lambda m,k:vcolorfn0(m)
elif props is None: vcfn = lambda m,k:vcolorfn0(m, weight=ws[k])
elif ws is None: vcfn = lambda m,k:vcolorfn0(m, property=props[k])
else: vcfn = lambda m,k:vcolorfn0(m, property=props[k], weight=ws[k])
return np.asarray([vcfn(r,kk) for (kk,r) in enumerate(tbl.rows)])
else:
return vcolorfn(m, **kwargs)
def angle_colors(*args, **kwargs):
'''
angle_colors(obj) yields an array of colors for the polar angle map of the given
property-bearing object (cortex, tesselation, mesh).
angle_colors(dict) yields an array of the color for the particular vertex property
mapping that is given as dict.
angle_colors() yields a functor version of angle_colors that can be called with one of the
above arguments; note that this is useful precisely because the returned function
preserves the arguments passed; e.g. angle_colors(weighted=False)(mesh) is equivalent to
angle_colors(mesh, weighted=False).
The following options are accepted:
* weighted (True) specifies whether to use weight as opacity.
* weight_min (0.2) specifies that below this weight value, the curvature (or null color)
should be plotted.
* property (Ellipsis) specifies the specific property that should be used as the
polar angle value; if Ellipsis, will attempt to auto-detect this value.
* weight (Ellipsis) specifies the specific property that should be used as the weight value.
* null_color ('curvature') specifies a color that should be used as the background.
'''
return retino_colors(vertex_angle_color, *args, **kwargs)
def eccen_colors(*args, **kwargs):
'''
eccen_colors(obj) yields an array of colors for the eccentricity map of the given
property-bearing object (cortex, tesselation, mesh).
eccen_colors(dict) yields an array of the color for the particular vertex property mapping
that is given as dict.
eccen_colors() yields a functor version of eccen_colors that can be called with one of the
above arguments; note that this is useful precisely because the returned function
preserves the arguments passed; e.g. eccen_colors(weighted=False)(mesh) is equivalent to
eccen_colors(mesh, weighted=False).
The following options are accepted:
* weighted (True) specifies whether to use weight as opacity.
* weight_min (0.2) specifies that below this weight value, the curvature (or null color)
should be plotted.
* property (Ellipsis) specifies the specific property that should be used as the
eccentricity value; if Ellipsis, will attempt to auto-detect this value.
* weight (Ellipsis) specifies the specific property that should be used as the weight value.
* null_color ('curvature') specifies a color that should be used as the background.
'''
return retino_colors(vertex_eccen_color, *args, **kwargs)
def sigma_colors(*args, **kwargs):
'''
sigma_colors(obj) yields an array of colors for the pRF-radius map of the given
property-bearing object (cortex, tesselation, mesh).
sigma_colors(dict) yields an array of the color for the particular vertex property mapping
that is given as dict.
sigma_colors() yields a functor version of sigma_colors that can be called with one of the
above arguments; note that this is useful precisely because the returned function
preserves the arguments passed; e.g. sigma_colors(weighted=False)(mesh) is equivalent to
sigma_colors(mesh, weighted=False).
Note: radius_colors() is an alias for sigma_colors().
The following options are accepted:
* weighted (True) specifies whether to use weight as opacity.
* weight_min (0.2) specifies that below this weight value, the curvature (or null color)
should be plotted.
* property (Ellipsis) specifies the specific property that should be used as the
eccentricity value; if Ellipsis, will attempt to auto-detect this value.
* weight (Ellipsis) specifies the specific property that should be used as the weight value.
* null_color ('curvature') specifies a color that should be used as the background.
'''
return retino_colors(vertex_sigma_color, *args, **kwargs)
radius_colors = sigma_colors
def varea_colors(*args, **kwargs):
'''
varea_colors(obj) yields an array of colors for the visual area map of the given
property-bearing object (cortex, tesselation, mesh).
varea_colors(dict) yields an array of the color for the particular vertex property mapping
that is given as dict.
varea_colors() yields a functor version of varea_colors that can be called with one of the
above arguments; note that this is useful precisely because the returned function
preserves the arguments passed; e.g. varea_colors(weighted=False)(mesh) is equivalent to
varea_colors(mesh, weighted=False).
The following options are accepted:
* weighted (True) specifies whether to use weight as opacity.
* weight_min (0.2) specifies that below this weight value, the curvature (or null color)
should be plotted.
* property (Ellipsis) specifies the specific property that should be used as the
eccentricity value; if Ellipsis, will attempt to auto-detect this value.
* weight (Ellipsis) specifies the specific property that should be used as the weight value.
* null_color ('curvature') specifies a color that should be used as the background.
'''
return retino_colors(vertex_varea_color, *args, **kwargs)
def colors_to_cmap(colors):
colors = np.asarray(colors)
if len(colors.shape) == 1: return colors_to_cmap([colors])[0]
if colors.shape[1] == 3:
colors = np.hstack((colors, np.ones((len(colors),1))))
steps = (0.5 + np.asarray(range(len(colors)-1), dtype=np.float))/(len(colors) - 1)
return matplotlib.colors.LinearSegmentedColormap(
'auto_cmap',
{clrname: ([(0, col[0], col[0])] +
[(step, c0, c1) for (step,c0,c1) in zip(steps, col[:-1], col[1:])] +
[(1, col[-1], col[-1])])
for (clridx,clrname) in enumerate(['red', 'green', 'blue', 'alpha'])
for col in [colors[:,clridx]]},
N=(len(colors)))
def guess_cortex_cmap(pname):
'''
guess_cortex_cmap(proptery_name) yields a tuple (cmap, (vmin, vmax)) of a cortical color map
appropriate to the given property name and the suggested value scaling for the cmap. If the
given property is not a string or is not recognized then the log_eccentricity axis is used
and the suggested vmin and vmax are None.
'''
import matplotlib as mpl
if isinstance(pname, mpl.colors.Colormap): pname = pname.name
if not pimms.is_str(pname): return ('eccenflat', cmap_eccenflat, (None, None), None)
if pname in colormaps: (cm,cmname) = (colormaps[pname],pname)
else:
# check each manually
cm = None
for (k,v) in six.iteritems(colormaps):
if pname.endswith(k):
(cmname,cm) = (k,v)
break
if cm is None:
for (k,v) in six.iteritems(colormaps):
if pname.startswith(k):
(cmname,cm) = (k,v)
break
# we prefer log-eccentricity when possible
if cm is None: return ('eccenflat', cmap_eccenflat, (None, None), None)
if ('log_'+cmname) in colormaps:
cmname = 'log_'+cmname
cm = colormaps[cmname]
return (cmname,) + (cm if len(cm) == 3 else cm + (None,))
def apply_cmap(zs, cmap, vmin=None, vmax=None, unit=None, logrescale=False):
'''
apply_cmap(z, cmap) applies the given cmap to the values in z; if vmin and/or vmax are passed,
they are used to scale z.
Note that this function can automatically rescale data into log-space if the colormap is a
neuropythy log-space colormap such as log_eccentricity. To enable this behaviour use the
optional argument logrescale=True.
'''
zs = pimms.mag(zs) if unit is None else pimms.mag(zs, unit)
zs = np.asarray(zs, dtype='float')
if pimms.is_str(cmap): cmap = matplotlib.cm.get_cmap(cmap)
if logrescale:
if vmin is None: vmin = np.log(np.nanmin(zs))
if vmax is None: vmax = np.log(np.nanmax(zs))
mn = np.exp(vmin)
u = zdivide(nanlog(zs + mn) - vmin, vmax - vmin, null=np.nan)
else:
if vmin is None: vmin = np.nanmin(zs)
if vmax is None: vmax = np.nanmax(zs)
u = zdivide(zs - vmin, vmax - vmin, null=np.nan)
u[np.isnan(u)] = -np.inf
return cmap(u)
def cortex_cmap_plot_2D(the_map, zs, cmap, vmin=None, vmax=None, axes=None, triangulation=None):
'''
cortex_cmap_plot_2D(map, zs, cmap, axes) plots the given cortical map values zs on the given
axes using the given given color map and yields the resulting polygon collection object.
cortex_cmap_plot_2D(map, zs, cmap) uses matplotlib.pyplot.gca() for the axes.
The following options may be passed:
* triangulation (None) may specify the triangularion object for the mesh if it has already
been created; otherwise it is generated fresh.
* axes (None) specify the axes on which to plot; if None, then matplotlib.pyplot.gca() is
used. If Ellipsis, then a tuple (triangulation, z, cmap) is returned; to recreate the plot,
one would call:
axes.tripcolor(triangulation, z, cmap, shading='gouraud', vmin=vmin, vmax=vmax)
* vmin (default: None) specifies the minimum value for scaling the property when one is passed
as the color option. None means to use the min value of the property.
* vmax (default: None) specifies the maximum value for scaling the property when one is passed
as the color option. None means to use the max value of the property.
'''
if triangulation is None:
triangulation = matplotlib.tri.Triangulation(the_map.coordinates[0], the_map.coordinates[1],
triangles=the_map.tess.indexed_faces.T)
if axes is Ellipsis: return (triangulation, zs, cmap)
return axes.tripcolor(triangulation, zs, cmap=cmap, shading='gouraud', vmin=vmin, vmax=vmax)
def cortex_rgba_plot_2D(the_map, rgba, axes=None, triangulation=None):
'''
cortex_rgba_plot_2D(map, rgba, axes) plots the given cortical map on the given axes using the
given (n x 4) matrix of vertex colors and yields the resulting polygon collection object.
cortex_rgba_plot_2D(map, rgba) uses matplotlib.pyplot.gca() for the axes.
The option triangulation may also be passed if the triangularion object has already been
created; otherwise it is generated fresh.
'''
cmap = colors_to_cmap(rgba)
zs = np.linspace(0.0, 1.0, the_map.vertex_count)
return cortex_cmap_plot_2D(the_map, zs, cmap, axes=axes, triangulation=triangulation)
def cortex_plot_colors(the_map,
color=None, cmap=None, vmin=None, vmax=None, alpha=None,
underlay='curvature', mask=None):
'''
cortex_plot_colors(mesh, opts...) yields the cortex colors as a matrix of RGBA rows for the
given mesh and options.
The following options are accepted:
* color (default: None) specifies the color to plot for each vertex; this argument may take a
number of forms:
* None, do not plot a color over the underlay (the default)
* a matrix of RGB or RGBA values, one per vertex
* a property vector or a string naming a property, in which case the cmap, vmin, and vmax
arguments are used to generate colors
* a function that, when passed a single argument, a dict of the properties of a single
vertex, yields an RGB or RGBA list for that vertex.
* cmap (default: 'eccenflat') specifies the colormap to use in plotting if the color
argument provided is a property.
* vmin (default: None) specifies the minimum value for scaling the property when one is passed
as the color option. None means to use the min value of the property.
* vmax (default: None) specifies the maximum value for scaling the property when one is passed
as the color option. None means to use the max value of the property.
* underlay (default: 'curvature') specifies the default underlay color to plot for the
cortical surface; it may be None, 'curvature', or a color.
* alpha (default None) specifies the alpha values to use for the color plot. If None, then
leaves the alpha values from color unchanged. If a single number, then all alpha values in
color are multiplied by that value. If a list of values, one per vertex, then this vector
is multiplied by the alpha values. Finally, any negative value is set instead of multiplied.
So, for example, if there were 3 vertices with:
* color = ((0,0,0,1), (0,0,1,0.5), (0,0,0.75,0,8))
* alpha = (-0.5, 1, 0.5)
then the resulting colors plotted will be ((0,0,0,0.5), (0,0,1,0.5), (0,0,0.75,0,4)).
* mask (default: None) specifies a mask to use for the mesh; this is passed through to_mask()
to figure out the masking. Those vertices not in the mask are not plotted (but they will be
plotted in the underlay if it is not None).
'''
# okay, let's interpret the color
if color is None:
color = np.full((the_map.vertex_count, 4), 0.5)
color[:,3] = 0
elif pimms.is_map(color) and len(color) == 1:
(k,v) = next(six.iteritems(color))
try: ktag = np.random.randint(sys.maxsize)
except Exception: ktag = np.random.randint(sys.maxint)
ktag = '%016x_%s' % (ktag, k)
the_map = the_map.with_prop({ktag:v})
return cortex_plot_colors(the_map, color=ktag,
cmap=cmap, vmin=vmin, vmax=vmax, alpha=alpha,
underlay=underlay, mask=mask)
try:
clr = matplotlib.colors.to_rgba(color)
# This is an rgb color to plot...
color = np.ones((the_map.vertex_count,4)) * matplotlib.colors.to_rgba(clr)
except Exception: pass
if pimms.is_vector(color) or pimms.is_str(color):
# it's a property that gets interpreted via the colormap
p = the_map.property(color)
# if the colormap is none, we can try to guess it
logtr = False
if cmap is None:
(cmapname,cmap,(vmn,vmx),unit) = guess_cortex_cmap(color)
logtr = cmapname.startswith('log_')
else:
cmapname = cmap if pimms.is_str(cmap) else cmap.name
if cmapname in colormaps:
q = colormaps[cmapname]
(cmap,(vmn,vmx),unit) = q if len(q) == 3 else (q + (None,))
else:
(vmn,vmx,unit) = (np.nanmin(p), np.nanmax(p), None)
p = pimms.mag(p) if unit is None else pimms.mag(p, unit)
if vmin is None: vmin = vmn
if vmax is None: vmax = vmx
# we use logrescale here because we assume that if color was 'eccentricity' or similar,
# we want to rescale according to that colormap:
color = apply_cmap(p, cmap, vmin=vmin, vmax=vmax, unit=unit, logrescale=logtr)
if not pimms.is_matrix(color):
# must be a function; let's try it...
color = to_rgba(the_map.map(color))
color = np.array(color)
if color.shape[1] != 4: color = np.hstack([color, np.ones([color.shape[0], 1])])
# okay, and the underlay...
if underlay is not None:
if pimms.is_str(underlay) and underlay.lower() in ['curvature', 'curv']:
underlay = apply_cmap(the_map.prop('curvature'), cmap_curvature, vmin=-1, vmax=1)
else:
try: underlay = np.ones((the_map.vertex_count, 4)) * to_rgba(underlay)
except Exception: raise ValueError('plot underlay failed: must be a color or curvature')
# okay, let's check on alpha...
if alpha is not None:
if pimms.is_number(alpha): alpha = np.full(color.shape[0], alpha)
else: alpha = the_map.property(alpha)
color[:,3] *= alpha
neg = (alpha < 0)
color[neg,3] = -alpha[neg]
alpha = color[:,3]
# and the mask...
if mask is not None:
ii = the_map.mask(mask, indices=True)
tmp = np.zeros(len(color))
tmp[ii] = color[ii,3]
color[:,3] = tmp
# then, blend with the underlay if need be
if underlay is not None:
color = color_overlap(underlay, color)
return color
def cortex_plot_2D(the_map,
color=None, cmap=None, vmin=None, vmax=None, alpha=None,
underlay='curvature', mask=None, axes=None, triangulation=None):
'''
cortex_plot_2D(map) yields a plot of the given 2D cortical mesh, map.
The following options are accepted:
* color (default: None) specifies the color to plot for each vertex; this argument may take a
number of forms:
* None, do not plot a color over the underlay (the default)
* a matrix of RGB or RGBA values, one per vertex
* a property vector or a string naming a property, in which case the cmap, vmin, and vmax
arguments are used to generate colors
* a function that, when passed a single argument, a dict of the properties of a single
vertex, yields an RGB or RGBA list for that vertex.
* cmap (default: 'eccenflat') specifies the colormap to use in plotting if the color
argument provided is a property.
* vmin (default: None) specifies the minimum value for scaling the property when one is passed
as the color option. None means to use the min value of the property.
* vmax (default: None) specifies the maximum value for scaling the property when one is passed
as the color option. None means to use the max value of the property.
* underlay (default: 'curvature') specifies the default underlay color to plot for the
cortical surface; it may be None, 'curvature', or a color.
* alpha (default None) specifies the alpha values to use for the color plot. If None, then
leaves the alpha values from color unchanged. If a single number, then all alpha values in
color are multiplied by that value. If a list of values, one per vertex, then this vector
is multiplied by the alpha values. Finally, any negative value is set instead of multiplied.
So, for example, if there were 3 vertices with:
* color = ((0,0,0,1), (0,0,1,0.5), (0,0,0.75,0,8))
* alpha = (-0.5, 1, 0.5)
then the resulting colors plotted will be ((0,0,0,0.5), (0,0,1,0.5), (0,0,0.75,0,4)).
* mask (default: None) specifies a mask to use for the mesh; thi sis passed through map.mask()
to figure out the masking. Those vertices not in the mask are not plotted (but they will be
plotted in the underlay if it is not None).
* axes (default: None) specifies a particular set of matplotlib pyplot axes that should be
used. If axes is Ellipsis, then instead of attempting to render the plot, a tuple of
(tri, zs, cmap) is returned; in this case, tri is a matplotlib.tri.Triangulation
object for the given map and zs and cmap are an array and colormap (respectively) that
will produce the correct colors. Without axes equal to Ellipsis, these would instead
be rendered as axes.tripcolor(tri, zs, cmap, shading='gouraud'). If axes is None, then
uses the current axes.
* triangulation (default: None) specifies the matplotlib triangulation object to use, if one
already exists; otherwise a new one is made.
'''
# parse the axes
if axes is None: axes = matplotlib.pyplot.gca()
# process the colors
color = cortex_plot_colors(the_map, color=color, cmap=cmap, vmin=vmin, vmax=vmax, alpha=alpha,
underlay=underlay, mask=mask)
# finally, we can make the plot!
return cortex_rgba_plot_2D(the_map, color, axes=axes, triangulation=triangulation)
# 3D Graphics ######################################################################################
# If we're using Python 2, we're compatible with pysurfer:
def _ipyvolume_load_error(*args, **kwargs):
raise RuntimeError('load failure: the requested object could not be loaded, probably ' +
'because you do not have ipyvolume installed correctly')
cortex_plot_3D = _ipyvolume_load_error
try:
import ipyvolume as ipv
def cortex_plot_3D(obj,
color=None, cmap=None, vmin=None, vmax=None, alpha=None,
underlay='curvature', mask=None, hemi=None, surface='inflated',
figure=Ellipsis, width=600, height=600, mesh_alpha=None,
view=None, camera_distance=100, camera_fov=None, camera_up=None):
'''
cortex_plot_3D(hemi) plots the inflated surface of the given cortex object hemi and returns the
ipyvolume figure object.
cortex_plot_3D(mesh) plots the given mesh.
cortex_plot_3D((hemi1, hemi2...)) plots all the given hemispheres or meshes.
cortex_plot_3D(mesh) yields a PySurfer Brain object for the given 3D cortical mesh.
cortex_plot_3D(subject) is equivalent to cortex_plot_3D((subject.lh, subject.rh)).
The following options are accepted:
* color (default: None) specifies the color to plot for each vertex; this argument may take a
number of forms:
* None, do not plot a color over the underlay (the default)
* a matrix of RGB or RGBA values, one per vertex
* a property vector or a string naming a property, in which case the cmap, vmin, and vmax
arguments are used to generate colors
* a function that, when passed a single argument, a dict of the properties of a single
vertex, yields an RGB or RGBA list for that vertex.
* cmap (default: 'eccenflat') specifies the colormap to use in plotting if the color
argument provided is a property.
* vmin (default: None) specifies the minimum value for scaling the property when one is passed
as the color option. None means to use the min value of the property.
* vmax (default: None) specifies the maximum value for scaling the property when one is passed
as the color option. None means to use the max value of the property.
* underlay (default: 'curvature') specifies the default underlay color to plot for the
cortical surface; it may be None, 'curvature', or a color.
* alpha (default None) specifies the alpha values to use for the color plot. If None, then
leaves the alpha values from color unchanged. If a single number, then all alpha values in
color are multiplied by that value. If a list of values, one per vertex, then this vector
is multiplied by the alpha values. Finally, any negative value is set instead of multiplied.
So, for example, if there were 3 vertices with:
* color = ((0,0,0,1), (0,0,1,0.5), (0,0,0.75,0,8))
* alpha = (-0.5, 1, 0.5)
then the resulting colors plotted will be ((0,0,0,0.5), (0,0,1,0.5), (0,0,0.75,0,4)).
* mesh_alpha (default: None) specifies the alpha to use for the mesh object itself. This may
be single value (0 for transparent, 1 for opaque) or a vector of values, one per vertex, or
a property name.
* mask (default: None) specifies a mask to use for the mesh; this is passed through to_mask()
to figure out the masking. Those vertices not in the mask are not plotted (but they will be
plotted in the underlay if it is not None).
* hemi (defaut: None) specifies the hemisphere to use. If the passed mesh object is actually a
subject or mesh pair then this specifies which hemisphere to use. If the passed object is a
mesh, then this overrides its chirality, if specified in meta_data. If two hemispheres are
given, then this may be 'both' or 'split' in accordinace with PySurfer's Brain() class.
* surface (default: 'white') specifies the surface to use if the mesh object passed is in fact
either a cortex or subject object.
* axes (default: None) specifies the ipyvolume axes that should be used; if None, then the
current axes are used.
* stylize (default: True) specifies whether or not neuropythy should apply additional standard
stylings to the ipyvolume figure; these include making sure that the plot range is
appropriate, that the axes and box are not plotted, and that the camera is in an
orthographic position looking at the mesh.
'''
# First, see if we've been passed a cortex or subject object...
mesh = []
for arg in (obj if pimms.is_vector(obj) else [obj]):
if geo.is_mesh(arg): mesh.append(arg)
elif geo.is_topo(arg): mesh.append(geo.to_mesh((arg, surface)))
elif mri.is_subject(arg):
for h in (hemi if pimms.is_vector(hemi) else [hemi]):
if geo.is_topo(h) or geo.is_mesh(h): hh = [h]
elif h in arg.hemis: hh = [arg.hemis[h]]
else:
h = to_hemi_str(h)
hh = [arg.lh, arg.rh] if h == 'lr' else [arg.hemis[h]]
for hhh in hh: mesh.append(geo.to_mesh((hhh, surface)))
# process the colors
rgba = np.concatenate(
[cortex_plot_colors(m, color=color, cmap=cmap, vmin=vmin, vmax=vmax,
alpha=alpha, underlay=underlay, mask=mask)
for m in mesh])
n = len(rgba)
# process the mesh_alpha parameter
if mesh_alpha is None:
mesh_alpha = [None for m in mesh]
elif pimms.is_scalar(mesh_alpha):
if mesh_alpha == 1: mesh_alpha = None
elif pimms.is_str(mesh_alpha): mesh_alpha = [m.prop(mesh_alpha) for m in mesh]
else: mesh_alpha = [mesh_alpha for m in mesh]
elif len(mesh_alpha) != len(mesh):
mesh_alpha = [mesh_alpha for m in mesh]
# Okay, setup the ipyv figure...
mns = np.full(3, np.inf)
mxs = np.full(3, -np.inf)
ms = ()
if figure is None: f = ipv.gcf()
elif figure is Ellipsis: f = ipv.figure(width=width, height=height)
else: f = figure
i0 = 0
for (m,ma) in zip(mesh, mesh_alpha):
(x,y,z) = m.coordinates
ii = slice(i0, i0 + m.vertex_count)
i0 += m.vertex_count
ipvm = ipv.plot_trisurf(x,y,z, m.tess.indexed_faces.T, color=rgba[ii,:3])
mns = np.nanmin([mns, np.nanmin([x,y,z], axis=1)], axis=0)
mxs = np.nanmax([mxs, np.nanmax([x,y,z], axis=1)], axis=0)
ms = ms + (ipvm,)
# handle mesh alpha, if given
if ma is not None:
if pimms.is_scalar(ma): ma = np.full([m.vertex_count, 1], ma)
else: ma = np.reshape(ma, [m.vertex_count, 1])
ipvm.color = np.hstack([ipvm.color, ma])
ipvm.material.transparent = True
# Figure out the bounding box...
szs = mxs - mns
mid = 0.5*(mxs + mns)
bsz = np.nanmax(szs)
# okay, set the plot limits
mxs = mid + 0.5*bsz
mns = mid - 0.5*bsz
ipv.pylab.xlim(mns[0],mxs[0])
ipv.pylab.ylim(mns[1],mxs[1])
ipv.pylab.zlim(mns[2],mxs[2])
# few other styling things
ipv.pylab.xlabel('')
ipv.pylab.ylabel('')
ipv.pylab.zlabel('')
ipv.style.box_off()
ipv.style.axes_off()
# figure out the view
d = camera_distance
up = None
if view is None: view = 'back'
if pimms.is_str(view):
if camera_distance is None: d = 10
view = view.lower()
(view, up) = (((0,-d,0), (0,0, 1)) if view in ['back','rear','posterior','p','-y'] else
((0, d,0), (0,0, 1)) if view in ['front','anterior','a','+y','y'] else
((0,0, d), (0, 1,0)) if view in ['top','superior','s','+z','z'] else
((0,0,-d), (0,-1,0)) if view in ['bottom','inferior','i','-z'] else
(( d,0,0), (0,0, 1)) if view in ['right','r','x','+x'] else
((-d,0,0), (0,0, 1)) if view in ['left','l','-x'] else
(view, None))
if pimms.is_str(view): raise ValueError('Unknown view: %s' % view)
if camera_up is not None: up = camera_up
if d is None: d = np.sqrt(np.sum(np.asarray(d)**2))
fov = camera_fov
if fov is None: fov = 180.0/np.pi * 2 * np.arctan(1.125/(2*d))
f.camera.position = tuple(view)
f.camera.up = tuple(up)
f.camera.fov = fov
f.camera.lookAt(tuple(mid))
warnings.warn('neuropythy: NOTE: due to a bug in ipyvolume, camera views cannot currently' +
' be set by neuropythy; however, if you click the reset (home) button in' +
' the upper-left corner of the figure, the requested view will be fixed.')
return f
except Exception: pass
def cortex_plot(mesh, *args, **opts):
'''
cortex_plot(mesh) calls either cortex_plot_2D or cortex_plot_3D depending on the dimensionality
of the given mesh, and yields the resulting graphics object. All optional arguments supported
by each is supported by cortex plot.
The following options are accepted:
* color (default: None) specifies the color to plot for each vertex; this argument may take a
number of forms:
* None, do not plot a color over the underlay (the default)
* a matrix of RGB or RGBA values, one per vertex
* a property vector or a string naming a property, in which case the cmap, vmin, and vmax
arguments are used to generate colors
* a function that, when passed a single argument, a dict of the properties of a single
vertex, yields an RGB or RGBA list for that vertex.
* cmap (default: 'eccenflat') specifies the colormap to use in plotting if the color
argument provided is a property.
* vmin (default: None) specifies the minimum value for scaling the property when one is passed
as the color option. None means to use the min value of the property.
* vmax (default: None) specifies the maximum value for scaling the property when one is passed
as the color option. None means to use the max value of the property.
* underlay (default: 'curvature') specifies the default underlay color to plot for the
cortical surface; it may be None, 'curvature', or a color.
* alpha (default None) specifies the alpha values to use for the color plot. If None, then
leaves the alpha values from color unchanged. If a single number, then all alpha values in
color are multiplied by that value. If a list of values, one per vertex, then this vector
is multiplied by the alpha values. Finally, any negative value is set instead of multiplied.
So, for example, if there were 3 vertices with:
* color = ((0,0,0,1), (0,0,1,0.5), (0,0,0.75,0,8))
* alpha = (-0.5, 1, 0.5)
then the resulting colors plotted will be ((0,0,0,0.5), (0,0,1,0.5), (0,0,0.75,0,4)).
* mask (default: None) specifies a mask to use for the mesh; thi sis passed through map.mask()
to figure out the masking. Those vertices not in the mask are not plotted (but they will be
plotted in the underlay if it is not None).
* hemi (defaut: None) specifies the hemisphere to use. If the passed mesh object is actually a
subject or mesh pair then this specifies which hemisphere to use. If the passed object is a
mesh, then this overrides its chirality, if specified in meta_data. If two hemispheres are
given, then this may be 'both' or 'split' in accordinace with PySurfer's Brain() class.
* surface (default: 'white') specifies the surface to use if the mesh object passed is in fact
either a cortex or subject object.
* axes (default: None) specifies a particular set of matplotlib pyplot axes that should be
used. If axes is Ellipsis, then instead of attempting to render the plot, a tuple of
(tri, zs, cmap) is returned; in this case, tri is a matplotlib.tri.Triangulation
object for the given map and zs and cmap are an array and colormap (respectively) that
will produce the correct colors. Without axes equal to Ellipsis, these would instead
be rendered as axes.tripcolor(tri, zs, cmap, shading='gouraud'). If axes is None, then
uses the current axes.
* triangulation (default: None) specifies the matplotlib triangulation object to use, if one
already exists; otherwise a new one is made.
'''
if not isinstance(mesh, geo.Mesh) or mesh.coordinates.shape[0] > 2:
# must be a 3D call
return cortex_plot_3D(mesh, *args, **opts)
else:
return cortex_plot_2D(mesh, *args, **opts)
class ROIDrawer:
'''
ROIDrawer(axes, mproj) creates a new ROIDrawer class that interacts with the plot on
the given axes to save the lines drawn on the plot as a neuropythy path trace, which
is available from the roiDrawer object as roiDrawer.trace once the user is finished.
ROIDrawer(axes, fmap) extracts the map projection from the meta-data in the given
flatmap; note that if the projection is not encoded in the flatmap's meta-data, an
error will be raised.
All arguments and keyword arguments following the first two are passed verbatim to the
axes.plot() function, which is used for plotting the drawn lines.
'''
def __init__(self, axes, mp, closed=True, event_handlers=None, plot_list=None, *args, **kw):
from neuropythy import (to_map_projection, is_map_projection, is_vset)
# Assumption: the axes have already been plotted, so there's no need to do any
# plotting here, aside from the lines we're about to draw
self.axes = axes
# get the map projection
if is_vset(mp) and 'projection' in mp.meta_data: mp = mp.meta_data['projection']
elif not is_map_projection(mp): mp = to_map_projection(mp)
self.map_projection = mp
# draw the initial lines
if len(args) == 0 and len(kw) == 0: args = ['k.-']
self.line = axes.plot([], [], *args, **kw)[0]
self.xs = list(self.line.get_xdata())
self.ys = list(self.line.get_ydata())
self.connections = [
self.line.figure.canvas.mpl_connect('button_press_event', self.on_button),
self.line.figure.canvas.mpl_connect('key_press_event', self.on_key),
self.line.figure.canvas.mpl_connect('close_event', self.on_close)]
# get rid of the mesh for path traces (because we don't want to export meshes on accident)
if mp.mesh is not None: mp = mp.copy(mesh=None)
self.trace = geo.PathTrace(mp, np.zeros((2,0)), closed=closed,
meta_data={'roi_drawer':self})
self.closed = closed
self.event_handlers = event_handlers
self.plot_list = plot_list
if plot_list is not None and len(plot_list) > 0:
for p in plot_list[1:]: p.set_visible(False)
plot_list[0].set_visible(True)
self.current_plot = 0
def end(self, success=True):
from neuropythy import path_trace
# we've finished; clean up and make the line
if success:
if len(self.xs) < 1: raise ValueError('Drawn line has no points')
pts = np.transpose([self.xs, self.ys])
# remove us from the trace meta-data if we're in it
rd = self.trace.meta_data.get('roi_drawer')
if rd is self: self.trace.meta_data = self.trace.meta_data.discard('roi_drawer')
self.trace.persist()
else: self.trace = None
if self.line:
for conn in self.connections:
self.line.figure.canvas.mpl_disconnect(conn)
# redraw the final version:
if self.closed:
self.xs.append(self.xs[0])
self.ys.append(self.ys[0])
self.line.set_data(self.xs, self.ys)
self.line.figure.canvas.draw()
matplotlib.pyplot.close(self.line.figure)
# clear everything
self.connection = None
self.line = None
self.xs = None
self.ys = None
def on_close(self, event):
self.end(success=False)
def on_key(self, event):
import matplotlib.collections
if event.inaxes != self.line.axes: return
if event.key == 'tab':
# we go to the next plot...
if not self.plot_list: return
self.plot_list[self.current_plot].set_visible(False)
tmp = self.current_plot
self.current_plot = (self.current_plot + 1) % len(self.plot_list)
a = self.plot_list[self.current_plot]
if isinstance(a, matplotlib.collections.TriMesh): a.set_visible(True)
else: self.plot_list[self.current_plot] = a(self.axes)
a.figure.canvas.draw()
def on_button(self, event):
endkeys = ['shift+control', 'shift+ctrl', 'control+shift', 'ctrl+shift']
if self.line is None: return
if event.inaxes != self.line.axes: return
# if shift is down, we delete the last point
if event.key == 'shift':
if len(self.xs) == 0: return
self.xs = self.xs[:-1]
self.ys = self.ys[:-1]
elif event.key in endkeys:
# we abort!
self.end(success=False)
return
elif (event.key is not None and self.event_handlers is not None and
event.key in self.event_handlers):
self.event_handlers[event.key](self, event)
return
else: # add the points
self.xs.append(event.xdata)
self.ys.append(event.ydata)
# redraw the line regardless
self.line.set_data(self.xs, self.ys)
self.line.figure.canvas.draw()
# and update the trace
self.trace.points = np.asarray([self.xs, self.ys])
if event.dblclick or event.key in ['control', 'ctrl']:
self.end(success=True)
def trace_roi(hemi, map_proj, axes, closed=True, event_handlers=None, plot_list=None, **kw):
'''
trace_roi(hemi, map_proj, axes) creates an ROIDrawer object that controls the tracing of lines
around an ROI in a 2D matplotlib plot and returns a not-yet-persistent immutable PathTrace
object with the ROIDrawer in its meta_data. The path trace is persisted as soon as the user
finished drawing their line; if the line is canceled, then the trace is never persisted.
ROI tracing is very simple: any point in the plot is appended to the path as it is clicked; in
order to eliminate the previous point, hold shift while clicking. To end the path, hold control
while clicking. To abort the path, hold both shift and control while clicking. (Double-clicking
should be equivalent to control-clicking, but this does not work in all setups.) In order to use
the ROI tracing, `%matplotlib notebook` is recommended.
The trace_roi() function accepts all options that can be passed to cortex_plot() as well as the
following options:
* closed (default: True) specifies whether the path-trace that is constructed should be closed
(True) or open (False).
* event_handlers (default: None) specifies additional event handlers (named by key) for the
ROIDrawer().
* plot_list (default: None) specifies a list of alternate TriMesh objects that can be plotted
cyclically when the user presses tab. TriMesh objects can be created by pyplot.triplot and
pyplot.tripcolor, which are used by the neuropythy cortex_plot function as well. If the
plot_list is not empty, then the first item of the list is immediately plotted on the axes.
Unlike in the ROIDrawer function itself, the plot_list may contain maps whose keys are
the various arguments (aside from the initial mesh argument) to cortex_plot.
'''
# okay, first off, if the plot_list has maps in it, we convert them using cortex_plot:
if plot_list is not None:
if geo.is_flatmap(hemi): fmap = hemi
elif geo.is_flatmap(map_proj): fmap = map_proj
elif not geo.is_map_projection(map_proj): fmap = geo.to_map_projection(map_proj)(hemi)
else: fmap = map_proj(hemi)
plot_list = [cortex_plot(fmap, axes=axes, **p) if pimms.is_map(p) else p
for p in plot_list]
# next, make the roi drawer
rd = ROIDrawer(axes, map_proj, closed=closed,
event_handlers=event_handlers, plot_list=plot_list)
return rd.trace
