Repository: mne-tools/mne-nirs (language: python)

Issue:
ENH: Add automated bad-channel marking
I recently migrated a pipeline from HOMER to MNE-Python. The one function I was missing was one that marked a channel in a recording if the optical density was over 3 standard deviations from the group mean (I think this was the criteria anyway, and it seems to match what they suggest offhand in [their publication](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2761652/#S11title)). Is it worth implementing something like this here? If so I have a bunch of recordings I can validate against. If not, I'm curious to know what you use instead @rob-luke !

It's not 100% clear to me, but I could see this being implemented as either 1) 3 std devs away from the mean *across all channels*, or 2) operating based on the std/mean of *each channel* independently -- probably the latter makes more sense. Then we'd want to use channel-specific annotations from MNE to mark segments as bad for segments of time during which it exceeded the threshold. Then users can choose to mark the entire channel bad for a given file/run (which the previous analysis I have access to did) or to have those segments marked bad in time. However, this latter method will be a bit tricky because in principle I think we'd want to exclude that interval from the GLM somehow -- which seems hard to do for the entire recording, let alone just for a single channel.

Mostly thinking out loud here @rob-luke, perhaps it makes sense to have a quick chat about this at some point!


Discussion / hints:
I've always found that the scalp coupling index identifies bad channels quite well, and so I didn't need a std based metric. Have you found cases where SCI doesn't mark the channels you want?

https://mne.tools/stable/generated/mne.preprocessing.nirs.scalp_coupling_index.html?highlight=coupling#mne.preprocessing.nirs.scalp_coupling_index

I have also been wanting to implement the peak spectral power metric for identifying bad segments of the recording (whereas sci finds bad channels for the entire recording). As SCI can incorrectly mark a channel as good if it has lots of artifacts. I think this might be a better method than the moving std approach

See https://github.com/lpollonini/qt-nirs and 

Hernandez, Samuel Montero, and Luca Pollonini. "NIRSplot: a tool for quality assessment of fNIRS scans." Optics and the Brain. Optical Society of America, 2020.
This is not to say we shouldn't add the std metric. I'm happy with that, more options the better. I just haven't had need for it. And would suggest the PSP method is a more considered approach.
After a quick chat @rob-luke and I settled on using peak spectral power (PSP) to do channel-wise annotations of bad segments, then use the scalp coupling index (SCI) computed using only good segments only. Sounds like @rob-luke has some existing MATLAB code that can be adapted with permission under our BSD 3-clause license

TL;DR: Not adding std-based bad channel marking, instead improving SCI using PSP
A few notes for myself...

The output should be able to generate a nice figure like: https://mne.tools/mne-bids-pipeline/examples/ds003392/sub-01_task-localizer_report.html#11

And the naming of arguments should be similar to https://github.com/mne-tools/mne-python/blob/f1e6097c36e7232514d67d4b04af568e36ac5758/mne/preprocessing/maxwell.py#L1927

Fix patch:
diff --git a/doc/api.rst b/doc/api.rst
index 1357995af..2983bcd15 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -24,6 +24,21 @@ Experimental Design and Analysis
    create_boxcar
 
 
+Preprocessing
+-------------
+
+.. currentmodule:: mne_nirs.preprocessing
+
+.. automodule:: mne_nirs.preprocessing
+   :no-members:
+   :no-inherited-members:
+
+.. autosummary::
+   :toctree: generated/
+
+   peak_power
+
+
 Signal Enhancement
 ------------------
 
@@ -74,6 +89,7 @@ Visualisation
    plot_glm_topo
    plot_glm_contrast_topo
    plot_glm_group_topo
+   plot_timechannel_quality_metric
 
 
 Simulation
diff --git a/environment.yml b/environment.yml
index 6721c0330..7c1282601 100644
--- a/environment.yml
+++ b/environment.yml
@@ -39,4 +39,4 @@ dependencies:
 - lxml
 - pip:
   - ipyvtklink
-  - https://github.com/nilearn/nilearn/archive/master.zip
+  - https://github.com/nilearn/nilearn/archive/main.zip
diff --git a/examples/general/plot_22_quality.py b/examples/general/plot_22_quality.py
new file mode 100644
index 000000000..07509dc99
--- /dev/null
+++ b/examples/general/plot_22_quality.py
@@ -0,0 +1,213 @@
+"""
+.. _ex-fnirs-quality:
+
+=========================
+Signal Quality Evaluation
+=========================
+
+This tutorial demonstrates how signal quality can be evaluated using MNE-NIRS.
+
+.. sidebar:: Relevant literature
+
+   Pollonini L et al., “PHOEBE: a method for real time mapping of
+   optodes-scalp coupling in functional near-infrared spectroscopy” in
+   Biomed. Opt. Express 7, 5104-5119 (2016).
+
+   Hernandez, Samuel Montero, and Luca Pollonini. "NIRSplot: a tool for
+   quality assessment of fNIRS scans." Optics and the Brain.
+   Optical Society of America, 2020.
+
+Ensuring your data is of high quality is essential to good scientific research.
+Evaluating data quality is an essential part of both data collection and
+analysis.
+For assessing data quality during data acquisition tools such as PHOEBE are
+highly recommended, these allow the researcher to address potential problems
+immediately and modify their experiment to improve the quality of recorded
+data.
+It is also important to assess data quality during data analysis.
+MNE-Python and MNE-NIRS provides several mechanisms to allow researchers
+to evaluate the quality of their data and to include this information in their
+downstream processing.
+Dedicated tools exist for quality evaluatiuon such as Hernandez (2020).
+This tutorial demonstrates methods in MNE-NIRS and MNE-Python for determining
+channels with poor signal
+quality, and methods for determining time segments of data that are of low
+quality in a subset of channels.
+
+Two methods are introduced in this tutorial. The scalp coupling index (SCI)
+and peak power (PP) metrics. Both of these methods examine the signal for
+the presence of a heart beat signal, which indicates that the sensors
+were in contact with the scalp. For further details see the papers listed
+in the relevant literature sidebar.
+
+"""
+# sphinx_gallery_thumbnail_number = 4
+
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import os
+import mne
+import numpy as np
+from itertools import compress
+import matplotlib.pyplot as plt
+
+from mne.preprocessing.nirs import optical_density
+from mne_nirs.preprocessing import peak_power
+from mne_nirs.visualisation import plot_timechannel_quality_metric
+
+###############################################################################
+# Import data
+# -----------
+#
+# Here we will work with the :ref:`fNIRS motor data <fnirs-motor-dataset>`.
+# We resample the data to make indexing exact times more convenient.
+# We then convert the data to optical density and plot the raw signal.
+
+fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
+fnirs_cw_amplitude_dir = os.path.join(fnirs_data_folder, 'Participant-1')
+raw_intensity = mne.io.read_raw_nirx(fnirs_cw_amplitude_dir, verbose=True)
+raw_intensity.load_data().resample(4.0, npad="auto")
+raw_od = optical_density(raw_intensity)
+raw_od.plot(n_channels=55, duration=4000, show_scrollbars=False, clipping=None)
+
+###############################################################################
+# From the above plot we observe that the data is relatively clean.
+# Later we will add some artificial bad sections to demonstrate the quality
+# evaluation metrics.
+
+
+###############################################################################
+# Scalp Coupling Index
+# --------------------
+#
+# The scalp coupling index (SCI) from Pollonini (2016) provides a measure of
+# the quality of the signal for a channel as determined using the entire
+# measurement duration. See Pollonini (2016) for further details of the
+# theory and implementation. Here we calculate the SCI for each channel
+# and view the distribution of values.
+
+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])
+
+###############################################################################
+# We observe that most of the channels have a good SCI of 1, but a few channels
+# have a poorer score. We can list the channels with an SCI below a threshold.
+# And we can mark these as bad in the MNE data. This way the functions
+# downstream will know that the data in these channels is bad.
+# For this example we set the threshold for a bad channel to be SCI < 0.8.
+# We then print a list of the bad channels and observe their are 10 channels
+# (five source-detector pairs) that are marked as bad.
+
+raw_od.info['bads'] = list(compress(raw_od.ch_names, sci < 0.7))
+print(raw_od.info['bads'])
+
+###############################################################################
+# We can plot the time course of the signal again and note that the bad
+# channels are now displayed in grey to indicate they are bad.
+
+raw_od.plot(n_channels=55, duration=4000, show_scrollbars=False, clipping=None)
+
+###############################################################################
+# Peak Power
+# ----------
+#
+# It may also be informative to view the quality of the signal at a finer
+# time resolution. The Peak Power metric provides a quality metric evalauted
+# over a 10 second window. This allows the user to view instances where
+# a subset of channels may be contaminated by artifacts for a short duration
+# of the recording.
+# This plot is based on the meg bad channel detection figures
+# available in mne-bids-pipeline.
+
+raw_od, scores, times = peak_power(raw_od)
+plot_timechannel_quality_metric(raw_od, scores, times, threshold=0.1)
+
+
+###############################################################################
+# Introduced Noise
+# ----------------
+#
+# The above data is quite clean, so here we add some noise to channels to
+# demonstrate that the algorithm is able to detect bad segments of data.
+# We add a step like signal at two time instances to mimic an artifact.
+# We add one artifact to channel S2-D4 at 1000 seconds
+# and one artifact to channel S5-D13 at 2000 seconds.
+
+# Add an artifact to channel S2-D4 at time 1000 seconds
+raw_od._data[12, 4000:4080] = np.linspace(0, -0.5, 80) + raw_od._data[12, 4000]
+
+# Add an artifact to channel S5-D13 at time 2000 seconds
+raw_od._data[34, 8000:8080] = np.linspace(0, 0.5, 80) + raw_od._data[34, 8000]
+
+###############################################################################
+# Next we plot just these channels to demonstrate that indeed an artifact
+# has been added.
+
+raw_od.copy().pick(picks = [12, 13, 34, 35]).\
+    plot(n_channels=55,duration=40000, show_scrollbars=False,
+         clipping=None, scalings={'fnirs_od': 0.2})
+
+
+###############################################################################
+# Peak Power Metric
+# -----------------
+#
+# To determine specific time and channel instances where data is of low
+# quality the peak power metric can be utilise (Pollonini, 2016).
+# Below the peak power is calculated and displayed for all channels.
+# The default threshold of 0.1 is applied to the quality scores to determine
+# bad segments of data.
+# The left facet illustrates the raw quality scores.
+# The right facet illustrates the channel and time segments that do not meet
+# the threshold criteria and are marked as bad.
+
+raw_od, scores, times = peak_power(raw_od)
+plot_timechannel_quality_metric(raw_od, scores, times, threshold=0.1)
+
+
+###############################################################################
+# Annotations
+# -----------
+#
+# Similar to how entire channels were marked as bad above, the peak power
+# function annotates the raw data structure to indicate where the bad
+# segments of data are. If we replot the subset of channels as above we note
+# that the bad segments are marked in red now, indicating this time section
+# contained bad data.
+# Note that although the red line appears as if the bad marking is present
+# for all channels, this is due to the plotting code. Internally the software
+# knows that only certain channels had a bad segment, and downstream processing
+# will only treat the specified bad channels as bad. I.e. when generating
+# epochs for time 1000 s, no epochs would be generated for channel S2-D4, but
+# channels would be generated for S5-D13 (as the artifact was only present
+# on S2-D4).
+
+raw_od.copy().pick(picks = [12, 13, 34, 35]).\
+    plot(n_channels=55,duration=40000, show_scrollbars=False,
+         clipping=None, scalings={'fnirs_od': 0.2})
+
+###############################################################################
+# These channel and time specific annotations are used by downstream
+# processing. For example, when extracting epochs if a specific channel
+# has an annotation indicating a bad segment in a specific time window, then
+# the epochs will be discarded for that specific channel.
+# Finally, we display the entire time signal after both channel and time
+# specific quality evaluation. Note that bad channels are now greyed out
+# and bad time segments are marked in red.
+
+raw_od.plot(n_channels=55, duration=4000, show_scrollbars=False, clipping=None)
+
+###############################################################################
+# Conclusion
+# ----------
+#
+# Two data quality metrics were presented and plotted.
+# One metric for determining a bad channel (scalp coupling index).
+# And one metric for determining specific time and channel segments that are
+# bad.
+# It was demonstrated how bad segments are visualised and how bad segments
+# are treated in downstream analysis.
diff --git a/mne_nirs/preprocessing/__init__.py b/mne_nirs/preprocessing/__init__.py
new file mode 100644
index 000000000..7c7bab24c
--- /dev/null
+++ b/mne_nirs/preprocessing/__init__.py
@@ -0,0 +1,7 @@
+"""Preprocessing."""
+
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+from ._peak_power import peak_power
diff --git a/mne_nirs/preprocessing/_peak_power.py b/mne_nirs/preprocessing/_peak_power.py
new file mode 100644
index 000000000..307dfc337
--- /dev/null
+++ b/mne_nirs/preprocessing/_peak_power.py
@@ -0,0 +1,106 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import numpy as np
+
+from scipy.signal import periodogram
+
+from mne import pick_types
+from mne.io import BaseRaw
+from mne.utils import _validate_type, verbose
+from mne.preprocessing.nirs import (_channel_frequencies,
+                                    _check_channels_ordered)
+from mne.filter import filter_data
+
+
+@verbose
+def peak_power(raw, time_window=10, threshold=0.1, l_freq=0.7, h_freq=1.5,
+               l_trans_bandwidth=0.3, h_trans_bandwidth=0.3,
+               verbose=False):
+    """
+    Compute peak spectral power metric from [1]_ and [2]_.
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        The haemoglobin data.
+    time_window : number
+        The duration of the window over which to calculate the metric.
+        Default is 10 seconds as in PHOEBE paper.
+    threshold : number
+        Values below this are marked as bad and annotated in the raw file.
+    %(l_freq)s
+    %(h_freq)s
+    %(l_trans_bandwidth)s
+    %(h_trans_bandwidth)s
+    %(verbose)s
+
+    Returns
+    -------
+    scores : array (n_nirs, n_windows)
+        Array of peak power values
+    times : list
+        List of the start and end times of each window used to compute the
+        peak spectral power.
+
+    References
+    ----------
+    .. [1] Pollonini L et al., “PHOEBE: a method for real time mapping of
+           optodes-scalp coupling in functional near-infrared spectroscopy” in
+           Biomed. Opt. Express 7, 5104-5119 (2016).
+    .. [2] Hernandez, Samuel Montero, and Luca Pollonini. "NIRSplot: a tool for
+           quality assessment of fNIRS scans." Optics and the Brain.
+           Optical Society of America, 2020.
+    """
+
+    raw = raw.copy().load_data()
+    _validate_type(raw, BaseRaw, 'raw')
+
+    if not len(pick_types(raw.info, fnirs='fnirs_od')):
+        raise RuntimeError('Scalp coupling index '
+                           'should be run on optical density data.')
+
+    freqs = np.unique(_channel_frequencies(raw.info))
+    picks = _check_channels_ordered(raw.info, freqs)
+
+    filtered_data = filter_data(raw._data, raw.info['sfreq'], l_freq, h_freq,
+                                picks=picks, verbose=verbose,
+                                l_trans_bandwidth=l_trans_bandwidth,
+                                h_trans_bandwidth=h_trans_bandwidth)
+
+    window_samples = int(np.ceil(time_window * raw.info['sfreq']))
+    n_windows = int(np.floor(len(raw) / window_samples))
+
+    scores = np.zeros((len(picks), n_windows))
+    times = []
+
+    for window in range(n_windows):
+
+        start_sample = int(window * window_samples)
+        end_sample = start_sample + window_samples
+
+        t_start = raw.times[start_sample]
+        t_stop = raw.times[end_sample]
+        times.append((t_start, t_stop))
+
+        for ii in picks[::2]:
+
+            c1 = filtered_data[ii][start_sample:end_sample]
+            c2 = filtered_data[ii + 1][start_sample:end_sample]
+
+            c1 = c1 / np.std(c1)
+            c2 = c2 / np.std(c2)
+
+            c = np.correlate(c1, c2, "full")
+            c = c / (window_samples)
+            [f, pxx] = periodogram(c, fs=raw.info['sfreq'], window='hamming')
+
+            scores[ii, window] = max(pxx)
+            scores[ii + 1, window] = max(pxx)
+
+            if (threshold is not None) & (max(pxx) < threshold):
+                raw.annotations.append(t_start, time_window, 'BAD_PeakPower',
+                                       ch_names=[raw.ch_names[ii:ii + 2]])
+
+    return raw, scores, times
diff --git a/mne_nirs/visualisation/__init__.py b/mne_nirs/visualisation/__init__.py
index 0e0c51eb7..936805e57 100644
--- a/mne_nirs/visualisation/__init__.py
+++ b/mne_nirs/visualisation/__init__.py
@@ -7,3 +7,4 @@
 from ._plot_nirs_source_detector import plot_nirs_source_detector
 from._plot_GLM_topo import (plot_glm_topo, plot_glm_contrast_topo,
                             plot_glm_group_topo)
