ENH: Allow passing source detector distances to beer lambert law

doc/changes/dev/13875.newfeature.rst

@@ -0,0 +1,6 @@
+- Fix :func:`mne.preprocessing.nirs.beer_lambert_law` to support explicitly
+  provided ``sd_distances`` during Beer-Lambert conversion. When valid
+  source-detector distances are already present in ``raw.info``, overriding them
+  now emits a warning, and when all inferred distances are invalid and no
+  explicit distances are provided, a ``ValueError`` is raised instead of
+  silently producing zero concentrations. (`#13862 <https://github.com/mne-tools/mne-python/pull/13862>`_ by `Kalle Makela`_)

mne/io/hitachi/tests/test_hitachi.py

@@ -226,7 +226,7 @@ def test_hitachi_basic(
     want = [np.nan] * (n_ch - 4)
     assert_allclose(distances, want, atol=0.0)
     raw_od_bad = optical_density(raw)
-    with pytest.warns(RuntimeWarning, match="will be zero"):
+    with pytest.raises(ValueError, match="all zero or NaN"):
         beer_lambert_law(raw_od_bad, ppf=6)
     # bad distances (too big)
     if versions[0] == "1.18" and len(fnames) == 1:

mne/preprocessing/nirs/__init__.py

@@ -5,6 +5,7 @@
 # Copyright the MNE-Python contributors.
 
 from .nirs import (
+    _has_source_detector_distances,
     short_channels,
     source_detector_distances,
     _check_channels_ordered,

mne/preprocessing/nirs/_beer_lambert_law.py

@@ -11,10 +11,15 @@
 from ..._fiff.constants import FIFF
 from ...io import BaseRaw
 from ...utils import _validate_type, pinv, warn
-from ..nirs import _channel_frequencies, _validate_nirs_info, source_detector_distances
+from ..nirs import (
+    _channel_frequencies,
+    _has_source_detector_distances,
+    _validate_nirs_info,
+    source_detector_distances,
+)
 
 
-def beer_lambert_law(raw, ppf=6.0):
+def beer_lambert_law(raw, ppf=6.0, *, sd_distances=None):
     r"""Convert NIRS optical density data to haemoglobin concentration.
 
     Parameters
@@ -26,6 +31,12 @@ def beer_lambert_law(raw, ppf=6.0):
 
         .. versionchanged:: 1.7
            Support for different factors for the two wavelengths.
+    sd_distances : array-like | float | None
+        Source-detector distances in meters. If ``None``, distances are read
+        from ``raw.info['chs']``. If array-like, the values must have a distance
+        for each channel, matching the order in ``info['chs']``.
+
+        .. versionadded:: 1.13
 
     Returns
     -------
@@ -70,9 +81,14 @@ def beer_lambert_law(raw, ppf=6.0):
         )
 
     abs_coef = _load_absorption(unique_freqs)  # shape (n_wavelengths, 2)
-    distances = source_detector_distances(raw.info, picks="all")
+    distances = _get_sd_distances(raw, sd_distances)
     bad = ~np.isfinite(distances[picks])
     bad |= distances[picks] <= 0
+    if bad.all():
+        raise ValueError(
+            "Source-detector distances are all zero or NaN. Consider setting a "
+            "montage with raw.set_montage or providing sd_distances."
+        )
     if bad.any():
         warn(
             "Source-detector distances are zero or NaN, some resulting "
@@ -129,6 +145,39 @@ def beer_lambert_law(raw, ppf=6.0):
     return raw
 
 
+def _get_sd_distances(raw, sd_distances):
+    """Get source-detector distances for each channel.
+
+    Returns
+    -------
+    dists : array of float
+        Array containing distances in meters.
+        Of shape equal to number of channels.
+    """
+    if sd_distances is None:
+        # picks="all" used here instead of picks s.t. distance indices match raw
+        return source_detector_distances(raw.info, picks="all")
+    elif _has_source_detector_distances(raw.info, picks="all"):
+        warn("Source-detector distances in raw.info[] will be overridden")
+    _validate_type(sd_distances, ("numeric", "array-like"), "sd_distances")
+    sd_distances = np.array(sd_distances, float)
+    n_channels = len(raw.info["chs"])
+    if sd_distances.ndim == 0:
+        return np.full(n_channels, sd_distances)
+    if sd_distances.ndim != 1:
+        raise ValueError(
+            "sd_distances must be a float or a 1D array-like, got "
+            f"shape {sd_distances.shape}"
+        )
+    if len(sd_distances) != n_channels:
+        raise ValueError(
+            "sd_distances must be a float or an array-like with length matching "
+            f"the len(raw.info['chs']) ({n_channels}), "
+            f"got length {len(sd_distances)}"
+        )
+    return sd_distances
+
+
 def _load_absorption(freqs):
     """Load molar extinction coefficients."""
     # Data from https://omlc.org/spectra/hemoglobin/summary.html

mne/preprocessing/nirs/nirs.py

@@ -62,6 +62,12 @@ def short_channels(info, threshold=0.01):
     return source_detector_distances(info) < threshold
 
 
+def _has_source_detector_distances(info, picks=None):
+    """Return True if source-detector distances can be computed."""
+    distances = source_detector_distances(info, picks=picks)
+    return len(distances) > 0 and np.all(distances > 0)
+
+
 def _channel_frequencies(info):
     """Return the light frequency for each channel."""
     # Only valid for fNIRS data before conversion to haemoglobin

mne/preprocessing/nirs/tests/test_beer_lambert_law.py

@@ -2,17 +2,22 @@
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
+
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose
 
+from mne import create_info
 from mne.datasets import testing
 from mne.datasets.testing import data_path
-from mne.io import BaseRaw, read_raw_fif, read_raw_nirx, read_raw_snirf
+from mne.io import BaseRaw, RawArray, read_raw_fif, read_raw_nirx, read_raw_snirf
 from mne.preprocessing.nirs import (
     _channel_frequencies,
     beer_lambert_law,
     optical_density,
+    source_detector_distances,
 )
+from mne.preprocessing.nirs._beer_lambert_law import _get_sd_distances
 from mne.utils import _validate_type
 
 testing_path = data_path(download=False)
@@ -112,3 +117,80 @@ def test_beer_lambert_v_matlab():
             + matlab_data["type"][idx]
         )
         assert raw.info["ch_names"][idx] == matlab_name
+
+
+def test_beer_lambert_sd_distances():
+    """Test Beer-Lambert conversion with explicit source-detector distances."""
+    data = np.array(
+        [[0.1, 0.2, 0.3], [0.15, 0.25, 0.35], [0.4, 0.5, 0.6], [0.45, 0.55, 0.65]]
+    )
+    # Ch names chosen to test reordered indices
+    ch_names = ["S1_D1 760", "S1_D1 850", "S10_D10 760", "S10_D10 850"]
+
+    # Case 1: valid locations, sd_distances=None
+    raw = RawArray(data, create_info(ch_names, sfreq=1.0, ch_types="fnirs_od"))
+    sd_distances = [0.03, 0.03, 0.03, 0.03]
+    for idx, (freq, distance) in enumerate(zip([760, 850, 760, 850], sd_distances)):
+        raw.info["chs"][idx]["loc"][3:6] = [0.0, 0.0, 0.0]
+        raw.info["chs"][idx]["loc"][6:9] = [distance, 0.0, 0.0]
+        raw.info["chs"][idx]["loc"][9] = freq
+    expected = beer_lambert_law(raw)
+
+    # Case 2: valid locations, sd_distances=<arr>
+    with pytest.warns(RuntimeWarning, match=r"(?i)will be overridden"):
+        actual = beer_lambert_law(raw, sd_distances=sd_distances)
+    assert actual.ch_names == expected.ch_names
+    assert_allclose(actual.get_data(), expected.get_data(), rtol=1e-12, atol=0)
+
+    # Case 3: no locations, sd_distances=None
+    for idx in range(len(raw.info["chs"])):
+        raw.info["chs"][idx]["loc"][3:9] = np.nan
+    assert np.isnan(source_detector_distances(raw.info)).all()
+    with pytest.raises(
+        ValueError, match=r"(?i)source-detector distances are all zero or NaN"
+    ):
+        beer_lambert_law(raw)
+
+    # Case 4: no locations, sd_distances=<arr>
+    actual = beer_lambert_law(raw, sd_distances=sd_distances)
+    assert actual.ch_names == expected.ch_names
+    assert_allclose(actual.get_data(), expected.get_data(), rtol=1e-12, atol=0)
+
+    # Case 5: no locations, sd_distances=<scalar>
+    actual = beer_lambert_law(raw, sd_distances=sd_distances[0])
+    assert actual.ch_names == expected.ch_names
+    assert_allclose(actual.get_data(), expected.get_data(), rtol=1e-12, atol=0)
+
+
+def test_get_sd_distances():
+    """Test source-detector distance selection and validation."""
+    raw = RawArray(
+        np.zeros((4, 3)),
+        create_info(
+            ["S1_D1 760", "S1_D1 850", "S2_D2 760", "S2_D2 850"], 1.0, "fnirs_od"
+        ),
+    )
+    expected = np.array([0.03, 0.03, 0.04, 0.04])
+    for idx, (freq, distance) in enumerate(zip([760, 850, 760, 850], expected)):
+        raw.info["chs"][idx]["loc"][3:6] = [0.0, 0.0, 0.0]
+        raw.info["chs"][idx]["loc"][6:9] = [distance, 0.0, 0.0]
+        raw.info["chs"][idx]["loc"][9] = freq
+
+    assert_allclose(_get_sd_distances(raw, None), expected, rtol=1e-12, atol=0)
+    with pytest.warns(RuntimeWarning, match=r"(?i)will be overridden"):
+        assert_allclose(_get_sd_distances(raw, expected), expected, rtol=1e-12, atol=0)
+    with pytest.warns(RuntimeWarning, match=r"(?i)will be overridden"):
+        assert_allclose(
+            _get_sd_distances(raw, 0.05), np.full(4, 0.05), rtol=1e-12, atol=0
+        )
+
+    for idx in range(len(raw.info["chs"])):
+        raw.info["chs"][idx]["loc"][3:9] = np.nan
+    assert_allclose(_get_sd_distances(raw, expected), expected, rtol=1e-12, atol=0)
+
+    with pytest.raises(ValueError, match=r"1D array-like"):
+        _get_sd_distances(raw, np.ones((2, 2)))
+    with pytest.raises(ValueError, match=r"length matching"):
+        _get_sd_distances(raw, [0.03, 0.03])
+    with pytest.raises(TypeError, match=r"sd_distances"):
+        _get_sd_distances(raw, "foo")

mne/preprocessing/nirs/tests/test_nirs.py

@@ -18,6 +18,7 @@
     _check_channels_ordered,
     _fnirs_optode_names,
     _fnirs_spread_bads,
+    _has_source_detector_distances,
     _optode_position,
     _validate_nirs_info,
     beer_lambert_law,
@@ -538,6 +539,22 @@ def test_optode_names():
     assert_array_equal(det_names, [f"D{n}" for n in ["1", "11", "17"]])
 
 
+def test_has_source_detector_distances():
+    """Ensure source-detector distance availability is detected."""
+    ch_names = ["S1_D1 760", "S1_D1 850", "S2_D1 760", "S2_D1 850"]
+    info = create_info(ch_names=ch_names, ch_types=np.repeat("fnirs_od", 4), sfreq=1.0)
+    assert not _has_source_detector_distances(info)
+    for idx in range(2):
+        info["chs"][idx]["loc"][3:6] = [0.0, 0.0, 0.0]
+        info["chs"][idx]["loc"][6:9] = [0.03, 0.0, 0.0]
+    assert _has_source_detector_distances(info, picks=[0, 1])
+    assert not _has_source_detector_distances(info)  # some chs missing
+    for idx in range(2, 4):
+        info["chs"][idx]["loc"][3:6] = [0.01, 0.0, 0.0]
+        info["chs"][idx]["loc"][6:9] = [0.04, 0.0, 0.0]
+    assert _has_source_detector_distances(info)
+
+
 @testing.requires_testing_data
 def test_optode_loc():
     """Ensure optode location extraction is correct."""