MRG, ENH: Add NIRSport support

doc/changes/latest.inc

@@ -19,12 +19,18 @@ Current (0.24.dev0)
 
 .. |New Contributor| replace:: **New Contributor**
 
+.. |David Julien| replace:: **David Julien**
+
+.. |Romain Derollepot| replace:: **Romain Derollepot**
+
 Enhancements
 ~~~~~~~~~~~~
 .. - Add something cool (:gh:`9192` **by new contributor** |New Contributor|_)
 
 - New function :func:`mne.chpi.get_chpi_info` to retrieve basic information about the cHPI system used when recording MEG data (:gh:`9369` by `Richard Höchenberger`_)
 
+- Add support for NIRSport devices to `mne.io.read_raw_nirx` (:gh:`9348` **by new contributor** |David Julien|_, **new contributor** |Romain Derollepot|_, `Robert Luke`_, and `Eric Larson`_)
+
 
 Bugs
 ~~~~

doc/changes/names.inc

@@ -389,3 +389,7 @@
 .. _Jack Zhang: https://github.com/jackz314
 
 .. _Felix Klotzsche: https://github.com/eioe
+
+.. _David Julien: https://github.com/Swy7ch
+
+.. _Romain Derollepot: https://github.com/rderollepot

doc/preprocessing.rst

@@ -71,6 +71,7 @@ Projections:
    annotate_flat
    annotate_movement
    annotate_muscle_zscore
+   annotate_nan
    compute_average_dev_head_t
    compute_current_source_density
    compute_fine_calibration

mne/datasets/utils.py

@@ -254,7 +254,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
-    releases = dict(testing='0.118', misc='0.9')
+    releases = dict(testing='0.119', misc='0.9')
     # And also update the "md5_hashes['testing']" variable below.
     # To update any other dataset, update the data archive itself (upload
     # an updated version) and update the md5 hash.
@@ -349,7 +349,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
         sample='12b75d1cb7df9dfb4ad73ed82f61094f',
         somato='32fd2f6c8c7eb0784a1de6435273c48b',
         spm='9f43f67150e3b694b523a21eb929ea75',
-        testing='0ca196c6dc4966570e14151cdf7ad4a5',
+        testing='2e7c60a055228928bd39f68892b3d488',
         multimodal='26ec847ae9ab80f58f204d09e2c08367',
         fnirs_motor='c4935d19ddab35422a69f3326a01fef8',
         opm='370ad1dcfd5c47e029e692c85358a374',

mne/io/nirx/nirx.py

@@ -17,19 +17,18 @@
 from ...annotations import Annotations
 from ...transforms import apply_trans, _get_trans
 from ...utils import (logger, verbose, fill_doc, warn, _check_fname,
-                      _validate_type)
+                      _validate_type, _check_option, _mask_to_onsets_offsets)
 
 
 @fill_doc
-def read_raw_nirx(fname, preload=False, verbose=None):
+def read_raw_nirx(fname, saturated='annotate', preload=False, verbose=None):
     """Reader for a NIRX fNIRS recording.
 
-    This function has only been tested with NIRScout devices.
-
     Parameters
     ----------
     fname : str
         Path to the NIRX data folder or header file.
+    %(saturated)s
     %(preload)s
     %(verbose)s
 
@@ -41,8 +40,12 @@ def read_raw_nirx(fname, preload=False, verbose=None):
     See Also
     --------
     mne.io.Raw : Documentation of attribute and methods.
+
+    Notes
+    -----
+    %(nirx_notes)s
     """
-    return RawNIRX(fname, preload, verbose)
+    return RawNIRX(fname, saturated, preload, verbose)
 
 
 def _open(fname):
@@ -57,20 +60,27 @@ class RawNIRX(BaseRaw):
     ----------
     fname : str
         Path to the NIRX data folder or header file.
+    %(saturated)s
     %(preload)s
     %(verbose)s
 
     See Also
     --------
     mne.io.Raw : Documentation of attribute and methods.
