Make sure motion parameters are reported in the same convention for tortoise and eddy (#1084)

* Export DWI head-motion parameters in RAS

Motion confounds were reported in backend-specific frames: eddy/SHORELine
via get_fsl_motion_params (FSL) and DIFFPREP via raw TORTOISE Okan params
(LPS). Relative to the RAS world frame this flips trans_x/rot_x (eddy) and
trans_x/trans_y/rot_x/rot_y (DIFFPREP), so the two backends disagreed and
neither matched RAS-defined ground truth.

Add get_ras_motion_params (per-volume ITK/LPS transform -> RAS, decomposed
about the image center) and use it in CombineMotions (eddy/SHORELine);
negate the x,y components of the DIFFPREP Okan params (LPS -> RAS). The
confounds trans_*/rot_* columns are now RAS for every HMC backend.

Add test_get_ras_motion_params_no_axis_flip, which verifies no axis flip on
a radiological grid.

* also make the ecc parameters appear

* update diffprep param test

---------

Co-authored-by: Taylor Salo <salot@pennmedicine.upenn.edu>

Author: mattcieslak

qsiprep/interfaces/confounds.py

@@ -35,6 +35,7 @@
 class GatherConfoundsInputSpec(BaseInterfaceInputSpec):
     fd = File(exists=True, desc='input framewise displacement')
     motion = File(exists=True, desc='input motion parameters')
+    ec = File(exists=True, desc='eddy-current field parameters (headed TSV)')
     sliceqc_file = File(exists=True, desc='output from sliceqc')
     original_files = traits.List(desc='original grouping of each volume')
     original_bvecs = InputMultiObject(File(exists=True), desc='original bvec files')
@@ -45,6 +46,7 @@ class GatherConfoundsInputSpec(BaseInterfaceInputSpec):
 class GatherConfoundsOutputSpec(TraitedSpec):
     confounds_file = File(exists=True, desc='output confounds file')
     confounds_list = traits.List(traits.Str, desc='list of headers')
+    confounds_metadata = traits.Dict(desc='per-column descriptions for the confounds JSON sidecar')
 
 
 class GatherConfounds(SimpleInterface):
@@ -61,6 +63,7 @@ def _run_interface(self, runtime):
             fdisp=self.inputs.fd,
             sliceqc_file=self.inputs.sliceqc_file,
             motion=self.inputs.motion,
+            ec=self.inputs.ec,
             original_files=self.inputs.original_files,
             original_bvals=concatenate_bvals(self.inputs.original_bvals, None),
             original_bvecs=concatenate_bvecs(self.inputs.original_bvecs),
@@ -69,12 +72,55 @@ def _run_interface(self, runtime):
         )
         self._results['confounds_file'] = combined_out
         self._results['confounds_list'] = confounds_list
+        columns = pd.read_csv(combined_out, sep='\t', nrows=0).columns.tolist()
+        self._results['confounds_metadata'] = _confounds_column_metadata(columns)
         return runtime
 
 
