ENH: Add annotation of head velocity (#757)

docs/source/settings/preprocessing/maxfilter.md

@@ -22,4 +22,6 @@ tags:
         - mf_mc_t_window
         - mf_mc_gof_limit
         - mf_mc_dist_limit
+        - mf_mc_rotation_velocity_limit
+        - mf_mc_translation_velocity_limit
         - mf_filter_chpi

docs/source/v1.5.md.inc

@@ -1,8 +1,8 @@
 ## v1.5.0 (unreleased)
 
-[//]: # (### :new: New features & enhancements)
+### :new: New features & enhancements
 
-[//]: # (- Whatever (#000 by @whoever))
+- Added support for annotating bad segments based on head movement velocity (#757 by @larsoner)
 
 [//]: # (### :warning: Behavior changes)
 

mne_bids_pipeline/_config.py

@@ -707,6 +707,18 @@
 Minimum distance (m) to accept for cHPI position fitting.
 """
 
+mf_mc_rotation_velocity_limit: Optional[float] = None
+"""
+The rotation velocity limit (degrees/second) to use when annotating
+movement-compensated data. If `None`, no annotations will be added.
+"""
+
+mf_mc_translation_velocity_limit: Optional[float] = None
+"""
+The translation velocity limit (meters/second) to use when annotating
+movement-compensated data. If `None`, no annotations will be added.
+"""
+
 mf_filter_chpi: Optional[bool] = None
 """
 Use mne.chpi.filter_chpi after Maxwell filtering. Can be None to use

mne_bids_pipeline/_report.py

@@ -1198,6 +1198,7 @@ def _add_raw(
     title: str,
     tags: tuple = (),
     raw: Optional[BaseRaw] = None,
+    extra_html: Optional[str] = None,
 ):
     if bids_path_in.run is not None:
         title += f", run {repr(bids_path_in.run)}"
@@ -1208,13 +1209,22 @@ def _add_raw(
         or bids_path_in.run in cfg.plot_psd_for_runs
         or bids_path_in.task in cfg.plot_psd_for_runs
     )
+    tags = ("raw", f"run-{bids_path_in.run}") + tags
     with mne.use_log_level("error"):
         report.add_raw(
             raw=raw or bids_path_in,
             title=title,
             butterfly=5,
             psd=plot_raw_psd,
-            tags=("raw", f"run-{bids_path_in.run}") + tags,
+            tags=tags,
             # caption=bids_path_in.basename,  # TODO upstream
             replace=True,
         )
+        if extra_html is not None:
+            report.add_html(
+                extra_html,
+                title=title,
+                tags=tags,
+                section=title,
+                replace=True,
+            )

mne_bids_pipeline/steps/preprocessing/_03_maxfilter.py

@@ -282,6 +282,42 @@ def run_maxwell_filter(
             t_window=cfg.mf_mc_t_window,
         )
 
+    if cfg.mf_mc and (
+        cfg.mf_mc_rotation_velocity_limit is not None
+        or cfg.mf_mc_translation_velocity_limit is not None
+    ):
+        movement_annot, _ = mne.preprocessing.annotate_movement(
+            raw_sss,
+            pos=head_pos,
+            rotation_velocity_limit=cfg.mf_mc_rotation_velocity_limit,
+            translation_velocity_limit=cfg.mf_mc_translation_velocity_limit,
+        )
+        perc_time = 100 / raw_sss.times[-1]
+        extra_html = list()
+        for kind, unit in (("translation", "m"), ("rotation", "°")):
+            limit = getattr(cfg, f"mf_mc_{kind}_velocity_limit")
+            if limit is None:
+                continue
+            desc = (f"BAD_mov_{kind[:5]}_vel",)
+            tot_time = np.sum(
+                movement_annot.duration[movement_annot.description == desc]
+            )
+            perc = perc_time * tot_time
+            logger_meth = logger.warning if perc > 20 else logger.info
+            msg = (
+                f"{kind.capitalize()} velocity exceeded {limit} {unit}/s "
+                f"limit for {tot_time:0.1f} s ({perc:0.1f}%)"
+            )
+            logger_meth(**gen_log_kwargs(message=msg))
+            extra_html.append(f"<li>{msg}</li>")
+        extra_html = (
+            "<p>The raw data were annotated with the following movement-related bad "
+            f"segment annotations:<ul>{''.join(extra_html)}</ul></p>"
+        )
+        raw_sss.set_annotations(raw_sss.annotations + movement_annot)
+    else:
+        movement_annot = extra_html = None
+
     out_files["sss_raw"] = bids_path_out
     msg = f"Writing {out_files['sss_raw'].fpath.relative_to(cfg.deriv_root)}"
     logger.info(**gen_log_kwargs(message=msg))
@@ -307,13 +343,15 @@ def run_maxwell_filter(
     ) as report:
         msg = "Adding Maxwell filtered raw data to report."
         logger.info(**gen_log_kwargs(message=msg))
+
         _add_raw(
             cfg=cfg,
             report=report,
             bids_path_in=out_files["sss_raw"],
             title="Raw (maxwell filtered)",
             tags=("sss",),
             raw=raw_sss,
+            extra_html=extra_html,
         )
 
     assert len(in_files) == 0, in_files.keys()
@@ -345,6 +383,8 @@ def get_config(
         mf_destination=config.mf_destination,
         mf_int_order=config.mf_int_order,
         mf_mc_t_window=config.mf_mc_t_window,
+        mf_mc_rotation_velocity_limit=config.mf_mc_rotation_velocity_limit,
+        mf_mc_translation_velocity_limit=config.mf_mc_translation_velocity_limit,
         **_import_data_kwargs(config=config, subject=subject),
     )
     return cfg

mne_bids_pipeline/tests/configs/config_ds004229.py

@@ -27,6 +27,8 @@
 mf_mc_t_step_min = 0.5  # just for speed!
 mf_mc_t_window = 0.2  # cleaner cHPI filtering on this dataset
 mf_filter_chpi = False  # for speed, not needed as we low-pass anyway
+mf_mc_rotation_velocity_limit = 30.0  # deg/s for annotations
+mf_mc_translation_velocity_limit = 20e-3  # m/s
 ch_types = ["meg"]
 
 l_freq = None