ENH: Enable sensor-specific OPM coregistration in mne coreg (#11405)

Author: larsoner

.circleci/config.yml

@@ -223,13 +223,16 @@ jobs:
             - data-cache-fsaverage
       - restore_cache:
           keys:
-            - data-cache-bst-phantom-ctf
+            - data-cache-bst-raw
       - restore_cache:
           keys:
-            - data-cache-bst-raw
+            - data-cache-bst-phantom-ctf
       - restore_cache:
           keys:
             - data-cache-bst-phantom-elekta
+      - restore_cache:
+          keys:
+            - data-cache-bst-phantom-kernel
       - restore_cache:
           keys:
             - data-cache-bst-auditory
@@ -368,18 +371,22 @@ jobs:
           key: data-cache-fsaverage
           paths:
             - ~/mne_data/MNE-fsaverage-data  # (762 M)
-      - save_cache:
-          key: data-cache-bst-phantom-ctf
-          paths:
-            - ~/mne_data/MNE-brainstorm-data/bst_phantom_ctf  # (177 M)
       - save_cache:
           key: data-cache-bst-raw
           paths:
             - ~/mne_data/MNE-brainstorm-data/bst_raw  # (830 M)
+      - save_cache:
+          key: data-cache-bst-phantom-ctf
+          paths:
+            - ~/mne_data/MNE-brainstorm-data/bst_phantom_ctf  # (177 M)
       - save_cache:
           key: data-cache-bst-phantom-elekta
           paths:
             - ~/mne_data/MNE-brainstorm-data/bst_phantom_elekta  # (1.4 G)
+      - save_cache:
+          key: data-cache-bst-phantom-kernel
+          paths:
+            - ~/mne_data/MNE-phantom-kernel-data  # (362 M)
       - save_cache:
           key: data-cache-bst-auditory
           paths:

doc/api/datasets.rst

@@ -41,6 +41,7 @@ Datasets
    ucl_opm_auditory.data_path
    visual_92_categories.data_path
    phantom_4dbti.data_path
+   phantom_kernel.data_path
    refmeg_noise.data_path
    ssvep.data_path
    erp_core.data_path

doc/documentation/datasets.rst

@@ -293,6 +293,19 @@ the MEG center in La Timone hospital in Marseille.
 
     * :ref:`tut-phantom-4Dbti`
 
+Kernel OPM phantom dataset
+==========================
+:func:`mne.datasets.phantom_kernel.data_path`.
+
+This dataset was obtained with a Neuromag phantom in a Kernel Flux (720-sensor)
+system at ILABS at the University of Washington. Only 7 out of 42 possible modules
+were active for testing purposes, yielding 121 channels of data with limited coverage
+(mostly occipital and parietal).
+
+.. topic:: Examples
+
+    * :ref:`ex-kernel-opm-phantom`
+
 OPM
 ===
 :func:`mne.datasets.opm.data_path`

examples/datasets/kernel_phantom.py

@@ -0,0 +1,103 @@
+"""
+.. _ex-kernel-opm-phantom:
+
+Kernel OPM phantom data
+=======================
+
+In this dataset, a Neuromag phantom was placed inside the Kernel OPM helmet and
+stimulated with 7 modules active (121 channels). Here we show some example traces.
+"""
+
+import numpy as np
+
+import mne
+
+data_path = mne.datasets.phantom_kernel.data_path()
+fname = data_path / "phantom_32_100nam_raw.fif"
+raw = mne.io.read_raw_fif(fname).load_data()
+events = mne.find_events(raw, stim_channel="STI101")
+
+# Bads identified by inspecting averages
+raw.info["bads"] = [
+    "RC2.bx.ave",
+    "LC3.bx.ave",
+    "RC2.by.7",
+    "RC2.bz.7",
+    "RC2.bx.4",
+    "RC2.by.4",
+    "LC3.bx.5",
+]
+# Drop the module-average channels
+raw.drop_channels([ch_name for ch_name in raw.ch_names if ".ave" in ch_name])
+# Add field correction projectors
+raw.add_proj(mne.preprocessing.compute_proj_hfc(raw.info, order=2))
+raw.pick("meg", exclude="bads")
+raw.filter(0.5, 40)
+epochs = mne.Epochs(
+    raw,
+    events,
+    tmin=-0.1,
+    tmax=0.25,
+    decim=5,
+    preload=True,
+    baseline=(None, 0),
+)
+evoked = epochs["17"].average()  # a high-SNR dipole for these data
+fig = evoked.plot()
+t_peak = 0.016  # based on visual inspection of evoked
+fig.axes[0].axvline(t_peak, color="k", ls=":", lw=3, zorder=2)
+
+# %%
+# The data covariance has an interesting structure because of densely packed sensors:
+cov = mne.compute_covariance(epochs, tmax=-0.01)
+mne.viz.plot_cov(cov, raw.info)
+
+# %%
+# So let's be careful and compute rank ahead of time and regularize:
+
+rank = mne.compute_rank(epochs, tol=1e-3, tol_kind="relative")
+cov = mne.compute_covariance(epochs, tmax=-0.01, rank=rank, method="shrunk")
+mne.viz.plot_cov(cov, raw.info)
+
+# %%
+# Look at our alignment:
+
+sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)
+trans = mne.transforms.Transform("head", "mri", np.eye(4))
+align_kwargs = dict(
+    trans=trans,
+    bem=sphere,
+    surfaces={"outer_skin": 0.2},
+    show_axes=True,
+)
+mne.viz.plot_alignment(
+    raw.info,
+    coord_frame="meg",
+    meg=dict(sensors=1.0, helmet=0.05),
+    **align_kwargs,
+)
+
+# %%
+# Let's do dipole fits, which are not great because the dev_head_t is approximate and
+# the sensor coverage is sparse:
+
+data = list()
+for ii in range(1, 33):
+    evoked = epochs[str(ii)][1:-1].average().crop(t_peak, t_peak)
+    data.append(evoked.data[:, 0])
+evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.0)
+del epochs
+dip, residual = mne.fit_dipole(evoked, cov, sphere, n_jobs=None)
+actual_pos, actual_ori = mne.dipole.get_phantom_dipoles()
+actual_amp = np.ones(len(dip))  # fake amp, needed to create Dipole instance
+actual_gof = np.ones(len(dip))  # fake GOF, needed to create Dipole instance
+dip_true = mne.Dipole(dip.times, actual_pos, actual_amp, actual_ori, actual_gof)
+
+fig = mne.viz.plot_alignment(
+    evoked.info, coord_frame="head", meg="sensors", **align_kwargs
+)
+mne.viz.plot_dipole_locations(
+    dipoles=dip_true, mode="arrow", color=(0.0, 0.0, 0.0), fig=fig
+)
+mne.viz.plot_dipole_locations(dipoles=dip, mode="arrow", color=(0.2, 1.0, 0.5), fig=fig)
+mne.viz.set_3d_view(figure=fig, azimuth=30, elevation=70, distance=0.4)

