def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: No input surface.')

        triangleFilter = vtk.vtkTriangleFilter()
        triangleFilter.SetInput(self.Surface)
        triangleFilter.Update()

        decimationFilter = vtk.vtkDecimatePro()
        decimationFilter.SetInput(triangleFilter.GetOutput())
        decimationFilter.SetTargetReduction(self.TargetReduction)
        decimationFilter.SetBoundaryVertexDeletion(self.BoundaryVertexDeletion)
        decimationFilter.PreserveTopologyOn()
        decimationFilter.Update()

        cleaner = vtk.vtkCleanPolyData()
        cleaner.SetInput(decimationFilter.GetOutput())
        cleaner.Update()

        triangleFilter = vtk.vtkTriangleFilter()
        triangleFilter.SetInput(cleaner.GetOutput())
        triangleFilter.Update()

        self.Surface = triangleFilter.GetOutput()

        if self.Surface.GetSource():
            self.Surface.GetSource().UnRegisterAllOutputs()

    def __init__(self, module_manager):
        # initialise our base class
        ModuleBase.__init__(self, module_manager)


        # we'll be playing around with some vtk objects, this could
        # be anything
        self._triangleFilter = vtk.vtkTriangleFilter()
        self._curvatures = vtk.vtkCurvatures()
        self._curvatures.SetCurvatureTypeToMaximum()
        self._curvatures.SetInput(self._triangleFilter.GetOutput())

        # initialise any mixins we might have
        NoConfigModuleMixin.__init__(self,
                {'Module (self)' : self,
                    'vtkTriangleFilter' : self._triangleFilter,
                    'vtkCurvatures' : self._curvatures})

        module_utils.setup_vtk_object_progress(self, self._triangleFilter,
                                           'Triangle filtering...')
        module_utils.setup_vtk_object_progress(self, self._curvatures,
                                           'Calculating curvatures...')
        
        
        self.sync_module_logic_with_config()

 def WriteTecplotSurfaceFile(self):
     if self.OutputFileName == "":
         self.PrintError("Error: no OutputFileName.")
     self.PrintLog("Writing Tecplot file.")
     triangleFilter = vtk.vtkTriangleFilter()
     triangleFilter.SetInput(self.Surface)
     triangleFilter.PassVertsOff()
     triangleFilter.PassLinesOff()
     triangleFilter.Update()
     self.Surface = triangleFilter.GetOutput()
     f = open(self.OutputFileName, "w")
     line = "VARIABLES = X,Y,Z"
     arrayNames = []
     for i in range(self.Surface.GetPointData().GetNumberOfArrays()):
         array = self.Surface.GetPointData().GetArray(i)
         arrayName = array.GetName()
         if arrayName == None:
             continue
         if arrayName[-1] == "_":
             continue
         arrayNames.append(arrayName)
         if array.GetNumberOfComponents() == 1:
             line = line + "," + arrayName
         else:
             for j in range(array.GetNumberOfComponents()):
                 line = line + "," + arrayName + str(j)
     line = line + "\n"
     f.write(line)
     line = (
         "ZONE "
         + "N="
         + str(self.Surface.GetNumberOfPoints())
         + ","
         + "E="
         + str(self.Surface.GetNumberOfCells())
         + ","
         + "F=FEPOINT"
         + ","
         + "ET=TRIANGLE"
         + "\n"
     )
     f.write(line)
     for i in range(self.Surface.GetNumberOfPoints()):
         point = self.Surface.GetPoint(i)
         line = str(point[0]) + " " + str(point[1]) + " " + str(point[2])
         for arrayName in arrayNames:
             array = self.Surface.GetPointData().GetArray(arrayName)
             for j in range(array.GetNumberOfComponents()):
                 line = line + " " + str(array.GetComponent(i, j))
         line = line + "\n"
         f.write(line)
     for i in range(self.Surface.GetNumberOfCells()):
         cellPointIds = self.Surface.GetCell(i).GetPointIds()
         line = ""
         for j in range(cellPointIds.GetNumberOfIds()):
             if j > 0:
                 line = line + " "
             line = line + str(cellPointIds.GetId(j) + 1)
         line = line + "\n"
         f.write(line)

