WIP: Add scalp coupling index for NIRS data

doc/changes/latest.inc

@@ -97,6 +97,8 @@ Changelog
 
 - Add function to convert events to annotations :func:`mne.annotations_from_events` by `Nicolas Barascud`_
 
+- Add function to calculate scalp coupling index for fNIRS data :func:`mne.preprocessing.nirs.scalp_coupling_index` by `Robert Luke`_
+
 - Add ``item`` argument to :meth:`mne.Epochs.get_data` for faster access to NumPy data arrays compared to :meth:`mne.Epochs.__getitem__` for frequent access on large :class:`mne.Epochs` objects, by `Eric Larson`_
 
 - More accurate coordinate system for Easycap montages in :func:`mne.channels.make_standard_montage` by `Christian Brodbeck`_

doc/python_reference.rst

@@ -383,6 +383,7 @@ Projections:
    beer_lambert_law
    source_detector_distances
    short_channels
+   scalp_coupling_index
 
 EEG referencing:
 

doc/references.bib

@@ -1172,6 +1172,16 @@ @article{Pham2001
   year = {2001}
 }
 
+@article{pollonini2014auditory,
+  title={Auditory cortex activation to natural speech and simulated cochlear implant speech measured with functional near-infrared spectroscopy},
+  author={Pollonini, Luca and Olds, Cristen and Abaya, Homer and Bortfeld, Heather and Beauchamp, Michael S and Oghalai, John S},
+  journal={Hearing research},
+  volume={309},
+  pages={84--93},
+  year={2014},
+  publisher={Elsevier}
+}
+
 @article{RidgwayEtAl2012,
   author = {Ridgway, Gerard R. and Litvak, Vladimir and Flandin, Guillaume and Friston, Karl J. and Penny, Will D.},
   doi = {10.1016/j.neuroimage.2011.10.027},

mne/preprocessing/nirs/__init__.py

@@ -6,6 +6,7 @@
 #
 # License: BSD (3-clause)
 
-from .nirs import short_channels, source_detector_distances
+from .nirs import short_channels, source_detector_distances, _check_channels_ordered, _channel_frequencies
 from ._optical_density import optical_density
 from ._beer_lambert_law import beer_lambert_law
+from ._scalp_coupling_index import scalp_coupling_index

mne/preprocessing/nirs/_beer_lambert_law.py

@@ -6,15 +6,14 @@
 
 import os.path as op
 
-import re as re
 import numpy as np
 from scipy import linalg
 
 from ...io import BaseRaw
-from ...io.pick import _picks_to_idx
 from ...io.constants import FIFF
 from ...utils import _validate_type
-from ..nirs import source_detector_distances
+from ..nirs import source_detector_distances, _channel_frequencies,\
+    _check_channels_ordered
 
 
 def beer_lambert_law(raw, ppf=0.1):
@@ -59,35 +58,6 @@ def beer_lambert_law(raw, ppf=0.1):
     return raw
 
 
-def _channel_frequencies(raw):
-    """Return the light frequency for each channel."""
-    picks = _picks_to_idx(raw.info, 'fnirs_od')
-    freqs = np.empty(picks.size, int)
-    for ii in picks:
-        freqs[ii] = raw.info['chs'][ii]['loc'][9]
-    return freqs
-
-
-def _check_channels_ordered(raw, freqs):
-    """Check channels followed expected fNIRS format."""
-    # Every second channel should be same SD pair
-    # and have the specified light frequencies.
-    picks = _picks_to_idx(raw.info, 'fnirs_od')
-    for ii in picks[::2]:
-        ch1_name_info = re.match(r'S(\d+)_D(\d+) (\d+)',
-                                 raw.info['chs'][ii]['ch_name'])
-        ch2_name_info = re.match(r'S(\d+)_D(\d+) (\d+)',
-                                 raw.info['chs'][ii + 1]['ch_name'])
-
-        if (ch1_name_info.groups()[0] != ch2_name_info.groups()[0]) or \
-           (ch1_name_info.groups()[1] != ch2_name_info.groups()[1]) or \
-           (int(ch1_name_info.groups()[2]) != freqs[0]) or \
-           (int(ch2_name_info.groups()[2]) != freqs[1]):
-            raise RuntimeError('NIRS channels not ordered correctly')
-
-    return picks
-
-
 def _load_absorption(freqs):
     """Load molar extinction coefficients."""
     # Data from https://omlc.org/spectra/hemoglobin/summary.html

mne/preprocessing/nirs/_scalp_coupling_index.py

@@ -0,0 +1,65 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#          Eric Larson <larson.eric.d@gmail.com>
+#          Alexandre Gramfort <alexandre.gramfort@inria.fr>
+#
+# License: BSD (3-clause)
+
+import numpy as np
+
+from ...io import BaseRaw
+from ...utils import _validate_type, verbose
+from ..nirs import _channel_frequencies, _check_channels_ordered
+from ...filter import filter_data
+
+
+@verbose
+def scalp_coupling_index(raw, l_freq=0.7, h_freq=1.5,
+                         l_trans_bandwidth=0.3, h_trans_bandwidth=0.3,
+                         verbose=False):
+    r"""Calculate scalp coupling index.
+
+    This function calculates the scalp coupling index
+    :footcite:`pollonini2014auditory`. This is a measure of the quality of the
+    connection between the optode and the scalp.
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        The raw data.
+    %(l_freq)s
+    %(h_freq)s
+    %(l_trans_bandwidth)s
+    %(h_trans_bandwidth)s
+    %(verbose)s
+
+    Returns
+    -------
+    sci : array of float
+        Array containing scalp coupling index for each channel.
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    raw = raw.copy().load_data()
+    _validate_type(raw, BaseRaw, 'raw')
+
+    freqs = np.unique(_channel_frequencies(raw))
+    picks = _check_channels_ordered(raw, freqs)
+
+    if not len(picks):
+        raise RuntimeError('Scalp coupling index '
+                           'should be run on optical density data.')
+
+    filtered_data = filter_data(raw._data, raw.info['sfreq'], l_freq, h_freq,
+                                picks=picks, verbose=verbose,
+                                l_trans_bandwidth=h_trans_bandwidth,
+                                h_trans_bandwidth=l_trans_bandwidth)
+
+    sci = np.zeros(picks.shape)
+    for ii in picks[::2]:
+        c = np.corrcoef(filtered_data[ii], filtered_data[ii + 1])[0][1]
+        sci[ii] = c
+        sci[ii + 1] = c
+
+    return sci

mne/preprocessing/nirs/nirs.py

@@ -4,6 +4,7 @@
 #
 # License: BSD (3-clause)
 
+import re
 import numpy as np
 from scipy import linalg
 
@@ -53,3 +54,32 @@ def short_channels(info, threshold=0.01):
         Of shape equal to number of channels.
     """
     return source_detector_distances(info) < threshold
+
+
+def _channel_frequencies(raw):
+    """Return the light frequency for each channel."""
+    picks = _picks_to_idx(raw.info, 'fnirs_od')
+    freqs = np.empty(picks.size, int)
+    for ii in picks:
+        freqs[ii] = raw.info['chs'][ii]['loc'][9]
+    return freqs
+
+
+def _check_channels_ordered(raw, freqs):
+    """Check channels followed expected fNIRS format."""
+    # Every second channel should be same SD pair
+    # and have the specified light frequencies.
+    picks = _picks_to_idx(raw.info, 'fnirs_od')
+    for ii in picks[::2]:
+        ch1_name_info = re.match(r'S(\d+)_D(\d+) (\d+)',
+                                 raw.info['chs'][ii]['ch_name'])
+        ch2_name_info = re.match(r'S(\d+)_D(\d+) (\d+)',
+                                 raw.info['chs'][ii + 1]['ch_name'])
+
+        if (ch1_name_info.groups()[0] != ch2_name_info.groups()[0]) or \
+           (ch1_name_info.groups()[1] != ch2_name_info.groups()[1]) or \
+           (int(ch1_name_info.groups()[2]) != freqs[0]) or \
+           (int(ch2_name_info.groups()[2]) != freqs[1]):
+            raise RuntimeError('NIRS channels not ordered correctly')
+
+    return picks

mne/preprocessing/tests/test_scalp_coupling_index.py

@@ -0,0 +1,60 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#          Eric Larson <larson.eric.d@gmail.com>
+#          Alexandre Gramfort <alexandre.gramfort@inria.fr>
+#
+# License: BSD (3-clause)
+
+import os.path as op
+
+import pytest
+import numpy as np
+from numpy.testing import assert_allclose, assert_array_less
+
+from mne.datasets.testing import data_path
+from mne.io import read_raw_nirx
+from mne.preprocessing.nirs import optical_density, scalp_coupling_index
+from mne.datasets import testing
+
+fname_nirx_15_0 = op.join(data_path(download=False),
+                          'NIRx', 'nirx_15_0_recording')
+fname_nirx_15_2 = op.join(data_path(download=False),
+                          'NIRx', 'nirx_15_2_recording')
+fname_nirx_15_2_short = op.join(data_path(download=False),
+                                'NIRx', 'nirx_15_2_recording_w_short')
+
+
+@testing.requires_testing_data
+@pytest.mark.parametrize('fname', ([fname_nirx_15_2_short, fname_nirx_15_2,
+                                    fname_nirx_15_0]))
+@pytest.mark.parametrize('fmt', ('nirx', 'fif'))
+def test_scalp_coupling_index(fname, fmt, tmpdir):
+    """Test converting NIRX files."""
+    assert fmt in ('nirx', 'fif')
+    raw = read_raw_nirx(fname)
+    raw = optical_density(raw)
+    sci = scalp_coupling_index(raw)
+
+    # All values should be between -1 and +1
+    assert_array_less(sci, 1.0)
+    assert_array_less(sci * -1.0, 1.0)
+
+    # Fill in some data with known correlation values
+    new_data = np.random.rand(raw._data[0].shape[0])
+    # Set first two channels to perfect correlation
+    raw._data[0] = new_data
+    raw._data[1] = new_data
+    # Set next two channels to perfect correlation
+    raw._data[2] = new_data
+    raw._data[3] = new_data * 0.3  # check scale invariance
+    # Set next two channels to anti correlation
+    raw._data[4] = new_data
+    raw._data[5] = new_data * -1.0
+    # Set next two channels to be uncorrelated
+    # TODO: this might be a bad idea as sometimes random noise might correlate
+    raw._data[6] = new_data
+    raw._data[7] = np.random.rand(raw._data[0].shape[0])
+    # Check values
+    sci = scalp_coupling_index(raw)
+    assert_allclose(sci[0:6], [1, 1, 1, 1, -1, -1], atol=0.01)
+    assert np.abs(sci[6]) < 0.5
+    assert np.abs(sci[7]) < 0.5

tutorials/preprocessing/plot_70_fnirs_processing.py

@@ -19,6 +19,7 @@
 import os
 import numpy as np
 import matplotlib.pyplot as plt
+from itertools import compress
 
 import mne
 
@@ -53,6 +54,31 @@
 raw_od.plot(n_channels=len(raw_od.ch_names), duration=500)
 
 
+###############################################################################
+# Evaluating the quality of the data
+# ----------------------------------
+#
+# At this stage we can quantify the quality of the coupling
+# between the scalp and the optodes using the scalp coupling index. This
+# method looks at the presence of a prominent synchronous signal in the
+# frequency range of cardiac signals across both photodetected signals.
+#
+# As this data is clean and the coupling is good for all channels we will
+# not mark any channels as bad based on the scalp coupling index.
+
+sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)
+fig, ax = plt.subplots()
+ax.hist(sci)
+ax.set(xlabel='Scalp Coupling Index', ylabel='Count', xlim=[0, 1])
+
+
+###############################################################################
+# In this example we will mark all channels with a SCI less than 0.5 as bad
+# (this dataset is quite clean, so no channels are marked as bad).
+
+raw_od.info['bads'] = list(compress(raw_od.ch_names, sci < 0.5))
+
+
 ###############################################################################
 # Converting from optical density to haemoglobin
 # ----------------------------------------------