def _read_polydata(self):
        """
        This private method reads the given `filename` and return a vtkPolyData
        object containing all informations about file; to avoid useless IO
        operation on the same file, it stores polydata of the last file parsed
        and if user ask for this file, the polydata previously stored is
        returned.

        :return: polydata containing information about file.
        :rtype: vtkPolyData
        """
        # Polydata from `filename` is allready loaded; return it
        if self._cached_data is not None:
            return self._cached_data

        if not os.path.isfile(self._filename):
            raise RuntimeError("{0!s} doesn't exist".format(
                os.path.abspath(self._filename)))

        reader = vtk.vtkDataSetReader()
        reader.SetFileName(self._filename)
        reader.Update()
        data = reader.GetOutput()

        self._cached_data = data

        return data

def get_tensor_array(filename):
    """Returns vtk tensor array object which can have 'GetTuple9(i)' called on it."""
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(filename)
    reader.Update()

    output = reader.GetOutput()
    npoints = output.GetNumberOfPoints()
    print 'npoints:', npoints
    ncells = output.GetNumberOfCells()
    print 'ncells:', ncells
    nscalars = reader.GetNumberOfScalarsInFile()
    print 'nscalars:', nscalars
    print reader.GetScalarsNameInFile(0)
    ntensors = reader.GetNumberOfTensorsInFile()
    print 'ntensors:', ntensors
    print reader.GetTensorsNameInFile(0)

    output = reader.GetOutput()
    point_data = output.GetPointData()
    scalar_array = point_data.GetArray('scalar')
    tensor_array = point_data.GetArray('tensor')
    if not tensor_array:
        tensor_array = point_data.GetArray('tensors')
    return tensor_array

def load_vtk(t, case="free"):
    f = "rvat-log-" + case + "/velocity/"+str(t)+".vtk"
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(f)
    reader.Update()
    data = reader.GetOutput()
    
    npoints = data.GetNumberOfPoints()
    point = data.GetPoint(0)
    d = data.GetPointData()
    
    array = d.GetArray("Velocity")
    
    u = np.zeros(npoints)
    v = np.zeros(npoints)
    w = np.zeros(npoints)
    x = np.zeros(npoints)
    y = np.zeros(npoints)
    z = np.zeros(npoints)
    
    for n in range(npoints):
        x[n], y[n], z[n] = data.GetPoint(n)
        u[n], v[n], w[n] = array.GetTuple(n)
        
    u = u[np.where(x==1.0)[0]]
    v = v[np.where(x==1.0)[0]]
    w = w[np.where(x==1.0)[0]]
    y = y[np.where(x==1.0)[0]]
    z = z[np.where(x==1.0)[0]]
        
    yarray, zarray, [uarray, varray, warray] = ts.build_plane_arrays(y, z,
                                                                    [u, v, w])
    return yarray, zarray, uarray, varray, warray

def readCellsFromFile(cells, path, iteration, rank):
  import vtk
  import os.path
  filename = path.replace('__ITERATION__', str(iteration)).replace('__RANK__', str(rank))
  
  if os.path.exists(filename):
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(filename)
    reader.SetReadAllScalars(True)
    reader.Update()

    grid = reader.GetOutput()
    
    numberOfCells = grid.GetNumberOfCells()
    cellData = grid.GetCellData()
    qs = cellData.GetScalars("q0")

    for cellId in xrange(numberOfCells):
      vtkCell = grid.GetCell(cellId)
      
      q = qs.GetTuple(cellId)[0]
      cells.append(Cell(vtkCell, vtkCell.GetBounds()[:], q))
      
    return numberOfCells
  else:
    return 0

    def write(self, mesh_points, filename):
        """
        Writes a vtk file, called filename, copying all the
        structures from self.filename but the coordinates.
        `mesh_points` is a matrix that contains the new coordinates
        to write in the vtk file.

        :param numpy.ndarray mesh_points: it is a `n_points`-by-3
            matrix containing the coordinates of the points of the
            mesh
        :param string filename: name of the output file.
        """
        self._check_filename_type(filename)
        self._check_extension(filename)
        self._check_infile_instantiation()

        self.outfile = filename

        reader = vtk.vtkDataSetReader()
        reader.SetFileName(self.infile)
        reader.ReadAllVectorsOn()
        reader.ReadAllScalarsOn()
        reader.Update()
        data = reader.GetOutput()

        points = vtk.vtkPoints()
        for i in range(data.GetNumberOfPoints()):
            points.InsertNextPoint(mesh_points[i, :])

        data.SetPoints(points)

        writer = vtk.vtkDataSetWriter()
        writer.SetFileName(self.outfile)
        writer.SetInputData(data)
        writer.Write()