def MakeTorus():
    '''
    Make a torus as the source.
    :return: vtkPolyData with normal and scalar data.
    '''
    source = vtk.vtkSuperquadricSource();
    source.SetCenter(0.0, 0.0, 0.0)
    source.SetScale(1.0, 1.0, 1.0)
    source.SetPhiResolution(64)
    source.SetThetaResolution(64)
    source.SetThetaRoundness(1)
    source.SetThickness(0.5)
    source.SetSize(10)
    source.SetToroidal(1)

    # The quadric is made of strips, so pass it through a triangle filter as
    # the curvature filter only operates on polys
    tri = vtk.vtkTriangleFilter()
    tri.SetInputConnection(source.GetOutputPort())

    # The quadric has nasty discontinuities from the way the edges are generated
    # so let's pass it though a CleanPolyDataFilter and merge any points which
    # are coincident, or very close
    cleaner = vtk.vtkCleanPolyData()
    cleaner.SetInputConnection(tri.GetOutputPort())
    cleaner.SetTolerance(0.005)
    cleaner.Update()
    return CalculateCurvatures(MakeElevations(cleaner.GetOutput()))

def ugrid2pdata(
        ugrid,
        only_trianlges=False,
        verbose=0):

    mypy.my_print(verbose, "*** ugrid2pdata ***")

    filter_geometry = vtk.vtkGeometryFilter()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        filter_geometry.SetInputData(ugrid)
    else:
        filter_geometry.SetInput(ugrid)
    filter_geometry.Update()
    pdata = filter_geometry.GetOutput()

    if (only_trianlges):
        filter_triangle = vtk.vtkTriangleFilter()
        if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
            filter_triangle.SetInputData(pdata)
        else:
            filter_triangle.SetInput(pdata)
        filter_triangle.Update()
        pdata = filter_triangle.GetOutput()

    return pdata

