Add KIT phantom dataset (mne-tools#12105)

Co-authored-by: Eric Larson <[email protected]>

.circleci/config.yml

@@ -266,6 +266,9 @@ jobs:
       - restore_cache:
           keys:
             - data-cache-ucl-opm-auditory
+      - restore_cache:
+          keys:
+            - data-cache-phantom-kit
       - run:
           name: Get data
           # This limit could be increased, but this is helpful for finding slow ones
@@ -431,6 +434,10 @@ jobs:
           key: data-cache-ucl-opm-auditory
           paths:
             - ~/mne_data/auditory_OPM_stationary  # (4 G)
+      - save_cache:
+          key: data-cache-phantom-kit
+          paths:
+            - ~/mne_data/MNE-phantom-KIT-data  # (1 G)
 
 
   linkcheck:

.pre-commit-config.yaml

@@ -22,7 +22,7 @@ repos:
         name: ruff tutorials and examples
         # D103: missing docstring in public function
         # D400: docstring first line must end with period
-        args: ["--ignore=D103,D400"]
+        args: ["--ignore=D103,D400", "--fix"]
         files: ^tutorials/|^examples/
 
   # Codespell

doc/api/datasets.rst

@@ -40,10 +40,11 @@ Datasets
    spm_face.data_path
    ucl_opm_auditory.data_path
    visual_92_categories.data_path
+   phantom_kit.data_path
    phantom_4dbti.data_path
    phantom_kernel.data_path
    refmeg_noise.data_path
    ssvep.data_path
    erp_core.data_path
    epilepsy_ecog.data_path
-   eyelink.data_path
+   eyelink.data_path

doc/changes/devel.rst

@@ -44,6 +44,7 @@ Enhancements
 - By default MNE-Python creates matplotlib figures with ``layout='constrained'`` rather than the default ``layout='tight'`` (:gh:`12050`, :gh:`12103` by `Mathieu Scheltienne`_ and `Eric Larson`_)
 - Enhance :func:`~mne.viz.plot_evoked_field` with a GUI that has controls for time, colormap, and contour lines (:gh:`11942` by `Marijn van Vliet`_)
 - Add :class:`mne.viz.ui_events.UIEvent` linking for interactive colorbars, allowing users to link figures and change the colormap and limits interactively. This supports :func:`~mne.viz.plot_evoked_topomap`, :func:`~mne.viz.plot_ica_components`, :func:`~mne.viz.plot_tfr_topomap`, :func:`~mne.viz.plot_projs_topomap`, :meth:`~mne.Evoked.plot_image`, and :meth:`~mne.Epochs.plot_image` (:gh:`12057` by `Santeri Ruuskanen`_)
+- Add example KIT phantom dataset in :func:`mne.datasets.phantom_kit.data_path` and :ref:`tut-phantom-kit` (:gh:`12105` by `Judy D Zhu`_ and `Eric Larson`_)
 - :func:`~mne.epochs.make_metadata` now accepts ``tmin=None`` and ``tmax=None``, which will bound the time window used for metadata generation by event names (instead of a fixed time). That way, you can now for example generate metadata spanning from one cue or fixation cross to the next, even if trial durations vary throughout the recording (:gh:`12118` by `Richard Höchenberger`_)
 - Add support for passing multiple labels to :func:`mne.minimum_norm.source_induced_power` (:gh:`12026` by `Erica Peterson`_, `Eric Larson`_, and `Daniel McCloy`_ )
 - Added documentation to :meth:`mne.io.Raw.set_montage` and :func:`mne.add_reference_channels` to specify that montages should be set after adding reference channels (:gh:`12160` by `Jacob Woessner`_)

doc/documentation/datasets.rst

@@ -282,6 +282,18 @@ are richly annotated, and can be used for e.g. multiple regression estimation
 of EEG correlates of printed word processing.
 
 
+KIT phantom dataset
+=============================
+:func:`mne.datasets.phantom_kit.data_path`.
+
+This dataset was obtained with a phantom on a KIT system at
+Macquarie University in Sydney, Australia.
+
+.. topic:: Examples
+
+    * :ref:`tut-phantom-KIT`
+
+
 4D Neuroimaging / BTi dataset
 =============================
 :func:`mne.datasets.phantom_4dbti.data_path`.

doc/references.bib

@@ -2449,3 +2449,16 @@ @article{TierneyEtAl2022
 	doi = {10.1016/j.neuroimage.2022.119338},
 	author = {Tierney, Tim M. and Mellor, Stephanie nd O'Neill, George C. and Timms, Ryan C. and Barnes, Gareth R.},
 }
+
+
+@article{OyamaEtAl2015,
+	title = {Dry phantom for magnetoencephalography —{Configuration}, calibration, and contribution},
+	volume = {251},
+	issn = {0165-0270},
+	doi = {10.1016/j.jneumeth.2015.05.004},
+	journal = {Journal of Neuroscience Methods},
+	author = {Oyama, Daisuke and Adachi, Yoshiaki and Yumoto, Masato and Hashimoto, Isao and Uehara, Gen},
+	month = aug,
+	year = {2015},
+	pages = {24--36},
+}

mne/datasets/__init__.pyi

@@ -24,6 +24,7 @@ __all__ = [
     "opm",
     "phantom_4dbti",
     "phantom_kernel",
+    "phantom_kit",
     "refmeg_noise",
     "sample",
     "sleep_physionet",
@@ -52,6 +53,7 @@ from . import (
     opm,
     phantom_4dbti,
     phantom_kernel,
+    phantom_kit,
     refmeg_noise,
     sample,
     sleep_physionet,

mne/datasets/config.py

@@ -174,6 +174,14 @@
     config_key="MNE_DATASETS_OPM_PATH",
 )
 
+MNE_DATASETS["phantom_kit"] = dict(
+    archive_name="MNE-phantom-KIT-24bit.zip",
+    hash="md5:CAF82EE978DD473C7DE6C1034D9CCD45",
+    url="https://osf.io/download/svnt3/",
+    folder_name="MNE-phantom-KIT-data",
+    config_key="MNE_DATASETS_PHANTOM_KIT_PATH",
+)
+
 MNE_DATASETS["phantom_4dbti"] = dict(
     archive_name="MNE-phantom-4DBTi.zip",
     hash="md5:938a601440f3ffa780d20a17bae039ff",

mne/datasets/phantom_kit/__init__.py

@@ -0,0 +1,3 @@
+"""KIT phantom dataset."""
+
+from .phantom_kit import data_path, get_version

mne/datasets/phantom_kit/phantom_kit.py

@@ -0,0 +1,28 @@
+from ...utils import verbose
+from ..utils import _data_path_doc, _download_mne_dataset, _get_version, _version_doc
+
+
+@verbose
+def data_path(
+    path=None, force_update=False, update_path=True, download=True, *, verbose=None
+):  # noqa: D103
+    return _download_mne_dataset(
+        name="phantom_kit",
+        processor="unzip",
+        path=path,
+        force_update=force_update,
+        update_path=update_path,
+        download=download,
+    )
+
+
+data_path.__doc__ = _data_path_doc.format(
+    name="phantom_kit", conf="MNE_DATASETS_PHANTOM_KIT_PATH"
+)
+
+
+def get_version():  # noqa: D103
+    return _get_version("phantom_kit")
+
+
+get_version.__doc__ = _version_doc.format(name="phantom_kit")

mne/datasets/utils.py

@@ -324,7 +324,7 @@ def _download_all_example_data(verbose=True):
     paths = dict()
     for kind in (
         "sample testing misc spm_face somato hf_sef multimodal "
-        "fnirs_motor opm mtrf fieldtrip_cmc kiloword phantom_4dbti "
+        "fnirs_motor opm mtrf fieldtrip_cmc kiloword phantom_kit phantom_4dbti "
         "refmeg_noise ssvep epilepsy_ecog ucl_opm_auditory eyelink "
         "erp_core brainstorm.bst_raw brainstorm.bst_auditory "
         "brainstorm.bst_resting brainstorm.bst_phantom_ctf "

mne/dipole.py

@@ -1760,6 +1760,36 @@ def fit_dipole(
     return dipoles, residual
 
 
+# Every other row of Table 3 from OyamaEtAl2015
+_OYAMA = """
+0.00 56.29 -27.50
+32.50 56.29 5.00
+0.00 65.00 5.00
+-32.50 56.29 5.00
+0.00 56.29 37.50
+0.00 32.50 61.29
+-56.29 0.00 -27.50
+-56.29 32.50 5.00
+-65.00 0.00 5.00
+-56.29 -32.50 5.00
+-56.29 0.00 37.50
+-32.50 0.00 61.29
+0.00 -56.29 -27.50
+-32.50 -56.29 5.00
+0.00 -65.00 5.00
+32.50 -56.29 5.00
+0.00 -56.29 37.50
+0.00 -32.50 61.29
+56.29 0.00 -27.50
+56.29 -32.50 5.00
+65.00 0.00 5.00
+56.29 32.50 5.00
+56.29 0.00 37.50
+32.50 0.00 61.29
+0.00 0.00 70.00
+"""
+
+
 def get_phantom_dipoles(kind="vectorview"):
     """Get standard phantom dipole locations and orientations.
 
@@ -1772,6 +1802,11 @@ def get_phantom_dipoles(kind="vectorview"):
               The Neuromag VectorView phantom.
             ``otaniemi``
               The older Neuromag phantom used at Otaniemi.
+            ``oyama``
+              The phantom from :footcite:`OyamaEtAl2015`.
+
+        .. versionchanged:: 1.6
+           Support added for ``'oyama'``.
 
     Returns
     -------
@@ -1788,8 +1823,13 @@ def get_phantom_dipoles(kind="vectorview"):
     -----
     The Elekta phantoms have a radius of 79.5mm, and HPI coil locations
     in the XY-plane at the axis extrema (e.g., (79.5, 0), (0, -79.5), ...).
+
+    References
+    ----------
+    .. footbibliography::
     """
-    _check_option("kind", kind, ["vectorview", "otaniemi"])
+    _validate_type(kind, str, "kind")
+    _check_option("kind", kind, ["vectorview", "otaniemi", "oyama"])
     if kind == "vectorview":
         # these values were pulled from a scanned image provided by
         # Elekta folks
@@ -1811,18 +1851,43 @@ def get_phantom_dipoles(kind="vectorview"):
         y = np.concatenate((c, c, -a, -b, c, c, b, a))
         z = np.concatenate((b, a, b, a, b, a, a, b))
         signs = [-1] * 8 + [1] * 16 + [-1] * 8
+    else:
+        assert kind == "oyama"
+        xyz = np.fromstring(_OYAMA.strip().replace("\n", " "), sep=" ").reshape(25, 3)
+        xyz = np.repeat(xyz, 2, axis=0)
+        x, y, z = xyz.T
+        signs = [1] * 50
     pos = np.vstack((x, y, z)).T / 1000.0
+    # For Neuromag-style phantoms,
     # Locs are always in XZ or YZ, and so are the oris. The oris are
     # also in the same plane and tangential, so it's easy to determine
     # the orientation.
+    # For Oyama, vectors are orthogonal to the position vector and oriented with one
+    # pointed toward the north pole (except for the topmost points, which are just xy).
     ori = list()
     for pi, this_pos in enumerate(pos):
         this_ori = np.zeros(3)
         idx = np.where(this_pos == 0)[0]
         # assert len(idx) == 1
+        if len(idx) == 0:  # oyama
+            idx = [np.argmin(this_pos)]
         idx = np.setdiff1d(np.arange(3), idx[0])
         this_ori[idx] = (this_pos[idx][::-1] / np.linalg.norm(this_pos[idx])) * [1, -1]
-        this_ori *= signs[pi]
+        if kind == "oyama":
+            # Ensure it's orthogonal to the position vector
+            pos_unit = this_pos / np.linalg.norm(this_pos)
+            this_ori -= pos_unit * np.dot(this_ori, pos_unit)
+            this_ori /= np.linalg.norm(this_ori)
+            # This was empirically determined by looking at the dipole fits
+            if np.abs(this_ori[2]) >= 1e-6:  # if it's not in the XY plane
+                this_ori *= -1 * np.sign(this_ori[2])  # point downward
+            elif np.abs(this_ori[0]) < 1e-6:  # in the XY plane (at the north pole)
+                this_ori *= -1 * np.sign(this_ori[1])  # point backward
+            # Odd ones create a RH coordinate system with their ori
+            if pi % 2:
+                this_ori = np.cross(pos_unit, this_ori)
+        else:
+            this_ori *= signs[pi]
         # Now we have this quality, which we could uncomment to
         # double-check:
         # np.testing.assert_allclose(np.dot(this_ori, this_pos) /

mne/event.py

@@ -25,6 +25,7 @@
     _check_option,
     _get_stim_channel,
     _on_missing,
+    _pl,
     _validate_type,
     check_fname,
     fill_doc,
@@ -482,6 +483,7 @@ def find_stim_steps(raw, pad_start=None, pad_stop=None, merge=0, stim_channel=No
 def _find_events(
     data,
     first_samp,
+    *,
     verbose=None,
     output="onset",
     consecutive="increasing",
@@ -490,6 +492,7 @@ def _find_events(
     uint_cast=False,
     mask_type="and",
     initial_event=False,
+    ch_name=None,
 ):
     """Help find events."""
     assert data.shape[0] == 1  # data should be only a row vector
@@ -520,9 +523,9 @@ def _find_events(
             events = np.insert(events, 0, [first_samp, 0, initial_value], axis=0)
         else:
             logger.info(
-                "Trigger channel has a non-zero initial value of {} "
-                "(consider using initial_event=True to detect this "
-                "event)".format(initial_value)
+                f"Trigger channel {ch_name} has a non-zero initial value of "
+                "{initial_value} (consider using initial_event=True to detect this "
+                "event)"
             )
 
     events = _mask_trigs(events, mask, mask_type)
@@ -555,22 +558,22 @@ def _find_events(
         logger.info("Removing orphaned onset at the end of the file.")
         onset_idx = np.delete(onset_idx, -1)
 
+    _check_option("output", output, ("onset", "step", "offset"))
     if output == "onset":
         events = events[onset_idx]
     elif output == "step":
         idx = np.union1d(onset_idx, offset_idx)
         events = events[idx]
-    elif output == "offset":
+    else:
+        assert output == "offset"
         event_id = events[onset_idx, 2]
         events = events[offset_idx]
         events[:, 1] = events[:, 2]
         events[:, 2] = event_id
         events[:, 0] -= 1
-    else:
-        raise ValueError("Invalid output parameter %r" % output)
 
-    logger.info("%s events found" % len(events))
-    logger.info("Event IDs: %s" % np.unique(events[:, 2]))
+    logger.info(f"{len(events)} event{_pl(events)} found on stim channel {ch_name}")
+    logger.info(f"Event IDs: {np.unique(events[:, 2])}")
 
     return events
 
@@ -772,7 +775,7 @@ def find_events(
     data, _ = raw[picks, :]
 
     events_list = []
-    for d in data:
+    for d, ch_name in zip(data, stim_channel):
         events = _find_events(
             d[np.newaxis, :],
             raw.first_samp,
@@ -784,6 +787,7 @@ def find_events(
             uint_cast=uint_cast,
             mask_type=mask_type,
             initial_event=initial_event,
+            ch_name=ch_name,
         )
         # add safety check for spurious events (for ex. from neuromag syst.) by
         # checking the number of low sample events

mne/tests/test_dipole.py

@@ -471,14 +471,25 @@ def _check_roundtrip_fixed(dip, tmp_path):
             assert_allclose(ch_1[key], ch_2[key], err_msg=key)
 
 
-def test_get_phantom_dipoles():
+@pytest.mark.parametrize(
+    "kind, count",
+    [
+        ("vectorview", 32),
+        ("otaniemi", 32),
+        ("oyama", 50),
+    ],
+)
+def test_get_phantom_dipoles(kind, count):
     """Test getting phantom dipole locations."""
-    pytest.raises(ValueError, get_phantom_dipoles, 0)
-    pytest.raises(ValueError, get_phantom_dipoles, "foo")
-    for kind in ("vectorview", "otaniemi"):
-        pos, ori = get_phantom_dipoles(kind)
-        assert pos.shape == (32, 3)
-        assert ori.shape == (32, 3)
+    with pytest.raises(TypeError, match="must be an instance of"):
+        get_phantom_dipoles(0)
+    with pytest.raises(ValueError, match="Invalid value for"):
+        get_phantom_dipoles("foo")
+    pos, ori = get_phantom_dipoles(kind)
+    assert pos.shape == (count, 3)
+    assert ori.shape == (count, 3)
+    # pos should be orthogonal to ori for all dipoles
+    assert_allclose(np.sum(pos * ori, axis=1), 0.0, atol=1e-7)
 
 
 @testing.requires_testing_data