def ConvertVTKMatlab(input_filename,output_filename):
  import vtk.util.numpy_support as vtkNumPy 
  import numpy
  import scipy.io as scipyio
 
  extension = input_filename.split('.').pop()
  vtkReader = None
  if extension == 'vtk':
    vtkReader = vtk.vtkDataSetReader() 
  elif extension == 'vti':
    vtkReader = vtk.vtkXMLImageDataReader() 
  else:
    raise RuntimeError('unknown file type %s ' % input_filename)
  vtkReader.SetFileName( "%s" % (input_filename) ) 
  vtkReader.Update()
  imageDataVTK = vtkReader.GetOutput()
  dimensions = imageDataVTK.GetDimensions()
  spacing = imageDataVTK.GetSpacing()
  origin  = imageDataVTK.GetOrigin()
  print spacing, origin, dimensions
  #fem.SetImagingDimensions( dimensions ,origin,spacing) 

  image_point_data = imageDataVTK.GetPointData() 
  image_data       = vtkNumPy.vtk_to_numpy( image_point_data.GetArray(0) ) 
  # write numpy to disk in matlab
  #  indexing is painful.... reshape to dimensions and transpose 2d dimensions only
  scipyio.savemat( output_filename, {'spacing':spacing, 'origin':origin,'image':image_data.reshape(dimensions,order='F').transpose(1,0,2)})

def show_vtk_file(path):
	# open vtk file from arglist
	reader = vtk.vtkDataSetReader()
	reader.SetFileName(path)
	reader.Update()
	# read out data and scalarrange
	output = reader.GetOutput()
	scalarrange = output.GetScalarRange()
	# generate Mapper and set DataSource
	mapper = vtk.vtkDataSetMapper()
	mapper.SetInput(output)
	mapper.SetScalarRange(scalarrange)
	# create actor
	actor = vtk.vtkActor()
	actor.SetMapper(mapper)
	# build renderer
	renderer = vtk.vtkRenderer()
	window = vtk.vtkRenderWindow()
	window.AddRenderer(renderer)
	# create interaction
	interaction = vtk.vtkRenderWindowInteractor()
	interaction.SetRenderWindow(window)
	# add actor
	renderer.AddActor(actor)
	# show window and start inteaction and renderer
	interaction.Initialize()
	window.Render()
	interaction.Start()

def get_tensor_array(filename):
    """Returns vtk tensor array object which can have 'GetTuple9(i)' called on it."""
    print "Reading " + filename
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(filename)
    reader.Update()

    output = reader.GetOutput()
    print 'npoints:', output.GetNumberOfPoints()
    print 'ncells:', output.GetNumberOfCells()
    print 'nscalars:', reader.GetNumberOfScalarsInFile()
    print 'ntensors:', reader.GetNumberOfTensorsInFile()
    print 'ScalarName:', reader.GetScalarsNameInFile(0)
    print 'TensorName:', reader.GetTensorsNameInFile(0)

    output = reader.GetOutput()
    pointdata = output.GetPointData()
    #scalar_array = pointdata.GetArray('scalar')
    tensor_array = pointdata.GetTensors()
    if not tensor_array:
        tensor_array = pointdata.GetArray('tensor')
    if not tensor_array:
        tensor_array = pointdata.GetArray('tensors')
    if not tensor_array:
        tensor_array = pointdata.GetArray('Tensors_')
    if not tensor_array:
        tensor_array = pointdata.GetArray('tensor1')
    if not tensor_array:
        print "Cannot find tensors in %s" % filename
        sys.exit(1)
    return tensor_array

