From b9d182ca2c0fc85ac3777238915367493ee4bd73 Mon Sep 17 00:00:00 2001 From: Leonardo Zaggia Date: Wed, 25 Feb 2026 18:37:47 +0100 Subject: [PATCH 1/5] ENH: add PCA-based motion correction for fNIRS (motion_correct_pca) --- doc/api/preprocessing.rst | 2 + doc/changes/dev/13691.newfeature.rst | 1 + doc/changes/names.inc | 1 + mne/preprocessing/nirs/__init__.py | 1 + mne/preprocessing/nirs/_pca.py | 179 +++++++++++++++++++++++ mne/preprocessing/nirs/tests/test_pca.py | 125 ++++++++++++++++ 6 files changed, 309 insertions(+) create mode 100644 doc/changes/dev/13691.newfeature.rst create mode 100644 mne/preprocessing/nirs/_pca.py create mode 100644 mne/preprocessing/nirs/tests/test_pca.py diff --git a/doc/api/preprocessing.rst b/doc/api/preprocessing.rst index f1d60cfdc3b..ccb2578aefc 100644 --- a/doc/api/preprocessing.rst +++ b/doc/api/preprocessing.rst @@ -137,6 +137,8 @@ Projections: short_channels scalp_coupling_index temporal_derivative_distribution_repair + motion_correct_pca + pca :py:mod:`mne.preprocessing.ieeg`: diff --git a/doc/changes/dev/13691.newfeature.rst b/doc/changes/dev/13691.newfeature.rst new file mode 100644 index 00000000000..104bf295243 --- /dev/null +++ b/doc/changes/dev/13691.newfeature.rst @@ -0,0 +1 @@ +Add :func:`mne.preprocessing.nirs.motion_correct_pca` (alias :func:`mne.preprocessing.nirs.pca`) for PCA-based motion correction of fNIRS data, based on Homer3 ``hmrR_MotionCorrectPCA``, by :newcontrib:`Leonardo Zaggia`. diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 79db14d44e3..03c87ccc16a 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -181,6 +181,7 @@ .. _Laurent Le Mentec: https://github.com/LaurentLM .. _Leonardo Barbosa: https://github.com/noreun .. _Leonardo Rochael Almeida: https://github.com/leorochael +.. _Leonardo Zaggia: https://github.com/leonardozaggia .. _Liberty Hamilton: https://github.com/libertyh .. _Lorenzo Desantis: https://github.com/lorenzo-desantis/ .. _Lukas Breuer: https://www.researchgate.net/profile/Lukas-Breuer-2 diff --git a/mne/preprocessing/nirs/__init__.py b/mne/preprocessing/nirs/__init__.py index 6e35b59dd15..a8f1a792f0d 100644 --- a/mne/preprocessing/nirs/__init__.py +++ b/mne/preprocessing/nirs/__init__.py @@ -20,3 +20,4 @@ from ._beer_lambert_law import beer_lambert_law from ._scalp_coupling_index import scalp_coupling_index from ._tddr import temporal_derivative_distribution_repair, tddr +from ._pca import motion_correct_pca, pca diff --git a/mne/preprocessing/nirs/_pca.py b/mne/preprocessing/nirs/_pca.py new file mode 100644 index 00000000000..d9c399d88f2 --- /dev/null +++ b/mne/preprocessing/nirs/_pca.py @@ -0,0 +1,179 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import numpy as np +from scipy.linalg import svd + +from ...io import BaseRaw +from ...utils import _validate_type, verbose +from ..nirs import _validate_nirs_info + + +@verbose +def motion_correct_pca(raw, tInc, nSV=0.97, *, verbose=None): + """Apply PCA-based motion correction to fNIRS data. + + Extracts the motion-artifact portions of the signal, performs a singular + value decomposition across all fNIRS channels, removes the dominant + principal components, and reinserts the cleaned segments. + + Based on Homer3 v1.80.2 ``hmrR_MotionCorrectPCA.m`` + :footcite:`HuppertEtAl2009`. + + Parameters + ---------- + raw : instance of Raw + The raw fNIRS data (optical density or hemoglobin). + tInc : array-like of bool, shape (n_times,) + Global motion-artifact mask. ``True`` = clean sample, + ``False`` = motion artifact. + nSV : float | int + Number of principal components to remove. + + * ``0 < nSV < 1`` – remove the fewest components whose cumulative + explained variance reaches ``nSV`` (e.g. ``0.97`` removes + components explaining the top 97 %% of variance). + * ``nSV >= 1`` – remove exactly ``int(nSV)`` components. + + Default is ``0.97``. + %(verbose)s + + Returns + ------- + raw : instance of Raw + Data with PCA motion correction applied (copy). + svs : ndarray + Normalised singular values (sum = 1) of the PCA decomposition. + nSV : int + Actual number of principal components removed. + + Notes + ----- + There is a shorter alias ``mne.preprocessing.nirs.pca`` that + can be used instead of this function. + + References + ---------- + .. footbibliography:: + """ + _validate_type(raw, BaseRaw, "raw") + raw = raw.copy().load_data() + picks = _validate_nirs_info(raw.info) + + if not len(picks): + raise RuntimeError( + "PCA motion correction should be run on optical density or hemoglobin data." + ) + + tInc = np.asarray(tInc, dtype=bool) + n_times = raw._data.shape[1] + + # (n_picks, n_times) + y = raw._data[picks, :] + + # Motion-artifact samples only, shape (n_motion_samples, n_picks) + y_motion = y[:, ~tInc].T + + if y_motion.shape[0] == 0: + raise ValueError( + "No motion-artifact samples found (tInc is all True). " + "Mark artifact regions with False before calling this function." + ) + + # Z-score across time within the motion segments + y_mean = y_motion.mean(axis=0) + y_std = y_motion.std(axis=0) + y_std[y_std == 0] = 1.0 # avoid divide-by-zero for flat channels + y_zscore = (y_motion - y_mean) / y_std + + # PCA via SVD of the covariance-like matrix + yo = y_zscore.copy() + c = np.dot(y_zscore.T, y_zscore) + V, St, _ = svd(c) + + svs = St / np.sum(St) + + # Cumulative variance for fractional nSV selection + svsc = svs.copy() + for i in range(1, len(svs)): + svsc[i] = svsc[i - 1] + svs[i] + + if 0 < nSV < 1: + mask_keep = svsc < nSV + nSV = int(np.where(~mask_keep)[0][0]) + else: + nSV = int(nSV) + + ev = np.zeros((len(svs), 1)) + ev[:nSV] = 1 + ev = np.diag(np.squeeze(ev)) + + # Remove top PCs from motion-artifact segments + yc = yo - np.dot(np.dot(yo, V), np.dot(ev, V.T)) + yc = (yc * y_std) + y_mean # back to original scale + + # Identify motion-segment boundaries (1=good → 0=bad transitions) + lst_ms = np.where(np.diff(tInc.astype(int)) == -1)[0] # motion starts + lst_mf = np.where(np.diff(tInc.astype(int)) == 1)[0] # motion ends + + if len(lst_ms) == 0: + lst_ms = np.asarray([0]) + if len(lst_mf) == 0: + lst_mf = np.asarray([n_times - 1]) + if lst_ms[0] > lst_mf[0]: + lst_ms = np.insert(lst_ms, 0, 0) + if lst_ms[-1] > lst_mf[-1]: + lst_mf = np.append(lst_mf, n_times - 1) + + # Cumulative lengths used to index into yc + lst_mb = lst_mf - lst_ms + for ii in range(1, len(lst_mb)): + lst_mb[ii] = lst_mb[ii - 1] + lst_mb[ii] + lst_mb = lst_mb - 1 + + n_picks = len(picks) + cleaned_ts = y.T.copy() # (n_times, n_picks) + orig_ts = y.T.copy() + + for jj in range(n_picks): + # First motion segment + lst = np.arange(lst_ms[0], lst_mf[0]) + if lst_ms[0] > 0: + cleaned_ts[lst, jj] = ( + yc[: lst_mb[0] + 1, jj] - yc[0, jj] + cleaned_ts[lst[0], jj] + ) + else: + cleaned_ts[lst, jj] = ( + yc[: lst_mb[0] + 1, jj] - yc[lst_mb[0], jj] + cleaned_ts[lst[-1], jj] + ) + + # Intermediate non-motion and motion segments + for kk in range(len(lst_mf) - 1): + # Non-motion gap between MA[kk] and MA[kk+1] + lst = np.arange(lst_mf[kk] - 1, lst_ms[kk + 1] + 1) + cleaned_ts[lst, jj] = ( + orig_ts[lst, jj] - orig_ts[lst[0], jj] + cleaned_ts[lst[0], jj] + ) + + # Next motion segment + lst = np.arange(lst_ms[kk + 1], lst_mf[kk + 1]) + cleaned_ts[lst, jj] = ( + yc[lst_mb[kk] + 1 : lst_mb[kk + 1] + 1, jj] + - yc[lst_mb[kk] + 1, jj] + + cleaned_ts[lst[0], jj] + ) + + # Trailing non-motion segment + if lst_mf[-1] < n_times - 1: + lst = np.arange(lst_mf[-1] - 1, n_times) + cleaned_ts[lst, jj] = ( + orig_ts[lst, jj] - orig_ts[lst[0], jj] + cleaned_ts[lst[0], jj] + ) + + raw._data[picks, :] = cleaned_ts.T + return raw, svs, nSV + + +# provide a short alias +pca = motion_correct_pca diff --git a/mne/preprocessing/nirs/tests/test_pca.py b/mne/preprocessing/nirs/tests/test_pca.py new file mode 100644 index 00000000000..bc2d7f6bab6 --- /dev/null +++ b/mne/preprocessing/nirs/tests/test_pca.py @@ -0,0 +1,125 @@ +# Authors: The MNE-Python contributors. +# License: BSD-3-Clause +# Copyright the MNE-Python contributors. + +import numpy as np +import pytest +from numpy.testing import assert_allclose + +from mne.datasets import testing +from mne.datasets.testing import data_path +from mne.io import read_raw_nirx +from mne.preprocessing.nirs import ( + _validate_nirs_info, + motion_correct_pca, + optical_density, + pca, +) + +fname_nirx_15_2 = ( + data_path(download=False) / "NIRx" / "nirscout" / "nirx_15_2_recording" +) + + +def _make_tinc(n_times, bad_start, bad_stop): + """Return a per-time boolean mask (True=clean, False=motion).""" + tInc = np.ones(n_times, dtype=bool) + tInc[bad_start:bad_stop] = False + return tInc + + +@testing.requires_testing_data +@pytest.mark.parametrize("fname", ([fname_nirx_15_2])) +def test_motion_correct_pca_removes_shared_artifact(fname): + """Test PCA correction reduces a correlated motion artefact.""" + raw = read_raw_nirx(fname) + raw_od = optical_density(raw) + picks = _validate_nirs_info(raw_od.info) + n_times = raw_od._data.shape[1] + + # Inject a correlated step shift in the middle of the recording + mid = n_times // 2 + bad_start, bad_stop = mid - 10, mid + 10 + + # Save clean signal for comparison + original = raw_od._data[picks[0]].copy() + + max_shift = np.max(np.abs(np.diff(raw_od._data[picks[0]]))) + shift_amp = 20 * max_shift + + for pick in picks: + raw_od._data[pick, bad_start:bad_stop] += shift_amp + + tInc = _make_tinc(n_times, bad_start, bad_stop) + raw_od_corr, svs, nSV_ret = motion_correct_pca(raw_od, tInc=tInc, nSV=0.97) + + # Corrected signal should be closer to original than corrupted signal + mse_before = np.mean((raw_od._data[picks[0]] - original) ** 2) + mse_after = np.mean((raw_od_corr._data[picks[0]] - original) ** 2) + assert mse_after < mse_before + + +@testing.requires_testing_data +@pytest.mark.parametrize("fname", ([fname_nirx_15_2])) +def test_motion_correct_pca_svs_sum_to_one(fname): + """Test returned singular values are normalised and sum to 1.""" + raw = read_raw_nirx(fname) + raw_od = optical_density(raw) + n_times = raw_od._data.shape[1] + + tInc = _make_tinc(n_times, n_times // 2 - 5, n_times // 2 + 5) + _, svs, _ = motion_correct_pca(raw_od, tInc=tInc, nSV=0.97) + + assert_allclose(np.sum(svs), 1.0, rtol=1e-6) + + +@testing.requires_testing_data +@pytest.mark.parametrize("fname", ([fname_nirx_15_2])) +def test_motion_correct_pca_integer_nsv(fname): + """Test integer nSV removes exactly that many components.""" + raw = read_raw_nirx(fname) + raw_od = optical_density(raw) + n_times = raw_od._data.shape[1] + + tInc = _make_tinc(n_times, n_times // 3, n_times // 3 + 20) + _, _, nSV_ret = motion_correct_pca(raw_od, tInc=tInc, nSV=2) + assert nSV_ret == 2 + + +@testing.requires_testing_data +@pytest.mark.parametrize("fname", ([fname_nirx_15_2])) +def test_motion_correct_pca_returns_copy(fname): + """Test PCA correction does not modify the input Raw in place.""" + raw = read_raw_nirx(fname) + raw_od = optical_density(raw) + picks = _validate_nirs_info(raw_od.info) + n_times = raw_od._data.shape[1] + original = raw_od._data[picks[0]].copy() + + tInc = _make_tinc(n_times, 10, 30) + _, _, _ = motion_correct_pca(raw_od, tInc=tInc) + assert_allclose(raw_od._data[picks[0]], original) + + +@testing.requires_testing_data +@pytest.mark.parametrize("fname", ([fname_nirx_15_2])) +def test_motion_correct_pca_all_good_raises(fname): + """Test PCA correction raises when tInc has no artefact samples.""" + raw = read_raw_nirx(fname) + raw_od = optical_density(raw) + n_times = raw_od._data.shape[1] + + tInc = np.ones(n_times, dtype=bool) # all clean – no motion to correct + with pytest.raises(ValueError, match="No motion-artifact samples"): + motion_correct_pca(raw_od, tInc=tInc) + + +def test_pca_alias(): + """Test pca is an alias for motion_correct_pca.""" + assert pca is motion_correct_pca + + +def test_motion_correct_pca_wrong_type(): + """Test passing a non-Raw object raises TypeError.""" + with pytest.raises(TypeError): + motion_correct_pca(np.zeros((10, 100)), tInc=np.ones(100, dtype=bool)) From 956fb096504bf84a4772e523da0651710354e266 Mon Sep 17 00:00:00 2001 From: Leonardo Zaggia Date: Wed, 25 Feb 2026 19:19:48 +0100 Subject: [PATCH 2/5] DOC: add HuppertEtAl2009 and MolaviDumont2012 to references.bib --- doc/references.bib | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/doc/references.bib b/doc/references.bib index 12d63458bd1..cad83673772 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -866,6 +866,17 @@ @article{HouckClaus2020 year = {2020} } +@article{HuppertEtAl2009, + author = {Huppert, Theodore J. and Diamond, Solomon G. and Franceschini, Maria A. and Boas, David A.}, + doi = {10.1364/AO.48.00D280}, + journal = {Applied Optics}, + number = {10}, + pages = {D280-D298}, + title = {{{HomER}}: a review of time-series analysis methods for near-infrared spectroscopy of the brain}, + volume = {48}, + year = {2009} +} + @article{Hyvarinen1999, author = {Hyvärinen, Aapo}, doi = {10.1109/72.761722}, @@ -1254,6 +1265,17 @@ @misc{Mills2016 year = {2016} } +@article{MolaviDumont2012, + author = {Molavi, Behnam and Dumont, Guy A.}, + doi = {10.1088/0967-3334/33/2/259}, + journal = {Physiological Measurement}, + number = {2}, + pages = {259-270}, + title = {Wavelet-based motion artifact removal for functional near-infrared spectroscopy}, + volume = {33}, + year = {2012} +} + @article{MolinsEtAl2008, author = {Molins A, and Stufflebeam S. M., and Brown E. N., and Hämäläinen M. S.}, doi = {10.1016/j.neuroimage.2008.05.064}, From 927d0f5dd0443fbc6034b05936518ac8966b302d Mon Sep 17 00:00:00 2001 From: Leonardo Zaggia Date: Wed, 25 Feb 2026 19:28:03 +0100 Subject: [PATCH 3/5] DOC: remove short alias from API reference (only long name needed) --- doc/api/preprocessing.rst | 1 - 1 file changed, 1 deletion(-) diff --git a/doc/api/preprocessing.rst b/doc/api/preprocessing.rst index ccb2578aefc..32a40c12992 100644 --- a/doc/api/preprocessing.rst +++ b/doc/api/preprocessing.rst @@ -138,7 +138,6 @@ Projections: scalp_coupling_index temporal_derivative_distribution_repair motion_correct_pca - pca :py:mod:`mne.preprocessing.ieeg`: From 04db60f510870997b0a4fe6e293aefd6ff7bced0 Mon Sep 17 00:00:00 2001 From: Leonardo Zaggia Date: Wed, 25 Feb 2026 20:08:37 +0100 Subject: [PATCH 4/5] DOC: fix :func: reference for short alias in changelog --- doc/changes/dev/13691.newfeature.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/changes/dev/13691.newfeature.rst b/doc/changes/dev/13691.newfeature.rst index 104bf295243..674d0a98c7d 100644 --- a/doc/changes/dev/13691.newfeature.rst +++ b/doc/changes/dev/13691.newfeature.rst @@ -1 +1 @@ -Add :func:`mne.preprocessing.nirs.motion_correct_pca` (alias :func:`mne.preprocessing.nirs.pca`) for PCA-based motion correction of fNIRS data, based on Homer3 ``hmrR_MotionCorrectPCA``, by :newcontrib:`Leonardo Zaggia`. +Add :func:`mne.preprocessing.nirs.motion_correct_pca` (alias ``mne.preprocessing.nirs.pca``) for PCA-based motion correction of fNIRS data, based on Homer3 ``hmrR_MotionCorrectPCA``, by :newcontrib:`Leonardo Zaggia`. From 8d5f42973c28ccadaefd94c8ebeb46a42f8bc8ca Mon Sep 17 00:00:00 2001 From: Leonardo Zaggia Date: Thu, 26 Feb 2026 14:00:46 +0100 Subject: [PATCH 5/5] add inline attribution to Cedalion and Homer3 --- mne/preprocessing/nirs/_pca.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/mne/preprocessing/nirs/_pca.py b/mne/preprocessing/nirs/_pca.py index d9c399d88f2..cbca7ad4a85 100644 --- a/mne/preprocessing/nirs/_pca.py +++ b/mne/preprocessing/nirs/_pca.py @@ -2,6 +2,10 @@ # License: BSD-3-Clause # Copyright the MNE-Python contributors. +# The core logic for this implementation was adapted from the Cedalion project +# (https://github.com/ibs-lab/cedalion), which is originally based on Homer3 +# (https://github.com/BUNPC/Homer3). + import numpy as np from scipy.linalg import svd