def vtk_ensure_trilist(polydata):
    try:
        import vtk
        from vtk.util.numpy_support import vtk_to_numpy

        trilist = vtk_to_numpy(polydata.GetPolys().GetData())

        # 5 is the triangle type - if we have another type we need to
        # use a vtkTriangleFilter
        c = vtk.vtkCellTypes()
        polydata.GetCellTypes(c)

        if c.GetNumberOfTypes() != 1 or polydata.GetCellType(0) != 5:
            warnings.warn('Non-triangular mesh connectivity was detected - '
                          'this is currently unsupported and thus the '
                          'connectivity is being coerced into a triangular '
                          'mesh. This may have unintended consequences.')
            t_filter = vtk.vtkTriangleFilter()
            t_filter.SetInputData(polydata)
            t_filter.Update()
            trilist = vtk_to_numpy(t_filter.GetOutput().GetPolys().GetData())

        return trilist.reshape([-1, 4])[:, 1:]
    except Exception as e:
        warnings.warn(str(e))
        return None

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkTriangleFilter(), 'Processing.',
         ('vtkPolyData',), ('vtkPolyData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

        def MakeText(primitive):

            tf.update({primitive: vtk.vtkTriangleFilter()})
            tf[primitive].SetInputConnection(primitive.GetOutputPort())

            mp.update({primitive: vtk.vtkMassProperties()})
            mp[primitive].SetInputConnection(tf[primitive].GetOutputPort())

            # here we capture stdout and write it to a variable for processing.
            summary = StringIO.StringIO()
            # save the original stdout
            old_stdout = sys.stdout
            sys.stdout = summary

            print mp[primitive]
            summary = summary.getvalue()

            startSum = summary.find("  VolumeX")
            endSum = len(summary)
            print summary[startSum:]
            # Restore stdout
            sys.stdout = old_stdout

            vt.update({primitive: vtk.vtkVectorText()})
            vt[primitive].SetText(summary[startSum:])

            pdm.update({primitive: vtk.vtkPolyDataMapper()})
            pdm[primitive].SetInputConnection(vt[primitive].GetOutputPort())

            ta.update({primitive: vtk.vtkActor()})
            ta[primitive].SetMapper(pdm[primitive])
            ta[primitive].SetScale(0.2, 0.2, 0.2)
            return ta[primitive]

def readMeshFile(filename, verbose=False):
    """Read mesh file.
    The input format is determined by file name extension. Degenerate data gets
    removed and polygons get split into triangles to support varios restrictive
    output formats."""

    informat = path.splitext(options.infilename)[1].strip('.')
    # set reader based on filename extension
    if informat=='stl':
        reader = vtk.vtkSTLReader()
    elif informat=='vtk':
        reader = vtk.vtkPolyDataReader()
    elif informat=='obj':
        reader = vtk.vtkMNIObjectReader()
    #elif informat=='tag':
    #    reader = vtk.vtkMNITagPointReader()
    else:
        raise ValueError('cannot read input format' + informat)
    reader.SetFileName(filename)

    # merge duplicate points, and/or remove unused points and/or remove degenerate cells
    clean = vtk.vtkCleanPolyData()
    clean.SetInputConnection(reader.GetOutputPort())

    # convert input polygons and strips to triangles
    triangles = vtk.vtkTriangleFilter()
    triangles.SetInputConnection(clean.GetOutputPort())

    #triangles = reader.GetOutputPort()  # skipping above 'cleaning' doesn't work
    if verbose:
        print "read", filename

    return triangles

def convertToImplicitPolyDataDistance(polyData):
    triFilter = vtk.vtkTriangleFilter() # triFilter to clean the polydata
    triFilter.SetInputData(polyData)
    triFilter.Update()
    ImpDist = vtk.vtkImplicitPolyDataDistance()
    ImpDist.SetTolerance(1e-1)
    ImpDist.SetInput(triFilter.GetOutput())
    return ImpDist

 def cleanerAndTriangleFilter(self, inputModel):
     cleanerPolydata = vtk.vtkCleanPolyData()
     cleanerPolydata.SetInputData(inputModel.GetPolyData())
     cleanerPolydata.Update()
     triangleFilter = vtk.vtkTriangleFilter()
     triangleFilter.SetInputData(cleanerPolydata.GetOutput())
     triangleFilter.Update()
     inputModel.SetAndObservePolyData(triangleFilter.GetOutput())

    def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: no Surface.')

        if not self.vmtkRenderer:
            self.vmtkRenderer = vmtkrenderer.vmtkRenderer()
            self.vmtkRenderer.Initialize()
            self.OwnRenderer = 1

        self.vmtkRenderer.RegisterScript(self) 

        triangleFilter = vtk.vtkTriangleFilter()
        triangleFilter.SetInputData(self.Surface)
        triangleFilter.Update()

        self.Surface = triangleFilter.GetOutput()

        contourScalars = vtk.vtkDoubleArray()
        contourScalars.SetNumberOfComponents(1)
        contourScalars.SetNumberOfTuples(self.Surface.GetNumberOfPoints())
        contourScalars.SetName(self.ContourScalarsArrayName)
        contourScalars.FillComponent(0,self.OutsideValue)

        self.Surface.GetPointData().AddArray(contourScalars)
        self.Surface.GetPointData().SetActiveScalars(self.ContourScalarsArrayName)

        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInputData(self.Surface)
        mapper.ScalarVisibilityOn()
        self.Actor = vtk.vtkActor()
        self.Actor.SetMapper(mapper)
        self.Actor.GetMapper().SetScalarRange(-1.0,0.0)
        self.vmtkRenderer.Renderer.AddActor(self.Actor)

        self.ContourWidget = vtk.vtkContourWidget()
        self.ContourWidget.SetInteractor(self.vmtkRenderer.RenderWindowInteractor)

        rep = vtk.vtkOrientedGlyphContourRepresentation.SafeDownCast(self.ContourWidget.GetRepresentation())
        rep.GetLinesProperty().SetColor(1, 0.2, 0)
        rep.GetLinesProperty().SetLineWidth(3.0)

        pointPlacer = vtk.vtkPolygonalSurfacePointPlacer()
        pointPlacer.AddProp(self.Actor)
        pointPlacer.GetPolys().AddItem(self.Surface)
        rep.SetPointPlacer(pointPlacer)

        self.Interpolator = vtk.vtkPolygonalSurfaceContourLineInterpolator()
        self.Interpolator.GetPolys().AddItem(self.Surface)
        rep.SetLineInterpolator(self.Interpolator)

        self.vmtkRenderer.AddKeyBinding('space','Generate scalars',self.ScalarsCallback)
        self.vmtkRenderer.AddKeyBinding('d','Delete contour',self.DeleteContourCallback)
        self.vmtkRenderer.AddKeyBinding('i','Start interaction',self.InteractCallback)
        self.Display()

        if self.OwnRenderer:
            self.vmtkRenderer.Deallocate()

    def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: No input surface.')

        cleaner = vtk.vtkCleanPolyData()
        cleaner.SetInput(self.Surface)
        cleaner.Update()

        triangleFilter = vtk.vtkTriangleFilter()
        triangleFilter.SetInput(cleaner.GetOutput())
        triangleFilter.Update()

        self.Surface = triangleFilter.GetOutput()

        if self.ElementSizeMode == 'edgelength':
            self.TargetArea = 0.25 * 3.0**0.5 * self.TargetEdgeLength**2
        elif self.ElementSizeMode == 'edgelengtharray':
            calculator = vtk.vtkArrayCalculator()
            calculator.SetInput(self.Surface)
            calculator.AddScalarArrayName(self.TargetEdgeLengthArrayName,0)
            calculator.SetFunction("%f^2 * 0.25 * sqrt(3) * %s^2" % (self.TargetEdgeLengthFactor,self.TargetEdgeLengthArrayName))
            calculator.SetResultArrayName(self.TargetAreaArrayName)
            calculator.Update()
            self.MaxArea = 0.25 * 3.0**0.5 * self.MaxEdgeLength**2
            self.MinArea = 0.25 * 3.0**0.5 * self.MinEdgeLength**2
            self.Surface = calculator.GetOutput()

        surfaceRemeshing = vtkvmtk.vtkvmtkPolyDataSurfaceRemeshing()
        surfaceRemeshing.SetInput(self.Surface)
        if self.CellEntityIdsArrayName:
            surfaceRemeshing.SetCellEntityIdsArrayName(self.CellEntityIdsArrayName)
        if self.ElementSizeMode in ['area','edgelength']:
            surfaceRemeshing.SetElementSizeModeToTargetArea()
        elif self.ElementSizeMode in ['areaarray','edgelengtharray']:
            surfaceRemeshing.SetElementSizeModeToTargetAreaArray()
            surfaceRemeshing.SetTargetAreaArrayName(self.TargetAreaArrayName)
        else:
            self.PrintError('Error: unsupported ElementSizeMode.')
        surfaceRemeshing.SetTargetArea(self.TargetArea)
        surfaceRemeshing.SetTargetAreaFactor(self.TargetAreaFactor)
        surfaceRemeshing.SetMaxArea(self.MaxArea)
        surfaceRemeshing.SetMinArea(self.MinArea)
        surfaceRemeshing.SetNumberOfIterations(self.NumberOfIterations)
        surfaceRemeshing.SetNumberOfConnectivityOptimizationIterations(self.NumberOfConnectivityOptimizationIterations)
        surfaceRemeshing.SetRelaxation(self.Relaxation)
        surfaceRemeshing.SetMinAreaFactor(self.MinAreaFactor)
        surfaceRemeshing.SetAspectRatioThreshold(self.AspectRatioThreshold)
        surfaceRemeshing.SetInternalAngleTolerance(self.InternalAngleTolerance)
        surfaceRemeshing.SetNormalAngleTolerance(self.NormalAngleTolerance)
        surfaceRemeshing.SetCollapseAngleThreshold(self.CollapseAngleThreshold)
        surfaceRemeshing.SetPreserveBoundaryEdges(self.PreserveBoundaryEdges)
        surfaceRemeshing.Update()

        self.Surface = surfaceRemeshing.GetOutput()

        if self.Surface.GetSource():
            self.Surface.GetSource().UnRegisterAllOutputs()

    def _readstls(self):
        """
        Reads in all STL files contained in directories indicated by **ref_dir** and **def_dir**.
        Also calls **_make3Dmesh()** to create 3-D tetrahedral meshes.

        Returns
        -------
        rsurfs, dsurfs
        """
        for fname in sorted(os.listdir(self._ref_dir)):
            if '.stl' in fname.lower():
                reader = vtk.vtkSTLReader()
                reader.SetFileName(
                    str(os.path.normpath(self._ref_dir + os.sep + fname)))
                reader.Update()
                triangles = vtk.vtkTriangleFilter()
                triangles.SetInputConnection(reader.GetOutputPort())
                triangles.Update()
                self.rsurfs.append(triangles.GetOutput())
                massProps = vtk.vtkMassProperties()
                massProps.SetInputData(self.rsurfs[-1])
                massProps.Update()
                print(("Generating tetrahedral mesh from {:s}".format(fname)))
                self._make3Dmesh(
                    str(os.path.normpath(self._ref_dir + os.sep + fname)),
                    'MATERIAL', massProps.GetVolume())

        for fname in sorted(os.listdir(self._def_dir)):
            if '.stl' in fname.lower():
                reader = vtk.vtkSTLReader()
                reader.SetFileName(
                    str(os.path.normpath(self._def_dir + os.sep + fname)))
                reader.Update()
                triangles = vtk.vtkTriangleFilter()
                triangles.SetInputConnection(reader.GetOutputPort())
                triangles.Update()
                self.dsurfs.append(triangles.GetOutput())
                massProps = vtk.vtkMassProperties()
                massProps.SetInputData(self.dsurfs[-1])
                massProps.Update()
                print(("Generating tetrahedral mesh from {:s}".format(fname)))
                self._make3Dmesh(
                    str(os.path.normpath(self._def_dir + os.sep + fname)),
                    'SPATIAL', massProps.GetVolume())

  def applyFilters(self, state):

    surface = None
    surface = state.inputModelNode.GetPolyDataConnection()

    if state.decimation:
      triangle = vtk.vtkTriangleFilter()
      triangle.SetInputConnection(surface)
      decimation = vtk.vtkDecimatePro()
      decimation.SetTargetReduction(state.reduction)
      decimation.SetBoundaryVertexDeletion(state.boundaryDeletion)
      decimation.PreserveTopologyOn()
      decimation.SetInputConnection(triangle.GetOutputPort())
      surface = decimation.GetOutputPort()

    if state.smoothing:
      if state.smoothingMethod == "Laplace":
        smoothing = vtk.vtkSmoothPolyDataFilter()
        smoothing.SetBoundarySmoothing(state.boundarySmoothing)
        smoothing.SetNumberOfIterations(state.laplaceIterations)
        smoothing.SetRelaxationFactor(state.laplaceRelaxation)
        smoothing.SetInputConnection(surface)
        surface = smoothing.GetOutputPort()
      elif state.smoothingMethod == "Taubin":
        smoothing = vtk.vtkWindowedSincPolyDataFilter()
        smoothing.SetBoundarySmoothing(state.boundarySmoothing)
        smoothing.SetNumberOfIterations(state.taubinIterations)
        smoothing.SetPassBand(state.taubinPassBand)
        smoothing.SetInputConnection(surface)
        surface = smoothing.GetOutputPort()

    if state.normals:
      normals = vtk.vtkPolyDataNormals()
      normals.AutoOrientNormalsOn()
      normals.SetFlipNormals(state.flipNormals)
      normals.SetSplitting(state.splitting)
      normals.SetFeatureAngle(state.featureAngle)
      normals.ConsistencyOn()
      normals.SetInputConnection(surface)
      surface = normals.GetOutputPort()

    if state.cleaner:
      cleaner = vtk.vtkCleanPolyData()
      cleaner.SetInputConnection(surface)
      surface = cleaner.GetOutputPort()

    if state.connectivity:
      connectivity = vtk.vtkPolyDataConnectivityFilter()
      connectivity.SetExtractionModeToLargestRegion()
      connectivity.SetInputConnection(surface)
      surface = connectivity.GetOutputPort()

    state.outputModelNode.SetPolyDataConnection(surface)
    return True

def write_stl(ugrid, filename):
    surface_filter = vtk.vtkDataSetSurfaceFilter()
    surface_filter.SetInputData(ugrid)

    triangle_filter = vtk.vtkTriangleFilter()
    triangle_filter.SetInputConnection(surface_filter.GetOutputPort())

    writer = vtk.vtkSTLWriter()
    writer.SetFileName(filename)
    writer.SetInputConnection(triangle_filter.GetOutputPort())
    writer.Write()

def Merge(polydata_list):
    append = vtk.vtkAppendPolyData()

    for polydata in polydata_list:
        triangle = vtk.vtkTriangleFilter()
        triangle.SetInputData(polydata)
        append.AddInputData(triangle.GetOutput())

    clean = vtk.vtkCleanPolyData()
    clean.SetInputData(append.GetOutput())

    return append.GetOutput()

    def Execute(self):
        from vmtk import vmtkscripts
        if self.Surface == None:
            self.PrintError('Error: no Surface.')

        if not self.vmtkRenderer:
            self.vmtkRenderer = vmtkrenderer.vmtkRenderer()
            self.vmtkRenderer.Initialize()
            self.OwnRenderer = 1

        self.vmtkRenderer.RegisterScript(self) 

        triangleFilter = vtk.vtkTriangleFilter()
        triangleFilter.SetInputData(self.Surface)
        triangleFilter.Update()

        self.Surface = triangleFilter.GetOutput()

        self.tagviewer = vmtkscripts.vmtkSurfaceViewer()
        self.tagviewer.Surface = self.Surface
        self.tagviewer.vmtkRenderer = self.vmtkRenderer
        self.tagviewer.Representation = 'edges'
        self.tagviewer.Opacity = self.Opacity
        self.tagviewer.Execute()

        self.ContourWidget = vtk.vtkContourWidget()
        self.ContourWidget.SetInteractor(self.vmtkRenderer.RenderWindowInteractor)

        rep = vtk.vtkOrientedGlyphContourRepresentation.SafeDownCast(self.ContourWidget.GetRepresentation())
        rep.GetLinesProperty().SetColor(1, 0, 0)
        rep.GetLinesProperty().SetLineWidth(4.0)

        pointPlacer = vtk.vtkPolygonalSurfacePointPlacer()
        pointPlacer.AddProp(self.tagviewer.Actor)
        pointPlacer.GetPolys().AddItem(self.Surface)
        rep.SetPointPlacer(pointPlacer)

        self.Interpolator = vtk.vtkPolygonalSurfaceContourLineInterpolator()
        self.Interpolator.GetPolys().AddItem(self.Surface)
        rep.SetLineInterpolator(self.Interpolator)

        self.InputInfo("Building loop ...\n")

        self.vmtkRenderer.AddKeyBinding('space','Generate loop',self.LoopCallback)
        self.vmtkRenderer.AddKeyBinding('d','Delete contour',self.DeleteContourCallback)
        self.vmtkRenderer.AddKeyBinding('i','Start/stop contour drawing',self.InteractCallback)


        self.Display()


        if self.OwnRenderer:
            self.vmtkRenderer.Deallocate()

    def prepareModel(self, polyData):
        '''
        '''
        # import the vmtk libraries
        try:
            import vtkvmtkComputationalGeometryPython as vtkvmtkComputationalGeometry
            import vtkvmtkMiscPython as vtkvmtkMisc
        except ImportError:
            logging.error("Unable to import the SlicerVmtk libraries")

        capDisplacement = 0.0

        surfaceCleaner = vtk.vtkCleanPolyData()
        surfaceCleaner.SetInputData(polyData)
        surfaceCleaner.Update()

        surfaceTriangulator = vtk.vtkTriangleFilter()
        surfaceTriangulator.SetInputData(surfaceCleaner.GetOutput())
        surfaceTriangulator.PassLinesOff()
        surfaceTriangulator.PassVertsOff()
        surfaceTriangulator.Update()

        # new steps for preparation to avoid problems because of slim models (f.e. at stenosis)
        subdiv = vtk.vtkLinearSubdivisionFilter()
        subdiv.SetInputData(surfaceTriangulator.GetOutput())
        subdiv.SetNumberOfSubdivisions(1)
        subdiv.Update()

        smooth = vtk.vtkWindowedSincPolyDataFilter()
        smooth.SetInputData(subdiv.GetOutput())
        smooth.SetNumberOfIterations(20)
        smooth.SetPassBand(0.1)
        smooth.SetBoundarySmoothing(1)
        smooth.Update()

        normals = vtk.vtkPolyDataNormals()
        normals.SetInputData(smooth.GetOutput())
        normals.SetAutoOrientNormals(1)
        normals.SetFlipNormals(0)
        normals.SetConsistency(1)
        normals.SplittingOff()
        normals.Update()

        surfaceCapper = vtkvmtkComputationalGeometry.vtkvmtkCapPolyData()
        surfaceCapper.SetInputData(normals.GetOutput())
        surfaceCapper.SetDisplacement(capDisplacement)
        surfaceCapper.SetInPlaneDisplacement(capDisplacement)
        surfaceCapper.Update()

        outPolyData = vtk.vtkPolyData()
        outPolyData.DeepCopy(surfaceCapper.GetOutput())

        return outPolyData

    def _readstls(self):
        for fname in sorted(os.listdir(self._ref_dir)):
            if '.stl' in fname.lower():
                reader = vtk.vtkSTLReader()
                reader.SetFileName(self._ref_dir+'/'+fname)
                reader.Update()
                triangles = vtk.vtkTriangleFilter()
                triangles.SetInputConnection(reader.GetOutputPort())
                triangles.Update()
                self.rsurfs.append(triangles.GetOutput())

                dl = vtk.vtkDelaunay3D()
                dl.SetInputConnection(triangles.GetOutputPort())
                dl.Update()
                self.rmeshes.append(dl)

                vol, cent, axes = self._getMassProps(self.rmeshes[-1])
                self.rvols.append(vol)
                self.rcentroids.append(cent)
                self.raxes.append(axes)

        for fname in sorted(os.listdir(self._def_dir)):
            if '.stl' in fname.lower():
                reader = vtk.vtkSTLReader()
                reader.SetFileName(self._def_dir+'/'+fname)
                reader.Update()
                triangles = vtk.vtkTriangleFilter()
                triangles.SetInputConnection(reader.GetOutputPort())
                triangles.Update()
                self.dsurfs.append(triangles.GetOutput())

                dl = vtk.vtkDelaunay3D()
                dl.SetInputConnection(triangles.GetOutputPort())
                dl.Update()
                self.dmeshes.append(dl)

                vol, cent, axes = self._getMassProps(self.dmeshes[-1])
                self.dvols.append(vol)
                self.dcentroids.append(cent)
                self.daxes.append(axes)