fix regression tests, add new MeshWarpMaths interface

Author: oesteban

CHANGES

@@ -1,6 +1,11 @@
 Next release
 ============
 
+* ENH: New mesh.MeshWarpMaths to operate on surface-defined warpings
+       (https://github.com/nipy/nipype/pull/1016)
+* FIX: Refactor P2PDistance, change name to ComputeMeshWarp, add regression tests,
+       fix bug in area weighted distance, and added optimizations
+       (https://github.com/nipy/nipype/pull/1016)
 * FIX: Amend create_tbss_non_fa() workflow to match FSL's tbss_non_fa command. (https://github.com/nipy/nipype/pull/1033)
 * FIX: remove unused mandatory flag from spm normalize (https://github.com/nipy/nipype/pull/1048)
 * ENH: Update ANTSCorticalThickness interface (https://github.com/nipy/nipype/pull/1013)

nipype/algorithms/mesh.py

@@ -13,8 +13,9 @@
 
 
 import numpy as np
+from numpy.linalg import norm
 import os.path as op
-from scipy.spatial.distance import euclidean
+from nipype.external import six
 
 from .. import logging
 
@@ -79,14 +80,16 @@ class ComputeMeshWarp(BaseInterface):
     output_spec = ComputeMeshWarpOutputSpec
 
     def _triangle_area(self, A, B, C):
-        ABxAC = euclidean(A, B) * euclidean(A, C)
-        prod = np.dot(np.array(B) - np.array(A), np.array(C) - np.array(A))
+        A = np.array(A)
+        B = np.array(B)
+        C = np.array(C)
+        ABxAC = norm(A - B) * norm(A - C)
+        prod = np.dot(B - A, C - A)
         angle = np.arccos(prod / ABxAC)
         area = 0.5 * ABxAC * np.sin(angle)
         return area
 
     def _run_interface(self, runtime):
-        from numpy import linalg as nla
         try:
             from tvtk.api import tvtk
         except ImportError:
@@ -116,7 +119,7 @@ def _run_interface(self, runtime):
         diff = points2 - points1
         weights = np.ones(len(diff))
 
-        errvector = nla.norm(diff, axis=1)
+        errvector = norm(diff, axis=1)
         if self.inputs.metric == 'sqeuclidean':
             errvector = errvector ** 2
 
@@ -147,7 +150,6 @@ def _run_interface(self, runtime):
             file_name=op.abspath(self.inputs.out_warp))
         writer.set_input_data(out_mesh)
         writer.write()
-        #write_data(out_mesh, op.abspath(self.inputs.out_warp))
 
         self._distance = np.average(errvector, weights=weights)
         return runtime
@@ -160,6 +162,144 @@ def _list_outputs(self):
         return outputs
 
 
+class MeshWarpMathsInputSpec(BaseInterfaceInputSpec):
+    in_surf = File(exists=True, mandatory=True,
+                   desc=('Input surface in vtk format, with associated warp '
+                         'field as point data (ie. from ComputeMeshWarp'))
+    float_trait = traits.Either(traits.Float(1.0), traits.Tuple(
+        traits.Float(1.0), traits.Float(1.0), traits.Float(1.0)))
+
+    operator = traits.Either(
+        float_trait, File(exists=True), default=1.0, mandatory=True,
+        desc=('image, float or tuple of floats to act as operator'))
+
+    operation = traits.Enum('sum', 'sub', 'mul', 'div', usedefault=True,
+                            desc=('operation to be performed'))
+
+    out_warp = File('warp_maths.vtk', usedefault=True,
+                    desc='vtk file based on in_surf and warpings mapping it '
+                    'to out_file')
+    out_file = File('warped_surf.vtk', usedefault=True,
+                    desc='vtk with surface warped')
+
+
+class MeshWarpMathsOutputSpec(TraitedSpec):
+    out_warp = File(exists=True, desc=('vtk file with the vertex-wise '
+                                       'mapping of surface1 to surface2'))
+    out_file = File(exists=True,
+                    desc='vtk with surface warped')
+
+
+class MeshWarpMaths(BaseInterface):
+
+    """
+    Performs the most basic mathematical operations on the warping field
+    defined at each vertex of the input surface. A surface with scalar
+    or vector data can be used as operator for non-uniform operations.
+
+    .. warning:
+
+      A point-to-point correspondence between surfaces is required
+
+
+    Example
+    -------
+
+    >>> import nipype.algorithms.mesh as m
+    >>> mmath = m.MeshWarpMaths()
+    >>> mmath.inputs.in_surf = 'surf1.vtk'
+    >>> mmath.inputs.operator = 'surf2.vtk'
+    >>> mmath.inputs.operation = 'mul'
+    >>> res = mmath.run() # doctest: +SKIP
+
+    """
+
+    input_spec = MeshWarpMathsInputSpec
+    output_spec = MeshWarpMathsOutputSpec
+
+    def _run_interface(self, runtime):
+        try:
+            from tvtk.api import tvtk
+        except ImportError:
+            raise ImportError('Interface ComputeMeshWarp requires tvtk')
+
+        try:
+            from enthought.etsconfig.api import ETSConfig
+            ETSConfig.toolkit = 'null'
+        except ImportError:
+            iflogger.warn(('ETS toolkit could not be imported'))
+            pass
+        except ValueError:
+            iflogger.warn(('ETS toolkit is already set'))
+            pass
+
+        r1 = tvtk.PolyDataReader(file_name=self.inputs.in_surf)
+        vtk1 = r1.output
+        r1.update()
+        points1 = np.array(vtk1.points)
+
+        if vtk1.point_data.vectors is None:
+            raise RuntimeError(('No warping field was found in in_surf'))
+
+        operator = self.inputs.operator
+        opfield = np.ones_like(points1)
+
+        if isinstance(operator, six.string_types):
+            r2 = tvtk.PolyDataReader(file_name=self.inputs.surface2)
+            vtk2 = r2.output
+            r2.update()
+            assert(len(points1) == len(vtk2.points))
+
+            opfield = vtk2.point_data.vectors
+
+            if opfield is None:
+                opfield = vtk2.point_data.scalars
+
+            if opfield is None:
+                raise RuntimeError(
+                    ('No operator values found in operator file'))
+
+            opfield = np.array(opfield)
+
+            if opfield.shape[1] < points1.shape[1]:
+                opfield = np.array([opfield.tolist()] * points1.shape[1]).T
+        else:
+            operator = np.atleast_1d(operator)
+            opfield *= operator
+
+        warping = np.array(vtk1.point_data.vectors)
+
+        if self.inputs.operation == 'sum':
+            warping += opfield
+        elif self.inputs.operation == 'sub':
+            warping -= opfield
+        elif self.inputs.operation == 'mul':
+            warping *= opfield
+        elif self.inputs.operation == 'div':
+            warping /= opfield
+
+        vtk1.point_data.vectors = warping
+        writer = tvtk.PolyDataWriter(
+            file_name=op.abspath(self.inputs.out_warp))
+        writer.set_input_data(vtk1)
+        writer.write()
+
+        vtk1.point_data.vectors = None
+        vtk1.points = points1 + warping
+        writer = tvtk.PolyDataWriter(
+            file_name=op.abspath(self.inputs.out_file))
+        writer.set_input_data(vtk1)
+        writer.write()
+
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs['out_file'] = op.abspath(self.inputs.out_file)
+        outputs['out_warp'] = op.abspath(self.inputs.out_warp)
+        return outputs
+
+
 class P2PDistance(ComputeMeshWarp):
 
     """

nipype/algorithms/tests/test_auto_MeshWarpMaths.py

@@ -0,0 +1,36 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from nipype.testing import assert_equal
+from nipype.algorithms.mesh import MeshWarpMaths
+
+def test_MeshWarpMaths_inputs():
+    input_map = dict(float_trait=dict(),
+    ignore_exception=dict(nohash=True,
+    usedefault=True,
+    ),
+    in_surf=dict(mandatory=True,
+    ),
+    operation=dict(usedefault=True,
+    ),
+    operator=dict(mandatory=True,
+    ),
+    out_file=dict(usedefault=True,
+    ),
+    out_warp=dict(usedefault=True,
+    ),
+    )
+    inputs = MeshWarpMaths.input_spec()
+
+    for key, metadata in input_map.items():
+        for metakey, value in metadata.items():
+            yield assert_equal, getattr(inputs.traits()[key], metakey), value
+
+def test_MeshWarpMaths_outputs():
+    output_map = dict(out_file=dict(),
+    out_warp=dict(),
+    )
+    outputs = MeshWarpMaths.output_spec()
+
+    for key, metadata in output_map.items():
+        for metakey, value in metadata.items():
+            yield assert_equal, getattr(outputs.traits()[key], metakey), value
+

nipype/algorithms/tests/test_mesh_ops.py

@@ -6,7 +6,7 @@
 from shutil import rmtree
 from tempfile import mkdtemp
 
-from nipype.testing import (assert_equal, assert_raises,
+from nipype.testing import (assert_equal, assert_raises, skipif,
                             assert_almost_equal, example_data)
 
 import numpy as np
@@ -34,11 +34,11 @@ def test_ident_distances():
     dist_ident.inputs.surface2 = in_surf
     dist_ident.inputs.out_file = os.path.join(tempdir, 'distance.npy')
     res = dist_ident.run()
-    yield assert_equal, res.outputs.distance == 0.0
+    yield assert_equal, res.outputs.distance, 0.0
 
     dist_ident.inputs.weighting = 'area'
     res = dist_ident.run()
-    yield assert_equal, res.outputs.distance == 0.0
+    yield assert_equal, res.outputs.distance, 0.0
 
     os.chdir(curdir)
     rmtree(tempdir)
@@ -68,10 +68,10 @@ def test_trans_distances():
     dist.inputs.surface2 = warped_surf
     dist.inputs.out_file = os.path.join(tempdir, 'distance.npy')
     res = dist.run()
-    yield assert_equal, np.isclose(res.outputs.distance, np.norm(inc))
-    dist_ident.inputs.weighting = 'area'
-    res = dist_ident.run()
-    yield assert_equal, np.isclose(res.outputs.distance, np.norm(inc))
+    yield assert_almost_equal, res.outputs.distance, np.linalg.norm(inc), 4
+    dist.inputs.weighting = 'area'
+    res = dist.run()
+    yield assert_almost_equal, res.outputs.distance, np.linalg.norm(inc), 4
 
     os.chdir(curdir)
     rmtree(tempdir)