Skip to content

Instantly share code, notes, and snippets.

@kingjr

kingjr/sources_3d.ipynb

Last active January 17, 2018 14:31
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save kingjr/cdadefbb5f522ab7450ead746dc0a5b2 to your computer and use it in GitHub Desktop.
Save kingjr/cdadefbb5f522ab7450ead746dc0a5b2 to your computer and use it in GitHub Desktop.
Download ZIP
Example continuing on Chris' p3volume example to plot source estimates
Display the source blob
Display the rendered blob
Raw
sources_3d.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.4.2'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import ipyvolume\n",
"ipyvolume.__version__"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/jrking/anaconda2/lib/python2.7/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
" from ._conv import register_converters as _register_converters\n"
]
}
],
"source": [
"import numpy as np\n",
"import ipyvolume.pylab as p3\n",
"import nibabel\n",
"import mne\n",
"import surfer\n",
"import matplotlib.pyplot as plt\n",
"\n",
"data_path = mne.datasets.sample.data_path()\n",
"\n",
"# Read Anatomy\n",
"subjects_dir = data_path + '/subjects'\n",
"subject = 'sample'\n",
"\n",
"fname_pial = '/'.join((subjects_dir, subject, 'surf/lh.pial'))\n",
"rr, tris = nibabel.freesurfer.read_geometry(fname_pial)\n",
"tris = tris.astype(np.uint32)\n",
"\n",
"# Curvature\n",
"fname_curv = '/'.join((subjects_dir, subject, 'surf/lh.curv'))\n",
"curv = nibabel.freesurfer.read_morph_data(fname_curv)\n",
"\n",
"# Normalize curvature to make gray cortex\n",
"curv = (curv > 0).astype(float)\n",
"curv = (curv - 0.5) / 3 + 0.5\n",
"curv = curv[:, np.newaxis] * [1, 1, 1]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Load source estimate\n",
"stc_fname = data_path + '/MEG/sample/sample_audvis-meg-lh.stc'\n",
"stc = mne.read_source_estimate(stc_fname)\n",
"\n",
"# smooth to get a color on all vertices\n",
"adj_mat = surfer.utils.mesh_edges(tris)\n",
"smooth_mat = surfer.utils.smoothing_matrix(stc.vertices[0], adj_mat, 20, verbose=False)\n",
"data = stc.data[:len(stc.vertices[0])]\n",
"data = smooth_mat.dot(data)\n",
"\n",
"# transform into color\n",
"data -= data.min()\n",
"data /= data.max()\n",
"data[data>1.] = 1.\n",
"cmap = plt.get_cmap('hot')\n",
"colors = cmap(data.ravel())[:, :3]\n",
"colors = colors.reshape(np.r_[data.shape, 3]).transpose(1, 0, 2)\n",
"# Combine activations and cortex color\n",
"brain_colors = colors * data.T[:, :, None] * 2 + curv[None, :, :]\n",
"brain_colors[brain_colors>1] = 1"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "c5a034b2ea4941d8bc4fc8704b6bc7a8",
"version_major": 2,
"version_minor": 0
},
"text/html": [
"<p>Failed to display Jupyter Widget of type <code>VBox</code>.</p>\n",
"<p>\n",
" If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
" that the widgets JavaScript is still loading. If this message persists, it\n",
" likely means that the widgets JavaScript library is either not installed or\n",
" not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
" Widgets Documentation</a> for setup instructions.\n",
"</p>\n",
"<p>\n",
" If you're reading this message in another frontend (for example, a static\n",
" rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
" it may mean that your frontend doesn't currently support widgets.\n",
"</p>\n"
],
"text/plain": [
"VBox(children=(Figure(animation=10.0, animation_exponent=1.0, camera_center=[0.0, 0.0, 0.0], height=500, matrix_projection=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], matrix_world=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], meshes=[Mesh(color=array([[[0.33478543, 0.33333333, 0.33333333],\n",
" [0.33476499, 0.33333333, 0.33333333],\n",
" [0.33477181, 0.33333333, 0.33333333],\n",
" ...,\n",
" [0.3384923 , 0.33333333, 0.33333333],\n",
" [0.34028782, 0.33333333, 0.33333333],\n",
" [0.34108852, 0.33333333, 0.33333333]],\n",
"\n",
" [[0.34493321, 0.33333333, 0.33333333],\n",
" [0.34490664, 0.33333333, 0.33333333],\n",
" [0.34469758, 0.33333333, 0.33333333],\n",
" ...,\n",
" [0.33733186, 0.33333333, 0.33333333],\n",
" [0.33878156, 0.33333333, 0.33333333],\n",
" [0.34018128, 0.33333333, 0.33333333]],\n",
"\n",
" [[0.33863127, 0.33333333, 0.33333333],\n",
" [0.33854908, 0.33333333, 0.33333333],\n",
" [0.33846188, 0.33333333, 0.33333333],\n",
" ...,\n",
" [0.33977613, 0.33333333, 0.33333333],\n",
" [0.34202992, 0.33333333, 0.33333333],\n",
" [0.34204094, 0.33333333, 0.33333333]],\n",
"\n",
" ...,\n",
"\n",
" [[0.3550426 , 0.33333333, 0.33333333],\n",
" [0.35495527, 0.33333333, 0.33333333],\n",
" [0.35460489, 0.33333333, 0.33333333],\n",
" ...,\n",
" [0.35149677, 0.33333333, 0.33333333],\n",
" [0.34490638, 0.33333333, 0.33333333],\n",
" [0.34277955, 0.33333333, 0.33333333]],\n",
"\n",
" [[0.34754734, 0.33333333, 0.33333333],\n",
" [0.3475103 , 0.33333333, 0.33333333],\n",
" [0.34744286, 0.33333333, 0.33333333],\n",
" ...,\n",
" [0.34756657, 0.33333333, 0.33333333],\n",
" [0.34765854, 0.33333333, 0.33333333],\n",
" [0.3495473 , 0.33333333, 0.33333333]],\n",
"\n",
" [[0.33592404, 0.33333333, 0.33333333],\n",
" [0.33593499, 0.33333333, 0.33333333],\n",
" [0.33598781, 0.33333333, 0.33333333],\n",
" ...,\n",
" [0.36925096, 0.33333333, 0.33333333],\n",
" [0.36250966, 0.33333333, 0.33333333],\n",
" [0.36595438, 0.33333333, 0.33333333]]]), texture=None, triangles=array([[ 0, 1, 3],\n",
" [ 4, 3, 1],\n",
" [ 0, 61, 62],\n",
" ...,\n",
" [153723, 153732, 153919],\n",
" [153900, 153907, 153901],\n",
" [154073, 154322, 154245]], dtype=uint32), x=array([-21.12763786, -21.47982788, -21.88400841, ..., -13.28218079,\n",
" -15.01982021, -16.47038078]), y=array([-88.4329071 , -88.38152313, -88.28314209, ..., 75.84895325,\n",
" 75.05326843, 74.3368454 ]), z=array([13.29888344, 13.2219944 , 13.2074194 , ..., 57.0162735 ,\n",
" 58.20228577, 56.61045456]))], style={'box': {'visible': True}, 'axes': {'color': 'black', 'visible': True, 'ticklabel': {'color': 'black'}, 'label': {'color': 'black'}}, 'background-color': 'white'}, tf=None, xlim=[-115.89412903785706, 51.105207204818726], ylim=[-88.54894256591797, 78.45039367675781], zlim=[-37.87920904159546, 129.12012720108032]), HBox(children=(Play(value=0, interval=10, max=24), FloatSlider(value=0.0, max=24.0, step=1.0)))))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = p3.figure(width=500, height=500)\n",
"brain = p3.plot_trisurf(rr[:, 0], rr[:, 1], rr[:, 2],\n",
" triangles=tris, color=brain_colors)\n",
"p3.squarelim()\n",
"p3.animation_control(brain, sequence_length=len(stc.times), interval=10)\n",
"p3.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment