Make sure motion parameters are reported in the same convention for tortoise and eddy

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

@@ -303,28 +303,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

@@ -1121,6 +1121,7 @@ class _DIFFPREPMotionParamsInputSpec(BaseInterfaceInputSpec):
 
 class _DIFFPREPMotionParamsOutputSpec(TraitedSpec):
     spm_motion_file = File(exists=True)
+    diffprep_ec_file = File(exists=True)
 
 
 class DIFFPREPMotionParams(SimpleInterface):
@@ -1129,11 +1130,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
@@ -1142,7 +1144,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',
@@ -1151,6 +1156,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/test_interfaces_gradients.py

@@ -8,7 +8,47 @@
 import pytest
 import SimpleITK as sitk
 
-from qsiprep.interfaces.gradients import get_fsl_motion_params
+from qsiprep.interfaces.gradients import get_fsl_motion_params, get_ras_motion_params
+
+
+def test_get_ras_motion_params_no_axis_flip(tmp_path):
+    """RAS export recovers applied motion with correct sign on a radiological grid.
+
+    On a grid with a negative x on the affine diagonal, the FSL/LPS conventions
+    flip x; ``get_ras_motion_params`` must report the applied RAS motion with no
+    flip.
+    """
+    # radiological reference (negative x-diagonal, det < 0)
+    affine = np.diag([-2.0, 2.0, 2.0, 1.0])
+    ref_file = os.path.join(tmp_path, 'ref.nii.gz')
+    nb.Nifti1Image(np.zeros((10, 10, 10), dtype=np.float32), affine=affine).to_filename(ref_file)
+
+    conv = np.diag([-1.0, -1.0, 1.0, 1.0])  # RAS <-> LPS
+    itk_file = os.path.join(tmp_path, 'xfm.mat')
+
+    def ras_to_itk(m_ras):
+        m_lps = conv @ m_ras @ conv
+        aff = sitk.AffineTransform(3)
+        aff.SetMatrix(m_lps[:3, :3].ravel().tolist())
+        aff.SetTranslation(m_lps[:3, 3].tolist())
+        sitk.WriteTransform(aff, itk_file)
+
+    # pure +3 mm translation along RAS x -> +3, not -3
+    m = np.eye(4)
+    m[0, 3] = 3.0
+    ras_to_itk(m)
+    params = get_ras_motion_params(itk_file, ref_file)
+    assert params.shape == (12,)
+    np.testing.assert_allclose(params[9:12], [3.0, 0.0, 0.0], atol=1e-6)
+
+    # pure +5 deg rotation about RAS x -> rotvec x = +5 deg, not -5
+    th = np.deg2rad(5.0)
+    rx = np.array([[1, 0, 0], [0, np.cos(th), -np.sin(th)], [0, np.sin(th), np.cos(th)]])
+    m = np.eye(4)
+    m[:3, :3] = rx
+    ras_to_itk(m)
+    params = get_ras_motion_params(itk_file, ref_file)
+    np.testing.assert_allclose(params[6:9], [th, 0.0, 0.0], atol=1e-6)
 
 
 def test_get_fsl_motion_params_identity_transform(tmp_path):

qsiprep/workflows/dwi/base.py

@@ -486,7 +486,10 @@ def init_dwi_preproc_wf(
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB,
     )
-    workflow.connect([(confounds_wf, ds_confounds, [('outputnode.confounds_file', 'in_file')])])
+    workflow.connect([(confounds_wf, ds_confounds, [
+        ('outputnode.confounds_file', 'in_file'),
+        ('outputnode.confounds_metadata', 'meta_dict'),
+    ])])  # fmt:skip
 
     # Carpetplot and confounds plot
     conf_plot = pe.Node(DMRISummary(), name='conf_plot', mem_gb=mem_gb['resampled'])
@@ -505,6 +508,7 @@ def init_dwi_preproc_wf(
         (hmc_wf, confounds_wf, [
             ('outputnode.slice_quality', 'inputnode.sliceqc_file'),
             ('outputnode.motion_params', 'inputnode.motion_params'),
+            ('outputnode.ec_file', 'inputnode.ec_file'),
         ]),
         (pre_hmc_wf, confounds_wf, [
             ('outputnode.denoising_confounds', 'inputnode.denoising_confounds'),

qsiprep/workflows/dwi/confounds.py

@@ -70,6 +70,7 @@ def init_dwi_confs_wf():
             fields=[
                 'sliceqc_file',
                 'motion_params',
+                'ec_file',
                 'bval_file',
                 'bvec_file',
                 'original_files',
@@ -79,7 +80,8 @@ def init_dwi_confs_wf():
         name='inputnode',
     )
     outputnode = pe.Node(
-        niu.IdentityInterface(fields=['confounds_file', 'imputed_images']), name='outputnode'
+        niu.IdentityInterface(fields=['confounds_file', 'imputed_images', 'confounds_metadata']),
+        name='outputnode',
     )
 
     # Frame displacement
@@ -102,13 +104,17 @@ def init_dwi_confs_wf():
         (add_motion_headers, concat, [('out_file', 'motion')]),
         (inputnode, concat, [
             ('sliceqc_file', 'sliceqc_file'),
+            ('ec_file', 'ec'),
             ('bval_file', 'original_bvals'),
             ('bvec_file', 'original_bvecs'),
             ('original_files', 'original_files'),
             ('denoising_confounds', 'denoising_confounds'),
         ]),
         # Set outputs
-        (concat, outputnode, [('confounds_file', 'confounds_file')]),
+        (concat, outputnode, [
+            ('confounds_file', 'confounds_file'),
+            ('confounds_metadata', 'confounds_metadata'),
+        ]),
     ])  # fmt:skip
 
     return workflow

qsiprep/workflows/dwi/diffprep.py

@@ -259,6 +259,7 @@ def init_diffprep_hmc_wf(
                 'sdc_method',
                 'slice_quality',
                 'motion_params',
+                'ec_file',
                 'cnr_map',
                 'bvec_files_to_transform',
                 'dwi_files_to_transform',
@@ -463,7 +464,10 @@ def init_diffprep_hmc_wf(
             ('b0_indices', 'b0_indices'),
             ('forward_transforms', 'to_dwi_ref_affines'),
         ]),
-        (motion_params, outputnode, [('spm_motion_file', 'motion_params')]),
+        (motion_params, outputnode, [
+            ('spm_motion_file', 'motion_params'),
+            ('diffprep_ec_file', 'ec_file'),
+        ]),
 
         # Pre-SDC enhancement (report)
         (corrected_node, extract_b0s, [('corrected_dwi_file', 'dwi_series')]),