def getCell2vtp(vtkObj,ind):
    """
        Function gets a cell by ind and constructs a polydata from it.

    """

    # Get the cell
    cE = vtkObj.GetCell(ind)
    # Make the polygon
    if cE.GetCellType() == 11:
        # Use a cubeSource, much faster
        cube = vtk.vtkCubeSource()
        cube.SetBounds(cE.GetBounds())
        cube.Update()
        vtpObj = cube.GetOutput()
    else:
        polygons = vtk.vtkCellArray()
        for iF in range(cE.GetNumberOfFaces()):
            f = cE.GetFace(iF)
            poly = vtk.vtkPolygon()
            poly.GetPointIds().SetNumberOfIds(f.GetNumberOfPoints())
            for nr in range(f.GetNumberOfPoints()):
                poly.GetPointIds().SetId(nr,f.GetPointId(nr))
            polygons.InsertNextCell(poly)
        # Build the polydata
        vtpObj = vtk.vtkPolyData()
        vtpObj.SetPoints(obj.GetPoints())
        vtpObj.SetPolys(polygons)

    return polydata.normFilter(polydata.triangulatePolyData(vtpObj))

def polygons(indices):
    """
    Maps a numpy ndarray to an vtkCellArray of vtkPolygons

    Args:
        indices (numpy.ndarray<int>): A numpy.ndarray of shape (n,m) of indices that define n polygons with m points each

    Returns:
        vtk_polygons (vtk.vtkCellArray): VTK representation of the polygons
    """
    if not isinstance(indices, numpy.ndarray):
        raise Numpy2VtkFormatException(
            'polygons needs numpy array as input'
        )
    if len(indices.shape) != 2:
        raise Numpy2VtkFormatException(
            'polygons needs a nxm ndarray as input'
        )
    if indices.dtype != numpy.int:
        raise Numpy2VtkFormatException(
            'polygons needs to be numpy array of type numpy.int'
        )

    number_of_polygons = indices.shape[0]
    poly_shape = indices.shape[1]
    vtk_polygons = vtk.vtkCellArray()
    for j in range(0, number_of_polygons):
        polygon = vtk.vtkPolygon()
        polygon.GetPointIds().SetNumberOfIds(poly_shape)
        for i in range(0, poly_shape):
            polygon.GetPointIds().SetId(i, indices[j, i])
        vtk_polygons.InsertNextCell(polygon)
    return vtk_polygons

    def _update_vtk_objects(self):
        """When n is changed the thus the number of coordinates this function is needed
        to update the vtk objects with the new number of points."""
        # self._vtk_points.SetNumberOfPoints(len(self._points))
        # for i, c in enumerate(self._points):
        #     self._vtk_points.InsertPoint(i, c[0], c[1], c[2])
        self._vtk_points = _vtk.vtkPoints()
        for coordinates in self._points:
            self._vtk_points.InsertNextPoint(coordinates[0], coordinates[1], coordinates[2])

        self._vtk_polygons = _vtk.vtkCellArray()
        for polygon in self._polygons:
            vtk_polygon = _vtk.vtkPolygon()
            vtk_polygon.GetPointIds().SetNumberOfIds(3)
            for local_index, global_index in enumerate(polygon):
                vtk_polygon.GetPointIds().SetId(local_index, global_index)
            self._vtk_polygons.InsertNextCell(vtk_polygon)

        self._vtk_poly_data.SetPoints(self._vtk_points)
        self._vtk_poly_data.SetPolys(self._vtk_polygons)

        self._vtk_scalars = _vtk.vtkFloatArray()
        self._vtk_scalars.SetNumberOfValues(self._vtk_poly_data.GetPoints().GetNumberOfPoints())
        for i in range(self._vtk_scalars.GetNumberOfTuples()):
            self._vtk_scalars.SetValue(i, 0.)

        self._vtk_poly_data.GetPointData().SetScalars(self._vtk_scalars)
        self._vtk_poly_data.Modified()

