ENH: Add movement compensation (mne-tools#747)

Author: larsoner

.circleci/config.yml

@@ -48,7 +48,7 @@ jobs:
           name: Get Python running
           command: |
             pip install --upgrade --progress-bar off pip setuptools
-            pip install --upgrade --progress-bar off "autoreject @ https://api.github.com/repos/autoreject/autoreject/zipball/master" "mne[hdf5] @ git+https://github.com/mne-tools/mne-python" "mne-bids[full] @ https://api.github.com/repos/mne-tools/mne-bids/zipball/main"
+            pip install --upgrade --progress-bar off "autoreject @ https://api.github.com/repos/autoreject/autoreject/zipball/master" "mne[hdf5] @ git+https://github.com/mne-tools/mne-python" "mne-bids[full] @ https://api.github.com/repos/mne-tools/mne-bids/zipball/main" numba
             pip install -ve .[tests]
             pip install PyQt6
       - run:
@@ -402,6 +402,7 @@ jobs:
 
   test_ds004107:
     <<: *imageconfig
+    resource_class: large  # memory
     steps:
       - attach_workspace:
           at: ~/
@@ -768,6 +769,7 @@ jobs:
 
   test_ds004229:
     <<: *imageconfig
+    resource_class: large  # head position estimation
     steps:
       - attach_workspace:
           at: ~/

docs/source/settings/preprocessing/maxfilter.md

@@ -10,7 +10,16 @@ tags:
       members:
         - use_maxwell_filter
         - mf_st_duration
+        - mf_st_correlation
         - mf_head_origin
+        - mf_destination
+        - mf_int_order
         - mf_reference_run
         - mf_cal_fname
         - mf_ctc_fname
+        - mf_mc
+        - mf_mc_t_step_min
+        - mf_mc_t_window
+        - mf_mc_gof_limit
+        - mf_mc_dist_limit
+        - mf_filter_chpi

docs/source/settings/preprocessing/ssp_ica.md

@@ -21,6 +21,7 @@ tags:
         - eog_proj_from_average
         - ssp_reject_eog
         - ssp_reject_ecg
+        - ssp_ecg_channel
         - ica_reject
         - ica_algorithm
         - ica_l_freq

docs/source/v1.4.md.inc

@@ -1,8 +1,9 @@
 ## v1.4.0 (unreleased)
 
-[//]: # (### :new: New features & enhancements)
+### :new: New features & enhancements
 
-[//]: # (- Whatever (#000 by @whoever))
+- Add movement compensation and cHPI filtering to the Maxwell filtering step, along with additional configuration options (#747 by @larsoner)
+- Add option to specify `ssp_ecg_channel` to override the default value (#747 by @larsoner)
 
 [//]: # (### :warning: Behavior changes)
 
@@ -17,3 +18,4 @@
 - Fix bug when `mf_reference_run != runs[0]` (#742 by @larsoner)
 - Fix bug with too many JSON files found during empty room matching (#743 by @allermat)
 - Fix bug with outdated info on ch_types config option (#745 by @allermat)
+- Fix bug where SSP projectors were not added to the report (#747 by @larsoner)

mne_bids_pipeline/_config.py

@@ -576,6 +576,16 @@
     ```
 """
 
+mf_st_correlation: float = 0.98
+"""
+The correlation limit for spatio-temporal SSS (tSSS).
+
+???+ example "Example"
+    ```python
+    st_correlation = 0.98
+    ```
+"""
+
 mf_head_origin: Union[Literal["auto"], ArrayLike] = "auto"
 """
 `mf_head_origin` : array-like, shape (3,) | 'auto'
@@ -591,17 +601,42 @@
     ```
 """
 
-mf_reference_run: Optional[str] = None
+mf_destination: Union[Literal["reference_run"], ArrayLike] = "reference_run"
 """
 Despite all possible care to avoid movements in the MEG, the participant
 will likely slowly drift down from the Dewar or slightly shift the head
 around in the course of the recording session. Hence, to take this into
-account, we are realigning all data to a single position. For this, you need
-to define a reference run (typically the one in the middle of
-the recording session).
+account, we are realigning all data to a single position. For this, you can:
 
+1. Choose a reference run. Often one from the middle of the recording session
+   is a good choice. Set `mf_destination = "reference_run" and then set
+   [`config.mf_reference_run`](mne_bids_pipeline._config.mf_reference_run).
+   This will result in a device-to-head transformation that differs between
+   subjects.
+2. Choose a standard position in the MEG coordinate frame. For this, pass
+   a 4x4 transformation matrix for the device-to-head
+   transform. This will result in a device-to-head transformation that is
+   the same across all subjects.
+
+   ???+ example "A Standardized Position"
+   ```python
+   from mne.transforms import translation
+   mf_destination = translation(z=0.04)
+   ```
+"""
+
+mf_int_order: int = 8
+"""
+Internal order for the Maxwell basis. Can be set to something lower (e.g., 6
+or higher for datasets where lower or higher spatial complexity, respectively,
+is expected.
+"""
+
+mf_reference_run: Optional[str] = None
+"""
 Which run to take as the reference for adjusting the head position of all
-runs. If `None`, pick the first run.
+runs when [`mf_destination="reference_run"`](mne_bids_pipeline._config.mf_destination).
+If `None`, pick the first run.
 
 ???+ example "Example"
     ```python
@@ -638,6 +673,39 @@
     ```
 """  # noqa : E501
 
+mf_mc: bool = False
+"""
+If True, perform movement compensation on the data.
+"""
+
+mf_mc_t_step_min: float = 0.01
+"""
+Minimum time step to use during cHPI coil amplitude estimation.
+"""
+
+mf_mc_t_window: Union[float, Literal["auto"]] = "auto"
+"""
+The window to use during cHPI coil amplitude estimation and in cHPI filtering.
+Can be "auto" to autodetect a reasonable value or a float (in seconds).
+"""
+
+mf_mc_gof_limit: float = 0.98
+"""
+Minimum goodness of fit to accept for each cHPI coil.
+"""
+
+mf_mc_dist_limit: float = 0.005
+"""
+Minimum distance (m) to accept for cHPI position fitting.
+"""
+
+mf_filter_chpi: Optional[bool] = None
+"""
+Use mne.chpi.filter_chpi after Maxwell filtering. Can be None to use
+the same value as [`mf_mc`][mne_bids_pipeline._config.mf_mc].
+Only used when [`use_maxwell_filter=True`][mne_bids_pipeline._config.use_maxwell_filter]
+"""  # noqa: E501
+
 ###############################################################################
 # STIMULATION ARTIFACT
 # --------------------
@@ -1125,6 +1193,12 @@
     ```
 """
 
+ssp_ecg_channel: Optional[str] = None
+"""
+Channel to use for ECG SSP. Can be useful when the autodetected ECG channel
+is not reliable.
+"""
+
 # Rejection based on ICA
 # ~~~~~~~~~~~~~~~~~~~~~~
 ica_reject: Optional[Dict[str, float]] = None

mne_bids_pipeline/_config_import.py

@@ -8,6 +8,7 @@
 from typing import Optional, List
 
 import matplotlib
+import numpy as np
 import mne
 from mne.utils import _check_option, _validate_type
 
@@ -343,6 +344,25 @@ def _check_config(config: SimpleNamespace) -> None:
         ("raise", "warn", "ignore"),
     )
 
+    _validate_type(
+        config.mf_destination,
+        (str, list, tuple, np.ndarray),
+        "config.mf_destination",
+    )
+    if isinstance(config.mf_destination, str):
+        _check_option(
+            "config.mf_destination",
+            config.mf_destination,
+            ("reference_run",),
+        )
+    else:
+        destination = np.array(config.mf_destination, float)
+        if destination.shape != (4, 4):
+            raise ValueError(
+                "config.mf_destination, if array-like, must have shape (4, 4) "
+                f"but got shape {destination.shape}"
+            )
+
 
 _REMOVED_NAMES = {
     "debug": dict(

mne_bids_pipeline/_config_utils.py

@@ -170,7 +170,10 @@ def get_intersect_run(config: SimpleNamespace) -> List[str]:
 
 
 def get_runs(
-    *, config: SimpleNamespace, subject: str, verbose: bool = False
+    *,
+    config: SimpleNamespace,
+    subject: str,
+    verbose: bool = False,
 ) -> Union[List[str], List[None]]:
     """Returns a list of runs in the BIDS input data.
 
@@ -218,17 +221,43 @@ def get_runs(
             logger.info(**gen_log_kwargs(message=msg))
 
     if runs == "all":
-        runs = valid_runs
+        runs = list(valid_runs)
 
     if not runs:
-        return [None]
+        runs = [None]
     else:
         inclusion = set(runs).issubset(set(valid_runs))
         if not inclusion:
             raise ValueError(
                 f"Invalid run. It can be a subset of {valid_runs} but " f"got {runs}"
             )
-        return runs
+    return runs
+
+
+def get_runs_tasks(
+    *,
+    config: SimpleNamespace,
+    subject: str,
+    session: Optional[str],
+) -> List[Tuple[str]]:
+    """Get (run, task) tuples for all runs plus (maybe) rest."""
+    from ._import_data import _get_bids_path_in
+
+    runs = get_runs(config=config, subject=subject)
+    tasks = [config.task] * len(runs)
+    if config.process_rest and not config.task_is_rest:
+        run, task = None, "rest"
+        bids_path_in = _get_bids_path_in(
+            cfg=config,
+            subject=subject,
+            session=session,
+            run=run,
+            task=task,
+        )
+        if bids_path_in.fpath.exists():
+            runs.append(None)
+            tasks.append("rest")
+    return tuple(zip(runs, tasks))
 
 
 def get_mf_reference_run(config: SimpleNamespace) -> str:

mne_bids_pipeline/_import_data.py

@@ -9,6 +9,7 @@
 from ._config_utils import (
     get_mf_reference_run,
     get_runs,
+    get_datatype,
     _bids_kwargs,
 )
 from ._io import _read_json, _empty_room_match_path
@@ -508,29 +509,29 @@ def _find_breaks_func(
     raw.set_annotations(raw.annotations + break_annots)  # add to existing
 
 
-def _get_raw_paths(
+def _get_bids_path_in(
     *,
     cfg: SimpleNamespace,
     subject: str,
     session: Optional[str],
     run: Optional[str],
-    kind: Literal["raw", "sss"],
-    add_bads: bool = True,
-    include_mf_ref: bool = True,
-) -> dict:
-    # Construct the basenames of the files we wish to load, and of the empty-
-    # room recording we wish to save.
-    # The basenames of the empty-room recording output file does not contain
-    # the "run" entity.
+    task: Optional[str],
+    kind: Literal["orig", "sss"] = "orig",
+) -> BIDSPath:
+    # b/c can be used before this is updated
+    if hasattr(cfg, "datatype"):
+        datatype = cfg.datatype
+    else:
+        datatype = get_datatype(config=cfg)
     path_kwargs = dict(
         subject=subject,
         run=run,
         session=session,
-        task=cfg.task,
+        task=task or cfg.task,
         acquisition=cfg.acq,
         recording=cfg.rec,
         space=cfg.space,
-        datatype=cfg.datatype,
+        datatype=datatype,
         check=False,
     )
     if kind == "sss":
@@ -545,7 +546,31 @@ def _get_raw_paths(
         path_kwargs["extension"] = None
         path_kwargs["processing"] = cfg.proc
     bids_path_in = BIDSPath(**path_kwargs)
+    return bids_path_in
 
+
+def _get_raw_paths(
+    *,
+    cfg: SimpleNamespace,
+    subject: str,
+    session: Optional[str],
+    run: Optional[str],
+    kind: Literal["orig", "sss"],
+    add_bads: bool = True,
+    include_mf_ref: bool = True,
+) -> dict:
+    # Construct the basenames of the files we wish to load, and of the empty-
+    # room recording we wish to save.
+    # The basenames of the empty-room recording output file does not contain
+    # the "run" entity.
+    bids_path_in = _get_bids_path_in(
+        cfg=cfg,
+        subject=subject,
+        session=session,
+        run=run,
+        task=None,
+        kind=kind,
+    )
     in_files = dict()
     key = f"raw_run-{run}"
     in_files[key] = bids_path_in
@@ -556,8 +581,6 @@ def _get_raw_paths(
             in_files=in_files,
             key=key,
         )
-    orig_key = key
-
     if run == cfg.runs[0]:
         do = dict(
             rest=cfg.process_rest and not cfg.task_is_rest,
@@ -600,9 +623,7 @@ def _get_raw_paths(
                 )
             if include_mf_ref and task == "noise":
                 key = "raw_ref_run"
-                in_files[key] = (
-                    in_files[orig_key].copy().update(run=cfg.mf_reference_run)
-                )
+                in_files[key] = bids_path_in.copy().update(run=cfg.mf_reference_run)
                 _update_for_splits(in_files, key, single=True, allow_missing=True)
                 if not in_files[key].fpath.exists():
                     in_files.pop(key)
@@ -612,7 +633,6 @@ def _get_raw_paths(
                         in_files=in_files,
                         key=key,
                     )
-
     return in_files