+def _confounds_column_metadata(columns):
+    """Per-column descriptions for the confounds JSON sidecar.
+
+    Motion columns are RAS+ (translation mm, rotation rad). The eddy-current
+    columns are the raw per-volume field coefficients each backend fits; they are
+    described at the block level (linear / quadratic / centre) with a model
+    reference rather than a per-column basis, since the exact ordering is defined
+    by FSL eddy (``--flm``) and TORTOISE (``OkanQuadraticTransform``).
+    """
+    motion = {
+        'trans_x': 'Translation along RAS+ x (mm)',
+        'trans_y': 'Translation along RAS+ y (mm)',
+        'trans_z': 'Translation along RAS+ z (mm)',
+        'rot_x': 'Rotation about RAS+ x (radians)',
+        'rot_y': 'Rotation about RAS+ y (radians)',
+        'rot_z': 'Rotation about RAS+ z (radians)',
+    }
+    eddy_block = (
+        'FSL eddy first-level-model eddy-current field coefficient (from '
+        '.eddy_parameters). For --flm=quadratic there are 10 per volume: ~3 linear '
+        '(x, y, z), ~6 quadratic/cross, and 1 spare/constant. The exact per-column '
+        'basis is defined by FSL eddy; see its documentation.'
+    )
+    okan_block = (
+        'TORTOISE DIFFPREP OkanQuadraticTransform eddy parameter (cols 6-23 of the '
+        '24-parameter transform): the quadratic eddy-current field (~3 linear x/y/z '
+        'plus quadratic) and the rotation/eddy centres. The exact per-column basis '
+        'is defined by TORTOISE.'
+    )
+    meta = {}
+    for col in columns:
+        if col in motion:
+            meta[col] = {'Description': motion[col]}
+        elif col.startswith('eddy_ec_'):
+            meta[col] = {'Description': eddy_block, 'Source': 'FSL eddy'}
+        elif col.startswith('diffprep_ec_'):
+            meta[col] = {'Description': okan_block, 'Source': 'TORTOISE DIFFPREP'}
+    return meta
+
+
 def _gather_confounds(
     fdisp=None,
     motion=None,
+    ec=None,
     sliceqc_file=None,
     newpath=None,
     original_files=None,
@@ -126,7 +172,11 @@ def _adjust_indices(left_df, right_df):
 
     all_files = []
     confounds_list = []
-    for confound, name in ((fdisp, 'Framewise displacement'), (motion, 'Motion parameters')):
+    for confound, name in (
+        (fdisp, 'Framewise displacement'),
+        (motion, 'Motion parameters'),
+        (ec, 'Eddy-current parameters'),
+    ):
         if confound is not None and isdefined(confound):
             confounds_list.append(name)
             if os.path.exists(confound) and os.stat(confound).st_size > 0:

qsiprep/interfaces/eddy.py

@@ -302,28 +302,65 @@ def _format_arg(self, name, spec, value):
         return super()._format_arg(name, spec, value)
 
 
+def _fsl_to_ras_axis_flip(ref_file):
+    """Per-axis ±1 to convert FSL rigid params to RAS+ for ``ref_file``.
+
+    FSL reports motion in its radiological voxel frame; the RAS direction of each
+    axis follows the sign of the affine diagonal, with FSL flipping x for a
+    neurological (positive-determinant) image. Applied to both translation and
+    (in-frame Euler) rotation components. Assumes an axis-aligned affine, which
+    holds for essentially all DWI acquisitions.
+    """
+    aff = nb.load(ref_file).affine[:3, :3]
+    flip = np.sign(np.diag(aff))
+    if np.linalg.det(aff) > 0:
+        flip[0] *= -1
+    return flip
+
+
 class Eddy2SPMMotionInputSpec(BaseInterfaceInputSpec):
     eddy_motion = File(exists=True)
+    ref_file = File(exists=True, desc='reference image defining the FSL<->RAS axis convention')
 
 
 class Eddy2SPMMotionOututSpec(TraitedSpec):
     spm_motion_file = File(exists=True)
+    eddy_ec_file = File(exists=True)
 
 
 class Eddy2SPMMotion(SimpleInterface):
     input_spec = Eddy2SPMMotionInputSpec
     output_spec = Eddy2SPMMotionOututSpec
 
     def _run_interface(self, runtime):
-        # Load the eddy motion params File
+        # eddy_parameters columns: 0-5 are rigid motion (3 translation mm, 3 rotation rad) in
+        # FSL's radiological frame; the remaining columns are the eddy-current field
+        # coefficients (10 for --flm=quadratic: ~3 linear x/y/z + ~6 quadratic + 1 spare).
         eddy_motion = np.loadtxt(self.inputs.eddy_motion)
-        spm_motion = eddy_motion[:, :6]
+        if eddy_motion.ndim == 1:
+            eddy_motion = eddy_motion[np.newaxis, :]
+
+        # Rigid motion, converted FSL -> RAS+ so it matches the SHORELine/DIFFPREP export.
+        spm_motion = eddy_motion[:, :6].astype(float)
+        if isdefined(self.inputs.ref_file):
+            flip = _fsl_to_ras_axis_flip(self.inputs.ref_file)
+            spm_motion[:, :3] *= flip
+            spm_motion[:, 3:6] *= flip
         spm_motion_file = fname_presuffix(
             self.inputs.eddy_motion, suffix='spm_rp.txt', use_ext=False, newpath=runtime.cwd
         )
         np.savetxt(spm_motion_file, spm_motion)
         self._results['spm_motion_file'] = spm_motion_file
 
+        # Eddy-current field coefficients -> headed TSV confounds columns (eddy_ec_NN).
+        ec = np.atleast_2d(eddy_motion[:, 6:])
+        ec_file = fname_presuffix(
+            self.inputs.eddy_motion, suffix='eddy_ec.tsv', use_ext=False, newpath=runtime.cwd
+        )
+        header = '\t'.join(f'eddy_ec_{i:02d}' for i in range(ec.shape[1]))
+        np.savetxt(ec_file, ec, delimiter='\t', header=header, comments='')
+        self._results['eddy_ec_file'] = ec_file
+
         return runtime
 
 

qsiprep/interfaces/gradients.py

@@ -243,7 +243,7 @@ def _run_interface(self, runtime):
         output_spm_fname = os.path.join(runtime.cwd, 'spm_movpar.txt')
         ref_file = self.inputs.ref_file
         for motion_file in self.inputs.transform_files:
-            collected_motion.append(get_fsl_motion_params(motion_file, ref_file, runtime.cwd))
+            collected_motion.append(get_ras_motion_params(motion_file, ref_file))
 
         final_motion = np.row_stack(collected_motion)
         cols = [
@@ -722,6 +722,48 @@ def get_trans_from_offset(image_center, rotmat):
     return np.concatenate([scale, shear, rotation, translation])
 
 
+def get_ras_motion_params(itk_file, ref_file):
+    """Decompose a per-volume rigid transform into **RAS+** motion parameters.
+
+    Mirror of :func:`get_fsl_motion_params`, but the translation (mm) and
+    rotation (rotation-vector radians) are expressed in the scanner-independent
+    RAS+ world frame instead of FSL's radiological frame. eddy/SHORELine report
+    FSL and DIFFPREP reports LPS; decomposing every backend's transform in RAS
+    removes those per-backend axis flips, so the exported motion is directly
+    comparable to RAS-defined ground truth regardless of grid orientation.
+    """
+    import SimpleITK as sitk
+
+    tfm = sitk.ReadTransform(itk_file)
+    try:
+        aff = sitk.AffineTransform(tfm)
+    except RuntimeError:  # composite with a single affine
+        aff = sitk.AffineTransform(sitk.CompositeTransform(tfm).GetNthTransform(0))
+    mat = np.array(aff.GetMatrix()).reshape(3, 3)
+    center = np.array(aff.GetCenter())
+    offset = np.array(aff.GetTranslation())
+    # ITK stores transforms in LPS: y = mat @ (x - center) + center + offset
+    m_lps = np.eye(4)
+    m_lps[:3, :3] = mat
+    m_lps[:3, 3] = offset + center - mat @ center
+    # LPS -> RAS flips x and y
+    conv = np.diag([-1.0, -1.0, 1.0, 1.0])
+    m_ras = conv @ m_lps @ conv
+
+    _, rotmat, scale, shear = decompose44(m_ras)
+    rotation = R.from_matrix(rotmat).as_rotvec()
+
+    src_img = nb.load(ref_file)
+    src_center = (np.array(src_img.shape[:3]) - 1) / 2
+    center_mm = nb.affines.apply_affine(src_img.affine, src_center) - src_img.affine[:3, 3]
+    translation = np.zeros(3)
+    for i in range(3):
+        translation[i] = (m_ras[i, 3] - center_mm[i]) + (
+            m_ras[i, 0] * center_mm[0] + m_ras[i, 1] * center_mm[1] + m_ras[i, 2] * center_mm[2]
+        )
+    return np.concatenate([scale, shear, rotation, translation])
+
+
 def match_transforms(dwi_files, transforms, b0_indices):
     original_b0_indices = np.array(b0_indices)
     num_dwis = len(dwi_files)

qsiprep/interfaces/tortoise.py

@@ -1131,6 +1131,7 @@ class _DIFFPREPMotionParamsInputSpec(BaseInterfaceInputSpec):
 
 class _DIFFPREPMotionParamsOutputSpec(TraitedSpec):
     spm_motion_file = File(exists=True)
+    diffprep_ec_file = File(exists=True)
 
 
 class DIFFPREPMotionParams(SimpleInterface):
@@ -1139,11 +1140,12 @@ class DIFFPREPMotionParams(SimpleInterface):
 
     The output columns are the leading 6 parameters of TORTOISE's
     ``OkanQuadraticTransform`` in SPM realignment-parameter order
-    (translation_x/y/z in mm of LPS physical coordinate, rotation_x/y/z as
-    Euler angles in radians). The remaining 18 Okan parameters encode the
-    eddy-current polynomial + rotation/eddy centre and are intentionally
-    dropped -- they are not rigid head motion. Units match the eddy and
-    SHORELine SPM motion files (translation mm, rotation radians).
+    (translation_x/y/z in mm, rotation_x/y/z as Euler angles in radians),
+    converted from TORTOISE's native LPS to **RAS+** so they match the
+    eddy/SHORELine motion files (which qsiprep now also exports in RAS via
+    :func:`~qsiprep.interfaces.gradients.get_ras_motion_params`). The remaining
+    18 Okan parameters encode the eddy-current polynomial + rotation/eddy
+    centre and are intentionally dropped -- they are not rigid head motion.
     """
 
     input_spec = _DIFFPREPMotionParamsInputSpec
@@ -1152,7 +1154,10 @@ class DIFFPREPMotionParams(SimpleInterface):
     def _run_interface(self, runtime):
         rows = _read_okan_transformations(self.inputs.transformations_file)
         params = np.asarray(rows, dtype=float)
-        spm_motion = params[:, :6]
+        # Okan params are LPS physical; LPS->RAS is a 180deg rotation about z,
+        # so negate the x and y components of both translation and rotation.
+        spm_motion = params[:, :6].copy()
+        spm_motion[:, [0, 1, 3, 4]] *= -1.0
         spm_motion_file = fname_presuffix(
             self.inputs.transformations_file,
             suffix='_spm_rp.txt',
@@ -1161,6 +1166,16 @@ def _run_interface(self, runtime):
         )
         np.savetxt(spm_motion_file, spm_motion)
         self._results['spm_motion_file'] = spm_motion_file
+
+        # Okan eddy-current + rotation/eddy-centre parameters (cols 6-23) -> headed TSV
+        # confounds columns (diffprep_ec_NN): ~3 linear x/y/z + quadratic + centres.
+        ec = params[:, 6:24]
+        ec_file = fname_presuffix(
+            self.inputs.transformations_file, suffix='_ec.tsv', use_ext=False, newpath=runtime.cwd
+        )
+        header = '\t'.join(f'diffprep_ec_{i:02d}' for i in range(ec.shape[1]))
+        np.savetxt(ec_file, ec, delimiter='\t', header=header, comments='')
+        self._results['diffprep_ec_file'] = ec_file
         return runtime
 
 

qsiprep/tests/data/diffprep_outputs.txt

@@ -18,6 +18,7 @@ sub-tester/anat/sub-tester_space-ACPC_desc-preproc_T1w.json
 sub-tester/anat/sub-tester_space-ACPC_desc-preproc_T1w.nii.gz
 sub-tester/anat/sub-tester_space-ACPC_dseg.nii.gz
 sub-tester/dwi
+sub-tester/dwi/sub-tester_acq-HASC55AP_desc-confounds_timeseries.json
 sub-tester/dwi/sub-tester_acq-HASC55AP_desc-confounds_timeseries.tsv
 sub-tester/dwi/sub-tester_acq-HASC55AP_space-ACPC_desc-brain_mask.nii.gz
 sub-tester/dwi/sub-tester_acq-HASC55AP_space-ACPC_desc-image_qc.tsv

qsiprep/tests/data/drbuddi_rpe_outputs.txt

@@ -7,6 +7,7 @@ logs/CITATION.tex
 sub-tinytensors
 sub-tinytensors.html
 sub-tinytensors/dwi
+sub-tinytensors/dwi/sub-tinytensors_desc-confounds_timeseries.json
 sub-tinytensors/dwi/sub-tinytensors_desc-confounds_timeseries.tsv
 sub-tinytensors/dwi/sub-tinytensors_space-ACPC_desc-brain_mask.nii.gz
 sub-tinytensors/dwi/sub-tinytensors_space-ACPC_desc-image_qc.tsv

qsiprep/tests/data/drbuddi_shoreline_epi_outputs.txt

@@ -7,6 +7,7 @@ logs/CITATION.tex
 sub-tinytensors
 sub-tinytensors.html
 sub-tinytensors/dwi
+sub-tinytensors/dwi/sub-tinytensors_dir-PA_desc-confounds_timeseries.json
 sub-tinytensors/dwi/sub-tinytensors_dir-PA_desc-confounds_timeseries.tsv
 sub-tinytensors/dwi/sub-tinytensors_dir-PA_space-ACPC_desc-brain_mask.nii.gz
 sub-tinytensors/dwi/sub-tinytensors_dir-PA_space-ACPC_desc-image_qc.tsv

qsiprep/tests/data/drbuddi_tensorline_epi_outputs.txt

@@ -7,6 +7,7 @@ logs/CITATION.tex
 sub-PNC
 sub-PNC.html
 sub-PNC/dwi
+sub-PNC/dwi/sub-PNC_acq-realistic_desc-confounds_timeseries.json
 sub-PNC/dwi/sub-PNC_acq-realistic_desc-confounds_timeseries.tsv
 sub-PNC/dwi/sub-PNC_acq-realistic_space-ACPC_desc-brain_mask.nii.gz
 sub-PNC/dwi/sub-PNC_acq-realistic_space-ACPC_desc-image_qc.tsv

qsiprep/tests/data/dscsdsi_outputs.txt

@@ -18,6 +18,7 @@ sub-tester/anat/sub-tester_space-ACPC_desc-preproc_T1w.json
 sub-tester/anat/sub-tester_space-ACPC_desc-preproc_T1w.nii.gz
 sub-tester/anat/sub-tester_space-ACPC_dseg.nii.gz
 sub-tester/dwi
+sub-tester/dwi/sub-tester_acq-HASC55AP_desc-confounds_timeseries.json
 sub-tester/dwi/sub-tester_acq-HASC55AP_desc-confounds_timeseries.tsv
 sub-tester/dwi/sub-tester_acq-HASC55AP_space-ACPC_desc-brain_mask.nii.gz
 sub-tester/dwi/sub-tester_acq-HASC55AP_space-ACPC_desc-image_qc.tsv

qsiprep/tests/data/dsdti_nofmap_outputs.txt

@@ -18,6 +18,7 @@ sub-PNC/anat/sub-PNC_space-ACPC_desc-preproc_T1w.json
 sub-PNC/anat/sub-PNC_space-ACPC_desc-preproc_T1w.nii.gz
 sub-PNC/anat/sub-PNC_space-ACPC_dseg.nii.gz
 sub-PNC/dwi
+sub-PNC/dwi/sub-PNC_acq-realistic_desc-confounds_timeseries.json
 sub-PNC/dwi/sub-PNC_acq-realistic_desc-confounds_timeseries.tsv
 sub-PNC/dwi/sub-PNC_acq-realistic_desc-pepolar_qc.tsv
 sub-PNC/dwi/sub-PNC_acq-realistic_space-ACPC_desc-brain_mask.nii.gz