diff --git a/.travis.yml b/.travis.yml index 0c88a5a..69116ba 100644 --- a/.travis.yml +++ b/.travis.yml @@ -6,9 +6,9 @@ matrix: env: - ENV_FILE=environment_py27.yml - ENV_NAME=img_pipe_py2 - - python: "3.5" + - python: "3.6" env: - - ENV_FILE=environment_py35.yml + - ENV_FILE=environment_py36.yml - ENV_NAME=img_pipe_py3 before_script: - "export DISPLAY=:99.0" diff --git a/environment_py27.yml b/environment_py27.yml index 1bcb4c1..83ee13b 100644 --- a/environment_py27.yml +++ b/environment_py27.yml @@ -59,5 +59,5 @@ dependencies: - traits==4.6.0 - traitsui==5.1.0 - wcwidth==0.1.7 + - nwbext_ecog==0.7.2 prefix: ~/anaconda/envs/img_pipe_py2 - diff --git a/environment_py35.yml b/environment_py35.yml index 5414302..ccc0ecd 100644 --- a/environment_py35.yml +++ b/environment_py35.yml @@ -63,4 +63,3 @@ dependencies: - sympy==1.1.1 - tqdm==4.15.0 prefix: ~/anaconda/envs/img_pipe_py3 - diff --git a/environment_py36.yml b/environment_py36.yml new file mode 100644 index 0000000..79f6449 --- /dev/null +++ b/environment_py36.yml @@ -0,0 +1,121 @@ +name: img_pipe_py36 +channels: + - defaults + - conda-forge + - statiskit +dependencies: + #- appnope=0.1.0 #=py36hf537a9a_0 + - apptools=4.4.0=py36_1 + - backcall=0.1.0=py36_0 + - backports=1.0=py36_1 + - backports.os=0.1.1=py36_0 + - blas=1.0=mkl + - bzip2=1.0.6 #=h1de35cc_5 + - ca-certificates=2019.1.23=0 + - certifi=2019.3.9=py36_0 + - configobj=5.0.6=py36_1 + - curl=7.64.1 #=ha441bb4_0 + - cython=0.29.7 #=py36h0a44026_0 + - dbus=1.13.6 #=h90a0687_0 + - decorator=4.4.0=py36_1 + - envisage=4.7.1=py_0 + - expat=2.2.6 #=h0a44026_0 + - freetype=2.9.1 #=hb4e5f40_0 + - future=0.17.1=py36_0 + - gettext=0.19.8.1 #=h15daf44_3 + - glib=2.56.2 #=hd9629dc_0 + - hdf4=4.2.13 #=h39711bb_2 + - hdf5=1.10.4 #=hfa1e0ec_0 + - icu=58.2 #=h4b95b61_1 + - importlib_metadata=0.9=py36_0 + - intel-openmp=2019.3=199 + - ipython=7.5.0 #=py36h39e3cac_0 + - ipython_genutils=0.2.0 #=py36h241746c_0 + - jedi=0.13.3=py36_0 + - jpeg=9b #=he5867d9_2 + - jsoncpp=1.8.4 #=h04f5b5a_0 + - krb5=1.16.1 #=hddcf347_7 + - libcurl=7.64.1 #=h051b688_0 + - libcxx #=4.0.1 #=hcfea43d_1 + - libcxxabi #=4.0.1 #=hcfea43d_1 + - libedit=3.1.20181209 #=hb402a30_0 + - libffi=3.2.1 #=h475c297_4 + - libgfortran #=3.0.1 #=h93005f0_2 + - libiconv=1.15 #=hdd342a3_7 + - libnetcdf=4.6.1 #=hd5207e6_2 + - libogg=1.3.2 #=h1de35cc_0 + - libpng=1.6.37 #=ha441bb4_0 + - libssh2=1.8.2 #=ha12b0ac_0 + - libtheora=1.1.1 #=hb4e5f40_1 + - libtiff=4.0.10 #=hcb84e12_2 + - libvorbis=1.3.6 #=h1de35cc_0 + - libxml2=2.9.9 #=hab757c2_0 + - lz4-c=1.8.1.2 #=h1de35cc_0 + - mkl=2019.3=199 + - mkl_fft=1.0.12 #=py36h5e564d8_0 + - mkl_random=1.0.2 #=py36h27c97d8_0 + - ncurses=6.1 #=h0a44026_1 + - numpy=1.16.3 #=py36hacdab7b_0 + - numpy-base=1.16.3 #=py36h6575580_0 + - openssl=1.1.1b #=h1de35cc_1 + - parso=0.4.0=py_0 + - path.py=12.0.1=py_0 + - pcre=8.43 #=h0a44026_0 + - pexpect=4.7.0=py36_0 + - pickleshare=0.7.5=py36_0 + - pip=19.1.1=py36_0 + - prompt_toolkit=2.0.9=py36_0 + - ptyprocess=0.6.0=py36_0 + - pyface=6.0.0=py36_0 + - pygments=2.4.0=py_0 + - pyqt=5.9.2 #=py36h655552a_2 + - python=3.6.8 #=haf84260_0 + - qt=5.9.7 #=h468cd18_1 + - readline=7.0 #=h1de35cc_5 + - scipy=1.2.1 #=py36h1410ff5_0 + - setuptools=41.0.1=py36_0 + - simplegeneric=0.8.1=py36_2 + - sip=4.19.8 #=py36h0a44026_0 + - six=1.12.0=py36_0 + - sqlite=3.28.0 #=ha441bb4_0 + - tbb=2019.4 #=h04f5b5a_0 + - tk=8.6.8 #=ha441bb4_0 + - traitlets=4.3.2 #=py36h65bd3ce_0 + - traits=5.1.1 #=py36h1de35cc_0 + - traitsui=6.0.0=py36_1 + - vtk=8.2.0 #=py36h9bafd54_200 + - wcwidth=0.1.7 #=py36h8c6ec74_0 + - wheel=0.33.4=py36_0 + - xz=5.2.4 #=h1de35cc_4 + - zipp=0.3.3=py36_1 + - zlib=1.2.11 #=h1de35cc_3 + - zstd=1.3.7 #=h5bba6e5_0 + - pip: + - chardet==3.0.4 + - configparser==3.7.4 + - cycler==0.10.0 + - h5py==2.9.0 + - hdmf==1.0.3 + - idna==2.8 + - img-pipe==2019.3.15.2 + - kiwisolver==1.1.0 + - matplotlib==3.0.3 + - mayavi==4.6.2 + - mne==0.17.2 + - mpmath==1.1.0 + - nibabel==2.4.0 + - nipy==0.4.2 + - nwbext-ecog==0.7.2 + - pandas==0.24.2 + - pymcubes==0.0.9 + - pynwb==1.0.2 + - pyparsing==2.4.0 + - python-dateutil==2.8.0 + - pytz==2019.1 + - pyvtk==0.5.18 + - requests==2.22.0 + - ruamel-yaml==0.15.96 + - sympy==1.4 + - tqdm==4.32.1 + - urllib3==1.25.2 +#prefix: /anaconda3/envs/img_pipe_py36 diff --git a/gallery/nwb_tutorial_0.png b/gallery/nwb_tutorial_0.png new file mode 100644 index 0000000..2dacdd0 Binary files /dev/null and b/gallery/nwb_tutorial_0.png differ diff --git a/gallery/nwb_tutorial_1.png b/gallery/nwb_tutorial_1.png new file mode 100644 index 0000000..6a57281 Binary files /dev/null and b/gallery/nwb_tutorial_1.png differ diff --git a/gallery/nwb_tutorial_2.png b/gallery/nwb_tutorial_2.png new file mode 100644 index 0000000..9a56ede Binary files /dev/null and b/gallery/nwb_tutorial_2.png differ diff --git a/gallery/nwb_tutorial_3.png b/gallery/nwb_tutorial_3.png new file mode 100644 index 0000000..ee28baa Binary files /dev/null and b/gallery/nwb_tutorial_3.png differ diff --git a/gallery/nwb_tutorial_4.png b/gallery/nwb_tutorial_4.png new file mode 100644 index 0000000..98d1ed8 Binary files /dev/null and b/gallery/nwb_tutorial_4.png differ diff --git a/img_pipe/plotting/__init__.py b/img_pipe/plotting/__init__.py index acde2dc..3d982a2 100644 --- a/img_pipe/plotting/__init__.py +++ b/img_pipe/plotting/__init__.py @@ -1,3 +1,5 @@ from .ctmr_brain_plot import * from .plot_brain import * -#from mlab_3D_to_2D import * \ No newline at end of file +#from mlab_3D_to_2D import * + +from .plot_from_nwb import * diff --git a/img_pipe/plotting/plot_from_nwb.py b/img_pipe/plotting/plot_from_nwb.py new file mode 100644 index 0000000..e2cfb4d --- /dev/null +++ b/img_pipe/plotting/plot_from_nwb.py @@ -0,0 +1,196 @@ +# plot_from_nwb.py +# author: Luiz Tauffer +# date: 20.04.2019 + +import matplotlib.pyplot as plt +import mayavi +import scipy.io +from .ctmr_brain_plot import ctmr_gauss_plot, el_add +from img_pipe.SupplementalFiles import FS_colorLUT as FS_colorLUT +import numpy as np +import os + +import pynwb +import nwbext_ecog + + +def plot_from_nwb(subj_file, + hem='stereo', + opacity=1.0, + roi = [], + electrodes=False, + elec_roi=[], + elec_space='original', + fpath_atlas=''): + ''' + This function will plot the mesh structures and electrode positions + contained in a NWB file. + + Parameters + ---------- + subj_file : str + subject file path + hem : str, optional + hemisphere: 'lh', 'rh' or 'stereo'. Default='stereo' + opacity : float, optional + opacity of the mesh, between 0.0 and 1.0. Default=1.0 + roi : list, optional + list of ROI (dictionaries) to be highlighted. Default=[] + electrodes : boolean, optional + if True plots electrodes. Default=False + elec_roi : list, optional + list of ROI of electrodes to plot. Default=[] + elec_space : str, optional + 'original' or 'warped'. Default='original' + fpath_atlas : string, necessary if elec_space=='warped' + path to atlas model file + + + Returns + ------- + mesh : mayavi brain mesh + mlab : mayavi mlab scene + None : if an error occurred + ''' + + # Check for file path + exists = os.path.isfile(subj_file) + if exists: + io = pynwb.NWBHDF5IO(subj_file,'r') + nwb = io.read() + else: + print('ERROR: File or path does not exist!') + print('Given path to file: ', subj_file) + print('The standard path is given by your local Freesurfer installation', + '$SUBJECTS_DIR variable. If this is a wrong path, put the correct', + 'path as an argument, e.g.:') + return None + + # Check for Brain Atlas + if elec_space=='warped': + # Check for file path + exists = os.path.isfile(fpath_atlas) + if exists: + io_atlas = pynwb.NWBHDF5IO(fpath_atlas,'r') + nwb_atlas = io_atlas.read() + else: + print('ERROR: Could not find Brain Atlas!') + print('Searched path: ', fpath_atlas) + print('The standard path is given by your local Freesurfer installation', + '$SUBJECTS_DIR variable. Please save the Atlas file in this folder') + return None + + # Subject ID + subject_id = nwb.subject.subject_id + + # Get mesh and plot the pial surfaces + if (hem=='lh') or (hem=='rh'): + if elec_space=='original': + pial_mesh = nwb.subject.cortical_surfaces.surfaces[hem+'_pial'] + elif elec_space=='warped': + pial_mesh = nwb_atlas.subject.cortical_surfaces.surfaces[hem+'_pial'] + else: + print('ERROR: Invalid elec_space name.') + return None + mesh, mlab = ctmr_gauss_plot(pial_mesh.faces[()], + pial_mesh.vertices[()], + opacity=opacity, + color=(0.8, 0.8, 0.8)) + + elif hem=='stereo': + if elec_space=='original': #subject's brain mesh + pial_mesh_l = nwb.subject.cortical_surfaces['lh_pial'] + pial_mesh_r = nwb.subject.cortical_surfaces['rh_pial'] + elif elec_space=='warped': #atlas' brain mesh + pial_mesh_l = nwb_atlas.subject.cortical_surfaces.surfaces['lh_pial'] + pial_mesh_r = nwb_atlas.subject.cortical_surfaces.surfaces['rh_pial'] + else: + print('ERROR: Invalid elec_space name.') + return None + + mesh, mlab = ctmr_gauss_plot(pial_mesh_l.faces[()], + pial_mesh_l.vertices[()], + opacity=opacity, + color=(0.8, 0.8, 0.8)) + + mesh, mlab = ctmr_gauss_plot(pial_mesh_r.faces[()], + pial_mesh_r.vertices[()], + opacity=opacity, + color=(0.8, 0.8, 0.8), + new_fig=False) + + # all anatomic regions + anatomic_regions = list(nwb.subject.cortical_surfaces.surfaces.keys()) + + # Paint chosen ROI in distinguished colors + for ind, rg in enumerate(roi): + full_name = rg['name']+'_pial' + if full_name in anatomic_regions: + rg_color = rg['color'] + rg_mesh = nwb.subject.cortical_surfaces.surfaces[full_name] + + mesh, mlab = ctmr_gauss_plot(rg_mesh.faces[()], + rg_mesh.vertices[()], + new_fig=False, + color=rg_color) + else: + print(full_name,' not present in the provided data file!') + print('Check list of available parcellated areas with:') + print(' list(nwb.subject.cortical_surfaces.surfaces.keys())') + + + + # Plot electrodes + if electrodes: + # Electrodes positions + if elec_space=='original': + x_pos = nwb.electrodes['x'][:] + y_pos = nwb.electrodes['y'][:] + z_pos = nwb.electrodes['z'][:] + elec_pos = np.zeros((x_pos.shape[0],3)) + elec_pos[:,0] = x_pos + elec_pos[:,1] = y_pos + elec_pos[:,2] = z_pos + elif elec_space=='warped': + x_pos = nwb.electrodes['x_warped'][:] + y_pos = nwb.electrodes['y_warped'][:] + z_pos = nwb.electrodes['z_warped'][:] + elec_pos = np.zeros((x_pos.shape[0],3)) + elec_pos[:,0] = x_pos + elec_pos[:,1] = y_pos + elec_pos[:,2] = z_pos + + region = nwb.electrodes.columns[4][:] + + # Import freesurfer color lookup table as a dictionary + cmap = FS_colorLUT.get_lut() + + # Make a list of electrode numbers + elec_numbers = np.arange(x_pos.shape[0])+1 + + # Find all the unique brain areas in this subject + region = nwb.electrodes.columns[4][:] + all_regions = list(set(region)) + + # If no ROIs were chosen for electrodes + if len(elec_roi)==0: + elec_roi = all_regions + + # Loop through electrodes and plot them in each brain area + for ind, pos in enumerate(elec_pos): + b = region[ind] + if b[0:3]!='ctx' and b[0:4] != 'Left' and b[0:5] != 'Right' and b[0:5] != 'Brain' and b != 'Unknown': + if x_pos[ind]<0: #left side + this_label = 'ctx-lh-%s'%(b) + else: #right side + this_label = 'ctx-rh-%s'%(b) + + if (this_label != '') and (b in elec_roi): + try: #check if 'this_label' exists within cmap + el_color = np.array(cmap[this_label])/255. + el_add(np.atleast_2d(pos), color=tuple(el_color)) + except: + el_add(np.atleast_2d(pos)) + + mlab.show() + return mesh, mlab diff --git a/setup.py b/setup.py index ce55100..c441e95 100644 --- a/setup.py +++ b/setup.py @@ -17,7 +17,7 @@ def run(self): packages = find_packages(), include_package_data = True, setup_requires=['cython','numpy','scipy'], - install_requires=['numpy','scipy','pyvtk','mayavi','pymcubes','mne','nibabel','nipy','matplotlib<2.0.0','configparser','tqdm'], + install_requires=['numpy','scipy','pyvtk','mayavi','pymcubes','mne','nibabel','nipy','matplotlib','configparser','tqdm'], dependency_links=['http://www.vtk.org/files/release/6.3/VTK-6.3.0.tar.gz','https://downloads.sourceforge.net/project/pyqt/PyQt4/PyQt-4.12/PyQt4_gpl_mac-4.12.tar.gz?r=&ts=1487961590&use_mirror=superb-sea2'], cmdclass={'install':MyInstall}, classifiers = [ diff --git a/tutorials/img_pipe_plot_from_nwb.ipynb b/tutorials/img_pipe_plot_from_nwb.ipynb new file mode 100644 index 0000000..bc3e41f --- /dev/null +++ b/tutorials/img_pipe_plot_from_nwb.ipynb @@ -0,0 +1,229 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plotting from `NWB` files\n", + "In this tutorial you'll learn to use the `plot_from_nwb()` function to plot the 3D mesh of brain structures and electrode positions contained in a [Neurodata Without Borders: Neurophysiology 2.0](https://www.nwb.org/) file format.\n", + "\n", + "Let's start importing the relevant modules:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import os\n", + "from img_pipe import plotting" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now indicate the location of your `nwb` file. The standard way is to place it inside the Freesurfers subjects path, defined by the environment variable `SUBJECTS_DIR`. You can also select which hemisphere to plot (the default is both, 'stereo').\n", + "\n", + "First let us take a look at the subjects's brain structure: " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tauffer/anaconda3/envs/img_pipe_py3t/lib/python3.5/site-packages/mayavi-4.5.0-py3.5-linux-x86_64.egg/tvtk/array_handler.py:268: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n", + " assert not numpy.issubdtype(z.dtype, complex), \\\n" + ] + } + ], + "source": [ + "subj_file = 'ben_subjects/EC125_B22.nwb'\n", + "\n", + "# Call function and generate mlab image\n", + "mesh, mlab = plotting.plot_from_nwb(subj_file)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You should see an interactive pop up window with an image similar to this:\n", + "![simple_brain](../gallery/nwb_tutorial_0.png)\n", + "\n", + "We can also emphasize any **Region of Interest** defined in the `NWB` file: " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tauffer/anaconda3/envs/img_pipe_py3t/lib/python3.5/site-packages/mayavi-4.5.0-py3.5-linux-x86_64.egg/tvtk/array_handler.py:268: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n", + " assert not numpy.issubdtype(z.dtype, complex), \\\n" + ] + } + ], + "source": [ + "lh_lOFC = {\"name\":'lh_lOFC', \"opacity\":1., \"color\":(.5,.3,.5)}\n", + "lh_LHinsula = {\"name\":'lh_LHinsula', \"opacity\":1., \"color\":(.1,.1,.7)}\n", + "roi = [lh_lOFC, lh_LHinsula]\n", + "\n", + "mesh, mlab = plotting.plot_from_nwb(subj_file, opacity=.3, roi=roi)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You should see an almost transparent brain with colored selected regions:\n", + "![rois](../gallery/nwb_tutorial_1.png)\n", + "\n", + "Next, let's plot the **electrodes** in their original coordinates." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tauffer/anaconda3/envs/img_pipe_py3t/lib/python3.5/site-packages/mayavi-4.5.0-py3.5-linux-x86_64.egg/tvtk/array_handler.py:268: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n", + " assert not numpy.issubdtype(z.dtype, complex), \\\n" + ] + } + ], + "source": [ + "mesh, mlab = plotting.plot_from_nwb(subj_file, hem='stereo', opacity=0.3, electrodes=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You should see an almost transparent brain with electrodes colored by region:\n", + "![elecs_original](../gallery/nwb_tutorial_2.png)\n", + "\n", + "We can select a specific group of **electrodes by location**, and plot only the relevant hemisphere:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tauffer/anaconda3/envs/img_pipe_py3t/lib/python3.5/site-packages/mayavi-4.5.0-py3.5-linux-x86_64.egg/tvtk/array_handler.py:268: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n", + " assert not numpy.issubdtype(z.dtype, complex), \\\n" + ] + } + ], + "source": [ + "elec_roi_list = ['rostralmiddlefrontal', 'postcentral']\n", + "\n", + "mesh, mlab = plotting.plot_from_nwb(subj_file, hem='rh', opacity=0.8, electrodes=True, elec_roi=elec_roi_list)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You should see just the right hemisphere with the electrodes of selected locations:\n", + "![elecs_loc](../gallery/nwb_tutorial_3.png)\n", + "\n", + "Finally, we can visualize the electrodes at their **warped coordinates**:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tauffer/anaconda3/envs/img_pipe_py3t/lib/python3.5/site-packages/mayavi-4.5.0-py3.5-linux-x86_64.egg/tvtk/array_handler.py:268: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n", + " assert not numpy.issubdtype(z.dtype, complex), \\\n" + ] + } + ], + "source": [ + "fs_dir = os.environ['SUBJECTS_DIR']\n", + "fpath_atlas = os.path.join(fs_dir, 'cvs_avg35_inMNI152.nwb')\n", + "\n", + "mesh, mlab = plotting.plot_from_nwb(subj_file, hem='rh', opacity=0.4, electrodes=True, \n", + " elec_space='warped', fpath_atlas=fpath_atlas)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You should see just the right hemisphere and the electrodes locations warped to the atlas:\n", + "![elecs_warped](../gallery/nwb_tutorial_4.png)\n", + "\n", + "\n", + "\n", + "**Useful references:**
\n", + "[Neurodata Without Borders: Neurophysiology 2.0](https://www.nwb.org/)\n", + "[NWB for Python](https://pynwb.readthedocs.io/en/stable/index.html)
\n", + "[nwbext_ecog: An NWB extension for ECoG data](https://github.com/bendichter/nwbext_ecog)
\n", + "[img_pipe: Image processing pipeline for ECoG data](https://github.com/ChangLabUcsf/img_pipe)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "hide_input": false, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + }, + "widgets": { + "state": {}, + "version": "1.1.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}