search

mne-tools / mne-python

MNE: Magnetoencephalography (MEG) and Electroencephalography (EEG) in Python
https://mne.tools
BSD 3-Clause "New" or "Revised" License
2.71k stars 1.32k forks source link

ENH: Measure head motion from cHPI #6607

Open jhouck opened 5 years ago

jhouck commented 5 years ago

It might be helpful to know how much head motion there was in a scan when cHPI was used, either for QA/reporting or to help decide whether to use movement compensation in maxfilter/maxwell_filter.

If we do something like what happens for fMRI, the steps would be (roughly) to 1) use chpi.read_head_pos or chpi._calculate_head positions to get the quaternions, 2) convert the translations to mm (as in viz.plot_head_positions), 3) convert the rotations to radians to degrees*, 4) concatenate translations and rotations into an array, 5) do a diff, 6) take the euclidean norm at each measurement (e.g., second), and then 7) take the mean and standard deviation of the norms to get motion per second.

*The logic of treating degrees of rotation as mm is borrowed from afni, see e.g. https://afni.nimh.nih.gov/afni/community/board/read.php?1,82641,82645#msg-82645 .

And I made a gist: https://gist.github.com/jhouck/b0e93a60f17c708c5a342dce354551e2

It might make sense to use some standard time with chpi._calculate_head positions, like setting both t_step_min and t_step_max to 1.0, but maybe that would be left to the user.

Thoughts?

larsoner commented 5 years ago

One thing I don't understand is what you mean by "rotations to radians to degrees". Rotations occur about an axis, or you can say there is some rotation angle between two rigid affines (/ quaternions). So is the idea that you would just compute the angle between veridical "up" / RAS on every time point? The problem is that this does not distinguish the angle of rotation between two positions -- you can be 45 degrees off from +Z as "up" in a lot of different orientations (opposites sides of which differ from one another by 90 degrees).

If you are just going to compute "angle" and then take a diff, it seems like what you really want to do is compute the angle difference between consecutive positions, which you can just do directly (we have an _angle_between_quats function that does it for example).

jhouck commented 5 years ago

I thought the Euler angles had a known orientation. "Rotations to radians to degrees" means I used head_pos_to_trans_rot_t to get the translation and rotation, then converted the rotation matrix to Euler angles, then converted the Euler angles from radians to degrees, before doing the diff on the combined translation (mm) and rotation (degrees that we pretend are mm) array.

Similar steps could be done through _angle_between_quats but that operates on quaternions, and after head_pos_to_trans_rot_t I've already got the rotation matrices for each time point. If I swap in rotation_angles for mat2euler I get similar, but not identical, results, so I suspect I've done something incorrect along the way.

For what it's worth I started this a couple of years ago (back when chpi was still part of io) and did not keep very good notes on the various decision points.

larsoner commented 5 years ago

I thought the Euler angles had a known orientation. "Rotations to radians to degrees" means I used head_pos_to_trans_rot_t to get the translation and rotation, then converted the rotation matrix to Euler angles, then converted the Euler angles from radians to degrees, before doing the diff on the combined translation (mm) and rotation (degrees that we pretend are mm) array.

Ahh yes Euler angles would work. There is a problem with degeneracy/gimbal lock but in practice I doubt it matters.

Similar steps could be done through _angle_between_quats but that operates on quaternions,

What you get from read_head_pos is time (:1), quaternions (1:4), translation (4:7), ..., so you actually have fewer operations/conversions if stick with quaternions and _angle_between_quats rather than convert quat->affine->euler then take a diff. Plus that's three diffs along the "euler angle" coordinate axes rather than a single one having to do with the total angle moved. The equivalent thing for position would be to compute the distance between the points rather than the distance along the three coordinate axes.

These are probably a bit nicer because these measures (one angle between the two orientations, one distance between the two translations) do not depend at all on your reference frame. In the other case, depending on what you do with the deltas along each dimension (rotx, roty, rotz, dx, dy, dz) your definition of X/Y/Z directions could matter. And I'm not sure there is much to be gained from keeping the terms separate like this.

jhouck commented 5 years ago

These are probably a bit nicer because these measures (one angle between the two orientations, one distance between the two translations) do not depend at all on your reference frame. In the other case, depending on what you do with the deltas along each dimension (rotx, roty, rotz, dx, dy, dz) your definition of X/Y/Z directions could matter. And I'm not sure there is much to be gained from keeping the terms separate like this.

Offhand, the only advantages of keeping the terms separate might be 1) it would make it easier to use the existing plotting code and 2) it gives us similar motion output to fMRI packages. The importance of either of these is unclear.

I'll be back in the office on Wednesday and can take another stab at it then.

larsoner commented 4 years ago

@jhouck we have some of these calculations done and used in annotate_movement now.

Maybe a simple and cool thing to do would be to add to plot_head_movement the translational and rotational velocities? They could even be on the same axes in a different color with a second (right) yaxis.