ENH WIP: Use geodesic distance in FramewiseDisplacement

nipype/algorithms/confounds.py

@@ -13,6 +13,7 @@
 
 import nibabel as nb
 import numpy as np
+
 from numpy.polynomial import Legendre
 from scipy import linalg
 
@@ -236,6 +237,15 @@ class FramewiseDisplacementInputSpec(BaseInterfaceInputSpec):
         "NIPY",
         desc="Source of movement parameters",
         mandatory=True)
+    metric = traits.Enum(
+        'L1',
+        'riemannian',
+        usedefault=True,
+        mandatory=True,
+        desc='Distance metric to apply: '
+        'L1 = Manhattan distance (original definition); '
+        'riemannian = Riemannian distance on the '
+        'Special Euclidean group in 3D (geodesic)')
     radius = traits.Float(
         50,
         usedefault=True,
@@ -307,9 +317,41 @@ def _run_interface(self, runtime):
             axis=1,
             arr=mpars,
             source=self.inputs.parameter_source)
-        diff = mpars[:-1, :6] - mpars[1:, :6]
-        diff[:, 3:6] *= self.inputs.radius
-        fd_res = np.abs(diff).sum(axis=1)
+
+        if self.inputs.metric == 'L1':
+            diff = mpars[:-1, :6] - mpars[1:, :6]
+            diff[:, 3:6] *= self.inputs.radius
+            fd_res = np.abs(diff).sum(axis=1)
+
+        elif self.inputs.metric == 'riemannian':
+            from geomstats.invariant_metric import InvariantMetric
+            from geomstats.special_euclidean_group import SpecialEuclideanGroup
+
+            SE3_GROUP = SpecialEuclideanGroup(n=3)
+            SO3_GROUP = SE3_GROUP.rotations
+            DIM_TRANSLATIONS = SE3_GROUP.translations.dimension
+            DIM_ROTATIONS = SE3_GROUP.rotations.dimension
+
+            # TODO(nina): Check SPM convention for Euler angles
+            # (Tait-Bryan angles)
+            so3pars = mpars[:, DIM_TRANSLATIONS:]
+            so3pars = SO3_GROUP.rotation_vector_from_tait_bryan_angles(
+                so3pars,
+                extrinsic_or_intrinsic='intrinsic',
+                order='xyz')
+
+            se3pars = np.hstack(
+                [so3pars, mpars[:, :DIM_TRANSLATIONS]])
+
+            diag_rotations = self.inputs.radius * np.ones(DIM_ROTATIONS)
+            diag_translations = np.ones(DIM_TRANSLATIONS)
+            diag = np.concatenate([diag_rotations, diag_translations])
+            inner_product = np.diag(diag)
+            metric = InvariantMetric(
+                       group=SE3_GROUP,
+                       inner_product_mat_at_identity=inner_product,
+                       left_or_right='left')
+            fd_res = metric.dist(se3pars[:-1], se3pars[1:])
 
         self._results = {
             'out_file': op.abspath(self.inputs.out_file),

nipype/algorithms/tests/test_auto_FramewiseDisplacement.py

@@ -12,6 +12,7 @@ def test_FramewiseDisplacement_inputs():
         out_figure=dict(usedefault=True, ),
         out_file=dict(usedefault=True, ),
         parameter_source=dict(mandatory=True, ),
+        metric=dict(usedefault=True, ),
         radius=dict(usedefault=True, ),
         save_plot=dict(usedefault=True, ),
         series_tr=dict(),
@@ -21,6 +22,8 @@ def test_FramewiseDisplacement_inputs():
     for key, metadata in list(input_map.items()):
         for metakey, value in list(metadata.items()):
             assert getattr(inputs.traits()[key], metakey) == value
+
+
 def test_FramewiseDisplacement_outputs():
     output_map = dict(
         fd_average=dict(),

nipype/algorithms/tests/test_confounds.py

@@ -17,10 +17,18 @@
 except ImportError:
     pass
 
+nogeomstats = True
+try:
+    import geomstats
+    nogeomstats = False
+except ImportError:
+    pass
+
 
 def test_fd(tmpdir):
     tempdir = tmpdir.strpath
     ground_truth = np.loadtxt(example_data('fsl_motion_outliers_fd.txt'))
+
     fdisplacement = FramewiseDisplacement(
         in_file=example_data('fsl_mcflirt_movpar.txt'),
         out_file=tempdir + '/fd.txt',
@@ -37,6 +45,29 @@ def test_fd(tmpdir):
     assert np.abs(ground_truth.mean() - res.outputs.fd_average) < 1e-2
 
 
+@pytest.mark.skipif(nonitime, reason="nitime is not installed")
+def test_fd_riemannian(tmpdir):
+    tempdir = tmpdir.strpath
+    # TODO(nina): Adapt ground_truth w. SPM Euler angles convention
+    ground_truth = np.loadtxt(example_data('fsl_motion_outliers_fd.txt'))
+
+    fdisplacement = FramewiseDisplacement(
+        in_file=example_data('fsl_mcflirt_movpar.txt'),
+        out_file=tempdir + '/fd.txt',
+        parameter_source="FSL",
+        metric='riemannian')
+    res = fdisplacement.run()
+
+    with open(res.outputs.out_file) as all_lines:
+        for line in all_lines:
+            assert 'FramewiseDisplacement' in line
+            break
+
+    assert np.allclose(
+        ground_truth, np.loadtxt(res.outputs.out_file, skiprows=1), atol=.16)
+    assert np.abs(ground_truth.mean() - res.outputs.fd_average) < 1e-2
+
+
 @pytest.mark.skipif(nonitime, reason="nitime is not installed")
 def test_dvars(tmpdir):
     ground_truth = np.loadtxt(example_data('ds003_sub-01_mc.DVARS'))

nipype/info.py

@@ -160,7 +160,8 @@ def get_nipype_gitversion():
     'doc': ['Sphinx>=1.4', 'numpydoc', 'matplotlib', 'pydotplus', 'pydot>=1.2.3'],
     'tests': TESTS_REQUIRES,
     'specs': ['yapf'],
-    'nipy': ['nitime', 'nilearn', 'dipy', 'nipy', 'matplotlib'],
+    'nipy': ['nitime', 'nilearn', 'dipy', 'nipy', 'matplotlib',
+             'geomstats>=1.1.0; python_version >= "3.5"'],
     'profiler': ['psutil>=5.0'],
     'duecredit': ['duecredit'],
     'xvfbwrapper': ['xvfbwrapper'],