BUG: Fix bug with optical_density

doc/changes/latest.inc

@@ -97,6 +97,8 @@ Bugs
 
 - Fix bug with `mne.io.read_raw_fieldtrip` and `mne.read_epochs_fieldtrip` where channel positions were not set properly (:gh:`9447` by `Eric Larson`_)
 
+- Fix bug with `mne.preprocessing.nirs.optical_density` where protection against zero values was not guaranteed (:gh:`9522` by `Eric Larson`_)
+
 - :func:`mne.concatenate_raws` now raises an exception if ``raw.info['dev_head_t']`` differs between files. This behavior can be controlled using the new ``on_mismatch`` parameter (:gh:`9438` by `Richard Höchenberger`_)
 
 - Fixed bug in :meth:`mne.Epochs.drop_bad` where subsequent rejections failed if they only specified thresholds for a subset of the channel types used in a previous rejection (:gh:`9485` by `Richard Höchenberger`_).
@@ -111,4 +113,4 @@ Bugs
 
 API changes
 ~~~~~~~~~~~
-- In `mne.compute_source_morph`, the ``niter_affine`` and ``niter_sdr`` parameters have been replaced by ``niter`` and ``pipeline`` parameters for more consistent and finer-grained control of registration/warping steps and iteration (:gh:`9505` by `Alex Rockhill`_ and `Eric Larson`_)
+- In `mne.compute_source_morph`, the ``niter_affine`` and ``niter_sdr`` parameters have been replaced by ``niter`` and ``pipeline`` parameters for more consistent and finer-grained control of registration/warping steps and iteration (:gh:`9505` by `Alex Rockhill`_ and `Eric Larson`_)

mne/preprocessing/nirs/_optical_density.py

@@ -8,18 +8,20 @@
 
 from ...io import BaseRaw
 from ...io.constants import FIFF
-from ...utils import _validate_type, warn
+from ...utils import _validate_type, warn, verbose
 from ...io.pick import _picks_to_idx
 from ..nirs import _channel_frequencies, _check_channels_ordered
 
 
-def optical_density(raw):
+@verbose
+def optical_density(raw, *, verbose=None):
     r"""Convert NIRS raw data to optical density.
 
     Parameters
     ----------
     raw : instance of Raw
         The raw data.
+    %(verbose)s
 
     Returns
     -------
@@ -32,20 +34,26 @@ def optical_density(raw):
         raw.info, np.unique(_channel_frequencies(raw.info, nominal=True)))
 
     picks = _picks_to_idx(raw.info, 'fnirs_cw_amplitude')
-    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 np.any(raw._data[picks] <= 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
+        min_ = np.inf
+        for pi in picks:
+            np.abs(raw._data[pi], out=raw._data[pi])
+            min_ = min(min_, raw._data[pi].min() or min_)
+        # avoid == 0
+        for pi in picks:
+            np.maximum(raw._data[pi], min_, out=raw._data[pi])
+
+    data_means = np.mean(raw.get_data(picks=picks), axis=1)
+    for ii, pi in enumerate(picks):
+        raw._data[pi] /= data_means[ii]
+        np.log(raw._data[pi], out=raw._data[pi])
+        raw._data[pi] *= -1
+        raw.info['chs'][pi]['coil_type'] = FIFF.FIFFV_COIL_FNIRS_OD
 
     return raw

mne/preprocessing/nirs/_scalp_coupling_index.py

@@ -42,25 +42,28 @@ def scalp_coupling_index(raw, l_freq=0.7, h_freq=1.5,
     ----------
     .. footbibliography::
     """
-    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.')
+    if 'fnirs_od' not in raw:
+        raise RuntimeError(
+            'Scalp coupling index should be run on optical density data.')
 
     freqs = np.unique(_channel_frequencies(raw.info, nominal=True))
     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)
+    raw = raw.copy().pick(picks).load_data()
+    zero_mask = np.std(raw._data, axis=-1) == 0
+    filtered_data = raw.filter(
+        l_freq, h_freq, l_trans_bandwidth=l_trans_bandwidth,
+        h_trans_bandwidth=h_trans_bandwidth).get_data()
 
     sci = np.zeros(picks.shape)
-    for ii in picks[::2]:
-        c = np.corrcoef(filtered_data[ii], filtered_data[ii + 1])[0][1]
+    for ii in range(0, len(picks), 2):
+        with np.errstate(invalid='ignore'):
+            c = np.corrcoef(filtered_data[ii], filtered_data[ii + 1])[0][1]
+        if not np.isfinite(c):  # someone had std=0
+            c = 0
         sci[ii] = c
         sci[ii + 1] = c
-
+    sci[zero_mask] = 0
     return sci

mne/preprocessing/nirs/_tddr.py

@@ -132,6 +132,8 @@ def _TDDR(signal, sample_rate):
         sigma = 1.4826 * np.median(dev)
 
         # Step 3d. Scale deviations by standard deviation and tuning parameter
+        if sigma == 0:
+            break
         r = dev / (sigma * tune)
 
         # Step 3e. Calculate new weights according to Tukey's biweight function

mne/preprocessing/nirs/tests/test_optical_density.py

@@ -37,8 +37,10 @@ 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)
+    raw._data[4, -1] = 0
+    with np.errstate(invalid='raise', divide='raise'):
+        with pytest.warns(RuntimeWarning, match='Negative'):
+            raw = optical_density(raw)
     assert 'fnirs_od' in raw
 
 

mne/preprocessing/nirs/tests/test_scalp_coupling_index.py

@@ -32,7 +32,7 @@
 def test_scalp_coupling_index(fname, fmt, tmpdir):
     """Test converting NIRX files."""
     assert fmt in ('nirx', 'fif')
-    raw = read_raw_nirx(fname).load_data()
+    raw = read_raw_nirx(fname)
     with pytest.raises(RuntimeError, match='Scalp'):
         scalp_coupling_index(raw)
 
@@ -58,11 +58,17 @@ def test_scalp_coupling_index(fname, fmt, tmpdir):
     # Set next two channels to be uncorrelated
     raw._data[6] = new_data
     raw._data[7] = rng.rand(raw._data[0].shape[0])
+    # Set next channel to have zero std
+    raw._data[8] = 0.
+    raw._data[9] = 1.
+    raw._data[10] = 2.
+    raw._data[11] = 3.
     # 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
+    assert_allclose(sci[8:12], 0, atol=1e-10)
 
     # Ensure function errors if wrong type is passed in
     raw = beer_lambert_law(raw)

mne/preprocessing/nirs/tests/test_temporal_derivative_distribution_repair.py

@@ -6,6 +6,7 @@
 
 import pytest
 import numpy as np
+from numpy.testing import assert_allclose
 
 from mne.datasets.testing import data_path
 from mne.io import read_raw_nirx
@@ -28,7 +29,12 @@ def test_temporal_derivative_distribution_repair(fname, tmpdir):
     max_shift = np.max(np.diff(raw._data[0]))
     shift_amp = 5 * max_shift
     raw._data[0, 0:30] = raw._data[0, 0:30] - (shift_amp)
+    # make one channel zero std
+    raw._data[1] = 0.
+    raw._data[2] = 1.
     assert np.max(np.diff(raw._data[0])) > shift_amp
     # Ensure that applying the algorithm reduces the step change
     raw = tddr(raw)
     assert np.max(np.diff(raw._data[0])) < shift_amp
+    assert_allclose(raw._data[1], 0.)  # unchanged
+    assert_allclose(raw._data[2], 1.)  # unchanged