ENH: Support for Multi-Wavelength (>2) NIRS/SNIRF Data Processing

[zEdS15B3GCwq]

Describe the new feature or enhancement

Goal: extend MNE-Python's NIRS processing to support an arbitrary number of wavelengths (≥2) instead of being limited (hard-coded in several places) to two.

The same idea has been raised in issue #9816 with the main difference that my issue is in preparation for a PR rather than a request. I'd like to discuss implementation details and to explore devs' opinion on how feasible such changes would be.

Please check the linked issue for additional information about why this would be a good idea. In short, (some) Shimadzu fNIRS devices work with 3 wavelengths, and in the future there might be other devices with more, but mne's SNIRF I/O and processing functions are coded to accept only 2 wavelengths.

From my perspective, these changes would be mostly beneficial as I'm using a Shimadzu device. While the device doesn't support .snirf output, there are ways to convert the data to that format. And while it's possible to preprocess and convert to Hb data first in some other tool(s) and then import to mne, it would be ideal to keep the entire processing pipeline to one tool.

Describe your proposed implementation

I've done some exploration in the codebase to identify relevant bits of code, both manually and with the help of Claude 4. Based on this, it seems that relevant parts of the code are in SNIRF I/O (mne.io.snirf) and NIRS preprocessing (mne.preprocessing.nirs). I'm sure some tests will need to be updated as well, but I haven't yet spent time on them.

I actually have a preliminary implementation of the described changes that you can see at my fork. That code is untested, but it should provide a general impression of the extent and complexity of required changes. Of course, tests will need to be added/updated in addition.

Major changes

• Channel data structure: with 2 wavelengths, channel data is stored as ABAB... . All functions handling OD and amplitude data will need to be updated not to rely on a fixed [::2], [1::2] sequence. I listed this as a major change due to the difficulty of locating every function that works this way. Claude thinks it's been able to find all of them (see list of files to be changed below). It still requires devs' experience to know if there are some I might not have been able to find.
• Beer-Lambert calculation (algorithmic change): The current implementation uses pseudo-inverse pinv() to calculate ΔOD = E × L × PPF × Δc. I think it's possible to update this with minimal changes, just changing the EL matrix to n x 2 where n is the number of wavelengths, instead of the current 2 x 2. After Beer-Lambert, excess channels will need to be dropped.
• Scalp coupling index (algorithmic change): I'm not aware of any reports where SCI was defined for >2 wavelengths. My idea here is to calculate all possible pairs of correlation, and use the minimum value. Intuitively, this should be correct as all wavelength channels should be correlated.

I'm hoping to receive feedback on the correctness of the above proposed algorithmic changes. For n=2 wavelengths, the new calculation methods wouldn't introduce any new changes (E in Beer-Lambert would still be 2x2, and there is only one correlation pair so the minimum is irrelevant), and if necessary, the current implementation could be preserved in a separate code path.

Minor changes

• SNIRF file reading: checks and validation steps need to be updated. While there's a number of validation checks, the overall complexity of the code changes isn't high.

List of affected files and functions:

• /mne/io/snirf/_snirf.py: class RawSNIRF(BaseRaw) has a simple check for 2 wavelengths and an error message.
• /mne/preprocessing/nirs/_beer_lambert_law.py: beer_lambert_law() logic needs to be updated to >=2 wavelengths, plus it will need to remove excess channels after assigning HbO/R data to the first 2.
• /mne/preprocessing/nirs/_beer_lambert_law.py: _load_absorption() is hard-coded to 2 wavelengths.
• /mne/preprocessing/nirs/_scalp_coupling_index.py: scalp_coupling_index() needs an algorithmic update
• /mne/preprocessing/nirs/nirs.py: _check_channels_ordered() has various validation checks and _fnirs_spread_bads() needs to spread bad markings over all wavelengths in group

Describe possible alternatives

The algorithmic changes (Beer-Lambert and SCI) could be implemented in another way. The Beer-Lambert pseudo-inverse should work with n>2 with minimal changes and the change seems mathematically correct. The SCI calculation using the minumum of all correlations appers to be the most conservative in terms of being able to spot any coupling issues, but without studies to rely on, it should be discussed.

