Derivation of motion parameters from HMC transforms

[effigies]

FSL motion parameters are not straightforwardly related to the transform matrices that they generate. They are calculated relative to the center of mass of the reference image, rather than the reference frame, which is what the affines encode.

This is the straightforward way to extract extract motion parameters from the affines:

import nitransforms as nt
import pandas as pd
from scipy.spatial.transform import Rotation

affines = nt.linear.load(motion_xfm)
trans = affines.matrix[:, :3, 3]
rot = Rotation.from_matrix(affine.matrix[:, :3, :3]).as_euler('XYZ')

params = pd.DataFrame(data=np.concatenate((trans, rot), axis=1), columns=['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z'])

Here are the original FSL (left) and affine-derived (right) parameters (rotations * 50 for comparability):

Here is the correlation matrix with the FSL parameters:

The correlations are reasonably good (>0.8 for all but trans_x).

Some code for generating similar plots:

from matplotlib import pyplot as plt
import seaborn as sns

fslparams = pd.read_csv(timeseries_tsv, delimiter='\t')[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']]

corr = pd.concat([fslparams, params], axis=1).corr().iloc[:6, 6:]
cmap = sns.diverging_palette(230, 20, as_cmap=True)
sns.heatmap(corr, cmap=cmap, vmax=1, annot=True)
plt.show()

ax1 = plt.subplot(2, 1, 1)
ax2 = plt.subplot(2, 1, 2)
sns.lineplot(fslparams * [1, 1, 1, 50, 50, 50], ax=ax1)
sns.lineplot(params * [1, 1, 1, 50, 50, 50], ax=ax2)
plt.show()

Related: #12

[effigies]

Okay, dug around the FSL source code and came up with the following:

import numpy as np
import nibabel as nb
import nitransforms as nt
import pandas as pd
from scipy.spatial import transform as sst
from scipy import ndimage as ndi

hmc = nt.linear.load(motion_xfm)
timeseries = pd.read_csv(timeseries_tsv, delimiter='\t')[['rot_x', 'rot_y', 'rot_z', 'trans_x', 'trans_y', 'trans_z']]
boldref = nb.load(boldref_nii)

center_of_gravity = np.diag(boldref.header.get_zooms()) @ ndi.center_of_mass(boldref.dataobj)

fsl_hmc = nt.io.fsl.FSLLinearTransformArray.from_ras(hmc.matrix, reference=boldref, moving=boldref)
fsl_matrix = np.stack([xfm['parameters'] for xfm in fsl_hmc.xforms])
mats = fsl_matrix[:, :3, :3].transpose(0, 2, 1)

rot_xyz = sst.Rotation.from_matrix(mats).as_euler('XYZ')
trans_xyz = fsl_matrix[:, :3, 3] - mats @ center_of_gravity + center_of_gravity

recon = pd.DataFrame(data=np.hstack((rot_xyz, trans_xyz)), columns=timeseries.columns)

assert np.abs((recon - timeseries).values).max() < 1e-3

[effigies]

@oesteban Is converting to/from motion parameters as calculated by the different tools (mcflirt, 3dvolreg, spm_realign) in-scope for nitransforms? If not, I'll implement this in fmriprep for now.