Quantify motion from head position file

mne_quantify_motion.py

#!/usr/bin/env python3
# Rationale for using degrees (pitch/roll/yaw) interchangeably with mm  is "borrowed" from afni/3dvolreg. Basically, at a radius of about
# 57 mm (57.2598), circumference=360 mm, and 1 degree of rotation == 1 mm shift on the spherical surface.
# See https://afni.nimh.nih.gov/afni/community/board/read.php?1,82641,82645#msg-82645
import os
import mne
from nibabel.eulerangles import mat2euler
from mne.fixes import einsum
import numpy as np

data_path = os.path.join(mne.datasets.testing.data_path(verbose=True), 'SSS')
quats = mne.chpi.read_head_pos(os.path.join(data_path, 'test_move_anon_raw.pos'))
   
trans, rot, t = mne.chpi.head_pos_to_trans_rot_t(quats)
use_trans = einsum('ijk,ik->ij', rot[:, :3, :3].transpose([0, 2, 1]),
                   -trans) * 1000 # Also convert from meters to mm
use_rot = rot.transpose([0, 2, 1])
rotations = np.zeros(use_trans.shape)    
for i, mat2 in enumerate(use_rot):
    rads2 = mat2euler(mat2) # rotations in radians around z, y, x axes
    yaw, roll, pitch = np.degrees(rads2)
    rotations[i,:] = pitch, roll, yaw 
mot_all = np.concatenate([np.expand_dims(t,1), use_trans, rotations], axis=1) # This is for output to a file, we don't need time in the array
mot_parms = np.diff(mot_all[:, 1:7], axis=0)   
# Get motion at each measurement (enorm of the 6 motion parms)
mot_norms = np.linalg.norm(mot_parms, axis=1)   
motion_mean = np.mean(mot_norms)
motion_std = motion_sd = np.std(mot_norms, ddof=1)