def getReaderWriter(file_name, out_dir=None):
    r = vtk.vtkDataSetReader()
    r.SetFileName(file_name)
    f_base = os.path.splitext(file_name)[0]
    r.Update()
    reader = None
    writer = None
    xmlsuffix = '.xml'
    map = {'StructuredPoints': '.vti', 'StructuredGrid': '.vts',
           'RectilinearGrid': '.vtr', 'UnstructuredGrid': '.vtu',
           'PolyData': '.vtp'}
    for i in ['StructuredPoints', 'StructuredGrid', 'RectilinearGrid',
              'UnstructuredGrid', 'PolyData']:
        if eval('r.IsFile%s()'%i):
            reader = eval('vtk.vtk%sReader()'%i)
            if i == 'StructuredPoints':
                writer = eval('vtk.vtkXMLImageDataWriter()')
            else:
                writer = eval('vtk.vtkXML%sWriter()'%i)
            xmlsuffix = map[i]
            break
    if not reader:
        return None, None

    reader.SetFileName(file_name)
    reader.Update()

    out_file = f_base + xmlsuffix
    if out_dir:
        out_file = os.path.join(out_dir,
                                os.path.basename(f_base) + xmlsuffix)
    writer.SetFileName(out_file)
    return reader, writer

def read_points(filename):
    """
    Load points of a VTK surface file.

    Parameters
    ----------
    filename : string
        path/filename of a VTK format file

    Returns
    -------
    points : list of lists of floats
        each element is a list of 3-D coordinates of a surface mesh vertex

    """
    import vtk

    Reader = vtk.vtkDataSetReader()
    Reader.SetFileName(filename)
    Reader.ReadAllScalarsOn()  # Activate the reading of all scalars
    Reader.Update()

    Data = Reader.GetOutput()
    points = [list(Data.GetPoint(point_id))
              for point_id in range(Data.GetNumberOfPoints())]

    return points

def read_vertices(Filename):
    """
    Load VERTICES segment from a VTK file (actually contains indices to vertices)

    Parameters
    ----------
    Filename : string
        The path/filename of a VTK format file.

    Returns
    -------
    indices : a list of integers
        Each element is an integer defined in VERTICES segment of the VTK file.
        The integer is an index referring to a point defined in POINTS segment of the VTK file.

    Notes ::

        We assume that VERTICES segment is organized as one line,
        the first column of which is the number of vertices.
        Vertices here are as vertices in VTK terminology. It may not be the vertices in your 3-D surface.

    """
    import vtk

    Reader = vtk.vtkDataSetReader()
    Reader.SetFileName(Filename)
    Reader.Update()

    Data = Reader.GetOutput()

    Vrts = Data.GetVerts()
    indices = [Vrts.GetData().GetValue(i) for i in range(1, Vrts.GetSize())]

    return indices

	def parse(self, filename):
		"""
		Method to parse the file `filename`. It returns a matrix with all the coordinates.

		:param string filename: name of the input file.

		:return: mesh_points: it is a `n_points`-by-3 matrix containing the coordinates of
			the points of the mesh
		:rtype: numpy.ndarray

		.. todo::

			- specify when it works
		"""
		self._check_filename_type(filename)
		self._check_extension(filename)

		self.infile = filename

		reader = vtk.vtkDataSetReader()
		reader.SetFileName(self.infile)
		reader.ReadAllVectorsOn()
		reader.ReadAllScalarsOn()
		reader.Update()
		data = reader.GetOutput()

		n_points = data.GetNumberOfPoints()
		mesh_points = np.zeros([n_points, 3])

		for i in range(n_points):
			mesh_points[i][0], mesh_points[i][1], mesh_points[i][2] = data.GetPoint(i)

		return mesh_points

    def _ReadFile(self,**kw):
        """Read in a VTK data structure from a vtk data file
        """


        fd = open( self.filepath, 'r' )
        
        # Instantiate the field object
        field = objects.field.Field()
        
        reader = vtk.vtkDataSetReader()
        reader.SetFileName(self.filepath)
        #field.vtkdata = reader.GetStructuredPointsOutput()

        data = reader.GetOutput()
        # Need to update so that the data gets pulled through.
        data.Update()
        if not data:
            raise Exception,"VTK_IO no data found while reading file: %s" % self.filepath

        if not data.GetClassName() == 'vtkStructuredPoints':
            raise Exception,"VTK_IO someone needs to tell me how to read something other than Structured Points!"


        field.dim = data.GetDimensions()
        #print data.GetDataDimension()
        origin =  data.GetOrigin()
        field.origin = objects.vector.Vector( origin )
        
        field.vtkdata = data
        field.title = self.name
        field.name = self.name
        self.fields.append(field)

