def get_topography_actor(x, y, topography, super_elevation, decimation,
                         smoothing=False):
    topography = topography.T
    topography *= super_elevation

    points, triangles, colors = setup_topography(x, y, topography,
                                                 decimation=decimation)
    # Create a polydata object
    trianglePolyData = vtk.vtkPolyData()

    # Add the geometry and topology to the polydata
    trianglePolyData.SetPoints(points)
    trianglePolyData.GetPointData().SetScalars(colors)
    trianglePolyData.SetPolys(triangles)

    # Clean the polydata so that the edges are shared !
    cleanPolyData = vtk.vtkCleanPolyData()
    cleanPolyData.SetInput(trianglePolyData)

    # Use a filter to smooth the data (will add triangles and smooth)
    # Use two different filters to show the difference
    mapper = vtk.vtkPolyDataMapper()

    if smoothing:
        smooth_loop = vtk.vtkLoopSubdivisionFilter()
        smooth_loop.SetNumberOfSubdivisions(2)
        smooth_loop.SetInputConnection(cleanPolyData.GetOutputPort())
        mapper.SetInputConnection(smooth_loop.GetOutputPort())

    else:
        mapper.SetInput(trianglePolyData)
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)

    return actor

    def ReadOBJFile(self, filename):
        '''
        @param filename: str 
        @rtype: None
        '''
        reader = vtk.vtkOBJReader()
        reader.SetFileName(filename)
        
        try:
            reader.Update()
            cleanFilter = vtk.vtkCleanPolyData()
            cleanFilter.SetInput(reader.GetOutput())
            cleanFilter.ConvertLinesToPointsOn()
            cleanFilter.ConvertPolysToLinesOn()
            cleanFilter.ConvertStripsToPolysOn()
            cleanFilter.PointMergingOn()
            cleanFilter.Update()
            surfaceFilter = vtk.vtkDataSetSurfaceFilter()
            surfaceFilter.SetInput(cleanFilter.GetOutput())
            surfaceFilter.Update()
            self.setDataSet(surfaceFilter.GetOutput())
            del cleanFilter
            del surfaceFilter
        except Exception, e:
#            del reader
            print e
            raise IOError, "Could not read the OBJ file! "

def cleanMesh(mesh, connectivityFilter=False):
    try:
        t = time.clock()
        connect = vtk.vtkPolyDataConnectivityFilter()
        clean = vtk.vtkCleanPolyData()

        if (connectivityFilter):
            if vtk.vtkVersion.GetVTKMajorVersion() >= 6:
                connect.SetInputData( mesh )
            else:
                connect.SetInput( mesh )
            connect.SetExtractionModeToLargestRegion()
            clean.SetInputConnection( connect.GetOutputPort() )
        else:
            if vtk.vtkVersion.GetVTKMajorVersion() >= 6:
                clean.SetInputData( mesh )
            else:
                clean.SetInput( mesh )

        clean.Update()
        print ("Surface cleaned")
        m2 = clean.GetOutput()
        print "    ", m2.GetNumberOfPolys(), "polygons"
        elapsedTime(t)
        clean = None
        connect = None
        return m2
    except:
        print "Surface cleaning failed"
        exc_type, exc_value, exc_traceback = sys.exc_info()
        traceback.print_exception(exc_type, exc_value, exc_traceback, limit=2, file=sys.stdout)
    return None

    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 create_mesh(M,P,T):
	n_naca_pts = 50
	write_test_file(M,P,T,-5.,5.,nsections=5)
	wing=create_wing('current_wing','output')
	n_sections=len(wing)
	n_section_pts = 2*n_naca_pts-1

	# build mesh
	vtk_model = vtk.vtkStructuredGrid()
	vtk_model.SetDimensions(n_section_pts,n_sections,1)
	# build points
	vtk_points = vtk.vtkPoints()

	for j in xrange(n_sections):
	    upper_pts = numpy.array([wing[j][1],wing[j][3]]).T
	    lower_pts = numpy.array([wing[j][2],wing[j][4]]).T
	    section_pts = numpy.concatenate((lower_pts[::-1],upper_pts[1:]))
	    for i in xrange(n_section_pts):
		vtk_points.InsertNextPoint(section_pts[i,0],wing[j][0],section_pts[i,1])
	# set points
	vtk_model.SetPoints(vtk_points)

	# convert to poly data    
	pdata_filter = vtk.vtkGeometryFilter()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    pdata_filter.SetInput(vtk_model)
	else:
	    pdata_filter.SetInputData(vtk_model)
	pdata_filter.Update()
	poly_data = pdata_filter.GetOutput()

	# compute normals
	norms = vtk.vtkPolyDataNormals()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    norms.SetInput(poly_data)
	else:
	    norms.SetInputData(poly_data)
	norms.ComputePointNormalsOff()
	norms.ComputeCellNormalsOn()
	norms.ConsistencyOn()
	norms.Update()

	# clean poly data
	clean_poly = vtk.vtkCleanPolyData()
	clean_poly.ToleranceIsAbsoluteOn()
	clean_poly.SetAbsoluteTolerance(1.e-6)
	if vtk.VTK_MAJOR_VERSION <= 5:
	    clean_poly.SetInput(norms.GetOutput())
	else:
	    clean_poly.SetInputData(norms.GetOutput())
	clean_poly.Update()

	# write output mesh
	writer = vtk.vtkXMLPolyDataWriter()
	if vtk.VTK_MAJOR_VERSION <= 5:
	    writer.SetInput(clean_poly.GetOutput())
	else:
	    writer.SetInputData(clean_poly.GetOutput())
	writer.SetFileName('output.vtp')
	writer.Write()