def gen_outline(pts,color,size):
	'''
	Returns an outline actor with specified pts, color and size. Incoming pnts should be ordered.
	'''
	if color[0]<=1 and color != None:
		color=(int(color[0]*255),int(color[1]*255),int(color[2]*255))
	if color[0]>=1 and color != None:
		color=(color[0]/float(255),color[1]/float(255),color[2]/float(255))
	points=vtk.vtkPoints()

	for i in pts:
		points.InsertNextPoint(i)

	lineseg=vtk.vtkPolygon()
	lineseg.GetPointIds().SetNumberOfIds(len(pts))
	for i in range(len(pts)):
		lineseg.GetPointIds().SetId(i,i)
	linesegcells=vtk.vtkCellArray()
	linesegcells.InsertNextCell(lineseg)
	outline=vtk.vtkPolyData()
	outline.SetPoints(points)
	outline.SetVerts(linesegcells)
	outline.SetLines(linesegcells)
	Omapper=vtk.vtkPolyDataMapper()
	Omapper.SetInputData(outline)
	outlineActor=vtk.vtkActor()
	outlineActor.SetMapper(Omapper)
	outlineActor.GetProperty().SetColor(color)
	outlineActor.GetProperty().SetPointSize(size)
	return outlineActor,outline #actor/polydata

def gocad2vtp(gcFile):
    """"
    Function to read gocad polystructure file and makes VTK Polydata object (vtp).

    Input:
        gcFile: gocadFile with polysturcture

    """
    print "Reading GOCAD ts file..."
    vrtx, trgl = read_GOCAD_ts(gcFile)
    # Adjust the index
    trgl = trgl - 1

    # Make vtk pts
    ptsvtk = vtk.vtkPoints()
    ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1))

    # Make the polygon connection
    polys = vtk.vtkCellArray()
    for face in trgl:
        poly = vtk.vtkPolygon()
        poly.GetPointIds().SetNumberOfIds(len(face))
        for nrv, vert in enumerate(face):
            poly.GetPointIds().SetId(nrv,vert)
        polys.InsertNextCell(poly)

    # Make the polydata, structure of connections and vrtx
    polyData = vtk.vtkPolyData()
    polyData.SetPoints(ptsvtk)
    polyData.SetPolys(polys)

    return polyData

    def draw(self, graphics):
        cell, pointnums = super(Polygon, self).draw(graphics)

        polygon = vtk.vtkPolygon()
        polygon.GetPointIds().SetNumberOfIds(len(pointnums))
        for i, p in enumerate(pointnums):
            polygon.GetPointIds().SetId(i, p)
        cell.InsertNextCell(polygon)
        graphics.data.SetPolys(cell)

 def _createpolygon(self, pointnumbers):
     if depth(pointnumbers) >= 3:
         for p in pointnumbers:
             self._createpolygon(p)
     else:
         polygon = vtk.vtkPolygon()
         polygon.GetPointIds().SetNumberOfIds(len(pointnumbers))
         i = 0
         for p in pointnumbers:
             polygon.GetPointIds().SetId(i, p)
             i += 1
         self.polygons.InsertNextCell(polygon)

    def generateVolumetricLine(self):
        """ Generates volumetric lines. """

        if self.centroid is None:
            return

        points = vtk.vtkPoints()
        poly_s = vtk.vtkPolygon()

        tangent_dir = self.calculateElectrodePolygon(points, poly_s)

        self.polygonActor = self.prepareActorWithShadersAndMapper(points, poly_s, tangent_dir)

def arrayToPolydata(verts, faces):
    pts = vtk.vtkPoints()
    tris = vtk.vtkCellArray()

    for v in verts:
        pts.InsertNextPoint( [v[0], v[1], v[2]] )
    for f in faces:
        poly = vtk.vtkPolygon()
        poly.GetPointIds().InsertId(0,f[0])
        poly.GetPointIds().InsertId(1,f[1])
        poly.GetPointIds().InsertId(2,f[2])
        tris.InsertNextCell(poly.GetPointIds())
        
    return pts, tris