mne/_fiff/_digitization.py

@@ -293,6 +293,7 @@ def _get_data_as_dict_from_dig(dig, exclude_ref_channel=True):
     # Split up the dig points by category
     hsp, hpi, elp = list(), list(), list()
     fids, dig_ch_pos_location = dict(), list()
+    dig = [] if dig is None else dig
 
     for d in dig:
         if d["kind"] == FIFF.FIFFV_POINT_CARDINAL:
@@ -307,6 +308,8 @@ def _get_data_as_dict_from_dig(dig, exclude_ref_channel=True):
                 dig_ch_pos_location.append(d["r"])
 
     dig_coord_frames = set([d["coord_frame"] for d in dig])
+    if len(dig_coord_frames) == 0:
+        dig_coord_frames = set([FIFF.FIFFV_COORD_HEAD])
     if len(dig_coord_frames) != 1:
         raise RuntimeError(
             "Only single coordinate frame in dig is supported, "

mne/bem.py

@@ -1751,13 +1751,15 @@ def _add_gamma_multipliers(bem):
     FIFF.FIFFV_BEM_SURF_ID_SKULL: "outer skull",
     FIFF.FIFFV_BEM_SURF_ID_HEAD: "outer skin ",
     FIFF.FIFFV_BEM_SURF_ID_UNKNOWN: "unknown    ",
+    FIFF.FIFFV_MNE_SURF_MEG_HELMET: "MEG helmet ",
 }
 _sm_surf_name = {
     FIFF.FIFFV_BEM_SURF_ID_BRAIN: "brain",
     FIFF.FIFFV_BEM_SURF_ID_CSF: "csf",
     FIFF.FIFFV_BEM_SURF_ID_SKULL: "outer skull",
     FIFF.FIFFV_BEM_SURF_ID_HEAD: "outer skin ",
     FIFF.FIFFV_BEM_SURF_ID_UNKNOWN: "unknown    ",
+    FIFF.FIFFV_MNE_SURF_MEG_HELMET: "helmet",
 }
 
 
@@ -1850,7 +1852,8 @@ def _write_bem_surfaces_block(fid, surfs):
     """Write bem surfaces to open file handle."""
     for surf in surfs:
         start_block(fid, FIFF.FIFFB_BEM_SURF)
-        write_float(fid, FIFF.FIFF_BEM_SIGMA, surf["sigma"])
+        if "sigma" in surf:
+            write_float(fid, FIFF.FIFF_BEM_SIGMA, surf["sigma"])
         write_int(fid, FIFF.FIFF_BEM_SURF_ID, surf["id"])
         write_int(fid, FIFF.FIFF_MNE_COORD_FRAME, surf["coord_frame"])
         write_int(fid, FIFF.FIFF_BEM_SURF_NNODE, surf["np"])

mne/channels/channels.py

@@ -104,6 +104,10 @@ def _get_meg_system(info):
                 system = "ARTEMIS123"
                 have_helmet = False
                 break
+            elif coil_type == FIFF.FIFFV_COIL_KERNEL_OPM_MAG_GEN1:
+                system = "Kernel_Flux"
+                have_helmet = True
+                break
     else:
         system = "unknown"
         have_helmet = False

mne/coreg.py

@@ -1550,7 +1550,9 @@ def _setup_bem(self):
         low_res_path = _find_head_bem(self._subject, self._subjects_dir, high_res=False)
         if high_res_path is None and low_res_path is None:
             raise RuntimeError(
-                "No standard head model was " f"found for subject {self._subject}"
+                "No standard head model was "
+                f"found for subject {self._subject} in "
+                f"{self._subjects_dir}"
             )
         if high_res_path is not None:
             self._bem_high_res = _read_surface(
@@ -1987,9 +1989,9 @@ def fit_fiducials(
         return self
 
     def _setup_icp(self, n_scale_params):
-        head_pts = list()
-        mri_pts = list()
-        weights = list()
+        head_pts = [np.zeros((0, 3))]
+        mri_pts = [np.zeros((0, 3))]
+        weights = [np.zeros(0)]
         if self._has_dig_data and self._hsp_weight > 0:  # should be true
             head_pts.append(self._filtered_extra_points)
             mri_pts.append(