+from ._plot_quality_metrics import plot_timechannel_quality_metric
\ No newline at end of file
diff --git a/mne_nirs/visualisation/_plot_quality_metrics.py b/mne_nirs/visualisation/_plot_quality_metrics.py
new file mode 100644
index 000000000..0292c3094
--- /dev/null
+++ b/mne_nirs/visualisation/_plot_quality_metrics.py
@@ -0,0 +1,83 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import seaborn as sns
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def plot_timechannel_quality_metric(raw, scores, times, threshold=0.1):
+    """
+    Plot time x channel based quality metrics.
+
+    The left figure shows the raw score per channel and time.
+
+
+    Parameters
+    ----------
+    raw : instance of Raw
+        The haemoglobin data.
+    scores : array
+        The quality metric scores.
+    times : list of pairs
+        Start and end time for each quality metric.
+    %(l_freq)s
+    %(h_freq)s
+    %(l_trans_bandwidth)s
+    %(h_trans_bandwidth)s
+    %(verbose)s
+
+    Returns
+    -------
+    fig : figure
+        Matplotlib figure displaying raw scores and thresholded scores.
+    """
+
+    ch_names = raw.ch_names
+    cols = [np.round(t[0]) for t in times]
+
+    data_to_plot = pd.DataFrame(data=scores,
+                                columns=pd.Index(cols, name='Time (s)'),
+                                index=pd.Index(ch_names, name='Channel'))
+
+    n_chans = len(ch_names)
+    vsize = 0.2 * n_chans
+
+    # First, plot the "raw" scores.
+    fig, ax = plt.subplots(1, 2, figsize=(20, vsize))
+    fig.suptitle('Automated noisy channel detection: fNIRS',
+                 fontsize=16, fontweight='bold')
+    sns.heatmap(data=data_to_plot, cmap='Reds_r', vmin=0, vmax=1,
+                cbar_kws=dict(label='Score'), ax=ax[0])
+    [ax[0].axvline(x, ls='dashed', lw=0.25, dashes=(25, 15), color='gray')
+        for x in range(1, len(times))]
+    ax[0].set_title('All Scores', fontweight='bold')
+    markbad(raw, ax[0])
+
+    # Now, adjust the color range to highlight segments that exceeded the
+    # limit.
+
+    data_to_plot = pd.DataFrame(data=scores > threshold,
+                                columns=pd.Index(cols, name='Time (s)'),
+                                index=pd.Index(ch_names, name='Channel'))
+    sns.heatmap(data=data_to_plot, vmin=0, vmax=1,
+                cmap='Reds_r', cbar_kws=dict(label='Score'), ax=ax[1])
+    [ax[1].axvline(x, ls='dashed', lw=0.25, dashes=(25, 15), color='gray')
+        for x in range(1, len(times))]
+    ax[1].set_title('Scores < Limit', fontweight='bold')
+    markbad(raw, ax[1])
+
+    # The figure title should not overlap with the subplots.
+    fig.tight_layout(rect=[0, 0.03, 1, 0.95])
+
+    return fig
+
+
+def markbad(raw, ax):
+
+    [ax.axhline(y + 0.5, ls='solid', lw=2, color='black')
+     for y in np.where([ch in raw.info['bads'] for ch in raw.ch_names])[0]]
+
+    return ax
diff --git a/requirements.txt b/requirements.txt
index 91479cf33..3c19809ad 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -23,7 +23,7 @@ vtk<8.2; platform_system == "Darwin"
 vtk; platform_system != "Darwin"
 mayavi
 PySurfer[save_movie]