def prepFillarea(renWin,ren,farea,cmap=None):
  n = prepPrimitive(farea)
  if n==0:
    return
  for i in range(n):
    x   = farea.x[i]
    y   = farea.y[i]
    c   = farea.color[i]
    st  = farea.style[i]
    idx = farea.index[i]
    N = max(len(x),len(y))
    for a in [x,y]:
      while len(a)<n:
        a.append(a[-1])
    #Create points
    pts = vtk.vtkPoints()
    for j in range(N):
      pts.InsertNextPoint(x[j],y[j],0.)
    #Create polygon out of these points
    polygon = vtk.vtkPolygon()
    pid = polygon.GetPointIds()
    pid.SetNumberOfIds(N)
    for j in range(N):
      pid.SetId(j,j)
    polygons = vtk.vtkCellArray()
    polygons.InsertNextCell(polygon)

    polygonPolyData = vtk.vtkPolyData()
    geo,pts = project(pts,farea.projection,farea.worldcoordinate)
    polygonPolyData.SetPoints(pts)
    polygonPolyData.SetPolys(polygons)

    a = vtk.vtkActor()
    m = vtk.vtkPolyDataMapper()
    m.SetInputData(polygonPolyData)
    a.SetMapper(m)
    p = a.GetProperty()

    if cmap is None:
      if farea.colormap is not None:
        cmap = farea.colormap
      else:
        cmap = 'default'
    if isinstance(cmap,str):
      cmap = vcs.elements["colormap"][cmap]
    color = cmap.index[c]
    p.SetColor([C/100. for C in color])
    ren.AddActor(a)
    fitToViewport(a,ren,farea.viewport,wc=farea.worldcoordinate,geo=geo)
  return

