From e8baea60027e10423e86f0cf0ec9eaa50238fa47 Mon Sep 17 00:00:00 2001 From: Rob Luke Date: Tue, 19 Sep 2023 17:37:07 -0300 Subject: [PATCH] ENH: Add support for Artinis SNIRF data (#11926) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Eric Larson Co-authored-by: Daniel McCloy --- azure-pipelines.yml | 4 +-- doc/changes/devel.rst | 1 + mne/io/snirf/_snirf.py | 33 +++++++++++++++----- mne/io/snirf/tests/test_snirf.py | 44 ++++++++++++++++++++++++++- requirements_doc.txt | 2 +- requirements_testing_extra.txt | 2 +- tools/azure_dependencies.sh | 2 +- tutorials/io/30_reading_fnirs_data.py | 2 +- 8 files changed, 76 insertions(+), 14 deletions(-) diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 3ea1cd008a4..396cfe956b2 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -257,9 +257,9 @@ stages: 3.9 pip: TEST_MODE: 'pip' PYTHON_VERSION: '3.9' - 3.10 pip pre: + 3.11 pip pre: TEST_MODE: 'pip-pre' - PYTHON_VERSION: '3.10' + PYTHON_VERSION: '3.11' steps: - task: UsePythonVersion@0 inputs: diff --git a/doc/changes/devel.rst b/doc/changes/devel.rst index 3210c367d4d..4a3ade8c439 100644 --- a/doc/changes/devel.rst +++ b/doc/changes/devel.rst @@ -29,6 +29,7 @@ Enhancements - Added public :func:`mne.io.write_info` to complement :func:`mne.io.read_info` (:gh:`11918` by `Eric Larson`_) - Added option ``remove_dc`` to to :meth:`Raw.compute_psd() `, :meth:`Epochs.compute_psd() `, and :meth:`Evoked.compute_psd() `, to allow skipping DC removal when computing Welch or multitaper spectra (:gh:`11769` by `Nikolai Chapochnikov`_) - Add the possibility to provide a float between 0 and 1 as ``n_grad``, ``n_mag`` and ``n_eeg`` in `~mne.compute_proj_raw`, `~mne.compute_proj_epochs` and `~mne.compute_proj_evoked` to select the number of vectors based on the cumulative explained variance (:gh:`11919` by `Mathieu Scheltienne`_) +- Added support for Artinis fNIRS data files to :func:`mne.io.read_raw_snirf` (:gh:`11926` by `Robert Luke`_) - Add helpful error messages when using methods on empty :class:`mne.Epochs`-objects (:gh:`11306` by `Martin Schulz`_) - Add inferring EEGLAB files' montage unit automatically based on estimated head radius using :func:`read_raw_eeglab(..., montage_units="auto") ` (:gh:`11925` by `Jack Zhang`_, :gh:`11951` by `Eric Larson`_) - Add :class:`~mne.time_frequency.EpochsSpectrumArray` and :class:`~mne.time_frequency.SpectrumArray` to support creating power spectra from :class:`NumPy array ` data (:gh:`11803` by `Alex Rockhill`_) diff --git a/mne/io/snirf/_snirf.py b/mne/io/snirf/_snirf.py index 1a57282a961..5da8642301f 100644 --- a/mne/io/snirf/_snirf.py +++ b/mne/io/snirf/_snirf.py @@ -526,6 +526,11 @@ def _get_lengthunit_scaling(length_unit): def _extract_sampling_rate(dat): """Extract the sample rate from the time field.""" + # This is a workaround to provide support for Artinis data. + # It allows for a 1% variation in the sampling times relative + # to the average sampling rate of the file. + MAXIMUM_ALLOWED_SAMPLING_JITTER_PERCENTAGE = 1.0 + time_data = np.array(dat.get("nirs/data1/time")) sampling_rate = 0 if len(time_data) == 2: @@ -533,16 +538,30 @@ def _extract_sampling_rate(dat): sampling_rate = 1.0 / (time_data[1] - time_data[0]) else: # specified as time points - fs_diff = np.around(np.diff(time_data), decimals=4) - if len(np.unique(fs_diff)) == 1: + periods = np.diff(time_data) + uniq_periods = np.unique(periods.round(decimals=4)) + if uniq_periods.size == 1: # Uniformly sampled data - sampling_rate = 1.0 / np.unique(fs_diff).item() + sampling_rate = 1.0 / uniq_periods.item() else: - warn( - "MNE does not currently support reading " - "SNIRF files with non-uniform sampled data." + # Hopefully uniformly sampled data with some precision issues. + # This is a workaround to provide support for Artinis data. + mean_period = np.mean(periods) + sampling_rate = 1.0 / mean_period + ideal_times = np.linspace(time_data[0], time_data[-1], time_data.size) + max_jitter = np.max(np.abs(time_data - ideal_times)) + percent_jitter = 100.0 * max_jitter / mean_period + msg = ( + f"Found jitter of {percent_jitter:3f}% in sample times. Sampling " + f"rate has been set to {sampling_rate:1f}." ) - + if percent_jitter > MAXIMUM_ALLOWED_SAMPLING_JITTER_PERCENTAGE: + warn( + f"{msg} Note that MNE-Python does not currently support SNIRF " + "files with non-uniformly-sampled data." + ) + else: + logger.info(msg) time_unit = _get_metadata_str(dat, "TimeUnit") time_unit_scaling = _get_timeunit_scaling(time_unit) sampling_rate *= time_unit_scaling diff --git a/mne/io/snirf/tests/test_snirf.py b/mne/io/snirf/tests/test_snirf.py index 04a29e35726..141e34d00c6 100644 --- a/mne/io/snirf/tests/test_snirf.py +++ b/mne/io/snirf/tests/test_snirf.py @@ -18,6 +18,7 @@ ) from mne.transforms import apply_trans, _get_trans from mne._fiff.constants import FIFF +from mne.utils import catch_logging testing_path = data_path(download=False) @@ -471,7 +472,11 @@ def test_annotation_duration_from_stim_groups(): def test_birthday(tmp_path, monkeypatch): """Test birthday parsing.""" - snirf = pytest.importorskip("snirf") + try: + snirf = pytest.importorskip("snirf") + except AttributeError as exc: + # Until https://github.com/BUNPC/pysnirf2/pull/43 is released + pytest.skip(f"snirf import error: {exc}") fname = tmp_path / "test.snirf" with snirf.Snirf(str(fname), "w") as a: a.nirs.appendGroup() @@ -503,3 +508,40 @@ def test_birthday(tmp_path, monkeypatch): raw = read_raw_snirf(fname) assert raw.info["subject_info"]["birthday"] == (1950, 1, 1) + + +@requires_testing_data +def test_sample_rate_jitter(tmp_path): + """Test handling of jittered sample times.""" + from shutil import copy2 + + # Create a clean copy and ensure it loads without error + new_file = tmp_path / "snirf_nirsport2_2019.snirf" + copy2(snirf_nirsport2_20219, new_file) + read_raw_snirf(new_file) + + # Edit the file and add jitter within tolerance (0.99%) + with h5py.File(new_file, "r+") as f: + orig_time = np.array(f.get("nirs/data1/time")) + acceptable_time_jitter = orig_time.copy() + average_time_diff = np.mean(np.diff(orig_time)) + acceptable_time_jitter[-1] += 0.0099 * average_time_diff + del f["nirs/data1/time"] + f.flush() + f.create_dataset("nirs/data1/time", data=acceptable_time_jitter) + with catch_logging("info") as log: + read_raw_snirf(new_file) + lines = "\n".join(line for line in log.getvalue().splitlines() if "jitter" in line) + assert "Found jitter of 0.9" in lines + + # Add jitter of 1.01%, which is greater than allowed tolerance + with h5py.File(new_file, "r+") as f: + unacceptable_time_jitter = orig_time + unacceptable_time_jitter[-1] = unacceptable_time_jitter[-1] + ( + 0.0101 * average_time_diff + ) + del f["nirs/data1/time"] + f.flush() + f.create_dataset("nirs/data1/time", data=unacceptable_time_jitter) + with pytest.warns(RuntimeWarning, match="non-uniformly-sampled data"): + read_raw_snirf(new_file, verbose=True) diff --git a/requirements_doc.txt b/requirements_doc.txt index 7e2fa0c7986..a3cc904f979 100644 --- a/requirements_doc.txt +++ b/requirements_doc.txt @@ -1,5 +1,5 @@ # requirements for building docs -sphinx!=4.1.0,<6 +sphinx>=6 numpydoc pydata_sphinx_theme==0.13.3 git+https://github.com/sphinx-gallery/sphinx-gallery@master diff --git a/requirements_testing_extra.txt b/requirements_testing_extra.txt index 3c8e8986c80..d0318d9c707 100644 --- a/requirements_testing_extra.txt +++ b/requirements_testing_extra.txt @@ -7,4 +7,4 @@ EDFlib-Python pybv imageio>=2.6.1 imageio-ffmpeg>=0.4.1 -pysnirf2 +snirf diff --git a/tools/azure_dependencies.sh b/tools/azure_dependencies.sh index f5e16b635f5..072665d9c3c 100755 --- a/tools/azure_dependencies.sh +++ b/tools/azure_dependencies.sh @@ -9,7 +9,7 @@ elif [ "${TEST_MODE}" == "pip-pre" ]; then python -m pip install $STD_ARGS pip setuptools wheel packaging setuptools_scm python -m pip install $STD_ARGS --only-binary ":all:" --extra-index-url "https://www.riverbankcomputing.com/pypi/simple" PyQt6 PyQt6-sip PyQt6-Qt6 echo "Numpy etc." - python -m pip install $STD_ARGS --only-binary ":all:" --extra-index-url "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple" "numpy>=2.0.0.dev0" scipy statsmodels pandas scikit-learn matplotlib + python -m pip install $STD_ARGS --only-binary ":all:" --extra-index-url "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple" "numpy>=2.0.0.dev0" "scipy>=1.12.0.dev0" statsmodels pandas scikit-learn matplotlib echo "dipy" python -m pip install $STD_ARGS --only-binary ":all:" --extra-index-url "https://pypi.anaconda.org/scipy-wheels-nightly/simple" dipy echo "h5py" diff --git a/tutorials/io/30_reading_fnirs_data.py b/tutorials/io/30_reading_fnirs_data.py index b1c7a45e630..b2208c85941 100644 --- a/tutorials/io/30_reading_fnirs_data.py +++ b/tutorials/io/30_reading_fnirs_data.py @@ -38,7 +38,7 @@ is designed by the fNIRS community in an effort to facilitate sharing and analysis of fNIRS data. And is the official format of the Society for functional near-infrared spectroscopy (SfNIRS). -The manufacturers Gowerlabs, NIRx, Kernel, and Cortivision +The manufacturers Gowerlabs, NIRx, Kernel, Artinis, and Cortivision export data in the SNIRF format, and these files can be imported in to MNE. SNIRF is the preferred format for reading data in to MNE-Python. Data stored in the SNIRF format can be read in