+
+    Notes
+    -----
+    %(nirx_notes)s
     """
 
     @verbose
-    def __init__(self, fname, preload=False, verbose=None):
+    def __init__(self, fname, saturated, preload=False, verbose=None):
         from ...externals.pymatreader import read_mat
         from ...coreg import get_mni_fiducials  # avoid circular import prob
         logger.info('Loading %s' % fname)
         _validate_type(fname, 'path-like', 'fname')
+        _validate_type(saturated, str, 'saturated')
+        _check_option('saturated', saturated, ('annotate', 'nan', 'ignore'))
         fname = str(fname)
         if fname.endswith('.hdr'):
             fname = op.dirname(op.abspath(fname))
@@ -81,23 +91,43 @@ def __init__(self, fname, preload=False, verbose=None):
         files = dict()
         keys = ('hdr', 'inf', 'set', 'tpl', 'wl1', 'wl2',
                 'config.txt', 'probeInfo.mat')
+        nan_mask = dict()
         for key in keys:
             files[key] = glob.glob('%s/*%s' % (fname, key))
+            fidx = 0
             if len(files[key]) != 1:
-                raise RuntimeError('Expect one %s file, got %d' %
-                                   (key, len(files[key]),))
-            files[key] = files[key][0]
+                if key not in ('wl1', 'wl2'):
+                    raise RuntimeError(
+                        f'Need one {key} file, got {len(files[key])}')
+                noidx = np.where(['nosatflags_' in op.basename(x)
+                                  for x in files[key]])[0]
+                if len(noidx) != 1 or len(files[key]) != 2:
+                    raise RuntimeError(
+                        f'Need one nosatflags and one standard {key} file, '
+                        f'got {len(files[key])}')
+                # Here two files have been found, one that is called
+                # no sat flags. The nosatflag file has no NaNs in it.
+                noidx = noidx[0]
+                if saturated == 'ignore':
+                    # Ignore NaN and return values
+                    fidx = noidx
+                elif saturated == 'nan':
+                    # Return NaN
+                    fidx = 0 if noidx == 1 else 1
+                else:
+                    assert saturated == 'annotate'  # guaranteed above
+                    fidx = noidx
+                    nan_mask[key] = files[key][0 if noidx == 1 else 1]
+            files[key] = files[key][fidx]
         if len(glob.glob('%s/*%s' % (fname, 'dat'))) != 1:
             warn("A single dat file was expected in the specified path, but "
                  "got %d. This may indicate that the file structure has been "
                  "modified since the measurement was saved." %
                  (len(glob.glob('%s/*%s' % (fname, 'dat')))))
 
         # Read number of rows/samples of wavelength data
-        last_sample = -1
         with _open(files['wl1']) as fid:
-            for line in fid:
-                last_sample += 1
+            last_sample = fid.read().count('\n') - 1
 
         # Read header file
         # The header file isn't compliant with the configparser. So all the
@@ -112,7 +142,8 @@ def __init__(self, fname, preload=False, verbose=None):
         if hdr['GeneralInfo']['NIRStar'] not in ['"15.0"', '"15.2"', '"15.3"']:
             raise RuntimeError('MNE does not support this NIRStar version'
                                ' (%s)' % (hdr['GeneralInfo']['NIRStar'],))
-        if "NIRScout" not in hdr['GeneralInfo']['Device']:
+        if "NIRScout" not in hdr['GeneralInfo']['Device'] \
+                and "NIRSport" not in hdr['GeneralInfo']['Device']:
             warn("Only import of data from NIRScout devices have been "
                  "thoroughly tested. You are using a %s device. " %
                  hdr['GeneralInfo']['Device'])
@@ -299,43 +330,74 @@ def prepend(li, str):
             'sd_index': req_ind,
             'files': files,
             'bounds': bounds,
+            'nan_mask': nan_mask,
         }
+        # Get our saturated mask
+        annot_mask = None
+        for ki, key in enumerate(('wl1', 'wl2')):
+            if nan_mask.get(key, None) is None:
+                continue
+            mask = np.isnan(_read_csv_rows_cols(
+                nan_mask[key], 0, last_sample + 1, req_ind, {0: 0, 1: None}).T)
+            if saturated == 'nan':
+                nan_mask[key] = mask
+            else:
+                assert saturated == 'annotate'
+                if annot_mask is None:
+                    annot_mask = np.zeros(
+                        (len(info['ch_names']) // 2, last_sample + 1), bool)
+                annot_mask |= mask
+                nan_mask[key] = None  # shouldn't need again
 
         super(RawNIRX, self).__init__(
             info, preload, filenames=[fname], last_samps=[last_sample],
             raw_extras=[raw_extras], verbose=verbose)
 
+        # make onset/duration/description
+        onset, duration, description, ch_names = list(), list(), list(), list()
+        if annot_mask is not None:
+            for ci, mask in enumerate(annot_mask):
+                on, dur = _mask_to_onsets_offsets(mask)
+                on = on / info['sfreq']
+                dur = dur / info['sfreq']
+                dur -= on
+                onset.extend(on)
+                duration.extend(dur)
+                description.extend(['BAD_SATURATED'] * len(on))
+                ch_names.extend([self.ch_names[2 * ci:2 * ci + 2]] * len(on))
+
         # Read triggers from event file
         if op.isfile(files['hdr'][:-3] + 'evt'):
             with _open(files['hdr'][:-3] + 'evt') as fid:
                 t = [re.findall(r'(\d+)', line) for line in fid]
-            onset = np.zeros(len(t), float)
-            duration = np.zeros(len(t), float)
-            description = [''] * len(t)
-            for t_idx in range(len(t)):
-                binary_value = ''.join(t[t_idx][1:])[::-1]
-                trigger_frame = float(t[t_idx][0])
-                onset[t_idx] = (trigger_frame) * (1.0 / samplingrate)
-                duration[t_idx] = 1.0  # No duration info stored in files
-                description[t_idx] = int(binary_value, 2) * 1.
-            annot = Annotations(onset, duration, description)
-            self.set_annotations(annot)
+            for t_ in t:
+                binary_value = ''.join(t_[1:])[::-1]
+                trigger_frame = float(t_[0])
+                onset.append(trigger_frame / samplingrate)
+                duration.append(1.)  # No duration info stored in files
+                description.append(float(int(binary_value, 2)))
+                ch_names.append(list())
+        annot = Annotations(onset, duration, description, ch_names=ch_names)
+        self.set_annotations(annot)
 
     def _read_segment_file(self, data, idx, fi, start, stop, cals, mult):
         """Read a segment of data from a file.
 
         The NIRX machine records raw data as two different wavelengths.
         The returned data interleaves the wavelengths.
         """
-        sdindex = self._raw_extras[fi]['sd_index']
+        sd_index = self._raw_extras[fi]['sd_index']
 
-        wls = [
-            _read_csv_rows_cols(
+        wls = list()
+        for key in ('wl1', 'wl2'):
+            d = _read_csv_rows_cols(
                 self._raw_extras[fi]['files'][key],
-                start, stop, sdindex,
+                start, stop, sd_index,
                 self._raw_extras[fi]['bounds'][key]).T
-            for key in ('wl1', 'wl2')
-        ]
+            nan_mask = self._raw_extras[fi]['nan_mask'].get(key, None)
+            if nan_mask is not None:
+                d[nan_mask[:, start:stop]] = np.nan
+            wls.append(d)
 
         # TODO: Make this more efficient by only indexing above what we need.
         # For now let's just construct the full data matrix and index.
@@ -350,7 +412,10 @@ def _read_segment_file(self, data, idx, fi, start, stop, cals, mult):
 def _read_csv_rows_cols(fname, start, stop, cols, bounds):
     with open(fname, 'rb') as fid:
         fid.seek(bounds[start])
-        data = fid.read(bounds[stop] - bounds[start]).decode('latin-1')
+        args = list()
+        if bounds[1] is not None:
+            args.append(bounds[stop] - bounds[start])
+        data = fid.read(*args).decode('latin-1')
         x = np.fromstring(data, float, sep=' ')
     x.shape = (stop - start, -1)
     x = x[:, cols]

mne/io/nirx/tests/test_nirx.py

@@ -7,6 +7,7 @@
 import shutil
 import os
 import datetime as dt
+import numpy as np
 
 import pytest
 from numpy.testing import assert_allclose, assert_array_equal
@@ -15,6 +16,7 @@
 from mne.datasets.testing import data_path, requires_testing_data
 from mne.io import read_raw_nirx
 from mne.io.tests.test_raw import _test_raw_reader
+from mne.preprocessing import annotate_nan
 from mne.transforms import apply_trans, _get_trans
 from mne.preprocessing.nirs import source_detector_distances,\
     short_channels
@@ -31,6 +33,95 @@
                                 'NIRx', 'nirscout', 'nirx_15_3_recording')
 
 
+# This file has no saturated sections
+nirsport1_wo_sat = op.join(data_path(download=False), 'NIRx', 'nirsport_v1',
+                           'nirx_15_3_recording_wo_saturation')
+# This file has saturation, but not on the optode pairing in montage
+nirsport1_w_sat = op.join(data_path(download=False), 'NIRx', 'nirsport_v1',
+                          'nirx_15_3_recording_w_saturation_'
+                          'not_on_montage_channels')
+# This file has saturation in channels of interest
+nirsport1_w_fullsat = op.join(data_path(download=False), 'NIRx', 'nirsport_v1',
+                              'nirx_15_3_recording_w_'
+                              'saturation_on_montage_channels')
+
+
+@requires_testing_data
+@pytest.mark.filterwarnings('ignore:.*Extraction of measurement.*:')
+def test_nirsport_v1_wo_sat():
+    """Test NIRSport1 file with no saturation."""
+    raw = read_raw_nirx(nirsport1_wo_sat, preload=True)
+
+    # Test data import
+    assert raw._data.shape == (26, 164)
+    assert raw.info['sfreq'] == 10.416667
+
+    # By default real data is returned
+    assert np.sum(np.isnan(raw.get_data())) == 0
+
+    raw = read_raw_nirx(nirsport1_wo_sat, preload=True, saturated='nan')
+    data = raw.get_data()
+    assert data.shape == (26, 164)
+    assert np.sum(np.isnan(data)) == 0
+
+    raw = read_raw_nirx(nirsport1_wo_sat, saturated='annotate')
+    data = raw.get_data()
+    assert data.shape == (26, 164)
+    assert np.sum(np.isnan(data)) == 0
+
+
+@pytest.mark.filterwarnings('ignore:.*Extraction of measurement.*:')
+@requires_testing_data
+def test_nirsport_v1_w_sat():
+    """Test NIRSport1 file with NaNs but not in channel of interest."""
+    raw = read_raw_nirx(nirsport1_w_sat)
+
+    # Test data import
+    data = raw.get_data()
+    assert data.shape == (26, 176)
+    assert raw.info['sfreq'] == 10.416667
+    assert np.sum(np.isnan(data)) == 0
+
+    raw = read_raw_nirx(nirsport1_w_sat, saturated='nan')
+    data = raw.get_data()
+    assert data.shape == (26, 176)
+    assert np.sum(np.isnan(data)) == 0
+
+    raw = read_raw_nirx(nirsport1_w_sat, saturated='annotate')
+    data = raw.get_data()
+    assert data.shape == (26, 176)
+    assert np.sum(np.isnan(data)) == 0
+
+
+@pytest.mark.filterwarnings('ignore:.*Extraction of measurement.*:')
+@requires_testing_data
+@pytest.mark.parametrize('preload', (True, False))
+def test_nirsport_v1_w_bad_sat(preload):
+    """Test NIRSport1 file with NaNs."""
+    fname = nirsport1_w_fullsat
+    raw = read_raw_nirx(fname, preload=preload)
+    data = raw.get_data()
+    assert not np.isnan(data).any()
+    assert len(raw.annotations) == 5
+    # annotated version and ignore should have same data but different annot
+    raw_ignore = read_raw_nirx(fname, saturated='ignore', preload=preload)
+    assert_allclose(raw_ignore.get_data(), data)
+    assert len(raw_ignore.annotations) == 2
+    assert not any('NAN' in d for d in raw_ignore.annotations.description)
+    # nan version should not have same data, but we can give it the same annot
+    raw_nan = read_raw_nirx(fname, saturated='nan', preload=preload)
+    data_nan = raw_nan.get_data()
+    assert np.isnan(data_nan).any()
+    assert not np.allclose(raw_nan.get_data(), data)
+    raw_nan_annot = raw_ignore.copy()
+    raw_nan_annot.set_annotations(annotate_nan(raw_nan))
+    use_mask = np.where(raw.annotations.description == 'BAD_SATURATED')
+    for key in ('onset', 'duration'):
+        a = getattr(raw_nan_annot.annotations, key)[::2]  # one ch in each
+        b = getattr(raw.annotations, key)[use_mask]  # two chs in each
+        assert_allclose(a, b)
+
+
 @requires_testing_data
 def test_nirx_hdr_load():
     """Test reading NIRX files using path to header file."""