def genPoly(coords,pts,filled=True):
  N = pts.GetNumberOfPoints()
  if filled:
    poly = vtk.vtkPolygon()
  else:
    poly = vtk.vtkPolyLine()
  pid = poly.GetPointIds()
  n = len(coords)
  pid.SetNumberOfIds(n)
  for j in range(n):
    c = list(coords[j])
    if len(c)==2:
      c.append(0)
    pts.InsertNextPoint(*c)
    pid.SetId(j,j+N)
  return poly

    def _square_slice(self):
        """Return a polydata of a square slice."""
        polygons = vtk.vtkCellArray()
        for i in range(self._side-1):
            for j in range(self._side-1):
                corners = [(i, j), (i+1, j), (i+1, j+1), (i, j+1)]
                polygon = vtk.vtkPolygon()
                polygon.GetPointIds().SetNumberOfIds(4)
                for index, corner in enumerate(corners):
                    polygon.GetPointIds().SetId(index, corner[0]*self._side+corner[1])
                polygons.InsertNextCell(polygon)

        template_poly_data = vtk.vtkPolyData()
        template_poly_data.SetPoints(self._points)
        template_poly_data.GetPointData().SetScalars(self._image_values)
        template_poly_data.SetPolys(polygons)
        return template_poly_data

    def _square_slice(self):
        """Call in the beginning. Precalculates the polydata object
        without rotation."""
        self._polygons.Initialize()
        for i in range(self._side-1):
            for j in range(self._side-1):
                corners = [(i, j), (i+1, j), (i+1, j+1), (i, j+1)]
                polygon = _vtk.vtkPolygon()
                polygon.GetPointIds().SetNumberOfIds(4)
                for index, corner in enumerate(corners):
                    polygon.GetPointIds().SetId(index,
                                                corner[0]*self._side+corner[1])
                self._polygons.InsertNextCell(polygon)

        self._template_poly_data = _vtk.vtkPolyData()
        self._template_poly_data.SetPoints(self._points)
        self._template_poly_data.GetPointData().SetScalars(self._image_values)
        self._template_poly_data.SetPolys(self._polygons)

    def _generate_vtk_objects(self):
        """Generate vtk_points, vtk_polygons, vtk_poly_data, vtk_scalars,
        vtk_mapper and vtk_actor.
        This function must be called after _generate_points_and_polys()"""
        self._vtk_points = _vtk.vtkPoints()
        for coordinates in self._points:
            self._vtk_points.InsertNextPoint(coordinates[0], coordinates[1], coordinates[2])

        self._vtk_polygons = _vtk.vtkCellArray()
        for polygon in self._polygons:
            vtk_polygon = _vtk.vtkPolygon()
            vtk_polygon.GetPointIds().SetNumberOfIds(3)
            for local_index, global_index in enumerate(polygon):
                vtk_polygon.GetPointIds().SetId(local_index, global_index)
            self._vtk_polygons.InsertNextCell(vtk_polygon)

        self._vtk_poly_data = _vtk.vtkPolyData()
        self._vtk_poly_data.SetPoints(self._vtk_points)
        self._vtk_poly_data.SetPolys(self._vtk_polygons)

        self._vtk_scalars = _vtk.vtkFloatArray()
        self._vtk_scalars.SetNumberOfValues(self._vtk_poly_data.GetPoints().GetNumberOfPoints())
        for i in range(self._vtk_scalars.GetNumberOfTuples()):
            self._vtk_scalars.SetValue(i, 0.)

        self._vtk_poly_data.GetPointData().SetScalars(self._vtk_scalars)
        self._vtk_poly_data.Modified()

        self._vtk_mapper = _vtk.vtkPolyDataMapper()
        if VTK_VERSION < 6:
            self._vtk_mapper.SetInput(self._vtk_poly_data)
        else:
            self._vtk_mapper.SetInputData(self._vtk_poly_data)
        self._vtk_mapper.InterpolateScalarsBeforeMappingOn()
        self._vtk_mapper.UseLookupTableScalarRangeOn()
        self._vtk_actor = _vtk.vtkActor()
        self._vtk_actor.SetMapper(self._vtk_mapper)

        normals = _vtk.vtkFloatArray()
        normals.SetNumberOfComponents(3)
        for point in self._points:
            normals.InsertNextTuple(point)
        self._vtk_poly_data.GetPointData().SetNormals(normals)

    def _circular_slice(self):
        """Return a cellarray of a square slice."""
        polygons = vtk.vtkCellArray()
        for i in range(self._side-1):
            for j in range(self._side-1):
                corners = [(i, j), (i+1, j), (i+1, j+1), (i, j+1)]
                radius = max([numpy.sqrt((c[0] - self._side/2. + 0.5)**2 + (c[1] - self._side/2. + 0.5)**2)
                              for c in corners])
                if radius < self._side/2.:
                    polygon = vtk.vtkPolygon()
                    polygon.GetPointIds().SetNumberOfIds(4)
                    for index, corner in enumerate(corners):
                        polygon.GetPointIds().SetId(index, corner[0]*self._side+corner[1])
                    polygons.InsertNextCell(polygon)

        template_poly_data = vtk.vtkPolyData()
        template_poly_data.SetPoints(self._points)
        template_poly_data.GetPointData().SetScalars(self._image_values)
        template_poly_data.SetPolys(polygons)
        return template_poly_data

    def createActor1(self):
        # Create source
        polyDataPoints = vtk.vtkPoints()
        polyDataPoints.InsertPoint(0, 1., 0., 0.)
        polyDataPoints.InsertPoint(1, 1., 0., 10.)
        polyDataPoints.InsertPoint(2, 3., 0., 10.)
        polyDataPoints.InsertPoint(3, 2., 0., 0.)

        polygon = vtk.vtkPolygon()
        polygon.GetPointIds().SetNumberOfIds(4)
        polygon.GetPointIds().SetId(0, 0)
        polygon.GetPointIds().SetId(1, 1)
        polygon.GetPointIds().SetId(2, 2)
        polygon.GetPointIds().SetId(3, 3)

        polygons = vtk.vtkCellArray()
        polygons.InsertNextCell(polygon)

        polyData = vtk.vtkPolyData()
        polyData.SetPoints(polyDataPoints)
        #polyData.SetLines(polygons)
        polyData.SetPolys(polygons)

        rotationalExtrusionFilter = vtk.vtkRotationalExtrusionFilter()
        rotationalExtrusionFilter.SetInput(polyData)
        rotationalExtrusionFilter.SetResolution(10)
        rotationalExtrusionFilter.SetAngle(240)
        rotationalExtrusionFilter.SetTranslation(0)
 
        # Create a mapper
        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInputConnection(rotationalExtrusionFilter.GetOutputPort())
 
        # Create an actor
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
        #actor.GetProperty().SetColor(1.0, 0, 0)

        return actor

