ENH: Add support for Artinis SNIRF data (mne-tools#11926)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric Larson <[email protected]>
Co-authored-by: Daniel McCloy <[email protected]>

Author: 4 people

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:

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() <mne.io.Raw.compute_psd>`, :meth:`Epochs.compute_psd() <mne.Epochs.compute_psd>`, and :meth:`Evoked.compute_psd() <mne.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") <mne.io.read_raw_eeglab>` (: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 <numpy.ndarray>` data (:gh:`11803` by `Alex Rockhill`_)

mne/io/snirf/_snirf.py

@@ -526,23 +526,42 @@ 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:
         # specified as onset, samplerate
         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

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)

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

requirements_testing_extra.txt

@@ -7,4 +7,4 @@ EDFlib-Python
 pybv
 imageio>=2.6.1
 imageio-ffmpeg>=0.4.1
-pysnirf2
+snirf

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"

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