diff --git a/doc/python_reference.rst b/doc/python_reference.rst index e0a6566365e..eceb55a148f 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -359,6 +359,7 @@ Projections: read_ica run_ica corrmap + optical_density EEG referencing: diff --git a/mne/defaults.py b/mne/defaults.py index 9d362a988ac..8f3f265e480 100644 --- a/mne/defaults.py +++ b/mne/defaults.py @@ -11,20 +11,22 @@ 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', - fnirs_raw='k'), + fnirs_raw='k', fnirs_od='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', fnirs_raw='V'), + ecog='uV', hbo='uM', hbr='uM', ref_meg='fT', fnirs_raw='V', + fnirs_od='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, fnirs_raw=1), + hbo=1e6, hbr=1e6, ref_meg=1e15, fnirs_raw=1.0, fnirs_od=1.0), # 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., fnirs_raw=2e-2), + hbr=10e-6, whitened=10., fnirs_raw=2e-2, + fnirs_od=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.), @@ -35,7 +37,8 @@ ecg='ECG', emg='EMG', misc='misc', seeg='sEEG', bio='BIO', dipole='Dipole', ecog='ECoG', hbo='Oxyhemoglobin', ref_meg='Reference Magnetometers', fnirs_raw='fNIRS', - hbr='Deoxyhemoglobin', gof='Goodness of fit'), + fnirs_od='fNIRS', hbr='Deoxyhemoglobin', + gof='Goodness of fit'), mask_params=dict(marker='o', markerfacecolor='w', markeredgecolor='k', diff --git a/mne/preprocessing/__init__.py b/mne/preprocessing/__init__.py index 0c59396a91d..583202476b4 100644 --- a/mne/preprocessing/__init__.py +++ b/mne/preprocessing/__init__.py @@ -21,3 +21,4 @@ from .stim import fix_stim_artifact from .maxwell import maxwell_filter from .xdawn import Xdawn +from ._optical_density import optical_density diff --git a/mne/preprocessing/_optical_density.py b/mne/preprocessing/_optical_density.py new file mode 100644 index 00000000000..ea96835ed6b --- /dev/null +++ b/mne/preprocessing/_optical_density.py @@ -0,0 +1,47 @@ +# Authors: Robert Luke +# Eric Larson +# +# License: BSD (3-clause) + +from ..io import BaseRaw +from ..io.constants import FIFF +from ..utils import _validate_type, warn +from ..io.pick import _picks_to_idx + +import numpy as np + + +def optical_density(raw): + r"""Convert NIRS raw data to optical density. + + Parameters + ---------- + raw : instance of Raw + The raw data. + + Returns + ------- + raw : instance of Raw + The modified raw instance. + + """ + raw = raw.copy().load_data() + _validate_type(raw, BaseRaw, 'raw') + picks = _picks_to_idx(raw.info, 'fnirs_raw') + data_means = np.mean(raw.get_data(), axis=1) + + # The devices measure light intensity. Negative light intensities should + # not occur. If they do it is likely due to hardware or movement issues. + # Set all negative values to abs(x), this also has the benefit of ensuring + # that the means are all greater than zero for the division below. + if len(np.where(raw._data[picks] <= 0)[0]) > 0: + warn("Negative intensities encountered. Setting to abs(x)") + raw._data[picks] = np.abs(raw._data[picks]) + + for ii in picks: + raw._data[ii] /= data_means[ii] + np.log(raw._data[ii], out=raw._data[ii]) + raw._data[ii] *= -1 + raw.info['chs'][ii]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD + + return raw diff --git a/mne/preprocessing/tests/test_optical_density.py b/mne/preprocessing/tests/test_optical_density.py new file mode 100644 index 00000000000..1883039e1b1 --- /dev/null +++ b/mne/preprocessing/tests/test_optical_density.py @@ -0,0 +1,59 @@ +# Authors: Robert Luke +# Eric Larson +# +# License: BSD (3-clause) + +import os.path as op + +from mne.datasets.testing import data_path +from mne.io import read_raw_nirx, BaseRaw +from mne.preprocessing import optical_density +from mne.utils import _validate_type +from mne.datasets import testing + +import pytest as pytest +import numpy as np +from numpy.testing import assert_allclose + + +fname_nirx = op.join(data_path(download=False), + 'NIRx', 'nirx_15_2_recording_w_short') + + +@testing.requires_testing_data +def test_optical_density(): + """Test return type for optical density.""" + raw = read_raw_nirx(fname_nirx, preload=False) + assert 'fnirs_raw' in raw + assert 'fnirs_od' not in raw + raw = optical_density(raw) + _validate_type(raw, BaseRaw, 'raw') + assert 'fnirs_raw' not in raw + assert 'fnirs_od' in raw + + +@testing.requires_testing_data +def test_optical_density_zeromean(): + """Test that optical density can process zero mean data.""" + raw = read_raw_nirx(fname_nirx, preload=True) + raw._data[4] -= np.mean(raw._data[4]) + with pytest.warns(RuntimeWarning, match='Negative'): + raw = optical_density(raw) + assert 'fnirs_od' in raw + + +@testing.requires_testing_data +def test_optical_density_manual(): + """Test optical density on known values.""" + test_tol = 0.01 + raw = read_raw_nirx(fname_nirx, preload=True) + # log(1) = 0 + raw._data[4] = np.ones((145)) + # log(0.5)/-1 = 0.69 + # log(1.5)/-1 = -0.40 + test_data = np.tile([0.5, 1.5], 73)[:145] + raw._data[5] = test_data + + od = optical_density(raw) + assert_allclose(od.get_data([4]), 0.) + assert_allclose(od.get_data([5])[0, :2], [0.69, -0.4], atol=test_tol) diff --git a/mne/viz/raw.py b/mne/viz/raw.py index c8707189422..7c484179f52 100644 --- a/mne/viz/raw.py +++ b/mne/viz/raw.py @@ -332,7 +332,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', 'fnirs_raw']: + for t in ['hbo', 'hbr', 'fnirs_raw', 'fnirs_od']: inds += [pick_types(info, meg=False, ref_meg=False, fnirs=t, exclude=[])] types += [t] * len(inds[-1])