def makePolygon(polygon,elevation=0,triangulate=False):
    """
    Function to make a 2D vtk polygon PolyData from points making up a polygon.

    Inputs
        polygon - array list of vertices in clockwise order
        elevation - float, elevation of the polygon
        triangulate - boolean, triangulate the surface

    Output
        polyData - returns polydata of the surface polygon
    """

    # Make a enclosed polygons
    nrBasepts = len(polygon)
    ptsArr = np.hstack((polygon,np.ones((polygon.shape[0],1))*elevation))
    ppts = vtk.vtkPoints()
    ppts.SetData(npsup.numpy_to_vtk(ptsArr,deep=1))
    # Make the cells array of polygons
    polys = vtk.vtkCellArray()
    poly = vtk.vtkPolygon()
    poly.GetPointIds().SetNumberOfIds(len(ptsArr))
    for ind in np.arange(len(ptsArr)):
        poly.GetPointIds().SetId(ind,ind)
    polys.InsertNextCell(poly)

    polyPolyData = vtk.vtkPolyData()
    polyPolyData.SetPoints(ppts)
    polyPolyData.SetPolys(polys)
    if triangulate:
        triFilt = vtk.vtkTriangleFilter()
        triFilt.SetInputData(polyPolyData)
        triFilt.Update()
        return triFilt.GetOutput()
    else:
        return polyPolyData

def getpolyactor(poly, h=0):

    points = vtk.vtkPoints()
    for p in poly:
        print 'getpolyactor', p
        if len(p) == 2:
            p = [p[0], p[1], h]

        points.InsertNextPoint(p)

    # Create the polygon
    polygon = vtk.vtkPolygon()
    polygon.GetPointIds().SetNumberOfIds(len(poly)-1)
    for i in range(len(poly)-1):
        polygon.GetPointIds().SetId(i, i)
    # polygon.GetPointIds().SetId(0, 0)

    # Add the polygon to a list of polygons
    polygons = vtk.vtkCellArray()
    polygons.InsertNextCell(polygon)

    # Create a PolyData
    polygonPolyData = vtk.vtkPolyData()
    polygonPolyData.SetPoints(points)
    polygonPolyData.SetPolys(polygons)

    # Create a mapper and actor
    mapper = vtk.vtkPolyDataMapper()
    if vtk.VTK_MAJOR_VERSION <= 5:
        mapper.SetInput(polygonPolyData)
    else:
        mapper.SetInputData(polygonPolyData)

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    return actor

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

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

        self.Surface = triangleFilter.GetOutput()

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

        select = vtk.vtkImplicitSelectionLoop()
        select.SetLoop(self.Loop.GetPoints())
        normal = [0.0,0.0,0.0]
        centroid = [0.0,0.0,0.0]
        vtk.vtkPolygon().ComputeNormal(self.Loop.GetPoints(),normal)

        #compute centroid and check normals
        p = [0.0,0.0,0.0]
        for i in range (self.Loop.GetNumberOfPoints()):
            p = self.Loop.GetPoint(i)
            centroid[0] += p[0]
            centroid[1] += p[1]
            centroid[2] += p[2]
        centroid[0] = centroid[0] / self.Loop.GetNumberOfPoints()
        centroid[1] = centroid[1] / self.Loop.GetNumberOfPoints()
        centroid[2] = centroid[2] / self.Loop.GetNumberOfPoints()
        print "loop centroid", centroid

        locator = vtk.vtkPointLocator()
        locator.SetDataSet(self.Surface)
        locator.AutomaticOn()
        locator.BuildLocator()
        idsurface = locator.FindClosestPoint(centroid)

        if (self.Surface.GetPointData().GetNormals() == None):
            normalsFilter = vmtkscripts.vmtkSurfaceNormals()
            normalsFilter.Surface = self.Surface
            normalsFilter.NormalsArrayName = 'Normals'
            normalsFilter.Execute()
            self.Surface = normalsFilter.Surface
        normalsurface = [0.0,0.0,0.0]
        self.Surface.GetPointData().GetNormals().GetTuple(idsurface,normalsurface)
        print "loop normal: ", normal
        print "surface normal inside the loop: ", normalsurface
        check = vtk.vtkMath.Dot(normalsurface,normal)
        if check < 0:
            normal[0] = - normal[0]
            normal[1] = - normal[1]
            normal[2] = - normal[2]

        #compute plane
        proj = float(vtk.vtkMath.Dot(self.Loop.GetPoint(0),normal))
        point = [0.0,0.0,0.0]
        self.Loop.GetPoint(0,point)
        for i in range (self.Loop.GetNumberOfPoints()):
            tmp = vtk.vtkMath.Dot(self.Loop.GetPoint(i),normal)
            if tmp < proj:
                proj = tmp
                self.Loop.GetPoint(i,point)
        origin = [0.0,0.0,0.0]
        origin[0] = point[0] #- normal[0]
        origin[1] = point[1] #- normal[1]
        origin[2] = point[2] #- normal[2]
        plane=vtk.vtkPlane()
        plane.SetNormal(normal[0],normal[1],normal[2])
        plane.SetOrigin(origin[0],origin[1],origin[2])

        #set bool
        Bool = vtk.vtkImplicitBoolean()
        Bool.SetOperationTypeToDifference()
        Bool.AddFunction(select)
        Bool.AddFunction(plane)

        clipper=vtk.vtkClipPolyData()
        clipper.SetInput(self.Surface)
        clipper.SetClipFunction(Bool)
        clipper.GenerateClippedOutputOn()
        clipper.InsideOutOff()
        clipper.Update()

        self.Surface = clipper.GetOutput()