-https://github.com/nilearn/nilearn/archive/master.zip
+https://github.com/nilearn/nilearn/archive/main.zip
 xlrd
 imageio>=2.6.1
 imageio-ffmpeg>=0.4.1
diff --git a/requirements_doc.txt b/requirements_doc.txt
index d8b4e5b29..71cdb0c18 100644
--- a/requirements_doc.txt
+++ b/requirements_doc.txt
@@ -24,3 +24,4 @@ lets-plot
 memory_profiler
 Ipython
 sphinx-fontawesome
+seaborn


Test patch:
diff --git a/mne_nirs/preprocessing/tests/__init__.py b/mne_nirs/preprocessing/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/mne_nirs/preprocessing/tests/test_psp.py b/mne_nirs/preprocessing/tests/test_psp.py
new file mode 100644
index 000000000..27da62ee2
--- /dev/null
+++ b/mne_nirs/preprocessing/tests/test_psp.py
@@ -0,0 +1,21 @@
+# Authors: Robert Luke <mail@robertluke.net>
+#
+# License: BSD (3-clause)
+
+import os
+import pytest
+import mne
+
+from mne_nirs.preprocessing import peak_power
+
+
+@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
+def test_peak_power():
+
+    fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
+    fnirs_raw_dir = os.path.join(fnirs_data_folder, 'Participant-1')
+    raw = mne.io.read_raw_nirx(fnirs_raw_dir, verbose=True).load_data()
+    raw = mne.preprocessing.nirs.optical_density(raw)
+
+    raw, scores, times = peak_power(raw)
+    assert len(scores) == len(raw.ch_names)
diff --git a/mne_nirs/visualisation/tests/test_quality.py b/mne_nirs/visualisation/tests/test_quality.py
new file mode 100644
index 000000000..bfd9a3643
--- /dev/null
+++ b/mne_nirs/visualisation/tests/test_quality.py
@@ -0,0 +1,21 @@
+import os
+import pytest
+import mne
+
+from mne_nirs.preprocessing import peak_power
+from mne_nirs.visualisation import plot_timechannel_quality_metric
+
+
+@pytest.mark.filterwarnings('ignore:.*nilearn.glm module is experimental.*:')
+def test_peak_power():
+
+    fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
+    fnirs_raw_dir = os.path.join(fnirs_data_folder, 'Participant-1')
+    raw = mne.io.read_raw_nirx(fnirs_raw_dir, verbose=True).load_data()
+    raw = mne.preprocessing.nirs.optical_density(raw)
+
+    raw, scores, times = peak_power(raw)
+    assert len(scores) == len(raw.ch_names)
+
+    fig = plot_timechannel_quality_metric(raw, scores, times)
+    assert len(fig.axes) == 4
diff --git a/requirements_testing.txt b/requirements_testing.txt
index 039bd4db0..a495d3a4f 100644
--- a/requirements_testing.txt
+++ b/requirements_testing.txt
@@ -16,3 +16,4 @@ wheel
 # redundancy with docs as all doc notebooks are tested
 https://codeload.github.com/rob-luke/mne-bids/zip/nirs
 lets-plot
+seaborn
\ No newline at end of file