Axial to planar gradiometer transformation (#13196)

Co-authored-by: Konstantinos Tsilimparis <konstantinos@artinis.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

Author: 4 people

doc/api/preprocessing.rst

@@ -49,6 +49,7 @@ Projections:
    read_dig_localite
    make_standard_montage
    read_custom_montage
+   read_meg_canonical_info
    transform_to_head
    compute_dev_head_t
    read_layout

doc/changes/dev/13196.newfeature.rst

@@ -0,0 +1 @@
+Add support for :meth:`mne.io.Raw.interpolate_to` and related methods (and :func:`mne.channels.read_meg_canonical_info`) for interpolating MEG data across systems, by :newcontrib:`Konstantinos Tsilimparis`.

doc/changes/names.inc

@@ -170,6 +170,7 @@
 .. _Katarina Slama: https://github.com/katarinaslama
 .. _Katia Al-Amir: https://github.com/katia-sentry
 .. _Keith Doelling: https://github.com/kdoelling1919
+.. _Konstantinos Tsilimparis: https://contsili.github.io/
 .. _Kostiantyn Maksymenko: https://github.com/makkostya
 .. _Kristijan Armeni: https://github.com/kristijanarmeni
 .. _Kyle Mathewson: https://github.com/kylemath

examples/preprocessing/interpolate_to.py

@@ -2,21 +2,30 @@
 .. _ex-interpolate-to-any-montage:
 
 ======================================================
-Interpolate EEG data to any montage
+Interpolate MEG or EEG data to any montage
 ======================================================
 
-This example demonstrates how to interpolate EEG channels to match a given montage.
-This can be useful for standardizing
+This example demonstrates both EEG montage interpolation and MEG system
+transformation.
+
+For EEG, this can be useful for standardizing
 EEG channel layouts across different datasets (see :footcite:`MellotEtAl2024`).
 
 - Using the field interpolation for EEG data.
 - Using the target montage "biosemi16".
+- Using the MNE interpolation for MEG data to transform from Neuromag
+  (planar gradiometers and magnetometers) to CTF (axial gradiometers).
+
 
-In this example, the data from the original EEG channels will be
+In the first example, the data from the original EEG channels will be
 interpolated onto the positions defined by the "biosemi16" montage.
+
+In the second example, we will interpolate MEG data from a 306-sensor Neuromag
+to 275-sensor CTF system.
 """
 
 # Authors: Antoine Collas <contact@antoinecollas.fr>
+#          Konstantinos Tsilimparis <konstantinos.tsilimparis@outlook.com>
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
@@ -30,6 +39,9 @@
 ylim = (-10, 10)
 
 # %%
+# Part 1: EEG System Transformation
+# ==================================
+
 # Load EEG data
 data_path = sample.data_path()
 eeg_file_path = data_path / "MEG" / "sample" / "sample_audvis-ave.fif"
@@ -75,6 +87,61 @@
 )
 axs[2].set_title("Interpolated to Standard 1020 Montage using MNE interpolation")
 
+# %%
+# Part 2: MEG System Transformation
+# ==================================
+# We demonstrate transforming MEG data from Neuromag (planar gradiometers
+# and magnetometers) to CTF (axial gradiometers) sensor configuration.
+
+# Load the full evoked data with MEG channels
+evoked_meg = mne.read_evokeds(
+    eeg_file_path, condition="Left Auditory", baseline=(None, 0)
+)
+evoked_meg.pick("meg")
+
+print("Original Neuromag system:")
+print(f"  Number of magnetometers: {len(mne.pick_types(evoked_meg.info, meg='mag'))}")
+print(f"  Number of gradiometers: {len(mne.pick_types(evoked_meg.info, meg='grad'))}")
+
+# %%
+# Transform to CTF sensor configuration
+# ======================================
+
+# Interpolate Neuromag to CTF
+evoked_ctf = evoked_meg.copy().interpolate_to("ctf275", mode="accurate")
+
+print("\nTransformed to CTF system:")
+print(f"  Number of MEG channels: {len(mne.pick_types(evoked_ctf.info, meg=True))}")
+print(f"  Bad channels in original: {evoked_meg.info['bads']}")
+
+# %%
+# Compare evoked responses: Original Neuromag vs Transformed CTF
+# The data should be similar but projected onto different sensor arrays
+
+# Set consistent y-limits for comparison
+ylim_meg = dict(grad=[-300, 300], mag=[-600, 600], meg=[-300, 300])
+
+fig, axes = plt.subplots(3, 1, figsize=(10, 8), layout="constrained")
+
+# Plot original Neuromag gradiometers
+evoked_meg.copy().pick("grad").plot(
+    axes=axes[0], show=False, spatial_colors=True, ylim=ylim_meg, time_unit="s"
+)
+axes[0].set_title("Original Neuromag Planar Gradiometers", fontsize=14)
+
+
+# Plot original Neuromag magnetometers
+evoked_meg.copy().pick("mag").plot(
+    axes=axes[1], show=False, spatial_colors=True, ylim=ylim_meg, time_unit="s"
+)
+axes[1].set_title("Original Neuromag Magnetometers", fontsize=14)
+
+# Plot transformed CTF gradiometers
+evoked_ctf.plot(
+    axes=axes[2], show=False, spatial_colors=True, ylim=ylim_meg, time_unit="s"
+)
+axes[2].set_title("Transformed to CTF275 Axial Gradiometers", fontsize=14)
+
 # %%
 # References
 # ----------

mne/channels/__init__.pyi

@@ -30,6 +30,7 @@ __all__ = [
     "read_dig_localite",
     "read_dig_polhemus_isotrak",
     "read_layout",
+    "read_meg_canonical_info",
     "read_polhemus_fastscan",
     "read_vectorview_selection",
     "rename_channels",
@@ -75,6 +76,7 @@ from .montage import (
     read_dig_hpts,
     read_dig_localite,
     read_dig_polhemus_isotrak,
+    read_meg_canonical_info,
     read_polhemus_fastscan,
     transform_to_head,
 )

mne/channels/channels.py

@@ -41,7 +41,7 @@
     pick_info,
     pick_types,
 )
-from .._fiff.proj import _has_eeg_average_ref_proj, setup_proj
+from .._fiff.proj import setup_proj
 from .._fiff.reference import add_reference_channels, set_eeg_reference
 from .._fiff.tag import _rename_list
 from ..bem import _check_origin
@@ -61,6 +61,8 @@
     warn,
 )
 
+_ALLOWED_INTERPOLATION_MODES = ("accurate", "fast")
+
 
 def _get_meg_system(info):
     """Educated guess for the helmet type based on channels."""
@@ -968,161 +970,102 @@ def interpolate_bads(
 
         return self
 
-    def interpolate_to(self, sensors, origin="auto", method="spline", reg=0.0):
-        """Interpolate EEG data onto a new montage.
+    def interpolate_to(
+        self, sensors, origin="auto", method=None, mode="accurate", reg=0.0
+    ):
+        """Interpolate data onto a new sensor configuration.
 
-        .. warning::
-            Be careful, only EEG channels are interpolated. Other channel types are
-            not interpolated.
+        This method can interpolate EEG data onto a new montage or transform
+        MEG data to a different sensor configuration (e.g., Neuromag to CTF).
 
         Parameters
         ----------
-        sensors : DigMontage
-            The target montage containing channel positions to interpolate onto.
+        sensors : DigMontage | str
+            For EEG: A DigMontage object containing target channel positions.
+            For MEG: A string specifying the target MEG system. Currently
+            supported: ``'neuromag'``, ``'ctf151'`` or ``'ctf275'``.
         origin : array-like, shape (3,) | str
             Origin of the sphere in the head coordinate frame and in meters.
             Can be ``'auto'`` (default), which means a head-digitization-based
-            origin fit.
-        method : str
-            Method to use for EEG channels.
-            Supported methods are 'spline' (default) and 'MNE'.
+            origin fit. Used for both EEG and MEG interpolation.
+        method : str | None
+            Interpolation method to use.
+            For EEG: ``'spline'`` (default, same as None) or ``'MNE'``.
+            For MEG:  ``'MNE'`` (default, same as None).
+
+            .. versionchanged:: 1.10.0
+               Added support for MEG interpolation.
+        mode : str
+            Either ``'accurate'`` (default) or ``'fast'``, determines the
+            quality of the Legendre polynomial expansion used for
+            interpolation of MEG channels using the minimum-norm method.
+            Only used for MEG interpolation.
         reg : float
-            The regularization parameter for the interpolation method
-            (only used when the method is 'spline').
+            The regularization parameter for the interpolation method.
+            Only used when ``method='spline'`` for EEG channels.
 
         Returns
         -------
         inst : instance of Raw, Epochs, or Evoked
-            The instance with updated channel locations and data.
+            A new instance with interpolated data and updated channel
+            information.
 
         Notes
         -----
-        This method is useful for standardizing EEG layouts across datasets.
-        However, some attributes may be lost after interpolation.
+        **For EEG data:**
 
-        .. versionadded:: 1.10.0
-        """
-        from ..epochs import BaseEpochs, EpochsArray
-        from ..evoked import Evoked, EvokedArray
-        from ..forward._field_interpolation import _map_meg_or_eeg_channels
-        from ..io import RawArray
-        from ..io.base import BaseRaw
-        from .interpolation import _make_interpolation_matrix
-        from .montage import DigMontage
+        This method interpolates EEG channels onto a new montage using
+        spherical splines or minimum-norm estimation. Non-EEG channels
+        are preserved without modification.
 
-        # Check that the method option is valid.
-        _check_option("method", method, ["spline", "MNE"])
-        _validate_type(sensors, DigMontage, "sensors")
+        **For MEG data:**
 
-        # Get target positions from the montage
-        ch_pos = sensors.get_positions().get("ch_pos", {})
-        target_ch_names = list(ch_pos.keys())
-        if not target_ch_names:
-            raise ValueError(
-                "The provided sensors configuration has no channel positions."
-            )
+        This method transforms MEG data to a different sensor configuration
+        using field interpolation.
 
-        # Get original channel order
-        orig_names = self.info["ch_names"]
-
-        # Identify EEG channel
-        picks_good_eeg = pick_types(self.info, meg=False, eeg=True, exclude="bads")
-        if len(picks_good_eeg) == 0:
-            raise ValueError("No good EEG channels available for interpolation.")
-        # Also get the full list of EEG channel indices (including bad channels)
-        picks_remove_eeg = pick_types(self.info, meg=False, eeg=True, exclude=[])
-        eeg_names_orig = [orig_names[i] for i in picks_remove_eeg]
-
-        # Identify non-EEG channels in original order
-        non_eeg_names_ordered = [ch for ch in orig_names if ch not in eeg_names_orig]
-
-        # Create destination info for new EEG channels
-        sfreq = self.info["sfreq"]
-        info_interp = create_info(
-            ch_names=target_ch_names,
-            sfreq=sfreq,
-            ch_types=["eeg"] * len(target_ch_names),
-        )
-        info_interp.set_montage(sensors)
-        info_interp["bads"] = [ch for ch in self.info["bads"] if ch in target_ch_names]
-        # Do not assign "projs" directly.
-
-        # Compute the interpolation mapping
-        if method == "spline":
-            origin_val = _check_origin(origin, self.info)
-            pos_from = self.info._get_channel_positions(picks_good_eeg) - origin_val
-            pos_to = np.stack(list(ch_pos.values()), axis=0)
-
-            def _check_pos_sphere(pos):
-                d = np.linalg.norm(pos, axis=-1)
-                d_norm = np.mean(d / np.mean(d))
-                if np.abs(1.0 - d_norm) > 0.1:
-                    warn("Your spherical fit is poor; interpolation may be inaccurate.")
-
-            _check_pos_sphere(pos_from)
-            _check_pos_sphere(pos_to)
-            mapping = _make_interpolation_matrix(pos_from, pos_to, alpha=reg)
+        Common use cases for MEG transformation:
 
-        else:
-            assert method == "MNE"
-            info_eeg = pick_info(self.info, picks_good_eeg)
-            # If the original info has an average EEG reference projector but
-            # the destination info does not,
-            # update info_interp via a temporary RawArray.
-            if _has_eeg_average_ref_proj(self.info) and not _has_eeg_average_ref_proj(
-                info_interp
-            ):
-                # Create dummy data: shape (n_channels, 1)
-                temp_data = np.zeros((len(info_interp["ch_names"]), 1))
-                temp_raw = RawArray(temp_data, info_interp, first_samp=0)
-                # Using the public API, add an average reference projector.
-                temp_raw.set_eeg_reference(
-                    ref_channels="average", projection=True, verbose=False
-                )
-                # Extract the updated info.
-                info_interp = temp_raw.info
-            mapping = _map_meg_or_eeg_channels(
-                info_eeg, info_interp, mode="accurate", origin=origin
-            )
+        - Transform Neuromag data to CTF sensor layout for comparison
+        - Transform CTF data to Neuromag sensor layout
+        - Simulate what data would look like on a different MEG system
 
-        # Interpolate EEG data
-        data_good = self.get_data(picks=picks_good_eeg)
-        data_interp = mapping @ data_good
+        .. warning::
+            MEG field interpolation assumes that the head position relative
+            to the sensors is similar between systems. Large differences in
+            head position may affect interpolation accuracy.
 
-        # Create a new instance for the interpolated EEG channels
-        # TODO: Creating a new instance leads to a loss of information.
-        #       We should consider updating the existing instance in the future
-        #       by 1) drop channels, 2) add channels, 3) re-order channels.
-        if isinstance(self, BaseRaw):
-            inst_interp = RawArray(data_interp, info_interp, first_samp=self.first_samp)
-        elif isinstance(self, BaseEpochs):
-            inst_interp = EpochsArray(data_interp, info_interp)
-        else:
-            assert isinstance(self, Evoked)
-            inst_interp = EvokedArray(data_interp, info_interp)
-
-        # Merge only if non-EEG channels exist
-        if not non_eeg_names_ordered:
-            return inst_interp
-
-        inst_non_eeg = self.copy().pick(non_eeg_names_ordered).load_data()
-        inst_out = inst_non_eeg.add_channels([inst_interp], force_update_info=True)
-
-        # Reorder channels
-        # Insert the entire new EEG block at the position of the first EEG channel.
-        orig_names_arr = np.array(orig_names)
-        mask_eeg = np.isin(orig_names_arr, eeg_names_orig)
-        if mask_eeg.any():
-            first_eeg_index = np.where(mask_eeg)[0][0]
-            pre = orig_names_arr[:first_eeg_index]
-            new_eeg = np.array(info_interp["ch_names"])
-            post = orig_names_arr[first_eeg_index:]
-            post = post[~np.isin(orig_names_arr[first_eeg_index:], eeg_names_orig)]
-            new_order = np.concatenate((pre, new_eeg, post)).tolist()
+        .. versionadded:: 1.10.0
+        .. versionchanged:: 1.12.0
+           Added support for MEG interpolation to canonical systems.
+        """
+        from .interpolation import _interpolate_to_eeg, _interpolate_to_meg
+        from .montage import DigMontage
+
+        _check_preload(self, "interpolation")
+
+        # Determine if we're doing EEG or MEG interpolation
+        _validate_type(sensors, (str, DigMontage), "sensors")
+        _validate_type(method, (str, None), "method")
+        is_meg_interpolation = isinstance(sensors, str)
+        is_eeg_interpolation = isinstance(sensors, DigMontage)
+
+        if is_eeg_interpolation:
+            valid_methods = ["spline", "MNE"]
+            func = partial(_interpolate_to_eeg, method=method, reg=reg)
+            kind = "eeg"
         else:
-            new_order = orig_names
-        inst_out.reorder_channels(new_order)
-        return inst_out
+            assert is_meg_interpolation
+            valid_methods = ["MNE"]
+            _check_option("sensors", sensors, ["neuromag", "ctf151", "ctf275"])
+            func = partial(_interpolate_to_meg, mode=mode)
+            kind = "meg"
+
+        if method is None:
+            method = _handle_default("interpolation_method")[kind]
+        _check_option("mode", mode, _ALLOWED_INTERPOLATION_MODES)
+        extra = f"when doing {kind.upper()} interpolation"
+        _check_option("method", method, valid_methods, extra=extra)
+        return func(self, sensors, origin)
 
 
 @verbose