def decimate_file(input_vtk, reduction=0.5, smooth_steps=100, output_vtk=''):
    """
    Decimate vtk triangular mesh file with vtk.vtkDecimatePro.

    Parameters
    ----------
    input_vtk : string
        input vtk file with triangular surface mesh
    reduction : float
        fraction of mesh faces to remove
    do_smooth : Boolean
        smooth after decimation?
    output_vtk : string
        output decimated vtk file

    Returns
    -------
    output_vtk : string
        output decimated vtk file

    Examples
    --------
    >>> import os
    >>> from mindboggle.utils.mesh import decimate_file
    >>> path = os.environ['MINDBOGGLE_DATA']
    >>> input_vtk = os.path.join(path, 'arno', 'labels', 'label22.vtk')
    >>> output_vtk = 'decimated_file.vtk'
    >>> reduction = 0.9
    >>> smooth_steps = 0
    >>> decimate_file(input_vtk, reduction=reduction,
    >>>               smooth_steps=smooth_steps, output_vtk=output_vtk)
    >>> # View:
    >>> os.system('mayavi2 -d ' + output_vtk + ' -m Surface &')

    """
    import os
    import vtk

    from mindboggle.utils.io_vtk import read_vtk

    # Read VTK surface mesh file:
    faces, u1, u2, points, u4, labels, u5, u6 = read_vtk(input_vtk)

    # vtk points:
    vtk_points = vtk.vtkPoints()
    [vtk_points.InsertPoint(i, x[0], x[1], x[2]) for i,x in enumerate(points)]

    # vtk faces:
    vtk_faces = vtk.vtkCellArray()
    for face in faces:
        vtk_face = vtk.vtkPolygon()
        vtk_face.GetPointIds().SetNumberOfIds(3)
        vtk_face.GetPointIds().SetId(0, face[0])
        vtk_face.GetPointIds().SetId(1, face[1])
        vtk_face.GetPointIds().SetId(2, face[2])
        vtk_faces.InsertNextCell(vtk_face)

    # vtkPolyData:
    polydata = vtk.vtkPolyData()
    polydata.SetPoints(vtk_points)
    polydata.SetPolys(vtk_faces)

    # We want to preserve topology (not let any cracks form).
    # This may limit the total reduction possible.
    decimate = vtk.vtkDecimatePro()
    decimate.SetInput(polydata)
    decimate.SetTargetReduction(reduction)
    decimate.PreserveTopologyOn()

    # Export output:
    if not output_vtk:
        output_vtk = os.path.join(os.getcwd(), 'decimated_file.vtk')
    exporter = vtk.vtkPolyDataWriter()
    if smooth_steps > 0:
        smoother = vtk.vtkSmoothPolyDataFilter()
        smoother.SetInputConnection(decimate.GetOutputPort())
        smoother.SetNumberOfIterations(smooth_steps)
        exporter.SetInput(smoother.GetOutput())
    else:
        exporter.SetInput(decimate.GetOutput())
    exporter.SetFileName(output_vtk)
    exporter.Write()

    return output_vtk