diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index adafdbb054d..5a35de9143c 100644 --- a/doc/changes/latest.inc +++ b/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 ~~~~ diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 3c2e705a324..b03532a5dba 100644 --- a/doc/changes/names.inc +++ b/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 diff --git a/doc/preprocessing.rst b/doc/preprocessing.rst index 56c7d740c97..a4e1dded46e 100644 --- a/doc/preprocessing.rst +++ b/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 diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 07e10357d58..9b824fbe589 100644 --- a/mne/datasets/utils.py +++ b/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', diff --git a/mne/io/nirx/nirx.py b/mne/io/nirx/nirx.py index f2d31f17a93..ec3d7f5de5b 100644 --- a/mne/io/nirx/nirx.py +++ b/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,12 +91,34 @@ 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 " @@ -94,10 +126,8 @@ def __init__(self, fname, preload=False, verbose=None): (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,27 +330,55 @@ 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. @@ -327,15 +386,18 @@ def _read_segment_file(self, data, idx, fi, start, stop, cals, mult): 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] diff --git a/mne/io/nirx/tests/test_nirx.py b/mne/io/nirx/tests/test_nirx.py index 4f74f191878..ed92afc80d3 100644 --- a/mne/io/nirx/tests/test_nirx.py +++ b/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.""" diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py index 419870136db..618b0a14976 100644 --- a/mne/preprocessing/__init__.py +++ b/mne/preprocessing/__init__.py @@ -29,4 +29,5 @@ from ._regress import regress_artifact from ._fine_cal import (compute_fine_calibration, read_fine_calibration, write_fine_calibration) +from .annotate_nan import annotate_nan from .interpolate import equalize_bads diff --git a/mne/preprocessing/annotate_nan.py b/mne/preprocessing/annotate_nan.py new file mode 100644 index 00000000000..d7d7ef2991b --- /dev/null +++ b/mne/preprocessing/annotate_nan.py @@ -0,0 +1,35 @@ +# Author: David Julien +# +# License: BSD (3-clause) + +import numpy as np + +from ..annotations import Annotations +from ..utils import verbose +from .artifact_detection import _annotations_from_mask + + +@verbose +def annotate_nan(raw, *, verbose=None): + """Detect segments with NaN and return a new Annotations instance. + + Parameters + ---------- + raw : instance of Raw + Data to find segments with NaN values. + %(verbose)s + + Returns + ------- + annot : instance of Annotations + New channel-specific annotations for the data. + """ + data, times = raw.get_data(return_times=True) + onsets, durations, ch_names = list(), list(), list() + for row, ch_name in zip(data, raw.ch_names): + annot = _annotations_from_mask(times, np.isnan(row), 'BAD_NAN') + onsets.extend(annot.onset) + durations.extend(annot.duration) + ch_names.extend([[ch_name]] * len(annot)) + annot = Annotations(onsets, durations, 'BAD_NAN', ch_names=ch_names) + return annot diff --git a/mne/utils/docs.py b/mne/utils/docs.py index 68c9c89c7cb..efe0a30cffa 100644 --- a/mne/utils/docs.py +++ b/mne/utils/docs.py @@ -58,6 +58,36 @@ .. versionadded:: 0.22 """ % (_on_missing_base,) +docdict['saturated'] = """\ +saturated : str + Replace saturated segments of data with NaNs, can be: + + ``"ignore"`` + The measured data is returned, even if it contains measurements + while the amplifier was saturated. + ``"nan"`` + The returned data will contain NaNs during time segments + when the amplifier was saturated. + ``"annotate"`` (default) + The returned data will contain annotations specifying + sections the saturate segments. + + This argument will only be used if there is no .nosatflags file + (only if a NIRSport device is used and saturation occurred). + + .. versionadded:: 0.24 +""" +docdict['nirx_notes'] = """\ +This function has only been tested with NIRScout and NIRSport1 devices. + +The NIRSport device can detect if the amplifier is saturated. +Starting from NIRStar 14.2, those saturated values are replaced by NaNs +in the standard .wlX files. +The raw unmodified measured values are stored in another file +called .nosatflags_wlX. As NaN values can cause unexpected behaviour with +mathematical functions the default behaviour is to return the +saturated data. +""" # Cropping docdict['include_tmax'] = """