def VertFacetoPoly(new_pt, new_fc):
    """ Creates a vtk polydata object given points and triangular faces """
    
    # Convert points to vtkfloat object
    points = np.vstack(new_pt)
    vtkArray = VN.numpy_to_vtk(np.ascontiguousarray(points), deep=True)#, deep=True) 
    points = vtk.vtkPoints()
    points.SetData(vtkArray) 
    
    # Convert faces to vtk cells object
    ints = np.ones(len(new_fc), 'int')*3
    cells = np.hstack((ints.reshape(-1, 1), np.vstack(new_fc)))
    cells = np.ascontiguousarray(np.hstack(cells).astype('int64'))
    vtkcells = vtk.vtkCellArray()
    vtkcells.SetCells(cells.shape[0], VN.numpy_to_vtkIdTypeArray(cells, deep=True))

    # Create polydata object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(points)
    pdata.SetPolys(vtkcells)
    
    # Remove duplicate verticies
    clean = vtk.vtkCleanPolyData()
    clean.ConvertPolysToLinesOff()
    clean.ConvertLinesToPointsOff()
    clean.ConvertStripsToPolysOff()
    if vtk.vtkVersion().GetVTKMajorVersion() > 5:
        clean.SetInputData(pdata)
    else:
        clean.SetInput(pdata) 
    clean.Update()

    return clean.GetOutput()

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 __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkCleanPolyData(), 'Processing.',
         ('vtkPolyData',), ('vtkPolyData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

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 Execute(self):

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

        if (self.GroupId != -1) & (self.GroupIdsArrayName!=''):
            self.Surface.GetPointData().SetActiveScalars(self.GroupIdsArrayName)

        connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
        connectivityFilter.SetInput(self.Surface)
        connectivityFilter.ColorRegionsOff()       
        connectivityFilter.SetExtractionModeToLargestRegion()
        if self.GroupId != -1:
            connectivityFilter.ScalarConnectivityOn()
            scalarRange = [self.GroupId,self .GroupId]
            connectivityFilter.SetScalarRange(scalarRange)
        connectivityFilter.Update()

        self.Surface = connectivityFilter.GetOutput()
	
      	if self.CleanOutput == 1:
      	    cleaner = vtk.vtkCleanPolyData()
      	    cleaner.SetInput(connectivityFilter.GetOutput())
      	    cleaner.Update()

            self.Surface = cleaner.GetOutput()

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

    def Execute(self):

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

        if self.WidgetType == "box":
            self.ClipFunction = vtk.vtkPlanes()
        elif self.WidgetType == "sphere":
            self.ClipFunction = vtk.vtkSphere()

        self.Clipper = vtk.vtkClipPolyData()
        self.Clipper.SetInput(self.Surface)
        self.Clipper.SetClipFunction(self.ClipFunction)
        self.Clipper.GenerateClippedOutputOn()
        self.Clipper.InsideOutOn()
        
        if not self.vmtkRenderer:
            self.vmtkRenderer = vmtkrenderer.vmtkRenderer()
            self.vmtkRenderer.Initialize()
            self.OwnRenderer = 1

        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInput(self.Surface)
        mapper.ScalarVisibilityOff()
        self.Actor = vtk.vtkActor()
        self.Actor.SetMapper(mapper)
        self.vmtkRenderer.Renderer.AddActor(self.Actor)

        if self.WidgetType == "box":
            self.ClipWidget = vtk.vtkBoxWidget()
            self.ClipWidget.GetFaceProperty().SetColor(0.6,0.6,0.2)
            self.ClipWidget.GetFaceProperty().SetOpacity(0.25)
        elif self.WidgetType == "sphere":
            self.ClipWidget = vtk.vtkSphereWidget()
            self.ClipWidget.GetSphereProperty().SetColor(0.6,0.6,0.2)
            self.ClipWidget.GetSphereProperty().SetOpacity(0.25)
            self.ClipWidget.GetSelectedSphereProperty().SetColor(0.6,0.0,0.0)
            self.ClipWidget.GetSelectedSphereProperty().SetOpacity(0.75)
            self.ClipWidget.SetRepresentationToSurface()
            self.ClipWidget.SetPhiResolution(20)
            self.ClipWidget.SetThetaResolution(20)

        self.ClipWidget.SetInteractor(self.vmtkRenderer.RenderWindowInteractor)
        self.Display()

        self.Transform = vtk.vtkTransform()
        self.ClipWidget.GetTransform(self.Transform)

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

        if self.CleanOutput == 1:
            cleaner = vtk.vtkCleanPolyData()
            cleaner.SetInput(self.Surface)
            cleaner.Update()
            self.Surface = cleaner.GetOutput()

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

 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 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 cleanPolyData(polyData):
    cleanFilt = vtk.vtkCleanPolyData()
    cleanFilt.SetConvertLinesToPoints(0)
    cleanFilt.SetConvertPolysToLines(0)
    cleanFilt.SetConvertStripsToPolys(0)
    cleanFilt.SetPointMerging(1)
    cleanFilt.SetInputData(polyData)
    cleanFilt.Update()
    return cleanFilt.GetOutput()

def computeComponents(vertices, faces):
    pts, tris = arrayToPolydata(vertices, faces)
    polys = vtk.vtkPolyData()
    polys.SetPoints(pts)
    polys.SetPolys(tris)
    polys.Update()
    
    cleanFilter = vtk.vtkCleanPolyData()
    #cleanFilter.PointMergingOff()
    cleanFilter.ConvertStripsToPolysOff()
    cleanFilter.ConvertPolysToLinesOff()
    cleanFilter.ConvertLinesToPointsOff()
    cleanFilter.SetInput(polys)
    cleanFilter.Update()
    
    connFilter = vtk.vtkPolyDataConnectivityFilter()
    connFilter.ScalarConnectivityOff()
    connFilter.SetExtractionModeToAllRegions()
    connFilter.SetInput(cleanFilter.GetOutput())
    connFilter.Update()
    
    nregions = connFilter.GetNumberOfExtractedRegions()
    print "Regions extracted: %d" % (nregions)
    
    components = []
    for i in range(nregions):
        connFilter = vtk.vtkPolyDataConnectivityFilter()
        connFilter.ScalarConnectivityOff()
        connFilter.SetExtractionModeToSpecifiedRegions()
        connFilter.AddSpecifiedRegion(i)
        connFilter.SetInput(cleanFilter.GetOutput())
        connFilter.Update()
    
        cFilter = vtk.vtkCleanPolyData()
#        cFilter.ConvertStripsToPolysOff()
#        cFilter.ConvertPolysToLinesOff()
#        cFilter.ConvertLinesToPointsOff()
        cFilter.SetInput(connFilter.GetOutput())
        cFilter.Update()
    
        v,f = polydataToArray(cFilter.GetOutput())
        components.append([v,f])
        
    return components

    def __init__(self, module_manager):
        # initialise our base class
        ModuleBase.__init__(self, module_manager)
        # initialise any mixins we might have
        NoConfigModuleMixin.__init__(self)

        mm = self._module_manager

        self._cleaner = vtk.vtkCleanPolyData()

        self._tf = vtk.vtkTriangleFilter()
        self._tf.SetInput(self._cleaner.GetOutput())

        self._wspdf = vtk.vtkWindowedSincPolyDataFilter()
        # self._wspdf.SetNumberOfIterations(50)
        self._wspdf.SetInput(self._tf.GetOutput())
        self._wspdf.SetProgressText("smoothing")
        self._wspdf.SetProgressMethod(lambda s=self, mm=mm: mm.vtk_progress_cb(s._wspdf))

        self._cleaner2 = vtk.vtkCleanPolyData()
        self._cleaner2.SetInput(self._wspdf.GetOutput())

        self._curvatures = vtk.vtkCurvatures()
        self._curvatures.SetCurvatureTypeToMean()
        self._curvatures.SetInput(self._cleaner2.GetOutput())

        self._tf.SetProgressText("triangulating")
        self._tf.SetProgressMethod(lambda s=self, mm=mm: mm.vtk_progress_cb(s._tf))

        self._curvatures.SetProgressText("calculating curvatures")
        self._curvatures.SetProgressMethod(lambda s=self, mm=mm: mm.vtk_progress_cb(s._curvatures))

        self._inputFilter = self._tf

        self._inputPoints = None
        self._inputPointsOID = None
        self._giaGlenoid = None
        self._outsidePoints = None
        self._outputPolyDataARB = vtk.vtkPolyData()
        self._outputPolyDataHM = vtk.vtkPolyData()

        self._createViewFrame("Test Module View", {"vtkTriangleFilter": self._tf, "vtkCurvatures": self._curvatures})

        self._viewFrame.Show(True)

    def CentreLinePolyData(self):
        """Compute centrelines based on the profile, reusing our
        memoed copy or reading from the cache file if possible.
        """
        if (os.path.exists(self.CentreLineFile ) and 
            os.path.getmtime(self.CentreLineFile ) > os.path.getmtime(self.StlFile) and
            os.path.getmtime(self.CentreLineFile ) > os.path.getmtime(self.FileName)):
            # Cached!
            reader = vtk.vtkXMLPolyDataReader()
            reader.SetFileName(self.CentreLineFile)
            reader.Update()
            return reader.GetOutput()

        # Have to compute it
        
        # Read the STL file
        reader = vtk.vtkSTLReader()
        reader.SetFileName(profile.StlFile)
        
        # Find the seed points for centreline calculation
        # Use points one iolet radius back along the normal.
        
        outletPts = []
        def scale(iolet):
            pt = (iolet.Centre - iolet.Radius * iolet.Normal)
            pt = pt / self.LengthUnit
            return pt.magnitude
        
        for iolet in self.Iolets:
            if isinstance(iolet._iolet, Inlet):
                inletPt = scale(iolet)
            else:
                outletPts.append(scale(iolet))
                pass
            continue
        
        srcPts, tgtPts = FindSeeds(reader, inletPt, outletPts)
        
        # Lazy import since it's so slow!
        from vmtk import vtkvmtk
        centreliner = vtkvmtk.vtkvmtkPolyDataCenterlines()
        centreliner.SetInputConnection(reader.GetOutputPort())
    
        centreliner.SetSourceSeedIds(srcPts)
        centreliner.SetTargetSeedIds(tgtPts)
        centreliner.SetRadiusArrayName("radius")
        centreliner.SetCostFunction("1/R")

        cleaner = vtk.vtkCleanPolyData()
        cleaner.SetInputConnection(centreliner.GetOutputPort())
        
        writer = vtk.vtkXMLPolyDataWriter()
        writer.SetInputConnection(cleaner.GetOutputPort())
        writer.SetFileName(self.CentreLineFile)
        writer.Write()
        return cleaner.GetOutput()

def trim_mesh_with_cone(mesh, cone_point, cone_normal, cone_radius):
    """
    returns a mesh that contains only the triangles that where inside of the cone
    """
    cone_end = cone_point + cone_normal*100.0
    sq_radius = cone_radius*cone_radius
    w = cone_point
    v = cone_end
    n = w - v
    v_w_sq_len = dist2(v, w)
   
    points = mesh.GetPoints().GetData()
    cell_array = mesh.GetPolys()
    polygons = cell_array.GetData()
    triangles = vtk.vtkCellArray()
    for i in xrange(0,  cell_array.GetNumberOfCells()):
        triangle = [polygons.GetValue(j) for j in xrange(i*4+1, i*4+4)]
        
        p = points.GetTuple(triangle[0])
        delta = p - (v + (((p - v).dot(n)) / v_w_sq_len) * n)
        if delta.dot(delta) > sq_radius:
            continue
        
        p = points.GetTuple(triangle[1])
        delta = p - (v + (((p - v).dot(n)) / v_w_sq_len) * n)
        if delta.dot(delta) > sq_radius:
            continue
        
        p = points.GetTuple(triangle[2])
        delta = p - (v + (((p - v).dot(n)) / v_w_sq_len) * n)
        if delta.dot(delta) > sq_radius:
            continue
        
        cell = vtk.vtkTriangle()
        pointIds = cell.GetPointIds()
        pointIds.SetId(0, triangle[0])
        pointIds.SetId(1, triangle[1])
        pointIds.SetId(2, triangle[2])
        triangles.InsertNextCell(cell)
            
    # Create a polydata object
    trianglePolyData = vtk.vtkPolyData()

    # Add the geometry and topology to the polydata
    trianglePolyData.SetPoints(mesh.GetPoints())
    trianglePolyData.SetPolys(triangles)
    trianglePolyData.Update()
    
    #run the clean function here to remove the points that are not used
    cleanPolyData = vtk.vtkCleanPolyData()
    cleanPolyData.SetInput(trianglePolyData)
    cleanPolyData.Update()
    trimmed_mesh = cleanPolyData.GetOutput()

    return trimmed_mesh

    def DeformSurface(self):
        # interpolate and sample the displacement norms over the surface
        rbf = vtkvmtkcontrib.vtkvmtkRBFInterpolation2()
        rbf.SetSource(self.SourceSpheres)
        rbf.SetRBFTypeToBiharmonic()
        rbf.ComputeCoefficients()
        sampler = vtkvmtkcontrib.vtkvmtkPolyDataSampleFunction()
        sampler.SetInput(self.Surface)
        sampler.SetImplicitFunction(rbf)
        sampler.SetSampleArrayName("DisplacementNorms")
        sampler.Update()

        sampArray = sampler.GetOutput().GetPointData().GetArray("DisplacementNorms")

        ##Clamp the negative values to 0 and the positive values to one in a weight array
        calculator = vtk.vtkArrayCalculator()
        calculator.SetInput(sampler.GetOutput())
        calculator.AddScalarArrayName("DisplacementNorms")
        calculator.SetFunction("if( DisplacementNorms > 0 , iHat, jHat)")
        calculator.SetResultArrayName("Weights")
        calculator.SetResultArrayType(vtk.VTK_FLOAT)
        calculator.Update()

        # Create the transform
        thinPlateSplineTransform = vtk.vtkThinPlateSplineTransform()
        thinPlateSplineTransform.SetBasisToR()
        thinPlateSplineTransform.SetSourceLandmarks(self.SourcePoints)
        thinPlateSplineTransform.SetTargetLandmarks(self.TargetPoints)

        transform = vtk.vtkTransform()
        transform.Identity()
        transform2 = vtk.vtkTransform()
        transform2.Identity()

        # Apply weighted transform
        transformFilter = vtk.vtkWeightedTransformFilter()
        transformFilter.SetInput(calculator.GetOutput())
        transformFilter.SetNumberOfTransforms(3)
        transformFilter.SetWeightArray("Weights")
        transformFilter.SetTransform(thinPlateSplineTransform, 0)
        transformFilter.SetTransform(transform, 1)
        transformFilter.SetTransform(transform2, 2)
        transformFilter.Update()

        normalsFilter = vtk.vtkPolyDataNormals()
        normalsFilter.SetInput(transformFilter.GetOutput())
        normalsFilter.Update()

        # FIXME: the normal filter apparently introduced some holes in some meshes (wtf?). This filter cleans the mesh
        cleanFilter = vtk.vtkCleanPolyData()
        cleanFilter.SetInput(normalsFilter.GetOutput())
        cleanFilter.Update()

        self.DeformedSurface = cleanFilter.GetOutput()

def concomp(opts, argv):
	p = nv.readVTK(argv[0])
	c = vtk.vtkPolyDataConnectivityFilter()
	c.SetInput(p)
	c.SetExtractionModeToLargestRegion()
	c.Update()
	d = vtk.vtkCleanPolyData()
	d.SetInput(c.GetOutput())
	d.Update()
	p = d.GetOutput()
	nv.writeVTK(argv[1],p)