ENH: Convert NIRS raw data to optical density

doc/python_reference.rst

@@ -359,6 +359,7 @@ Projections:
    read_ica
    run_ica
    corrmap
+   optical_density
 
 EEG referencing:
 

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',

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

mne/preprocessing/_optical_density.py

@@ -0,0 +1,47 @@
+# Authors: Robert Luke <[email protected]>
+#          Eric Larson <[email protected]>
+#
+# 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

mne/preprocessing/tests/test_optical_density.py

@@ -0,0 +1,59 @@
+# Authors: Robert Luke <[email protected]>
+#          Eric Larson <[email protected]>
+#
+# 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)

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])