def Test3(datadir):
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(datadir + "/Data/blow.vtk")
    reader.UpdateInformation();
    reader.ReadAllScalarsOn()
    reader.ReadAllVectorsOn()

    dssf = vtk.vtkDataSetSurfaceFilter()
    dssf.SetInputConnection(reader.GetOutputPort())

    stripper = vtk.vtkStripper()
    stripper.SetInputConnection(dssf.GetOutputPort())

    f = vtk.vtkIntegrateAttributes()
    f.SetInputConnection(stripper.GetOutputPort())
    f.Update()

    result = f.GetOutputDataObject(0)
    val = result.GetPointData().GetArray("displacement1").GetValue(0)
    assert (val > 463.64 and val < 463.642)

    val = result.GetPointData().GetArray("thickness3").GetValue(0)
    assert (val > 874.61 and val < 874.618)

    val = result.GetCellData().GetArray("Area").GetValue(0)
    assert (val > 1145.405 and val < 1145.415)

	def plot(self, plot_file=None, save_fig=False):
		"""
		Method to plot a vtk file. If `plot_file` is not given it plots `self.infile`.

		:param string plot_file: the vtk filename you want to plot.
		:param bool save_fig: a flag to save the figure in png or not. If True the
			plot is not shown.
			
		:return: figure: matlplotlib structure for the figure of the chosen geometry
		:rtype: matplotlib.pyplot.figure
		"""
		if plot_file is None:
			plot_file = self.infile
		else:
			self._check_filename_type(plot_file)

		# Read the source file.		
		reader = vtk.vtkDataSetReader()
		reader.SetFileName(plot_file)
		reader.Update()

		data = reader.GetOutput()
		points = data.GetPoints()
		ncells = data.GetNumberOfCells()

		# for each cell it contains the indeces of the points that define the cell
		figure = plt.figure()
		axes = a3.Axes3D(figure)
		vtx = np.zeros((ncells, 3, 3))
		for i in range(0, ncells):
			for j in range(0, 3):
				cell = data.GetCell(i).GetPointId(j)
				vtx[i][j][0], vtx[i][j][1], vtx[i][j][2] = points.GetPoint(int(cell))
			tri = a3.art3d.Poly3DCollection([vtx[i]])
			tri.set_color('b')
			tri.set_edgecolor('k')
			axes.add_collection3d(tri)
		
		## Get the limits of the axis and center the geometry
		max_dim = np.array([np.max(vtx[:,:,0]), \
						np.max(vtx[:,:,1]), \
						np.max(vtx[:,:,2])])
		min_dim = np.array([np.min(vtx[:,:,0]), \
						np.min(vtx[:,:,1]), \
						np.min(vtx[:,:,2])])
		
		max_lenght = np.max(max_dim - min_dim)
		axes.set_xlim(-.6*max_lenght + (max_dim[0]+min_dim[0])/2, .6*max_lenght + (max_dim[0]+min_dim[0])/2)
		axes.set_ylim(-.6*max_lenght + (max_dim[1]+min_dim[1])/2, .6*max_lenght + (max_dim[1]+min_dim[1])/2)
		axes.set_zlim(-.6*max_lenght + (max_dim[2]+min_dim[2])/2, .6*max_lenght + (max_dim[2]+min_dim[2])/2)

		# Show the plot to the screen
		if not save_fig:
			plt.show()
		else:
			figure.savefig(plot_file.split('.')[0] + '.png')
			
		return figure

    def testGlyphs(self):
        '''Test if the glyphs are created nicely.'''
        reader = vtk.vtkDataSetReader()

        data_file = os.path.join(Testing.VTK_DATA_ROOT, "Data", "tensors.vtk")

        reader.SetFileName(data_file)

        g1 = SimpleGlyph(reader)
        g1.glyph.ColorGlyphsOff()
        g1.Update()

        g2 = SimpleGlyph(reader)
        g2.glyph.ExtractEigenvaluesOff()
        g2.Update()
        g2.SetPosition((2.0, 0.0, 0.0))

        g3 = SimpleGlyph(reader)
        g3.glyph.SetColorModeToEigenvalues()
        g3.glyph.ThreeGlyphsOn()
        g3.Update()
        g3.SetPosition((0.0, 2.0, 0.0))

        g4 = SimpleGlyph(reader)
        g4.glyph.SetColorModeToEigenvalues()
        g4.glyph.ThreeGlyphsOn()
        g4.glyph.SymmetricOn()
        g4.Update()
        g4.SetPosition((2.0, 2.0, 0.0))

        # 6Components symetric tensor
        g5 = SimpleGlyph(reader)
        g5.glyph.SetInputArrayToProcess(0, 0, 0, 0, "symTensors1")
        g5.SetPosition((4.0, 2.0, 0.0))
        g5.Update()

        ren = vtk.vtkRenderer()
        for i in (g1, g2, g3, g4, g5):
            for j in i.GetActors():
                ren.AddActor(j)

        ren.ResetCamera();

        cam = ren.GetActiveCamera()
        cam.Azimuth(-20)
        cam.Elevation(20)
        cam.Zoom(1.1)

        ren.SetBackground(0.5, 0.5, 0.5)

        renWin = vtk.vtkRenderWindow()
        renWin.AddRenderer(ren)
        renWin.Render()

        img_file = "TestTensorGlyph.png"
        Testing.compareImage(renWin, Testing.getAbsImagePath(img_file))
        Testing.interact()

def read(filename):
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(filename)
    reader.ReadAllScalarsOn()
    reader.Update()
    data = reader.GetOutput()
    header = reader.GetHeader()
    no_nodes = data.GetNumberOfPoints()
    no_connections = data.GetNumberOfCells()
    return (data, header, no_nodes, no_connections)

    def _create_vtk_reader(self,url,type):
        if type==1:
            self.reader=vtk.vtkDataSetReader()
            self.reader.SetFileName(url)
            self.reader.Update()
            print "zrobilem vtk readera"
        if type==0:
			print "convert_dat_file_to_vtk_reader(url)"
			
			self.reader=self.convert_dat_file_to_vtk_reader(url)

def get_num_fibers(filename):
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(filename)
    reader.Update()

    output = reader.GetOutput()
    #import pudb
    #pudb.set_trace()
    nlines = reader.GetOutput().GetNumberOfLines()
    return nlines

def load_vtk_cloud(vtk_cloud_file):
    reader = vtk.vtkDataSetReader()
    reader.SetFileName(vtk_cloud_file)
    reader.Update()

    actor = vtk.vtkActor()
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(reader.GetOutputPort())
    actor.SetMapper(mapper)

    return actor