[MRG] Adds iEEG Coordinates Writing That Is BIDS Compatible (#335)

* Adding updated bids root. * Fixing readme to reflect correct change. * Fixing examples to show the bids_root. * Initial commit. * Adding mne-python to git ignore. * Adding updates to incorporate stephan's changes. * BIDS compatible ieeg coordinates writing. * Fixing gitignore. * Adding example using mne data. * Adding fixes that updated to the correct version of write.py. * Adding some plotting at the end of the example. * Removing mayavi.' * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Alexandre Gramfort <alexandre.gramfort@m4x.org> * Adding updated mne dir. * Fix dataset path. * Fix examples from rebase. * Fixing issues that happened in rebase. * Fixing issues that happened in rebase. * Fixing issues that happened in rebase. * Fixing issues that happened in rebase. * Adding unit test for running ecog montage loading. * Adding unit test for reading in the ieeg data. * Fix coverage. * Fixing issues from agram. * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Alexandre Gramfort <alexandre.gramfort@m4x.org> * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Alexandre Gramfort <alexandre.gramfort@m4x.org> * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Alexandre Gramfort <alexandre.gramfort@m4x.org> * Fix order of imports. * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Merging remote changes. * Merging remote changes. * Trying to fix example. * Trying to fix example. * Fixing flake. * Adding updated example file. * Adding updated description. * Changing whatsnew, contributing and read function. * Fixing suggestions by Jas. * Fixing bug in calling verbose. * Fixing flake8. * Adding an updated example. * Update CONTRIBUTING.md Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Update CONTRIBUTING.md Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Rebase. * Fixing one line change of import. * Update mne_bids/write.py Co-Authored-By: Stefan Appelhoff <stefan.appelhoff@mailbox.org> * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Stefan Appelhoff <stefan.appelhoff@mailbox.org> * Adding side car jsons. * Adding additional info to side cars. * Increase test coverage. * Increase test coverage. * Increase test coverage. * Update mne_bids/write.py Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Adding ieeg example fixes based on agram. * Fixing kind comment from jas. * Update examples/convert_ieeg_to_bids.py Co-Authored-By: Mainak Jas <jasmainak@users.noreply.github.com> * Merging. * Fixing the formatting in instructions. Co-authored-by: Alexandre Gramfort <alexandre.gramfort@m4x.org> Co-authored-by: Mainak Jas <jasmainak@users.noreply.github.com> Co-authored-by: Stefan Appelhoff <stefan.appelhoff@mailbox.org>
mne-tools · Feb 3, 2020 · d90b59c · d90b59c
1 parent d3d7ef2
commit d90b59c
Show file tree

Hide file tree

Showing 8 changed files with 533 additions and 60 deletions.
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
@@ -53,3 +53,14 @@ The documentation can be built using sphinx. For that, please additionally
 install the following:
 
     $ pip install matplotlib nilearn sphinx numpydoc sphinx-gallery sphinx_bootstrap_theme pillow
+
+To build the documentation locally, one can run:
+
+    $ cd doc/
+    $ make html
+
+or
+
+    $ make html-noplot
+
+if you don't want to run the examples to build the documentation. This will result in a faster build but produce no plots in the examples.
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
@@ -63,6 +63,7 @@ Changelog
 - :func:`write_anat` allows writing T1 weighted MRI scans for subjects and optionally creating a T1w.json sidecar from a supplied :code:`trans` object, by `Stefan Appelhoff`_ (`#211 <https://github.com/mne-tools/mne-bids/pull/211>`_)
 - :func:`read_raw_bids` will return the the raw object with :code:`raw.info['bads']` already populated, whenever a :code:`channels.tsv` file is present, by `Stefan Appelhoff`_ (`#209 <https://github.com/mne-tools/mne-bids/pull/209>`_)
 - :func:`read_raw_bids` is now more likely to find event and channel sidecar json files, by `Marijn van Vliet`_ (`#233 <https://github.com/mne-tools/mne-bids/pull/233>`_)
+- Enhanced :func:`read_raw_bids` and :func:`write_raw_bids` for iEEG coordinates along with example and unit test, by `Adam Li`_ (`#335 <https://github.com/mne-tools/mne-bids/pull/335/>`_)
 
 Bug
 ~~~

diff --git a/examples/convert_eeg_to_bids.py b/examples/convert_eeg_to_bids.py
@@ -47,11 +47,8 @@
 # Conveniently, there is already a data loading function available with
 # MNE-Python:
 
-# Make a directory to save the data to
-home = os.path.expanduser('~')
-mne_dir = os.path.join(home, 'mne_data')
-if not os.path.exists(mne_dir):
-    os.makedirs(mne_dir)
+# get MNE directory w/ example data
+mne_dir = mne.get_config('MNE_DATASETS_SAMPLE_PATH')
 
 # Define which tasks we want to download.
 tasks = [2,  # This is 2 minutes of eyes closed rest
@@ -120,7 +117,7 @@
 subject_id = '001'  # zero padding to account for >100 subjects in this dataset
 task = 'resteyesclosed'
 raw_file = raw
-bids_root = os.path.join(home, 'mne_data', 'eegmmidb_bids')
+bids_root = os.path.join(mne_dir, 'eegmmidb_bids')
 
 ###############################################################################
 # Now we just need to specify a few more EEG details to get something sensible:

diff --git a/examples/convert_ieeg_to_bids.py b/examples/convert_ieeg_to_bids.py
@@ -0,0 +1,208 @@
+"""
+====================================
+07. Convert iEEG data to BIDS format
+====================================
+
+In this example, we use MNE-BIDS to create a BIDS-compatible directory of iEEG
+data. Specifically, we will follow these steps:
+
+1. Download some iEEG data from the
+   `PhysioBank database <https://physionet.org/physiobank/database>`_
+   and the
+   `MNE-ECoG ex <https://mne.tools/stable/auto_tutorials/misc/plot_ecog>`_.
+
+2. Load the data, extract information, and save in a new BIDS directory.
+
+3. Check the result and compare it with the standard.
+
+4. Confirm that written iEEG coordinates are the
+   same before :func:`write_raw_bids` was called.
+
+The iEEG data will be written by :func:`write_raw_bids` with
+the addition of extra metadata elements in the following files:
+
+    * sidecar.json
+    * electrodes.tsv
+    * coord_system.json
+    * events.tsv
+    * channels.tsv
+
+Compared to EEG data, the main differences are within the
+coord_system and electrodes files.
+For more information on these files,
+refer to the iEEG-BIDS specification.
+"""
+
+# Authors: Adam Li <adam2392@gmail.com>
+# License: BSD (3-clause)
+
+import os
+import tempfile
+from pprint import pprint
+
+import numpy as np
+from scipy.io import loadmat
+
+import mne
+from mne_bids import write_raw_bids, make_bids_basename, read_raw_bids
+from mne_bids.utils import print_dir_tree
+
+###############################################################################
+# Step 1: Download the data
+# -------------------------
+#
+# First, we need some data to work with. We will use the
+# data downloaded via MNE-Python's API.
+#
+# `<https://mne.tools/stable/generated/mne.datasets.misc.data_path>`_.
+#
+# Conveniently, there is already a data loading function available with
+# MNE-Python:
+
+# get MNE directory w/ example data
+mne_dir = mne.get_config('MNE_DATASETS_SAMPLE_PATH')
+
+# The electrode coords data are in the Matlab format: '.mat'.
+# This is easy to read in with :func:`scipy.io.loadmat` function.
+mat = loadmat(mne.datasets.misc.data_path() + '/ecog/sample_ecog.mat')
+ch_names = mat['ch_names'].tolist()
+ch_names = [x.strip() for x in ch_names]
+elec = mat['elec']  # electrode positions given in meters
+
+###############################################################################
+# Now we make a montage stating that the iEEG contacts are in MRI
+# coordinate system.
+montage = mne.channels.make_dig_montage(ch_pos=dict(zip(ch_names, elec)),
+                                        coord_frame='mri')
+print('Created %s channel positions' % len(ch_names))
+print(dict(zip(ch_names, elec)))
+
+###############################################################################
+# We also need to get some sample iEEG data,
+# so we will just generate random data from white noise.
+# Here is where you would use your own data if you had it.
+ieegdata = np.random.rand(len(ch_names), 1000)
+
+# We will create a :func:`mne.io.RawArray` object and
+# use the montage we created.
+info = mne.create_info(ch_names, 1000., 'ecog')
+raw = mne.io.RawArray(ieegdata, info)
+raw.set_montage(montage)
+
+###############################################################################
+# Let us confirm what our channel coordinates look like.
+
+# make a plot of the sensors in 2D plane
+raw.plot_sensors(ch_type='ecog')
+
+# Get the first 5 channels and show their locations.
+picks = mne.pick_types(raw.info, ecog=True)
+dig = [raw.info['dig'][pick] for pick in picks]
+chs = [raw.info['chs'][pick] for pick in picks]
+pos = np.array([ch['r'] for ch in dig[:5]])
+ch_names = np.array([ch['ch_name'] for ch in chs[:5]])
+print("The channel coordinates before writing into BIDS: ")
+pprint([x for x in zip(ch_names, pos)])
+
+###############################################################################
+# Step 2: Formatting as BIDS
+# --------------------------
+#
+# Now, let us format the `Raw` object into BIDS.
+
+###############################################################################
+# With this step, we have everything to start a new BIDS directory using
+# our data. To do that, we can use :func:`write_raw_bids`
+# Generally, :func:`write_raw_bids` tries to extract as much
+# meta data as possible from the raw data and then formats it in a BIDS
+# compatible way. :func:`write_raw_bids` takes a bunch of inputs, most of
+# which are however optional. The required inputs are:
+#
+# * :code:`raw`
+# * :code:`bids_basename`
+# * :code:`bids_root`
+#
+# ... as you can see in the docstring:
+print(write_raw_bids.__doc__)
+
+###############################################################################
+# Let us initialize some of the necessary data for the subject
+# There is a subject, and specific task for the dataset.
+subject_id = '001'  # zero padding to account for >100 subjects in this dataset
+task = 'testresteyes'
+
+# There is the root directory for where we will write our data.
+bids_root = os.path.join(mne_dir, 'ieegmmidb_bids')
+
+###############################################################################
+# Now we just need to specify a few iEEG details to make things work:
+# We need the basename of the dataset. In addition, write_raw_bids
+# requires a `filenames` of the Raw object to be non-empty, so since we
+# initialized the dataset from an array, we need to do a hack where we
+# temporarily save the data to disc before reading it back in.
+
+# Now convert our data to be in a new BIDS dataset.
+bids_basename = make_bids_basename(subject=subject_id,
+                                   task=task,
+                                   acquisition="ecog")
+
+# need to set the filenames if we are initializing data from array
+with tempfile.TemporaryDirectory() as tmp_root:
+    tmp_fpath = os.path.join(tmp_root, "test_raw.fif")
+    raw.save(tmp_fpath)
+    raw = mne.io.read_raw_fif(tmp_fpath)
+
+    # write `raw` to BIDS and anonymize it into BrainVision format
+    write_raw_bids(raw, bids_basename, bids_root=bids_root,
+                   anonymize=dict(daysback=30000), overwrite=True)
+
+###############################################################################
+# Step 3: Check and compare with standard
+# ---------------------------------------
+
+# Now we have written our BIDS directory.
+print_dir_tree(bids_root)
+
+###############################################################################
+# MNE-BIDS has created a suitable directory structure for us, and among other
+# meta data files, it started an `events.tsv` and `channels.tsv` and made an
+# initial `dataset_description.json` on top!
+#
+# Now it's time to manually check the BIDS directory and the meta files to add
+# all the information that MNE-BIDS could not infer. For instance, you must
+# describe iEEGReference and iEEGGround yourself. It's easy to find these by
+# searching for "n/a" in the sidecar files.
+#
+# `$ grep -i 'n/a' <bids_root>`
+#
+# Remember that there is a convenient javascript tool to validate all your BIDS
+# directories called the "BIDS-validator", available as a web version and a
+# command line tool:
+#
+# Web version: https://bids-standard.github.io/bids-validator/
+#
+# Command line tool: https://www.npmjs.com/package/bids-validator
+
+###############################################################################
+# Step 4: Plot output channels and check that they match!
+# -------------------------------------------------------
+#
+# Now we have written our BIDS directory. We can use
+# :func:`read_raw_bids` to read in the data.
+
+# read in the BIDS dataset and plot the coordinates
+bids_fname = bids_basename + "_ieeg.vhdr"
+raw = read_raw_bids(bids_fname, bids_root=bids_root)
+
+# get the first 5 channels and show their locations
+# this should match what was printed earlier.
+picks = mne.pick_types(raw.info, ecog=True)
+dig = [raw.info['dig'][pick] for pick in picks]
+chs = [raw.info['chs'][pick] for pick in picks]
+pos = np.array([ch['r'] for ch in dig[:5]])
+ch_names = np.array([ch['ch_name'] for ch in chs[:5]])
+print("The channel coordinates after writing into BIDS: ")
+pprint([x for x in zip(ch_names, pos)])
+
+# make a plot of the sensors in 2D plane
+raw.plot_sensors(ch_type='ecog')
diff --git a/mne_bids/read.py b/mne_bids/read.py
@@ -157,6 +157,65 @@ def _handle_events_reading(events_fname, raw):
     return raw
 
 
+def _handle_electrodes_reading(electrodes_fname, coord_frame, raw, verbose):
+    """Read associated electrodes.tsv and populate raw.
+
+    Handle xyz coordinates and coordinate frame of each channel.
+    Assumes units of coordinates are in 'm'.
+    """
+    logger.info('Reading electrode '
+                'coords from {}.'.format(electrodes_fname))
+    electrodes_dict = _from_tsv(electrodes_fname)
+    # First, make sure that ordering of names in channels.tsv matches the
+    # ordering of names in the raw data. The "name" column is mandatory in BIDS
+    ch_names_raw = list(raw.ch_names)
+    ch_names_tsv = electrodes_dict['name']
+
+    if ch_names_raw != ch_names_tsv:
+        msg = ('Channels do not correspond between raw data and the '
+               'channels.tsv file. For MNE-BIDS, the channel names in the '
+               'tsv MUST be equal and in the same order as the channels in '
+               'the raw data.\n\n'
+               '{} channels in tsv file: "{}"\n\n --> {}\n\n'
+               '{} channels in raw file: "{}"\n\n --> {}\n\n'
+               .format(len(ch_names_tsv), electrodes_fname, ch_names_tsv,
+                       len(ch_names_raw), raw.filenames, ch_names_raw)
+               )
+
+        # XXX: this could be due to MNE inserting a 'STI 014' channel as the
+        # last channel: In that case, we can work. --> Can be removed soon,
+        # because MNE will stop the synthesis of stim channels in the near
+        # future
+        if not (ch_names_raw[-1] == 'STI 014' and
+                ch_names_raw[:-1] == ch_names_tsv):
+            raise RuntimeError(msg)
+
+    if verbose:
+        print("The read in electrodes file is: \n", electrodes_dict)
+
+    # convert coordinates to float and create list of tuples
+    ch_names_raw = [x for i, x in enumerate(ch_names_raw)
+                    if electrodes_dict['x'][i] != "n/a"]
+    electrodes_dict['x'] = [float(x) for x in electrodes_dict['x']
+                            if x != "n/a"]
+    electrodes_dict['y'] = [float(x) for x in electrodes_dict['y']
+                            if x != "n/a"]
+    electrodes_dict['z'] = [float(x) for x in electrodes_dict['z']
+                            if x != "n/a"]
+
+    ch_locs = list(zip(electrodes_dict['x'],
+                       electrodes_dict['y'],
+                       electrodes_dict['z']))
+    ch_pos = dict(zip(ch_names_raw, ch_locs))
+
+    # create mne.DigMontage
+    montage = mne.channels.make_dig_montage(ch_pos=ch_pos,
+                                            coord_frame=coord_frame)
+    raw.set_montage(montage)
+
+    return raw
+
+
 def _handle_channels_reading(channels_fname, bids_fname, raw):
     """Read associated channels.tsv and populate raw.
 
@@ -296,6 +355,37 @@ def read_raw_bids(bids_fname, bids_root, extra_params=None,
     if channels_fname is not None:
         raw = _handle_channels_reading(channels_fname, bids_fname, raw)
 
+    # Try to find an associated electrodes.tsv and coordsystem.json
+    # to get information about the status and type of present channels
+    electrodes_fname = _find_matching_sidecar(bids_fname, bids_root,
+                                              'electrodes.tsv',
+                                              allow_fail=True)
+    coordsystem_fname = _find_matching_sidecar(bids_fname, bids_root,
+                                               'coordsystem.json',
+                                               allow_fail=True)
+    if electrodes_fname is not None:
+        if coordsystem_fname is None:
+            raise RuntimeError("BIDS mandates that the coordsystem.json "
+                               "should exist if electrodes.tsv does. "
+                               "Please create coordsystem.json for"
+                               "{}".format(bids_basename))
+        # Get MRI landmarks from the JSON sidecar
+        with open(coordsystem_fname, 'r') as fin:
+            coordsystem_json = json.load(fin)
+
+        # Get coordinate frames that electrode coordinates are in
+        if kind == "meg":
+            coord_frame = coordsystem_json['MEGCoordinateSystem']
+        elif kind == "ieeg":
+            coord_frame = coordsystem_json['iEEGCoordinateSystem']
+        else:  # noqa
+            raise RuntimeError("Kind {} not supported yet for "
+                               "coordsystem.json and "
+                               "electrodes.tsv.".format(kind))
+        # read in electrode coordinates and attach to raw
+        raw = _handle_electrodes_reading(electrodes_fname, coord_frame, raw,
+                                         verbose)
+
     # Try to find an associated sidecar.json to get information about the
     # recording snapshot
     sidecar_fname = _find_matching_sidecar(bids_fname, bids_root,