Additional context

The project mne-python will also require some changes, which will I will discuss in another issue opened there.

[welcome]

Hello! 👋 Thanks for opening your first issue here! ❤️ We will try to get back to you soon. 🚴

[larsoner]

I'm sure some tests will need to be updated as well, but I haven't yet spent time on them.

Yep! And just to orient the thinking here, in an ideal world, your PR could leave all existing tests as-is. They test 2-wavelength functionality, which should be unchanged by the PR. If your PR added new tests with multi-wavelength support, that would be best. That way we know we have kept the existing 2-wavelength functionality as-is while adding support for >= 2.

I actually have a preliminary implementation of the described changes that you can see at my fork.

The easiest way to view these is actually as a PR, can you open one? You can open it and draft mode and say it's just for discussion. It's what I would do to understand what you had to change, so might as well open it for others to see as well!

It still requires devs' experience to know if there are some I might not have been able to find.

Have you used git grep yourself? Try that first. Something like this has some false alarms but worth looking through:

git grep "\[1::2\]"

mne/channels/channels.py:    left = np.sort(np.concatenate([left, lobe[medians[1::2]]]))
mne/decoding/_mod_ged.py:        ix[1::2] = i[: len(i) // 2]
mne/io/besa/besa.py:    name_val_pairs = zip(parts[::2], parts[1::2])
mne/io/edf/tests/test_edf.py:    annot[1::2] = [a * 256 for a in annot[1::2]]
mne/io/edf/tests/test_edf.py:        list(map(sum, zip(annot[0::2], annot[1::2]))), dtype=np.int64
mne/io/nirx/tests/test_nirx.py:    assert prefixes[::2] == prefixes[1::2]
mne/io/tests/test_raw.py:            (picks[1::2], picks[0::2], np.setdiff1d(all_picks, picks))
mne/preprocessing/nirs/_beer_lambert_law.py:    for ii, jj in zip(picks[::2], picks[1::2]):
mne/preprocessing/nirs/nirs.py:    for ii, jj in zip(picks[::2], picks[1::2]):
mne/preprocessing/nirs/nirs.py:        for ii, jj in zip(picks[::2], picks[1::2]):
mne/preprocessing/nirs/nirs.py:    for ii, jj in zip(picks[::2], picks[1::2]):
mne/preprocessing/nirs/tests/test_nirs.py:    got_second = set(raw_reversed.ch_names[pick].split()[1] for pick in picks[1::2])
mne/tests/test_annotations.py:        annot[1::2],
mne/tests/test_event.py:    assert_array_equal(events[1::2], events1)
mne/tests/test_evoked.py:    aves1 = read_evokeds(fname)[1::2]
mne/viz/_brain/tests/test_brain.py:    fp = (fp[1::2] + fp[::2]) * 0.5
mne/viz/backends/_pyvista.py:        focalpoint = (bounds[1::2] + bounds[::2]) * 0.5
mne/viz/backends/_pyvista.py:        distance = max(bounds[1::2] - bounds[::2]) * 2.0
mne/viz/evoked.py:    tick_locs = np.linspace(0, 1, 2 * len(color_vals) + 1)[1::2]
mne/viz/topomap.py:                mask[np.ix_(picks[::2], time_idx)] | mask[np.ix_(picks[1::2], time_idx)]

I think we should add support here then figure out mne-nirs.

[zEdS15B3GCwq]

Thanks for the reply!

• Re tests: indeed, the current changes shouldn't affect n=2, I was thinking about adding test cases for n>2.
• Re PR: I'll open a PR, and I agree about code change comparison being easier there, but the contributing guide suggests opening an issue first, hence no PR yet.
• Re git grep: haven't tried that one. Besides the parts I've already noticed and some likely-relevant tests, I can see probably affected code (picks[1::2] and similar in test_raw.py, and topomap.py. Based on the above results, which files would you say are suspect?