diff --git a/doc/_includes/data_formats.rst b/doc/_includes/data_formats.rst index 3263f00b2ad..ede81a317c8 100644 --- a/doc/_includes/data_formats.rst +++ b/doc/_includes/data_formats.rst @@ -67,6 +67,8 @@ EEG :ref:`eXimia ` .nxe :func:`mn EEG :ref:`General data format ` .gdf :func:`mne.io.read_raw_gdf` EEG :ref:`Nicolet ` .data :func:`mne.io.read_raw_nicolet` + +NIRS :ref:`NIRx ` directory :func:`mne.io.read_raw_nirx` ============ ============================================ ========= =================================== More details are provided in the tutorials in the :ref:`tut-data-formats` diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 98f7e43933d..a1a31e992ea 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -12,6 +12,9 @@ Current (0.20.dev0) ------------------- +- Add reader for NIRx data in :func:`mne.io.read_raw_nirx` by `Robert Luke`_ + + Changelog ~~~~~~~~~ diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 6fb6e4a7b84..45b1915e9b5 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -260,3 +260,6 @@ .. _Abram Hindle: http://softwareprocess.es +.. _Robert Luke: https://github.com/rob-luke + + diff --git a/doc/conf.py b/doc/conf.py index fd07a9ea87e..b10ec037c37 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -595,7 +595,7 @@ def reset_warnings(gallery_conf, fname): 'n_segments', # Undocumented (on purpose) 'RawKIT', 'RawEximia', 'RawEGI', 'RawEEGLAB', 'RawEDF', 'RawCTF', 'RawBTi', - 'RawBrainVision', 'RawCurry', + 'RawBrainVision', 'RawCurry', 'RawNIRX', # sklearn subclasses 'mapping', 'to', 'any', # unlinkable diff --git a/doc/python_reference.rst b/doc/python_reference.rst index 40b947a1ca8..cd75080eaac 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -64,6 +64,7 @@ Reading raw data read_raw_gdf read_raw_kit read_raw_nicolet + read_raw_nirx read_raw_eeglab read_raw_brainvision read_raw_egi diff --git a/mne/channels/channels.py b/mne/channels/channels.py index 668a25ac7d3..85728bc5f2f 100644 --- a/mne/channels/channels.py +++ b/mne/channels/channels.py @@ -182,6 +182,8 @@ def compensation_grade(self): 'syst': FIFF.FIFFV_SYST_CH, 'bio': FIFF.FIFFV_BIO_CH, 'ecog': FIFF.FIFFV_ECOG_CH, + 'fnirs_raw': FIFF.FIFFV_FNIRS_CH, + 'fnirs_od': FIFF.FIFFV_FNIRS_CH, 'hbo': FIFF.FIFFV_FNIRS_CH, 'hbr': FIFF.FIFFV_FNIRS_CH} _human2unit = {'ecg': FIFF.FIFF_UNIT_V, @@ -197,6 +199,8 @@ def compensation_grade(self): 'syst': FIFF.FIFF_UNIT_NONE, 'bio': FIFF.FIFF_UNIT_V, 'ecog': FIFF.FIFF_UNIT_V, + 'fnirs_raw': FIFF.FIFF_UNIT_V, + 'fnirs_od': FIFF.FIFF_UNIT_NONE, 'hbo': FIFF.FIFF_UNIT_MOL, 'hbr': FIFF.FIFF_UNIT_MOL} _unit2human = {FIFF.FIFF_UNIT_V: 'V', @@ -423,6 +427,10 @@ def set_channel_types(self, mapping): coil_type = FIFF.FIFFV_COIL_FNIRS_HBO elif ch_type == 'hbr': coil_type = FIFF.FIFFV_COIL_FNIRS_HBR + elif ch_type == 'fnirs_raw': + coil_type = FIFF.FIFFV_COIL_FNIRS_RAW + elif ch_type == 'fnirs_od': + coil_type = FIFF.FIFFV_COIL_FNIRS_OD else: coil_type = FIFF.FIFFV_COIL_NONE self.info['chs'][c_ind]['coil_type'] = coil_type diff --git a/mne/cov.py b/mne/cov.py index bebd825903a..9a08326aeda 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -1136,8 +1136,8 @@ class _RegCovariance(BaseEstimator): """Aux class.""" def __init__(self, info, grad=0.1, mag=0.1, eeg=0.1, seeg=0.1, ecog=0.1, - hbo=0.1, hbr=0.1, store_precision=False, - assume_centered=False): + hbo=0.1, hbr=0.1, fnirs_raw=0.1, fnirs_od=0.1, + store_precision=False, assume_centered=False): self.info = info # For sklearn compat, these cannot (easily?) be combined into # a single dictionary @@ -1148,6 +1148,8 @@ def __init__(self, info, grad=0.1, mag=0.1, eeg=0.1, seeg=0.1, ecog=0.1, self.ecog = ecog self.hbo = hbo self.hbr = hbr + self.fnirs_raw = fnirs_raw + self.fnirs_od = fnirs_od self.store_precision = store_precision self.assume_centered = assume_centered @@ -1430,6 +1432,7 @@ def _smart_eigh(C, info, rank, scalings=None, projs=None, @verbose def regularize(cov, info, mag=0.1, grad=0.1, eeg=0.1, exclude='bads', proj=True, seeg=0.1, ecog=0.1, hbo=0.1, hbr=0.1, + fnirs_raw=0.1, fnirs_od=0.1, rank=None, scalings=None, verbose=None): """Regularize noise covariance matrix. @@ -1470,6 +1473,10 @@ def regularize(cov, info, mag=0.1, grad=0.1, eeg=0.1, exclude='bads', Regularization factor for HBO signals. hbr : float (default 0.1) Regularization factor for HBR signals. + fnirs_raw : float (default 0.1) + Regularization factor for fNIRS raw signals. + fnirs_od : float (default 0.1) + Regularization factor for fNIRS optical density signals. %(rank_None)s .. versionadded:: 0.17 diff --git a/mne/defaults.py b/mne/defaults.py index 8088e484059..83c1c3b3955 100644 --- a/mne/defaults.py +++ b/mne/defaults.py @@ -10,20 +10,21 @@ color=dict(mag='darkblue', grad='b', eeg='k', eog='k', ecg='m', emg='k', ref_meg='steelblue', misc='k', stim='k', resp='k', chpi='k', exci='k', ias='k', syst='k', seeg='saddlebrown', dipole='k', - gof='k', bio='k', ecog='k', hbo='darkblue', hbr='b'), + gof='k', bio='k', ecog='k', hbo='darkblue', hbr='b', + fnirs_raw='k'), units=dict(mag='fT', grad='fT/cm', eeg='uV', eog='uV', ecg='uV', emg='uV', misc='AU', seeg='mV', dipole='nAm', gof='GOF', bio='uV', - ecog='uV', hbo='uM', hbr='uM', ref_meg='fT'), + ecog='uV', hbo='uM', hbr='uM', ref_meg='fT', fnirs_raw='V'), # scalings for the units scalings=dict(mag=1e15, grad=1e13, eeg=1e6, eog=1e6, emg=1e6, ecg=1e6, misc=1.0, seeg=1e3, dipole=1e9, gof=1.0, bio=1e6, ecog=1e6, - hbo=1e6, hbr=1e6, ref_meg=1e15), + hbo=1e6, hbr=1e6, ref_meg=1e15, fnirs_raw=1), # rough guess for a good plot scalings_plot_raw=dict(mag=1e-12, grad=4e-11, eeg=20e-6, eog=150e-6, ecg=5e-4, emg=1e-3, ref_meg=1e-12, misc='auto', stim=1, resp=1, chpi=1e-4, exci=1, ias=1, syst=1, seeg=1e-4, bio=1e-6, ecog=1e-4, hbo=10e-6, - hbr=10e-6, whitened=10.), + hbr=10e-6, whitened=10., fnirs_raw=2e-2), scalings_cov_rank=dict(mag=1e12, grad=1e11, eeg=1e5, # ~100x scalings seeg=1e1, ecog=1e4, hbo=1e4, hbr=1e4), ylim=dict(mag=(-600., 600.), grad=(-200., 200.), eeg=(-200., 200.), @@ -33,7 +34,7 @@ titles=dict(mag='Magnetometers', grad='Gradiometers', eeg='EEG', eog='EOG', ecg='ECG', emg='EMG', misc='misc', seeg='sEEG', bio='BIO', dipole='Dipole', ecog='ECoG', hbo='Oxyhemoglobin', - ref_meg='Reference Magnetometers', + ref_meg='Reference Magnetometers', fnirs_raw='fNIRS', hbr='Deoxyhemoglobin', gof='Goodness of fit'), mask_params=dict(marker='o', markerfacecolor='w', diff --git a/mne/evoked.py b/mne/evoked.py index ebb262f6ee3..195592d7693 100644 --- a/mne/evoked.py +++ b/mne/evoked.py @@ -568,7 +568,7 @@ def get_peak(self, ch_type=None, tmin=None, tmax=None, .. versionadded:: 0.16 """ # noqa: E501 supported = ('mag', 'grad', 'eeg', 'seeg', 'ecog', 'misc', 'hbo', - 'hbr', 'None') + 'hbr', 'None', 'fnirs_raw', 'fnirs_od') data_picks = _pick_data_channels(self.info, with_ref_meg=False) types_used = {channel_type(self.info, idx) for idx in data_picks} @@ -602,7 +602,7 @@ def get_peak(self, ch_type=None, tmin=None, tmax=None, seeg = True elif ch_type == 'ecog': ecog = True - elif ch_type in ('hbo', 'hbr'): + elif ch_type in ('hbo', 'hbr', 'fnirs_raw', 'fnirs_od'): fnirs = ch_type if ch_type is not None: diff --git a/mne/io/__init__.py b/mne/io/__init__.py index fc14e29d599..1c9bd0f1e68 100644 --- a/mne/io/__init__.py +++ b/mne/io/__init__.py @@ -26,6 +26,7 @@ from . import fiff from . import kit from . import nicolet +from . import nirx from . import eeglab from . import pick @@ -43,6 +44,7 @@ from .artemis123 import read_raw_artemis123 from .eeglab import read_raw_eeglab, read_epochs_eeglab from .eximia import read_raw_eximia +from .nirx import read_raw_nirx from .fieldtrip import (read_raw_fieldtrip, read_epochs_fieldtrip, read_evoked_fieldtrip) @@ -66,5 +68,5 @@ read_raw_curry, read_raw_edf, read_raw_eeglab, read_raw_egi, read_raw_eximia, read_raw_fieldtrip, read_raw_fif, read_raw_gdf, read_raw_kit, read_raw_nicolet, set_bipolar_reference, set_eeg_reference, - show_fiff, write_fiducials, write_info, + show_fiff, write_fiducials, write_info, read_raw_nirx, ] diff --git a/mne/io/meas_info.py b/mne/io/meas_info.py index b820213990f..0da56b16e6d 100644 --- a/mne/io/meas_info.py +++ b/mne/io/meas_info.py @@ -53,6 +53,10 @@ seeg=(FIFF.FIFFV_SEEG_CH, FIFF.FIFFV_COIL_EEG, FIFF.FIFF_UNIT_V), bio=(FIFF.FIFFV_BIO_CH, FIFF.FIFFV_COIL_NONE, FIFF.FIFF_UNIT_V), ecog=(FIFF.FIFFV_ECOG_CH, FIFF.FIFFV_COIL_EEG, FIFF.FIFF_UNIT_V), + fnirs_raw=(FIFF.FIFFV_FNIRS_CH, FIFF.FIFFV_COIL_FNIRS_RAW, + FIFF.FIFF_UNIT_V), + fnirs_od=(FIFF.FIFFV_FNIRS_CH, FIFF.FIFFV_COIL_FNIRS_OD, + FIFF.FIFF_UNIT_NONE), hbo=(FIFF.FIFFV_FNIRS_CH, FIFF.FIFFV_COIL_FNIRS_HBO, FIFF.FIFF_UNIT_MOL), hbr=(FIFF.FIFFV_FNIRS_CH, FIFF.FIFFV_COIL_FNIRS_HBR, FIFF.FIFF_UNIT_MOL) ) diff --git a/mne/io/nirx/__init__.py b/mne/io/nirx/__init__.py new file mode 100644 index 00000000000..87c42b9c12b --- /dev/null +++ b/mne/io/nirx/__init__.py @@ -0,0 +1,7 @@ +"""fNIRS module for conversion to FIF.""" + +# Author: Robert Luke +# +# License: BSD (3-clause) + +from .nirx import read_raw_nirx diff --git a/mne/io/nirx/nirx.py b/mne/io/nirx/nirx.py new file mode 100644 index 00000000000..11fd13300ed --- /dev/null +++ b/mne/io/nirx/nirx.py @@ -0,0 +1,283 @@ +# Authors: Robert Luke +# +# License: BSD (3-clause) + +from configparser import ConfigParser, RawConfigParser +import glob as glob +import re as re + +import numpy as np + +from ..base import BaseRaw +from ..meas_info import create_info +from ...utils import logger, verbose, fill_doc +from ..constants import FIFF +from ...annotations import Annotations + + +@fill_doc +def read_raw_nirx(fname, preload=False, verbose=None): + """Reader for a NIRX fNIRS recording. + + Parameters + ---------- + fname : str + Path to the NIRX data folder. + %(preload)s + %(verbose)s + + Returns + ------- + raw : instance of RawNIRX + A Raw object containing NIRX data. + + See Also + -------- + mne.io.Raw : Documentation of attribute and methods. + """ + return RawNIRX(fname, preload, verbose) + + +@fill_doc +class RawNIRX(BaseRaw): + """Raw object from a NIRX fNIRS file. + + Parameters + ---------- + fname : str + Path to the NIRX data folder. + %(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 import read_mat + logger.info('Loading %s' % fname) + + # Check if required files exist and store names for later use + files = dict() + keys = ('dat', 'evt', 'hdr', 'inf', 'set', 'tpl', 'wl1', 'wl2', + 'config.txt', 'probeInfo.mat') + for key in keys: + files[key] = glob.glob('%s/*%s' % (fname, key)) + if len(files[key]) != 1: + raise RuntimeError('Expect one %s file, got %d' % + (key, len(files[key]),)) + files[key] = files[key][0] + + # Read number of rows/samples of wavelength data + last_sample = -1 + for line in open(files['wl1']): + last_sample += 1 + + # Read participant information file + inf = ConfigParser(allow_no_value=True) + inf.read(files['inf']) + inf = inf._sections['Subject Demographics'] + + # Store subject information from inf file in mne format + # Note: NIRX also records "Study Type", "Experiment History", + # "Additional Notes", "Contact Information" and this information + # is currently discarded + subject_info = {} + names = inf['name'].split() + if len(names) > 0: + subject_info['first_name'] = \ + inf['name'].split()[0].replace("\"", "") + if len(names) > 1: + subject_info['last_name'] = \ + inf['name'].split()[-1].replace("\"", "") + if len(names) > 2: + subject_info['middle_name'] = \ + inf['name'].split()[-2].replace("\"", "") + # subject_info['birthday'] = inf['age'] # TODO: not formatted properly + subject_info['sex'] = inf['gender'].replace("\"", "") + # Recode values + if subject_info['sex'] in {'M', 'Male', '1'}: + subject_info['sex'] = FIFF.FIFFV_SUBJ_SEX_MALE + elif subject_info['sex'] in {'F', 'Female', '2'}: + subject_info['sex'] = FIFF.FIFFV_SUBJ_SEX_FEMALE + # NIRStar does not record an id, or handedness by default + + # Read header file + # The header file isn't compliant with the configparser. So all the + # text between comments must be removed before passing to parser + with open(files['hdr']) as f: + hdr_str = f.read() + hdr_str = re.sub('#.*?#', '', hdr_str, flags=re.DOTALL) + hdr = RawConfigParser() + hdr.read_string(hdr_str) + + # Check that the file format version is supported + if hdr['GeneralInfo']['NIRStar'] != "\"15.2\"": + raise RuntimeError('Only NIRStar version 15.2 is supported') + + # Parse required header fields + + # Extract frequencies of light used by machine + fnirs_wavelengths = [int(s) for s in + re.findall(r'(\d+)', + hdr['ImagingParameters']['Wavelengths'])] + + # Extract source-detectors + sources = np.asarray([int(s) for s in re.findall(r'(\d+)-\d+:\d+', + hdr['DataStructure']['S-D-Key'])], int) + detectors = np.asarray([int(s) for s in re.findall(r'\d+-(\d+):\d+', + hdr['DataStructure']['S-D-Key'])], int) + + # Determine if short channels are present and on which detectors + has_short = np.array(hdr['ImagingParameters']['ShortBundles'], int) + short_det = [int(s) for s in + re.findall(r'(\d+)', + hdr['ImagingParameters']['ShortDetIndex'])] + short_det = np.array(short_det, int) + + # Extract sampling rate + samplingrate = float(hdr['ImagingParameters']['SamplingRate']) + + # Read information about probe/montage/optodes + # A word on terminology used here: + # Sources produce light + # Detectors measure light + # Sources and detectors are both called optodes + # Each source - detector pair produces a channel + # Channels are defined as the midpoint between source and detector + mat_data = read_mat(files['probeInfo.mat'], uint16_codec=None) + requested_channels = mat_data['probeInfo']['probes']['index_c'] + src_locs = mat_data['probeInfo']['probes']['coords_s3'] * 10 + det_locs = mat_data['probeInfo']['probes']['coords_d3'] * 10 + ch_locs = mat_data['probeInfo']['probes']['coords_c3'] * 10 + + # Determine requested channel indices + # The wl1 and wl2 files include all possible source - detector pairs. + # But most of these are not relevant. We want to extract only the + # subset requested in the probe file + req_ind = np.array([], int) + for req_idx in range(requested_channels.shape[0]): + sd_idx = np.where((sources == requested_channels[req_idx][0]) & + (detectors == requested_channels[req_idx][1])) + req_ind = np.concatenate((req_ind, sd_idx[0])) + req_ind = req_ind.astype(int) + + # Generate meaningful channel names + def prepend(list, str): + str += '{0}' + list = [str.format(i) for i in list] + return(list) + snames = prepend(sources[req_ind], 'S') + dnames = prepend(detectors[req_ind], '-D') + sdnames = [m + str(n) for m, n in zip(snames, dnames)] + sd1 = [s + ' ' + str(fnirs_wavelengths[0]) for s in sdnames] + sd2 = [s + ' ' + str(fnirs_wavelengths[1]) for s in sdnames] + chnames = [val for pair in zip(sd1, sd2) for val in pair] + + # Create mne structure + info = create_info(chnames, + samplingrate, + ch_types='fnirs_raw') + info.update({'subject_info': subject_info}) + + # Store channel, source, and detector locations + # The channel location is stored in the first 3 entries of loc. + # The source location is stored in the second 3 entries of loc. + # The detector location is stored in the third 3 entries of loc. + # NIRx NIRSite uses MNI coordinates. + for ch_idx2 in range(requested_channels.shape[0]): + # Find source and store location + src = int(requested_channels[ch_idx2, 0]) - 1 + info['chs'][ch_idx2 * 2]['loc'][3:6] = src_locs[src, :] + info['chs'][ch_idx2 * 2 + 1]['loc'][3:6] = src_locs[src, :] + # Find detector and store location + det = int(requested_channels[ch_idx2, 1]) - 1 + info['chs'][ch_idx2 * 2]['loc'][6:9] = det_locs[det, :] + info['chs'][ch_idx2 * 2 + 1]['loc'][6:9] = det_locs[det, :] + # Store channel location + # Channel locations for short channels are bodged, + # for short channels use the source location and add small offset + if (has_short > 0) & (len(np.where(short_det == det + 1)[0]) > 0): + info['chs'][ch_idx2 * 2]['loc'][:3] = src_locs[src, :] + info['chs'][ch_idx2 * 2 + 1]['loc'][:3] = src_locs[src, :] + info['chs'][ch_idx2 * 2]['loc'][0] += 0.8 + info['chs'][ch_idx2 * 2 + 1]['loc'][0] += 0.8 + else: + info['chs'][ch_idx2 * 2]['loc'][:3] = ch_locs[ch_idx2, :] + info['chs'][ch_idx2 * 2 + 1]['loc'][:3] = ch_locs[ch_idx2, :] + raw_extras = {"sd_index": req_ind, 'files': files} + + super(RawNIRX, self).__init__( + info, preload, filenames=[fname], last_samps=[last_sample], + raw_extras=[raw_extras], verbose=verbose) + + # Read triggers from event file + t = [re.findall(r'(\d+)', line) for line in open(files['evt'])] + 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) + + 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'] + + wls = [ + _read_csv_rows_cols( + self._raw_extras[fi]['files'][key], + start, stop, sdindex, len(self.ch_names) // 2).T + for key in ('wl1', 'wl2') + ] + + # TODO: Make this more efficient by only indexing above what we need. + # For now let's just construct the full data matrix and index. + # Interleave wavelength 1 and 2 to match channel names: + this_data = np.zeros((len(wls[0]) * 2, stop - start)) + this_data[0::2, :] = wls[0] + this_data[1::2, :] = wls[1] + data[:] = this_data[idx] + + return data + + def _probe_distances(self): + """Return the distance between each source-detector pair.""" + dist = [np.linalg.norm(ch['loc'][3:6] - ch['loc'][6:9]) + for ch in self.info['chs']] + return np.array(dist, float) + + def _short_channels(self, threshold=10.0): + """Return a vector indicating which channels are short. + + Channels with distance less than `threshold` are reported as short. + """ + return self._probe_distances() < threshold + + +def _read_csv_rows_cols(fname, start, stop, cols, n_cols): + # The following is equivalent to: + # x = pandas.read_csv(fname, header=None, usecols=cols, skiprows=start, + # nrows=stop - start, delimiter=' ') + # But does not require Pandas, and is hopefully fast enough, as the + # reading should be done in C (CPython), as should the conversion to float + # (NumPy). + x = np.zeros((stop - start, n_cols)) + with open(fname, 'r') as fid: + for li, line in enumerate(fid): + if li >= start: + if li >= stop: + break + x[li - start] = np.array(line.split(), float)[cols] + return x diff --git a/mne/io/nirx/tests/__init__.py b/mne/io/nirx/tests/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/mne/io/nirx/tests/test_nirx.py b/mne/io/nirx/tests/test_nirx.py new file mode 100644 index 00000000000..c0ea4c4c298 --- /dev/null +++ b/mne/io/nirx/tests/test_nirx.py @@ -0,0 +1,89 @@ +# -*- coding: utf-8 -*- +# Authors: Robert Luke +# simplified BSD-3 license + +import os.path as op + +import numpy as np + +from mne.io import read_raw_nirx +from mne.io.tests.test_raw import _test_raw_reader +from mne.utils import run_tests_if_main +from mne.datasets.testing import data_path, requires_testing_data + +fname_nirx = op.join(data_path(download=False), + 'NIRx', 'nirx_15_2_recording_w_short') + + +@requires_testing_data +def test_nirx(): + """Test reading NIRX files.""" + raw = read_raw_nirx(fname_nirx, 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'][0] == "S1-D1 760" + assert raw.info['ch_names'][1] == "S1-D1 850" + assert raw.info['ch_names'][2] == "S1-D9 760" + assert raw.info['ch_names'][3] == "S1-D9 850" + assert raw.info['ch_names'][24] == "S5-D13 760" + assert raw.info['ch_names'][25] == "S5-D13 850" + + # Test info import + assert raw.info['subject_info']['sex'] == 1 + assert raw.info['subject_info']['first_name'] == "MNE" + assert raw.info['subject_info']['middle_name'] == "Test" + assert raw.info['subject_info']['last_name'] == "Recording" + + # Test distance between optodes matches values from + # nirsite https://github.com/mne-tools/mne-testing-data/pull/51 + # step 4 figure 2 + allowed_distance_error = 0.2 + distances = raw._probe_distances() + assert abs(distances[0] - 30.4) < allowed_distance_error + assert abs(distances[2] - 7.8) < allowed_distance_error + assert abs(distances[4] - 31.0) < allowed_distance_error + assert abs(distances[6] - 8.6) < allowed_distance_error + assert abs(distances[8] - 41.6) < allowed_distance_error + assert abs(distances[10] - 7.2) < allowed_distance_error + assert abs(distances[12] - 38.9) < allowed_distance_error + assert abs(distances[14] - 7.5) < allowed_distance_error + assert abs(distances[16] - 55.8) < allowed_distance_error + assert abs(distances[18] - 56.2) < allowed_distance_error + assert abs(distances[20] - 56.1) < allowed_distance_error + assert abs(distances[22] - 56.5) < allowed_distance_error + assert abs(distances[24] - 7.7) < allowed_distance_error + + # Test which channels are short + # These are the ones marked as red at + # https://github.com/mne-tools/mne-testing-data/pull/51 step 4 figure 2 + short_channels = raw._short_channels() + assert short_channels[0] is np.False_ + assert short_channels[2] is np.True_ + assert short_channels[4] is np.False_ + assert short_channels[6] is np.True_ + assert short_channels[8] is np.False_ + short_channels = raw._short_channels(threshold=3) + assert short_channels[0] is np.False_ + assert short_channels[2] is np.False_ + short_channels = raw._short_channels(threshold=50) + assert short_channels[0] is np.True_ + assert short_channels[2] is np.True_ + + # Test trigger events + assert raw.annotations[0]['description'] == '3.0' + assert raw.annotations[1]['description'] == '2.0' + assert raw.annotations[2]['description'] == '1.0' + + +@requires_testing_data +def test_nirx_standard(): + """Test standard operations.""" + _test_raw_reader(read_raw_nirx, fname=fname_nirx, + boundary_decimal=1) # low fs + + +run_tests_if_main() diff --git a/mne/io/pick.py b/mne/io/pick.py index 72e9b0bc39c..4df7e3c3a31 100644 --- a/mne/io/pick.py +++ b/mne/io/pick.py @@ -49,6 +49,10 @@ def get_channel_types(): dipole=dict(kind=FIFF.FIFFV_DIPOLE_WAVE), gof=dict(kind=FIFF.FIFFV_GOODNESS_FIT), ecog=dict(kind=FIFF.FIFFV_ECOG_CH), + fnirs_raw=dict(kind=FIFF.FIFFV_FNIRS_CH, + coil_type=FIFF.FIFFV_COIL_FNIRS_RAW), + fnirs_od=dict(kind=FIFF.FIFFV_FNIRS_CH, + coil_type=FIFF.FIFFV_COIL_FNIRS_OD), hbo=dict(kind=FIFF.FIFFV_FNIRS_CH, coil_type=FIFF.FIFFV_COIL_FNIRS_HBO), hbr=dict(kind=FIFF.FIFFV_FNIRS_CH, @@ -90,7 +94,10 @@ def get_channel_types(): 'meg': ('unit', {FIFF.FIFF_UNIT_T_M: 'grad', FIFF.FIFF_UNIT_T: 'mag'}), 'fnirs': ('coil_type', {FIFF.FIFFV_COIL_FNIRS_HBO: 'hbo', - FIFF.FIFFV_COIL_FNIRS_HBR: 'hbr'}), + FIFF.FIFFV_COIL_FNIRS_HBR: 'hbr', + FIFF.FIFFV_COIL_FNIRS_RAW: 'fnirs_raw', + FIFF.FIFFV_COIL_FNIRS_OD: 'fnirs_od', + }), } @@ -254,6 +261,10 @@ def _triage_fnirs_pick(ch, fnirs): return True elif ch['coil_type'] == FIFF.FIFFV_COIL_FNIRS_HBR and fnirs == 'hbr': return True + elif ch['coil_type'] == FIFF.FIFFV_COIL_FNIRS_RAW and fnirs == 'fnirs_raw': + return True + elif ch['coil_type'] == FIFF.FIFFV_COIL_FNIRS_OD and fnirs == 'fnirs_od': + return True return False @@ -382,14 +393,15 @@ def pick_types(info, meg=True, eeg=False, stim=False, eog=False, ecg=False, for key in ('grad', 'mag'): param_dict[key] = meg if isinstance(fnirs, bool): - for key in ('hbo', 'hbr'): + for key in ('hbo', 'hbr', 'fnirs_raw', 'fnirs_od'): param_dict[key] = fnirs for k in range(nchan): ch_type = channel_type(info, k) try: pick[k] = param_dict[ch_type] except KeyError: # not so simple - assert ch_type in ('grad', 'mag', 'hbo', 'hbr', 'ref_meg') + assert ch_type in ('grad', 'mag', 'hbo', 'hbr', 'ref_meg', + 'fnirs_raw', 'fnirs_od') if ch_type in ('grad', 'mag'): pick[k] = _triage_meg_pick(info['chs'][k], meg) elif ch_type == 'ref_meg': @@ -679,7 +691,8 @@ def channel_indices_by_type(info, picks=None): """ idx_by_type = {key: list() for key in _PICK_TYPES_KEYS if key not in ('meg', 'fnirs')} - idx_by_type.update(mag=list(), grad=list(), hbo=list(), hbr=list()) + idx_by_type.update(mag=list(), grad=list(), hbo=list(), hbr=list(), + fnirs_raw=list(), fnirs_od=list()) picks = _picks_to_idx(info, picks, none='all', exclude=(), allow_empty=True) for k in picks: @@ -746,7 +759,7 @@ def _contains_ch_type(info, ch_type): _validate_type(ch_type, 'str', "ch_type") meg_extras = ['mag', 'grad', 'planar1', 'planar2'] - fnirs_extras = ['hbo', 'hbr'] + fnirs_extras = ['hbo', 'hbr', 'fnirs_raw', 'fnirs_od'] valid_channel_types = sorted([key for key in _PICK_TYPES_KEYS if key != 'meg'] + meg_extras + fnirs_extras) _check_option('ch_type', ch_type, valid_channel_types) @@ -850,19 +863,21 @@ def _check_excludes_includes(chs, info=None, allow_bads=False): misc=False, resp=False, chpi=False, exci=False, ias=False, syst=False, seeg=True, dipole=False, gof=False, bio=False, ecog=True, fnirs=True) _PICK_TYPES_KEYS = tuple(list(_PICK_TYPES_DATA_DICT.keys()) + ['ref_meg']) -_DATA_CH_TYPES_SPLIT = ('mag', 'grad', 'eeg', 'seeg', 'ecog', 'hbo', 'hbr') +_DATA_CH_TYPES_SPLIT = ('mag', 'grad', 'eeg', 'seeg', 'ecog', 'hbo', 'hbr', + 'fnirs_raw', 'fnirs_od') _DATA_CH_TYPES_ORDER_DEFAULT = ('mag', 'grad', 'eeg', 'eog', 'ecg', 'emg', 'ref_meg', 'misc', 'stim', 'resp', 'chpi', 'exci', 'ias', 'syst', 'seeg', 'bio', - 'ecog', 'hbo', 'hbr', 'whitened') + 'ecog', 'hbo', 'hbr', 'fnirs_raw', 'fnirs_od', + 'whitened') # Valid data types, ordered for consistency, used in viz/evoked. _VALID_CHANNEL_TYPES = ('eeg', 'grad', 'mag', 'seeg', 'eog', 'ecg', 'emg', 'dipole', 'gof', 'bio', 'ecog', 'hbo', 'hbr', - 'misc') + 'fnirs_raw', 'fnirs_od', 'misc') _MEG_CH_TYPES_SPLIT = ('mag', 'grad', 'planar1', 'planar2') -_FNIRS_CH_TYPES_SPLIT = ('hbo', 'hbr') +_FNIRS_CH_TYPES_SPLIT = ('hbo', 'hbr', 'fnirs_raw', 'fnirs_od') def _pick_data_channels(info, exclude='bads', with_ref_meg=True): diff --git a/mne/io/tests/test_constants.py b/mne/io/tests/test_constants.py index bb239ce0c07..f561463fc29 100644 --- a/mne/io/tests/test_constants.py +++ b/mne/io/tests/test_constants.py @@ -26,8 +26,7 @@ 'has_key', 'iteritems', 'iterkeys', 'itervalues', # Py2 'viewitems', 'viewkeys', 'viewvalues', # Py2 ) -_tag_ignore_names = ( # for fiff-constants pending updates -) +_tag_ignore_names = () # for fiff-constants pending updates _ignore_incomplete_enums = ( # XXX eventually we could complete these 'bem_surf_id', 'cardinal_point_cardiac', 'cond_model', 'coord', 'dacq_system', 'diffusion_param', 'gantry_type', 'map_surf', diff --git a/mne/io/tests/test_pick.py b/mne/io/tests/test_pick.py index 367702a8af9..8ba2bdbeda9 100644 --- a/mne/io/tests/test_pick.py +++ b/mne/io/tests/test_pick.py @@ -259,11 +259,12 @@ def test_pick_bio(): def test_pick_fnirs(): """Test picking fNIRS channels.""" - names = 'A1 A2 Fz O hbo1 hbo2 hbr1'.split() - types = 'mag mag eeg eeg hbo hbo hbr'.split() + names = 'A1 A2 Fz O hbo1 hbo2 hbr1 fnirsRaw1 fnirsRaw2 fnirsOD1'.split() + types = 'mag mag eeg eeg hbo hbo hbr fnirs_raw fnirs_raw fnirs_od'.split() info = create_info(names, 1024., types) picks_by_type = [('mag', [0, 1]), ('eeg', [2, 3]), - ('hbo', [4, 5]), ('hbr', [6])] + ('hbo', [4, 5]), ('hbr', [6]), + ('fnirs_raw', [7, 8]), ('fnirs_od', [9])] assert_indexing(info, picks_by_type) @@ -492,6 +493,18 @@ def test_picks_to_idx(): assert_array_equal(np.arange(len(info['ch_names'])), _picks_to_idx(info, 'all')) assert_array_equal([0], _picks_to_idx(info, 'data')) + info = create_info(['a', 'b'], 1000., ['fnirs_raw', 'fnirs_od']) + assert_array_equal(np.arange(2), _picks_to_idx(info, 'fnirs')) + assert_array_equal([0], _picks_to_idx(info, 'fnirs_raw')) + assert_array_equal([1], _picks_to_idx(info, 'fnirs_od')) + info = create_info(['a', 'b'], 1000., ['fnirs_raw', 'misc']) + assert_array_equal(np.arange(len(info['ch_names'])), + _picks_to_idx(info, 'all')) + assert_array_equal([0], _picks_to_idx(info, 'data')) + info = create_info(['a', 'b'], 1000., ['fnirs_od', 'misc']) + assert_array_equal(np.arange(len(info['ch_names'])), + _picks_to_idx(info, 'all')) + assert_array_equal([0], _picks_to_idx(info, 'data')) run_tests_if_main() diff --git a/mne/io/tests/test_raw.py b/mne/io/tests/test_raw.py index f5d6504c109..a66b3cb66cf 100644 --- a/mne/io/tests/test_raw.py +++ b/mne/io/tests/test_raw.py @@ -41,7 +41,8 @@ def test_orig_units(): BaseRaw(info, last_samps=[1], orig_units=True) -def _test_raw_reader(reader, test_preloading=True, test_kwargs=True, **kwargs): +def _test_raw_reader(reader, test_preloading=True, test_kwargs=True, + boundary_decimal=2, **kwargs): """Test reading, writing and slicing of raw classes. Parameters @@ -51,6 +52,10 @@ def _test_raw_reader(reader, test_preloading=True, test_kwargs=True, **kwargs): test_preloading : bool Whether not preloading is implemented for the reader. If True, both cases and memory mapping to file are tested. + test_kwargs : dict + Test _init_kwargs support. + boundary_decimal : int + Number of decimals up to which the boundary should match. **kwargs : Arguments for the reader. Note: Do not use preload as kwarg. Use ``test_preloading`` instead. @@ -140,7 +145,7 @@ def _test_raw_reader(reader, test_preloading=True, test_kwargs=True, **kwargs): assert_array_almost_equal(concat_raw.annotations.onset[idx], expected_bad_boundary_onset, - decimal=2) + decimal=boundary_decimal) if raw.info['meas_id'] is not None: for key in ['secs', 'usecs', 'version']: diff --git a/mne/viz/raw.py b/mne/viz/raw.py index bd4d4644d1e..86852fdf5e2 100644 --- a/mne/viz/raw.py +++ b/mne/viz/raw.py @@ -331,7 +331,7 @@ def plot_raw(raw, events=None, duration=10.0, start=0.0, n_channels=20, for t in ['grad', 'mag']: inds += [pick_types(info, meg=t, ref_meg=False, exclude=[])] types += [t] * len(inds[-1]) - for t in ['hbo', 'hbr']: + for t in ['hbo', 'hbr', 'fnirs_raw']: inds += [pick_types(info, meg=False, ref_meg=False, fnirs=t, exclude=[])] types += [t] * len(inds[-1]) @@ -1056,7 +1056,6 @@ def _setup_browser_selection(raw, kind, selector=True): topo_ax = plt.subplot2grid((6, 1), (0, 0), rowspan=2, colspan=1) keys = np.concatenate([keys, ['Custom']]) order.update({'Custom': list()}) # custom selection with lasso - plot_sensors(raw.info, kind='select', ch_type='all', axes=topo_ax, ch_groups=kind, title='', show=False) fig_selection.radio = RadioButtons(rax, [key for key in keys diff --git a/mne/viz/utils.py b/mne/viz/utils.py index 9027910c709..71e66ec15b0 100644 --- a/mne/viz/utils.py +++ b/mne/viz/utils.py @@ -2950,16 +2950,17 @@ def center_cmap(cmap, vmin, vmax, name="cmap_centered"): def _set_psd_plot_params(info, proj, picks, ax, area_mode): """Set PSD plot params.""" import matplotlib.pyplot as plt - _data_types = ('mag', 'grad', 'eeg', 'seeg', 'ecog') + _data_types = ('mag', 'grad', 'eeg', 'seeg', 'ecog', 'fnirs_raw') _check_option('area_mode', area_mode, [None, 'std', 'range']) picks = _picks_to_idx(info, picks) # XXX this could be refactored more with e.g., plot_evoked # XXX when it's refactored, Report._render_raw will need to be updated - megs = ['mag', 'grad', False, False, False] - eegs = [False, False, True, False, False] - seegs = [False, False, False, True, False] - ecogs = [False, False, False, False, True] + megs = ['mag', 'grad', False, False, False, False] + eegs = [False, False, True, False, False, False] + seegs = [False, False, False, True, False, False] + ecogs = [False, False, False, False, True, False] + fnirss = [False, False, False, False, False, 'fnirs_raw'] titles = _handle_default('titles', None) units = _handle_default('units', None) scalings = _handle_default('scalings', None) @@ -2967,10 +2968,10 @@ def _set_psd_plot_params(info, proj, picks, ax, area_mode): titles_list = list() units_list = list() scalings_list = list() - for meg, eeg, seeg, ecog, name in zip(megs, eegs, seegs, ecogs, - _data_types): + for meg, eeg, seeg, ecog, fnirs, name in zip(megs, eegs, seegs, ecogs, + fnirss, _data_types): these_picks = pick_types(info, meg=meg, eeg=eeg, seeg=seeg, ecog=ecog, - ref_meg=False, exclude=[]) + ref_meg=False, fnirs=fnirs, exclude=[]) these_picks = np.intersect1d(these_picks, picks) if len(these_picks) > 0: picks_list.append(these_picks) @@ -3114,7 +3115,7 @@ def _plot_psd(inst, fig, freqs, psd_list, picks_list, titles_list, info['chs'] = [inst.info['chs'][p] for p in picks] valid_channel_types = ['mag', 'grad', 'eeg', 'seeg', 'eog', 'ecg', 'emg', 'dipole', 'gof', 'bio', 'ecog', 'hbo', - 'hbr', 'misc'] + 'hbr', 'misc', 'fnirs_raw', 'fnirs_od'] ch_types_used = list() for this_type in valid_channel_types: if this_type in types: diff --git a/tutorials/io/plot_30_reading_fnirs_data.py b/tutorials/io/plot_30_reading_fnirs_data.py new file mode 100644 index 00000000000..0e9470f7a9e --- /dev/null +++ b/tutorials/io/plot_30_reading_fnirs_data.py @@ -0,0 +1,38 @@ +# -*- coding: utf-8 -*- +r""" +.. _tut-importing-fnirs-data: + +================================= +Importing data from fNIRS devices +================================= + +MNE includes various functions and utilities for reading NIRS +data and optode locations. + +.. contents:: Page contents + :local: + :depth: 2 + + +.. _import-nirx: + +NIRx (directory) +================================ + +NIRx recordings can be read in using :func:`mne.io.read_raw_nirx`. +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. + + +Storing of optode locations +=========================== + +NIRs devices consist of light sources and light detectors. +A channel is formed by source-detector pairs. +MNE stores the location of the channels, sources, and detectors. + +""" # noqa:E501