A bunch of bugfixes and improvements

Author: oesteban

nipype/algorithms/mesh.py

@@ -13,6 +13,7 @@
 
 
 import numpy as np
+import os.path as op
 from scipy.spatial.distance import euclidean
 
 from .. import logging
@@ -29,14 +30,19 @@ class P2PDistanceInputSpec(BaseInterfaceInputSpec):
     surface2 = File(exists=True, mandatory=True,
                     desc=("Test surface (vtk format) from which compute "
                           "distance."))
-    weighting = traits.Enum("none", "surface", usedefault=True,
-                            desc=('"none": no weighting is performed, '
-                                  '"surface": edge distance is weighted by the '
-                                  'corresponding surface area'))
+    weighting = traits.Enum(
+        "none", "area", usedefault=True,
+        desc=('"none": no weighting is performed, "area": vertex distances are'
+              'weighted by the total area of faces corresponding to the '
+              'vertex'))
+    out_file = File('distance.npy', usedefault=True,
+                    desc='numpy file keeping computed distances and weights')
 
 
 class P2PDistanceOutputSpec(TraitedSpec):
     distance = traits.Float(desc="computed distance")
+    out_file = File(exists=True,
+                    desc='numpy file keeping computed distances and weights')
 
 
 class P2PDistance(BaseInterface):
@@ -78,6 +84,9 @@ def _run_interface(self, runtime):
         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.surface1)
         r2 = tvtk.PolyDataReader(file_name=self.inputs.surface2)
@@ -86,32 +95,36 @@ def _run_interface(self, runtime):
         r1.update()
         r2.update()
         assert(len(vtk1.points) == len(vtk2.points))
-        d = 0.0
-        totalWeight = 0.0
 
-        points = vtk1.points
-        faces = vtk1.polys.to_array().reshape(-1, 4).astype(int)[:, 1:]
+        points1 = np.array(vtk1.points)
+        points2 = np.array(vtk2.points)
+
+        diff = np.linalg.norm(points1 - points2, axis=1)
+        weights = np.ones(len(diff))
+
+        if (self.inputs.weighting == 'area'):
+            faces = vtk1.polys.to_array().reshape(-1, 4).astype(int)[:, 1:]
 
-        for p1, p2 in zip(points, vtk2.points):
-            weight = 1.0
-            if (self.inputs.weighting == 'surface'):
+            for i, p1 in enumerate(points2):
                 # compute surfaces, set in weight
-                weight = 0.0
-                point_faces = faces[(faces[:, :] == 0).any(axis=1)]
+                w = 0.0
+                point_faces = faces[(faces[:, :] == i).any(axis=1)]
 
                 for idset in point_faces:
-                    p1 = points[int(idset[0])]
-                    p2 = points[int(idset[1])]
-                    p3 = points[int(idset[2])]
-                    weight = weight + self._triangle_area(p1, p2, p3)
+                    fp1 = points1[int(idset[0])]
+                    fp2 = points1[int(idset[1])]
+                    fp3 = points1[int(idset[2])]
+                    w += self._triangle_area(fp1, fp2, fp3)
+                weights[i] = w
 
-            d += weight * euclidean(p1, p2)
-            totalWeight = totalWeight + weight
+        result = np.vstack([diff, weights])
+        np.save(op.abspath(self.inputs.out_file), result.transpose())
 
-        self._distance = d / totalWeight
+        self._distance = np.average(diff, weights=weights)
         return runtime
 
     def _list_outputs(self):
         outputs = self._outputs().get()
+        outputs['out_file'] = op.abspath(self.inputs.out_file)
         outputs['distance'] = self._distance
         return outputs