|
| 1 | +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- |
| 2 | +# vi: set ft=python sts=4 ts=4 sw=4 et: |
| 3 | +''' |
| 4 | +Miscellaneous algorithms for 2D contours and 3D triangularized meshes handling |
| 5 | +
|
| 6 | + Change directory to provide relative paths for doctests |
| 7 | + >>> import os |
| 8 | + >>> filepath = os.path.dirname( os.path.realpath( __file__ ) ) |
| 9 | + >>> datadir = os.path.realpath(os.path.join(filepath, '../testing/data')) |
| 10 | + >>> os.chdir(datadir) |
| 11 | +
|
| 12 | +''' |
| 13 | + |
| 14 | + |
| 15 | +import os |
| 16 | +import os.path as op |
| 17 | +import warnings |
| 18 | +import nibabel as nb |
| 19 | +import numpy as np |
| 20 | +from scipy.spatial.distance import euclidean |
| 21 | + |
| 22 | +try: |
| 23 | + from tvtk.api import tvtk |
| 24 | +except ImportError: |
| 25 | + warnings.warn('tvtk or vtk not installed') |
| 26 | + raise |
| 27 | + |
| 28 | +from .. import logging |
| 29 | + |
| 30 | +from ..interfaces.base import (BaseInterface, traits, TraitedSpec, File, |
| 31 | + InputMultiPath, OutputMultiPath, |
| 32 | + BaseInterfaceInputSpec, isdefined) |
| 33 | +from ..utils.filemanip import fname_presuffix, split_filename |
| 34 | +iflogger = logging.getLogger('interface') |
| 35 | + |
| 36 | + |
| 37 | +class P2PDistanceInputSpec(BaseInterfaceInputSpec): |
| 38 | + surface1 = File(exists=True, mandatory=True, |
| 39 | + desc="Reference surface (vtk format) to which compute distance.") |
| 40 | + surface2 = File(exists=True, mandatory=True, |
| 41 | + desc="Test surface (vtk format) from which compute distance.") |
| 42 | + weighting = traits.Enum("none", "surface", desc='""none": no weighting is performed\ |
| 43 | + "surface": edge distance is weighted by the corresponding surface area',usedefault=True) |
| 44 | + |
| 45 | +class P2PDistanceOutputSpec(TraitedSpec): |
| 46 | + distance = traits.Float(desc="computed distance") |
| 47 | + |
| 48 | +class P2PDistance(BaseInterface): |
| 49 | + """ |
| 50 | + Calculates a point-to-point (p2p) distance between two corresponding meshes or contours. |
| 51 | + Therefore, a point-to-point correspondence between nodes is required |
| 52 | +
|
| 53 | + Example |
| 54 | + ------- |
| 55 | +
|
| 56 | + >>> import nipype.algorithms.mesh as mesh |
| 57 | + >>> dist = mesh.P2PDistance() |
| 58 | + >>> dist.inputs.surface1 = 'surf1.vtk' |
| 59 | + >>> dist.inputs.surface2 = 'surf2.vtk' |
| 60 | + >>> res = dist.run() # doctest: +SKIP |
| 61 | + """ |
| 62 | + |
| 63 | + def _run_interface(self, runtime): |
| 64 | + r1 = tvtk.PolyDataReader( file_name=self.inputs.surface1 ) |
| 65 | + r2 = tvtk.PolyDataReader( file_name=self.inputs.surface2 ) |
| 66 | + vtk1 = r1.output |
| 67 | + vtk2 = r2.output |
| 68 | + r1.update() |
| 69 | + r2.update() |
| 70 | + |
| 71 | + d = 0.0 |
| 72 | + totalWeight = 0.0 |
| 73 | + for p1,p2 in zip( vtk1.points, vtk2.points ): |
| 74 | + weight = 1.0 |
| 75 | + if (self.inputs.method == 'surface'): |
| 76 | + #compute surfaces, set in weight |
| 77 | + print "Error, not implemented" |
| 78 | + |
| 79 | + d+= euclidean( p1, p2 ) |
| 80 | + totalWeight = totalWeight + weight |
| 81 | + |
| 82 | + self._distance = d / weight |
| 83 | + return runtime |
| 84 | + |
| 85 | + def _list_outputs(self): |
| 86 | + outputs = self._outputs().get() |
| 87 | + outputs['distance'] = self._distance |
| 88 | + return outputs |
| 89 | + |
0 commit comments