MRG: Add cHPI info to MEG sidecars (#794)

* Add cHPI info to MEG sidecars

Fixes #792

* Use warnings filter

* Fix tests for MNE stable

* Remove accidentally added version constraint

* Call BIDS validator

doc/requirements.txt

@@ -6,3 +6,4 @@ numpydoc
 matplotlib
 pillow
 pandas
+setuptools

doc/whats_new.rst

@@ -40,6 +40,7 @@ Enhancements
 - Add writing simultaneous EEG-iEEG recordings via :func:`mne_bids.write_raw_bids`. The desired output datatype must be specified in the :class:`mne_bids.BIDSPath` object, by `Richard Köhler`_ (:gh:`774`)
 - :func:`mne_bids.write_raw_bids` gained a new keyword argument ``symlink``, which allows to create symbolic links to the original data files instead of copying them over. Currently works for ``FIFF`` files on macOS and Linux, by `Richard Höchenberger`_ (:gh:`778`)
 - :class:`mne_bids.BIDSPath` now has property getter and setter methods for all BIDS entities, i.e., you can now do things like ``bids_path.subject = 'foo'`` and don't have to resort to ``bids_path.update()``. This also ensures you'll get proper completion suggestions from your favorite Python IDE, by `Richard Höchenberger`_ (:gh:`786`)
+- :func:`mne_bids.write_raw_bids` now stores information about continuous head localization measurements (e.g., Elekta/Neuromag cHPI) in the MEG sidecar file, by `Richard Höchenberger`_ (:gh:`794`)
 
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
@@ -53,6 +54,7 @@ Requirements
 ^^^^^^^^^^^^
 
 - For downloading `OpenNeuro <https://openneuro.org>`_ datasets, ``openneuro-py`` is now required to run the examples and build the documentation, by `Alex Rockhill`_ (:gh:`753`)
+- MNE-BIDS now depends on `setuptools <https://setuptools.readthedocs.io>`_. This package is normally installed by your Python distribution automatically, so we don't expect any users to be affected by this change, by `Richard Höchenberger`_ (:gh:`794`)
 
 Bug fixes
 ^^^^^^^^^

mne_bids/read.py

@@ -241,7 +241,7 @@ def _handle_scans_reading(scans_fname, raw, bids_path, verbose=False):
 
 
 def _handle_info_reading(sidecar_fname, raw, verbose=None):
-    """Read associated sidecar.json and populate raw.
+    """Read associated sidecar JSON and populate raw.
 
     Handle PowerLineFrequency of recording.
     """
@@ -268,6 +268,29 @@ def _handle_info_reading(sidecar_fname, raw, verbose=None):
                              "Sidecar JSON is -> {} ".format(line_freq))
 
     raw.info["line_freq"] = line_freq
+
+    # get cHPI info
+    chpi = sidecar_json.get('ContinuousHeadLocalization')
+    if chpi is None:
+        # no cHPI info in the sidecar – leave raw.info unchanged
+        pass
+    elif chpi is True:
+        hpi_freqs = sidecar_json['HeadCoilFrequency']
+        hpi_freqs_data, _, _ = mne.chpi.get_chpi_info(raw.info)
+        if not np.allclose(hpi_freqs, hpi_freqs_data):
+            raise ValueError(
+                f'The cHPI coil frequencies in the sidecar file '
+                f'{sidecar_fname}:\n    {hpi_freqs}\ndiffer from what is '
+                f'stored in the raw data:\n    {hpi_freqs_data}\n'
+                f'Cannot proceed.'
+            )
+    else:
+        if raw.info['hpi_subsystem']:
+            logger.info('Dropping cHPI information stored in raw data, '
+                        'following specification in sidecar file')
+        raw.info['hpi_subsystem'] = None
+        raw.info['hpi_meas'] = []
+
     return raw
 
 

mne_bids/tests/test_read.py

@@ -22,7 +22,7 @@
                             category=ImportWarning)
     import mne
 from mne.io.constants import FIFF
-from mne.utils import requires_nibabel, object_diff
+from mne.utils import requires_nibabel, object_diff, requires_version
 from mne.utils import assert_dig_allclose
 from mne.datasets import testing, somato
 
@@ -59,11 +59,15 @@
 somato_raw_fname = op.join(somato_path, 'sub-01', 'meg',
                            'sub-01_task-somato_meg.fif')
 
+# Data with cHPI info
+raw_fname_chpi = op.join(data_path, 'SSS', 'test_move_anon_raw.fif')
+
 warning_str = dict(
     channel_unit_changed='ignore:The unit for chann*.:RuntimeWarning:mne',
     meas_date_set_to_none="ignore:.*'meas_date' set to None:RuntimeWarning:"
                           "mne",
     nasion_not_found='ignore:.*nasion not found:RuntimeWarning:mne',
+    maxshield='ignore:.*Internal Active Shielding:RuntimeWarning:mne'
 )
 
 
@@ -418,7 +422,7 @@ def test_handle_scans_reading(tmpdir):
 
 @pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
 def test_handle_info_reading(tmpdir):
-    """Test reading information from a BIDS sidecar.json file."""
+    """Test reading information from a BIDS sidecar JSON file."""
     # read in USA dataset, so it should find 50 Hz
     raw = _read_raw_fif(raw_fname)
 
@@ -502,6 +506,46 @@ def test_handle_info_reading(tmpdir):
         assert raw.info['line_freq'] == 55
 
 
+@requires_version('mne', '0.24.dev0')
+@pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
+@pytest.mark.filterwarnings(warning_str['maxshield'])
+def test_handle_chpi_reading(tmpdir):
+    """Test reading of cHPI information."""
+    raw = _read_raw_fif(raw_fname_chpi, allow_maxshield=True)
+    root = tmpdir.mkdir('chpi')
+    bids_path = BIDSPath(subject='01', session='01',
+                         task='audiovisual', run='01',
+                         root=root, datatype='meg')
+    bids_path = write_raw_bids(raw, bids_path)
+
+    raw_read = read_raw_bids(bids_path)
+    assert raw_read.info['hpi_subsystem'] is not None
+
+    # cause conflicts between cHPI info in sidecar and raw data
+    meg_json_path = bids_path.copy().update(suffix='meg', extension='.json')
+    with open(meg_json_path, 'r', encoding='utf-8') as f:
+        meg_json_data = json.load(f)
+
+    # cHPI frequency mismatch
+    meg_json_data_freq_mismatch = meg_json_data.copy()
+    meg_json_data_freq_mismatch['HeadCoilFrequency'][0] = 123
+    with open(meg_json_path, 'w', encoding='utf-8') as f:
+        json.dump(meg_json_data_freq_mismatch, f)
+
+    with pytest.raises(ValueError, match='cHPI coil frequencies'):
+        raw_read = read_raw_bids(bids_path)
+
+    # cHPI "off" according to sidecar, but present in the data
+    meg_json_data_chpi_mismatch = meg_json_data.copy()
+    meg_json_data_chpi_mismatch['ContinuousHeadLocalization'] = False
+    with open(meg_json_path, 'w', encoding='utf-8') as f:
+        json.dump(meg_json_data_chpi_mismatch, f)
+
+    raw_read = read_raw_bids(bids_path)
+    assert raw_read.info['hpi_subsystem'] is None
+    assert raw_read.info['hpi_meas'] == []
+
+
 @pytest.mark.filterwarnings(warning_str['nasion_not_found'])
 @pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
 def test_handle_eeg_coords_reading(tmpdir):

mne_bids/tests/test_write.py

@@ -20,7 +20,8 @@
 import json
 from pathlib import Path
 import codecs
-from distutils.version import LooseVersion
+
+from pkg_resources import parse_version
 
 import numpy as np
 from numpy.testing import assert_array_equal, assert_array_almost_equal
@@ -74,7 +75,8 @@
                          'pytest.PytestUnraisableExceptionWarning',
     encountered_data_in='ignore:Encountered data in*.:RuntimeWarning:mne',
     edf_warning=r'ignore:^EDF\/EDF\+\/BDF files contain two fields .*'
-                r':RuntimeWarning:mne'
+                r':RuntimeWarning:mne',
+    maxshield='ignore:.*Internal Active Shielding:RuntimeWarning:mne'
 )
 
 
@@ -374,6 +376,7 @@ def test_line_freq(line_freq, _bids_validate, tmpdir):
 
 @requires_version('pybv', '0.4')
 @pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
+@pytest.mark.filterwarnings(warning_str['maxshield'])
 def test_fif(_bids_validate, tmpdir):
     """Test functionality of the write_raw_bids conversion for fif."""
     bids_root = tmpdir.mkdir('bids1')
@@ -606,15 +609,21 @@ def test_fif(_bids_validate, tmpdir):
     write_raw_bids(raw, bids_path, events_data=events, event_id=event_id,
                    overwrite=True)
 
-    # test whether extra points in raw.info['dig'] are correctly used
-    # to set DigitizedHeadShape in the json sidecar
     meg_json = _find_matching_sidecar(
         bids_path.copy().update(root=bids_root),
         suffix='meg', extension='.json')
+
     with open(meg_json, 'r') as fin:
         meg_json_data = json.load(fin)
-        # unchanged sample data includes Extra points
-        assert meg_json_data['DigitizedHeadPoints'] is True
+
+    # no cHPI info is contained in the sample data
+    assert meg_json_data['ContinuousHeadLocalization'] is False
+    assert meg_json_data['HeadCoilFrequency'] == []
+
+    # test whether extra points in raw.info['dig'] are correctly used
+    # to set DigitizedHeadShape in the json sidecar
+    # unchanged sample data includes Extra points
+    assert meg_json_data['DigitizedHeadPoints'] is True
 
     # drop extra points from raw.info['dig'] and write again
     raw_no_extra_points = raw.copy()
@@ -636,6 +645,28 @@ def test_fif(_bids_validate, tmpdir):
         # DigitizedHeadPoints should be false
         assert meg_json_data['DigitizedHeadPoints'] is False
 
+    # test data with cHPI info
+    data_path = testing.data_path()
+    raw_fname = op.join(data_path, 'SSS', 'test_move_anon_raw.fif')
+    raw = _read_raw_fif(raw_fname, allow_maxshield=True)
+
+    root = tmpdir.mkdir('chpi')
+    bids_path = bids_path.copy().update(root=root, datatype='meg')
+    bids_path = write_raw_bids(raw, bids_path)
+    _bids_validate(bids_path.root)
+
+    meg_json = bids_path.copy().update(suffix='meg', extension='.json')
+    with open(meg_json, 'r') as fin:
+        meg_json_data = json.load(fin)
+
+    if parse_version(mne.__version__) > parse_version('0.23'):
+        assert meg_json_data['ContinuousHeadLocalization'] is True
+        assert_array_almost_equal(meg_json_data['HeadCoilFrequency'],
+                                  [83., 143., 203., 263., 323.])
+    else:
+        assert meg_json_data['ContinuousHeadLocalization'] is False
+        assert meg_json_data['HeadCoilFrequency'] == []
+
 
 @pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
 def test_fif_dtype(_bids_validate, tmpdir):
@@ -1374,7 +1405,7 @@ def test_bdf(_bids_validate, tmpdir):
         write_raw_bids(mne.concatenate_raws([raw.copy(), raw]), bids_path,
                        overwrite=True)
 
-    if LooseVersion(mne.__version__) >= LooseVersion('0.23'):
+    if parse_version(mne.__version__) >= parse_version('0.23'):
         raw.info['sfreq'] -= 10  # changes raw.times, but retains its dimension
     else:
         raw._times = raw._times / 5

mne_bids/write.py

@@ -16,6 +16,8 @@
 import shutil
 from collections import defaultdict, OrderedDict
 
+from pkg_resources import parse_version
+
 import numpy as np
 from scipy import linalg
 import mne
@@ -654,6 +656,16 @@ def _sidecar_json(raw, task, manufacturer, fname, datatype, overwrite=False,
                     raw.filenames[0].endswith('.fif'):
                 digitized_head_points = True
 
+    # Compile cHPI information, if any.
+    chpi = False
+    hpi_freqs = np.array([])
+    if (datatype == 'meg' and
+            parse_version(mne.__version__) > parse_version('0.23')):
+        hpi_freqs, _, _ = mne.chpi.get_chpi_info(info=raw.info,
+                                                 on_missing='ignore')
+        if hpi_freqs.size > 0:
+            chpi = True
+
     # Define datatype-specific JSON dictionaries
     ch_info_json_common = [
         ('TaskName', task),
@@ -668,7 +680,9 @@ def _sidecar_json(raw, task, manufacturer, fname, datatype, overwrite=False,
         ('DigitizedLandmarks', digitized_landmark),
         ('DigitizedHeadPoints', digitized_head_points),
         ('MEGChannelCount', n_megchan),
-        ('MEGREFChannelCount', n_megrefchan)]
+        ('MEGREFChannelCount', n_megrefchan),
+        ('ContinuousHeadLocalization', chpi),
+        ('HeadCoilFrequency', list(hpi_freqs))]
     ch_info_json_eeg = [
         ('EEGReference', 'n/a'),
         ('EEGGround', 'n/a'),
@@ -1111,8 +1125,8 @@ def write_raw_bids(raw, bids_path, events_data=None,
         events = mne.find_events(raw, min_duration=0.002)
         write_raw_bids(..., events_data=events)
 
-    See the documentation of `mne.find_events` for more information on event
-    extraction from ``STIM`` channels.
+    See the documentation of :func:`mne.find_events` for more information on
+    event extraction from ``STIM`` channels.
 
     When anonymizing ``.edf`` files, then the file format for EDF limits
     how far back we can set the recording date. Therefore, all anonymized
@@ -1213,14 +1227,14 @@ def write_raw_bids(raw, bids_path, events_data=None,
         raise ValueError(msg)
 
     datatype = _handle_datatype(raw, bids_path.datatype, verbose)
-    bids_path = bids_path.copy()
-    bids_path = bids_path.update(
-        datatype=datatype, suffix=datatype, extension=ext)
+    bids_path = (bids_path.copy()
+                 .update(datatype=datatype, suffix=datatype, extension=ext))
 
     # check whether the info provided indicates that the data is emptyroom
     # data
     emptyroom = False
-    if bids_path.subject == 'emptyroom' and bids_path.task == 'noise':
+    if (bids_path.datatype == 'meg' and bids_path.subject == 'emptyroom' and
+            bids_path.task == 'noise'):
         emptyroom = True
         # check the session date provided is consistent with the value in raw
         meas_date = raw.info.get('meas_date', None)

setup.cfg

@@ -39,6 +39,7 @@ install_requires =
   mne >= 0.21.2
   numpy>=1.15.4
   scipy>=1.1.0
+  setuptools
 packages = find:
 include_package_data = True