Time-varying position / head-movement correction would be the next natural step after #10780 I think.

IIRC the procedure is/could be here:

1. Extract time-varying rigid-body position + rotation of the system, typically using optical tracking: new mne.chpi.read_head_pos_optitrack, or mne.chpi.calculate_head_pos_optitrack, depending on the complexity of the problem (if optitrack directly gives pos+rot, or we actually need to do some calculation). It's a bit of a misnomer, but mne.chpi is where all our head-position related stuff lives, so the optitrack outputting pos+rot as a function of time is natural there. We could put it in a new mne.optitrack, or mne.preprocessing.opm or something if we want, but optical tracking seems pretty generic.
2. Clean up position+rotation parameters: new mne.chpi.clean_head_positions or so? Or maybe this should be built into calculate_head_pos_optitrack. TBD.
3. Add head positions to the raw data (as new channels) -- this is what maxwell_filter does, so it's natural for us to do it as well -- maybe mne.chpi.add_head_positions(raw, pos) and you can choose how to interpolate. (Under the hood this probably boils down to a mne.add_channels((raw, RawArray(...))) or so but we'll need to interpolate/upsample the positions+rotations, and setting the info correctly will be annoying -- so not totally trivial, thus worth wrapping.)
4. Regressing can then be done immediately with regress_artifact I think. OPM-specific options we need, if any, should hopefully be useful for other types of regression and thus reasonably implemented there.

@neurofractal did you ever play around/try regressing quaternions rather than Euler angles for rotation angle (by my bad memory you use Euler for regression)? We can convert back and forth but the .pos files we use for position + rotation information in space store in quaternions, so if those work just as well or better than Euler angles it's a slightly more natural fit here.