def ConvertVTItoMatlab(work_dir):
  import vtk
  import vtk.util.numpy_support as vtkNumPy 
  import numpy
  import scipy.io as scipyio
  import os
  # echo vtk version info
  print "using vtk version", vtk.vtkVersion.GetVTKVersion()

  directoryList = os.listdir(work_dir)
  files = filter(lambda x:os.path.isfile("%s/%s" % (work_dir,x) ),directoryList)
  vtiFiles = filter(lambda x: x.split(".").pop() == "vti" ,files)
  filenames = filter(lambda x: x.split(".").pop(0) == "updatefemROI5" ,vtiFiles)
  
  for idlist, fileID in enumerate( filenames):
    print "reading %s/%s" % (work_dir,fileID) 
    vtkReader = vtk.vtkXMLImageDataReader() 
    vtkReader.SetFileName( "%s/%s" % (work_dir,fileID) ) 
    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
    scipyio.savemat("%s/%s.mat" % (work_dir,fileID), {'spacing':spacing, 'origin':origin,'image':image_data})
    print "wrote %d of %d" % (idlist, len( filenames))

def load(file):
    datareader = v.vtkXMLImageDataReader()
    datareader.SetFileName(file)
    datareader.Update()
    data = datareader.GetOutput()

    return data

    def _load_slice(self, q, index, desc):

        if not self._volume:
            import vtk
            dirname = os.path.dirname(self.__dbfilename)
            vol_file = os.path.join(dirname, "cinema.vti")
            vr = vtk.vtkXMLImageDataReader()
            vr.SetFileName(vol_file)
            vr.Update()
            volume = vr.GetOutput()
            self._volume = volume
            self._vol_file = vol_file
        else:
            volume = self._volume

        ext = volume.GetExtent()
        width = ext[1]-ext[0]
        height = ext[3]-ext[2]

        from vtk.numpy_interface import dataset_adapter as dsa
        image = dsa.WrapDataObject(volume)

        oid = volume.ComputePointId([0, 0, index])
        nparray = image.PointData[0]
        imageslice = np.reshape(nparray[oid:oid+width*height], (width,height,3))

        doc = Document(desc, imageslice)
        doc.attributes = None
        return doc

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkXMLImageDataReader(), 'Reading vtkXMLImageData.',
         (), ('vtkXMLImageData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

def main(argv):
   #Just get something working for testing...
    try:
        opts, args = getopt.getopt(argv,"hi:", ["ifile="])
    except getopt.GetoptError as err:
        print 'tviewer.py -i <inputfile.vtk>'
        print (str(err))
    for opt, arg in opts:
        if opt == '-h':
            print 'tviewer.py -i <inputfile.vtk>'
            sys.exit()
        elif opt in ("-i", "--ifile"):
            inputfile = arg
    print("Going to load and generate png from ", inputfile)
    #Read data
    reader = vtk.vtkXMLImageDataReader()
    reader.SetFileName(inputfile)
    reader.Update()

    image = reader.GetOutput()
    #image.SetSpacing(1,1,1)
    #image.GetPointData().SetScalars(image.GetPointData().GetVectors())
    #Compute Q Criterion for texture mapping
    
    #Now put this in a png file
    castFilter = vtk.vtkImageCast()
    castFilter.SetInputData(image)
    castFilter.Update()
    
    w = vtk.vtkPNGWriter()
    w.SetInputData(castFilter.GetOutput())
    w.SetFileName("xyslice.png")
    w.Write()

    def __init__(self,offset,get=1,file='data',type='vtu',mirrorPlane=None,silent=0):
        self.offset=offset
        self.filenameout = file
        self.type = type
        self.isLoaded = False
        self.mirrorPlane=mirrorPlane
        self.silent = silent

        if type == 'vtu':
            self.filename=''.join([self.filenameout,repr(offset).zfill(4),'.vtu'])
            self.datareader = v.vtkXMLUnstructuredGridReader()
        elif type == 'pvtu':
            self.filename=''.join([self.filenameout,repr(offset).zfill(4),'.pvtu'])
            self.datareader = v.vtkXMLPUnstructuredGridReader()
        elif type == 'vti':
            self.filename=''.join([self.filenameout,repr(offset).zfill(4),'.vti'])
            self.datareader = v.vtkXMLImageDataReader()
        else:
            print 'Unknown filetype'
            
        if (self.silent == 0): print '========================================'
        if (self.silent == 0): print  'loading file %s' % (self.filename)

        if get != None:
            self.getAll()

def vti2mat(fileIn, fileOut):
    """Convert voxel-array from VTI to MAT (MATLAB(R)) format.
    
    Parameters
    ----------
    fileIn : str
        Path for input VTI file.
        
    fileOut : str
        Path for output MAT file.
        
    """
    
    import numpy as np
    import vtk
    import scipy.io as sio
    from vtk.util import numpy_support as nps
    from math_utils import lcmm
    
    reader = vtk.vtkXMLImageDataReader()
    reader.SetFileName(fileIn)
    reader.Update()
    vtkImageData = reader.GetOutput()
    dim = vtkImageData.GetDimensions()
    flatV = nps.vtk_to_numpy(vtkImageData.GetPointData().GetScalars())
    V = flatV.reshape(dim[::-1])
    spacing = np.array(vtkImageData.GetSpacing())[::-1]
    estimatedFactors = lcmm(*spacing) / spacing
    estimatedVoxelSize = 1. / estimatedFactors
    sio.savemat(fileOut, {'volume':V, 'spacing': spacing, 'estimated_voxel_size': estimatedVoxelSize})

 def getReader(self):
     reader = None;
     if self.srcFormat == "vti":
         reader = vtk.vtkXMLImageDataReader();
     reader.SetFileName(self.src);
     reader.Update();
     return reader;

 def loadVisualisations(self):
   for tuple in self.d_tuples:
     reader = vtk.vtkXMLImageDataReader()
     reader.SetFileName(tuple[0])
     reader.Update()
     self.d_window.open_vtk_data(reader.GetOutput())
     self.d_window.load_mm(tuple[1])

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 vti_to_numpy(filename, fieldName='density'):
    
    reader = vtk.vtkXMLImageDataReader()
    reader.SetFileName(filename)
    reader.Update()

    # get vtkImageData
    imageData = reader.GetOutput()

    sr = imageData.GetScalarRange()
    print("scalar range {}".format(sr))

    # get dimensions tuple
    dims = imageData.GetDimensions()
    print("dims {}".format(dims))
    
    # get vtk data
    vtk_data = imageData.GetPointData().GetArray(fieldName)

    # convert to numpy array
    numpy_data = numpy_support.vtk_to_numpy(vtk_data)

    #numpy_data = numpy_data.reshape(dims[0], dims[1], dims[2])
    numpy_data = numpy_data.reshape(dims[1], dims[0])
    print("shape in reader {}".format(numpy_data.shape))

    #numpy_data = numpy_data.transpose(2,1,0)
    
    return numpy_data

 def ReadVTKXMLImageFile(self):
     if (self.InputFileName == ''):
         self.PrintError('Error: no InputFileName.')
     self.PrintLog('Reading VTK XML image file.')
     reader = vtk.vtkXMLImageDataReader()
     reader.SetFileName(self.InputFileName)
     reader.Update()
     self.Image = reader.GetOutput()

def Import(filename):
    reader = vtk.vtkXMLImageDataReader()
    reader.SetFileName(filename)
    # TODO: Check if the code bellow is necessary
    reader.WholeSlicesOn()
    reader.Update()

    return reader.GetOutput()

 def __init__(self, vtifname, parallel=False):
     self.reader = vtk.vtkXMLImageDataReader()
     self.reader.SetFileName(vtifname)
     self.reader.Update()
     self.data = self.reader.GetOutput()  
     self.dim = self.data.GetDimensions()   
     self.s_scal = [self.dim[1]-1, self.dim[0]-1]
     self.s_vec = [self.dim[1]-1, self.dim[0]-1,3]

	def GetImageDataForBaseAndExtension(self, fileName, extension):
		"""
		:type fileName: basestring
		:type extension: basestring
		:rtype: vtkImageData
		"""
		if extension == DataReader.TypeMHA or extension == DataReader.TypeMHD:
			# Use a vktMetaImageReader
			imageReader = vtkMetaImageReader()
			imageReader.SetFileName(fileName)
			imageReader.Update()
			return imageReader.GetOutput()
		elif extension == DataReader.TypeDICOM:
			# Use a dicom reader
			dirName = os.path.dirname(fileName)
			return self.GetImageDataFromDirectory(dirName)
		elif extension == DataReader.TypeDAT:
			raise Exception("Support for .dat files is not implemented.")
			# Read in the .dat file byte by byte
			imageData = None

			import numpy as np
			with open(fileName, "rb") as f:
				dimensions = np.fromfile(f, np.int16, count=3)
				imageData = vtkImageData()
				imageData.SetDimensions(int(dimensions[0]), int(dimensions[1]), int(dimensions[2]))
				imageData.SetScalarTypeToFloat()
				imageData.SetNumberOfScalarComponents(1)
				imageData.AllocateScalars()
				imageData.Update()
				imageData.PrepareForNewData()

				fileData = np.fromfile(f, np.int16)

				dataIndex = 0
				for z in range(int(dimensions[2])):
					for y in range(int(dimensions[1])):
						for x in range(int(dimensions[0])):
							imageData.SetScalarComponentFromFloat(x, y, z, 0, float(fileData[dataIndex]))
							dataIndex += 1

			return imageData
		elif extension == DataReader.TypeVTI:
			# Use a XMLImageReader
			imageReader = vtkXMLImageDataReader()
			imageReader.SetFileName(fileName)
			imageReader.Update()
			return imageReader.GetOutput()
		elif extension == DataReader.TypeNRRD:
			# Use a NrrdReader
			imageReader = vtkNrrdReader()
			imageReader.SetFileName(fileName)
			imageReader.Update()
			return imageReader.GetOutput()
		else:
			assert False

def load_image_data(handle):
    """
    Load vtkImageData from XML file.

    :param handle: file name
    :type handle: str
    """
    reader = vtk.vtkXMLImageDataReader()
    reader.SetFileName(handle)
    return reader.GetOutput()

def multi_con_h2(filename):
    # Creating contour for h2 fileld
    
    data = vtk.vtkXMLImageDataReader()
    data.SetFileName(filename)#"/Users/y1275963/Documents/homework/multifield.0060.vti")
    data.Update()
    #Getting highest value
    [low,high] = data.GetOutput().GetPointData().GetArray("H2").GetRange()

    data_h2 = vtk.vtkAssignAttribute()
    data_h2.SetInputConnection(data.GetOutputPort())
    data_h2.Assign('H2','SCALARS','POINT_DATA')
    
    
    
    # Create a filter to extract an isosurface from
    # the dataset.  Connect the input of this
    # filter to the output of the reader.  Set the
    # value for the (one) isosurface that we want
    # to extract, and tell the filter to compute
    # normal vectors for each triangle on the
    # output isosurface.
    
    iso = vtk.vtkContourFilter()
    iso.SetInputConnection(data_h2.GetOutputPort())
    iso.ComputeNormalsOn()
    iso.SetNumberOfContours(10)
    for i in range(10):
        iso.SetValue(i, (1.0-0.1*i)*high) #Chnge here to change the isovalue
    
    
    # Create a mapper to traverse the polydata and
    # generate graphics commands for drawing the
    # polygons onto the output device.
    
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(iso.GetOutputPort())
    mapper.ScalarVisibilityOff()
    
    # Output primitives (i.e. the triangles making
    # up the isosurface) are managed as an actor;
    # we specify that all outputs are coloured
    # red.
    
    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(0, 0, 1)
    actor.GetProperty().SetOpacity(0.3)
    
    # Create a renderer which will output the
    # graphics commands onto a drawing surface within
    # a window.  By default the renderer will fill
    # all of the window.  We set the background
    # colour of the window to white.
    return actor

def vtkKWHeaderAnnotationEditorEntryPoint(parent, win):

    app = parent.GetApplication()
    
    # -----------------------------------------------------------------------
    
    # Create a render widget
    # Set the header annotation visibility and set some text
    
    hae_renderwidget = vtkKWRenderWidget()
    hae_renderwidget.SetParent(parent)
    hae_renderwidget.Create()
    
    hae_renderwidget.HeaderAnnotationVisibilityOn()
    hae_renderwidget.SetHeaderAnnotationText("Hello, World!")
    
    app.Script("pack %s -side right -fill both -expand y -padx 0 -pady 0",
        hae_renderwidget.GetWidgetName())
    
    # -----------------------------------------------------------------------
    
    # Create a volume reader
    
    hae_reader = vtkXMLImageDataReader()
    hae_reader.SetFileName(os.path.join(
        os.path.dirname(os.path.abspath(__file__)),
        "..", "..", "..", "..", "Data", "head100x100x47.vti"))
    
    # Create an image viewer
    # Use the render window and renderer of the renderwidget
    
    hae_viewer = vtkImageViewer2()
    hae_viewer.SetRenderWindow(hae_renderwidget.GetRenderWindow())
    hae_viewer.SetRenderer(hae_renderwidget.GetRenderer())
    hae_viewer.SetInput(hae_reader.GetOutput())
    hae_viewer.SetupInteractor(
        hae_renderwidget.GetRenderWindow().GetInteractor())
    
    hae_renderwidget.ResetCamera()
    
    # -----------------------------------------------------------------------
    
    # Create a header annotation editor
    # Connect it to the render widget
    
    hae_anno_editor = vtkKWHeaderAnnotationEditor()
    hae_anno_editor.SetParent(parent)
    hae_anno_editor.Create()
    hae_anno_editor.SetRenderWidget(hae_renderwidget)
    
    app.Script("pack %s -side left -anchor nw -expand n -padx 2 -pady 2",
        hae_anno_editor.GetWidgetName())
    
    return "TypeVTK"

def GetNumPyData(filename):
  print filename
  vtiReader = vtk.vtkXMLImageDataReader()
  vtiReader.SetFileName(filename)
  vtiReader.Update()
  vtiData = vtk.vtkImageCast()
  vtiData.SetOutputScalarTypeToDouble()
  vtiData.SetInput( vtiReader.GetOutput() )
  vtiData.Update( )
  vti_image = vtiData.GetOutput().GetPointData() 
  vti_array = vtkNumPy.vtk_to_numpy(vti_image.GetArray(0)) 
  return vti_array 

def vti_to_nparray(fname):
    r = vtk.vtkXMLImageDataReader()
    r.SetFileName(fname)
    r.Update()

    o = r.GetOutput()
    x, y, z = o.GetDimensions()
    
    m = numpy_support.vtk_to_numpy(o.GetPointData().GetScalars())
    m.shape = (z, y, x)
    m[:] = (m > 127) * 1
    
    return m, o.GetSpacing()