FIX: Use nominal wavelengths in BLL and SCI calculations

FIX: mne/io/hitachi tests with actual wavelength data now succeed

* Beer-Lambert Law (BLL) and Scalp Coupling Index (SCI) calculations
  used the info structure to determine the number of wavelengths, but
  that lead to errors as the info structure can contain arbitrary data
  (e.g. different wavelengths for each channel).
* BLL now uses the nominal wavelengths for the BLL calculation for all
  channels. Previously, it used the actual wavelengths of the first
  channel pair for all channels, which was incorrect if the channels
  have different actual wavelengths.
* In the future, there may be an option in BLL to use the actual
  freq values for each channel, to improve calculation accuracy.
* The Hitachi tests have data with actual wavelengths for each channel
  stored in the info structure. This caused issues with counting the
  number of wavelengths, which caused tests to fail.

Author: zEdS15B3GCwq

mne/preprocessing/nirs/_beer_lambert_law.py

@@ -11,7 +11,7 @@
 from ..._fiff.constants import FIFF
 from ...io import BaseRaw
 from ...utils import _validate_type, pinv, warn
-from ..nirs import _validate_nirs_info, source_detector_distances
+from ..nirs import _channel_frequencies, _validate_nirs_info, source_detector_distances
 
 
 def beer_lambert_law(raw, ppf=6.0):
@@ -37,8 +37,20 @@ def beer_lambert_law(raw, ppf=6.0):
     _validate_type(ppf, ("numeric", "array-like"), "ppf")
     ppf = np.array(ppf, float)
     picks = _validate_nirs_info(raw.info, fnirs="od", which="Beer-lambert")
-    # This is the one place we *really* need the actual/accurate frequencies
-    freqs = np.array([raw.info["chs"][pick]["loc"][9] for pick in picks], float)
+
+    # Use nominal channel frequencies
+    #
+    # Notes on implementation:
+    # 1. Frequencies are calculated the same way as in nirs._validate_nirs_info().
+    # 2. Wavelength values in the info structure may contain actual frequencies,
+    #    which may be used for more accurate calculation in the future.
+    # 3. nirs._channel_frequencies uses both cw_amplitude and OD data to determine
+    #    frequencies, whereas we only need those from OD here. Is there any chance
+    #    that they're different?
+    # 4. If actual frequencies were used, using np.unique() like below will lead to
+    #    errors. Instead, absorption coefficients will need to be calculated for
+    #    each individual frequency.
+    freqs = _channel_frequencies(raw.info)
 
     # Get unique wavelengths and determine number of wavelengths
     unique_freqs = np.unique(freqs)

mne/preprocessing/nirs/_scalp_coupling_index.py

@@ -6,7 +6,7 @@
 
 from ...io import BaseRaw
 from ...utils import _validate_type, verbose
-from ..nirs import _validate_nirs_info
+from ..nirs import _channel_frequencies, _validate_nirs_info
 
 
 @verbose
@@ -57,8 +57,9 @@ def scalp_coupling_index(
     ).get_data()
 
     # Determine number of wavelengths per source-detector pair
-    ch_wavelengths = [c["loc"][9] for c in raw.info["chs"]]
-    n_wavelengths = len(set(ch_wavelengths))
+    # We use nominal wavelengths as the info structure may contain arbitrary data.
+    freqs = _channel_frequencies(raw.info)
+    n_wavelengths = len(np.unique(freqs))
 
     # freqs = np.array([raw.info["chs"][pick]["loc"][9] for pick in picks], float)
     # n_wavelengths = len(set(unique_freqs))
@@ -67,6 +68,7 @@ def scalp_coupling_index(
 
     # Calculate all pairwise correlations within each group and use the minimum as SCI
     pair_indices = np.triu_indices(n_wavelengths, k=1)
+
     for gg in range(0, len(picks), n_wavelengths):
         group_data = filtered_data[gg : gg + n_wavelengths]
 

mne/preprocessing/nirs/nirs.py

@@ -155,6 +155,8 @@ def _check_channels_ordered(info, pair_vals, *, throw_errors=True, check_bads=Tr
         )
 
     # Ensure wavelength info exists for waveform data
+    # Note: currently, the only requirement for the wavelength field in info is
+    # that it cannot be NaN. It depends on the data readers what is stored in it.
     if len(picks_wave) > 0:
         all_freqs = [info["chs"][ii]["loc"][9] for ii in picks_wave]
         # test for nan values first as those mess up the output of set()
@@ -164,13 +166,6 @@ def _check_channels_ordered(info, pair_vals, *, throw_errors=True, check_bads=Tr
                 f'info["chs"] structure. The encoded wavelengths are {all_freqs}.',
                 throw_errors,
             )
-        # test if the info structure has the same number of freqs as pair_vals
-        if len(pair_vals) != len(set(all_freqs)):
-            picks = _throw_or_return_empty(
-                f"The {error_word} in info must match the number of wavelengths, "
-                f"but the data contains {len(set(all_freqs))} wavelengths instead.",
-                throw_errors,
-            )
 
     # Validate the channel naming scheme
     for pick in picks:

mne/preprocessing/nirs/tests/test_beer_lambert_law.py

@@ -73,7 +73,7 @@ def test_beer_lambert_unordered_errors():
         }
     )
     assert raw_od.ch_names[0] == "S1_D1 770" and raw_od.ch_names[1] == "S1_D1 840"
-    with pytest.raises(ValueError, match="must match the number of wavelengths"):
+    with pytest.raises(ValueError, match="NIRS channels not ordered correctly."):
         beer_lambert_law(raw_od)