diff --git a/doc/conf.py b/doc/conf.py index 930c2cf41d9..adb4ab238bb 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -631,7 +631,7 @@ def reset_warnings(gallery_conf, fname): 'n_dipoles_fwd', # Undocumented (on purpose) 'RawKIT', 'RawEximia', 'RawEGI', 'RawEEGLAB', 'RawEDF', 'RawCTF', 'RawBTi', - 'RawBrainVision', 'RawCurry', 'RawNIRX', 'RawGDF', + 'RawBrainVision', 'RawCurry', 'RawNIRX', 'RawGDF', 'RawSNIRF', # sklearn subclasses 'mapping', 'to', 'any', # unlinkable diff --git a/doc/python_reference.rst b/doc/python_reference.rst index c74edbdd512..4a11754e659 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -65,6 +65,7 @@ Reading raw data read_raw_kit read_raw_nicolet read_raw_nirx + read_raw_snirf read_raw_eeglab read_raw_brainvision read_raw_egi diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 1c923897ec3..00eee7994fa 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -239,7 +239,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.94', misc='0.6') + releases = dict(testing='0.95', misc='0.6') # And also update the "md5_hashes['testing']" variable below. # To update any other dataset, update the data archive itself (upload @@ -326,7 +326,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, sample='12b75d1cb7df9dfb4ad73ed82f61094f', somato='ea825966c0a1e9b2f84e3826c5500161', spm='9f43f67150e3b694b523a21eb929ea75', - testing='42b7d22edf5c4dbd0b47c59534b8815b', + testing='542b00d7335681c9818860b094d91205', multimodal='26ec847ae9ab80f58f204d09e2c08367', fnirs_motor='c4935d19ddab35422a69f3326a01fef8', opm='370ad1dcfd5c47e029e692c85358a374', diff --git a/mne/io/__init__.py b/mne/io/__init__.py index 99c0281abe0..8a723be59b1 100644 --- a/mne/io/__init__.py +++ b/mne/io/__init__.py @@ -45,6 +45,7 @@ from .eeglab import read_raw_eeglab, read_epochs_eeglab from .eximia import read_raw_eximia from .nirx import read_raw_nirx +from .snirf import read_raw_snirf from .fieldtrip import (read_raw_fieldtrip, read_epochs_fieldtrip, read_evoked_fieldtrip) from ._read_raw import read_raw diff --git a/mne/io/snirf/__init__.py b/mne/io/snirf/__init__.py new file mode 100644 index 00000000000..88810ac35d2 --- /dev/null +++ b/mne/io/snirf/__init__.py @@ -0,0 +1,7 @@ +"""SNIRF module for conversion to FIF.""" + +# Author: Robert Luke +# +# License: BSD (3-clause) + +from ._snirf import read_raw_snirf diff --git a/mne/io/snirf/_snirf.py b/mne/io/snirf/_snirf.py new file mode 100644 index 00000000000..c5bfd2a9c14 --- /dev/null +++ b/mne/io/snirf/_snirf.py @@ -0,0 +1,198 @@ +# Authors: Robert Luke +# +# License: BSD (3-clause) + +import re as re +import numpy as np + +from ..base import BaseRaw +from ..meas_info import create_info +from ...annotations import Annotations +from ...utils import logger, verbose, fill_doc, warn, requires_h5py + + +@requires_h5py +@fill_doc +def read_raw_snirf(fname, preload=False, verbose=None): + """Reader for a SNIRF fNIRS recording. + + This reader supports the .snirf file type only. Not the .jnirs version. + Only support for continuous wave SNIRF data is currently implemented. + + Parameters + ---------- + fname : str + Path to the SNIRF data file. + %(preload)s + %(verbose)s + + Returns + ------- + raw : instance of RawSNIRF + A Raw object containing SNIRF data. + + See Also + -------- + mne.io.Raw : Documentation of attribute and methods. + """ + return RawSNIRF(fname, preload, verbose) + + +def _open(fname): + return open(fname, 'r', encoding='latin-1') + + +@requires_h5py +@fill_doc +class RawSNIRF(BaseRaw): + """Raw object from a SNIRF fNIRS file. + + Parameters + ---------- + fname : str + Path to the SNIRF data file. + %(preload)s + %(verbose)s + + See Also + -------- + mne.io.Raw : Documentation of attribute and methods. + """ + + @verbose + def __init__(self, fname, preload=False, verbose=None): + from ...externals.pymatreader.utils import _import_h5py + h5py = _import_h5py() + logger.info('Loading %s' % fname) + + dat = h5py.File(fname, 'r') + + if 'data2' in dat['nirs'].keys(): + print("File contains multiple recordings. " + "This is not supported.") + + data = np.array(dat.get('/nirs/data1/dataTimeSeries')).T + last_samps = data.shape[1] - 1 + + samplingrate_raw = np.array(dat.get('nirs/data1/time')) + sampling_rate = 0 + if samplingrate_raw.shape == (2, 1): + # specified as onset/samplerate + warn("Onset sample rate SNIRF not yet supported") + else: + # specified as time points + fs_diff = np.around(np.diff(samplingrate_raw), decimals=4) + if len(np.unique(fs_diff)) == 1: + # Uniformly sampled data + sampling_rate = 1. / np.unique(fs_diff) + else: + # print(np.unique(fs_diff)) + warn("Non uniform sampled data not supported") + if sampling_rate == 0: + warn("Unable to extract sample rate") + + sources = np.array(dat.get('nirs/probe/sourceLabels')) + detectors = np.array(dat.get('nirs/probe/detectorLabels')) + sources = [s.decode('UTF-8') for s in sources] + detectors = [d.decode('UTF-8') for d in detectors] + + # Extract source and detector locations + detPos3D = np.array(dat.get('nirs/probe/detectorPos3D')) + srcPos3D = np.array(dat.get('nirs/probe/sourcePos3D')) + + assert len(sources) == srcPos3D.shape[0] + assert len(detectors) == detPos3D.shape[0] + + if detPos3D.shape[1] == 2: + detPos3D = np.hstack((detPos3D, np.zeros((len(detectors), 1)))) + srcPos3D = np.hstack((srcPos3D, np.zeros((len(sources), 1)))) + + # Extract wavelengths + fnirs_wavelengths = np.array(dat.get('nirs/probe/wavelengths')) + fnirs_wavelengths = [int(w) for w in fnirs_wavelengths] + + # Extract channels + def atoi(text): + return int(text) if text.isdigit() else text + + def natural_keys(text): + return [atoi(c) for c in re.split(r'(\d+)', text)] + + channels = np.array([name for name in dat['nirs']['data1'].keys()]) + channels_idx = np.array(['measurementList' in n for n in channels]) + channels = channels[channels_idx] + channels = sorted(channels, key=natural_keys) + + chnames = [] + for chan in channels: + src_idx = int(np.array(dat.get('nirs/data1/' + + chan + '/sourceIndex'))[0]) + det_idx = int(np.array(dat.get('nirs/data1/' + + chan + '/detectorIndex'))[0]) + wve_idx = int(np.array(dat.get('nirs/data1/' + + chan + '/wavelengthIndex'))[0]) + ch_name = sources[src_idx - 1] + '_' + detectors[det_idx - 1] + \ + ' ' + str(fnirs_wavelengths[wve_idx - 1]) + chnames.append(ch_name) + + # Create mne structure + info = create_info(chnames, + sampling_rate, + ch_types='fnirs_cw_amplitude') + + subject_info = {} + names = np.array(dat.get('nirs/metaDataTags/SubjectID')) + subject_info['first_name'] = names[0].decode('UTF-8') + info.update(subject_info=subject_info) + + LengthUnit = np.array(dat.get('/nirs/metaDataTags/LengthUnit')) + LengthUnit = LengthUnit[0].decode('UTF-8') + loc_scal = 1 + if "cm" in LengthUnit: + loc_scal = 100 + elif "mm" in LengthUnit: + loc_scal = 1000 + + for idx, chan in enumerate(channels): + src_idx = int(np.array(dat.get('nirs/data1/' + + chan + '/sourceIndex'))[0]) + det_idx = int(np.array(dat.get('nirs/data1/' + + chan + '/detectorIndex'))[0]) + wve_idx = int(np.array(dat.get('nirs/data1/' + + chan + '/wavelengthIndex'))[0]) + info['chs'][idx]['loc'][0:3] = srcPos3D[src_idx - 1, :] / loc_scal + info['chs'][idx]['loc'][3:6] = srcPos3D[src_idx - 1, :] / loc_scal + info['chs'][idx]['loc'][6:9] = detPos3D[det_idx - 1, :] / loc_scal + info['chs'][idx]['loc'][9] = fnirs_wavelengths[wve_idx - 1] + + super(RawSNIRF, self).__init__(info, preload, filenames=[fname], + last_samps=[last_samps], + verbose=verbose) + + # Extract annotations + annot = Annotations([], [], []) + for key in dat['nirs'].keys(): + if 'stim' in key: + data = np.array(dat.get('/nirs/' + key + '/data')) + if not data.shape == (): + annot.append(data[:, 0], 1.0, key[4:]) + self.set_annotations(annot) + + # Reorder channels to match expected ordering in MNE + num_chans = len(self.ch_names) + chans = [] + for idx in range(int(num_chans / 2)): + chans.append(idx) + chans.append(idx + int(num_chans / 2)) + self.pick(picks=chans) + + def _read_segment_file(self, data, idx, fi, start, stop, cals, mult): + """Read a segment of data from a file.""" + from ...externals.pymatreader.utils import _import_h5py + h5py = _import_h5py() + dat = h5py.File(self._filenames[0], 'r') + dat = np.array(dat.get('/nirs/data1/dataTimeSeries')).T + + data[:] = dat[idx, start:stop] + + return data diff --git a/mne/io/snirf/tests/__init__.py b/mne/io/snirf/tests/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/mne/io/snirf/tests/test_snirf.py b/mne/io/snirf/tests/test_snirf.py new file mode 100644 index 00000000000..b8c31e40544 --- /dev/null +++ b/mne/io/snirf/tests/test_snirf.py @@ -0,0 +1,80 @@ +# -*- coding: utf-8 -*- +# Authors: Robert Luke +# simplified BSD-3 license + +import os.path as op +from numpy.testing import assert_allclose + +from mne.datasets.testing import data_path, requires_testing_data +from mne.utils import run_tests_if_main, requires_h5py +from mne.io import read_raw_snirf, read_raw_nirx +from mne.io.tests.test_raw import _test_raw_reader + +fname_snirf_15_2_short = op.join(data_path(download=False), + 'SNIRF', + 'snirf_1_3_nirx_15_2_recording_w_short.snirf') + +fname_original = op.join(data_path(download=False), + 'NIRx', 'nirx_15_2_recording_w_short') + + +@requires_testing_data +@requires_h5py +def test_snirf_basic(): + """Test reading NIRX files.""" + raw = read_raw_snirf(fname_snirf_15_2_short, preload=True) + + # Test data import + assert raw._data.shape == (26, 145) + assert raw.info['sfreq'] == 12.5 + + # Test channel naming + assert raw.info['ch_names'][:4] == ["S1_D1 760", "S1_D1 850", + "S1_D9 760", "S1_D9 850"] + assert raw.info['ch_names'][24:26] == ["S5_D13 760", "S5_D13 850"] + + # Test frequency encoding + assert raw.info['chs'][0]['loc'][9] == 760 + assert raw.info['chs'][1]['loc'][9] == 850 + + assert_allclose([-8.6765 * 1e-2, 0.0049 * 1e-2, -2.6167 * 1e-2], + raw.info['chs'][1]['loc'][3:6], rtol=0.02) + + assert 'fnirs_cw_amplitude' in raw + + +@requires_testing_data +@requires_h5py +def test_snirf_against_nirx(): + """Test against file snirf was created from.""" + raw = read_raw_snirf(fname_snirf_15_2_short, preload=True) + raw_orig = read_raw_nirx(fname_original, preload=True) + + # Check annotations are the same + assert_allclose(raw.annotations.onset, raw_orig.annotations.onset) + assert_allclose([float(d) for d in raw.annotations.description], + [float(d) for d in raw_orig.annotations.description]) + assert_allclose(raw.annotations.duration, raw_orig.annotations.duration) + + # Check names are the same + assert raw.info['ch_names'] == raw_orig.info['ch_names'] + + # Check frequencies are the same + num_chans = len(raw.ch_names) + new_chs = raw.info['chs'] + ori_chs = raw_orig.info['chs'] + assert_allclose([new_chs[idx]['loc'][9] for idx in range(num_chans)], + [ori_chs[idx]['loc'][9] for idx in range(num_chans)]) + + assert_allclose(raw.get_data(), raw_orig.get_data()) + + +@requires_testing_data +@requires_h5py +def test_snirf_standard(): + """Test standard operations.""" + _test_raw_reader(read_raw_snirf, fname=fname_snirf_15_2_short, + boundary_decimal=0) # low fs + + +run_tests_if_main() diff --git a/tutorials/io/plot_30_reading_fnirs_data.py b/tutorials/io/plot_30_reading_fnirs_data.py index 0e9470f7a9e..3e80b06f144 100644 --- a/tutorials/io/plot_30_reading_fnirs_data.py +++ b/tutorials/io/plot_30_reading_fnirs_data.py @@ -23,9 +23,20 @@ The NIRx device stores data directly to a directory with multiple file types, MNE extracts the appropriate information from each file. -.. warning:: Information about device light wavelength is stored in - channel names. Manual modification of channel names is not - recommended. + +.. _import-snirf: + +SNIRF (.snirf) +================================ + +Data stored in the SNIRF format can be read in +using :func:`mne.io.read_raw_snirf`. +Details of the SNIRF format can be found on +The Society for functional near-infrared spectroscopy (SfNIRS) website. + +.. warning:: The SNIRF format has provisions for many different types of NIRS + recordings. MNE currently only supports continuous wave data + stored in the SNIRF format. Storing of optode locations @@ -35,4 +46,9 @@ A channel is formed by source-detector pairs. MNE stores the location of the channels, sources, and detectors. + +.. warning:: Information about device light wavelength is stored in + channel names. Manual modification of channel names is not + recommended. + """ # noqa:E501