def main(argv):
  if len(argv) < 2:
    print "usage:",argv[0]," data.vtk"
    exit(1)
  data_fn = argv[1]
  reader = vtk.vtkStructuredPointsReader()
  reader.SetFileName(data_fn)
  reader.Update()
  data = reader.GetOutput()
  updateColorOpacity()
  # composite function (using ray tracing)
  compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
  volumeMapper = vtk.vtkVolumeRayCastMapper()
  volumeMapper.SetVolumeRayCastFunction(compositeFunction)
  volumeMapper.SetInput(data)
  # make the volume
  #volume = vtk.vtkVolume()
  global volume
  volume.SetMapper(volumeMapper)
  volume.SetProperty(volumeProperty)
  # renderer
  renderer = vtk.vtkRenderer()
  renderWin = vtk.vtkRenderWindow()
  renderWin.AddRenderer(renderer)
  renderInteractor = vtk.vtkRenderWindowInteractor()
  renderInteractor.SetRenderWindow(renderWin)
  renderInteractor.AddObserver( vtk.vtkCommand.KeyPressEvent, keyPressed )
  renderer.AddVolume(volume)
  renderer.SetBackground(0,0,0)
  renderWin.SetSize(400, 400)
  renderInteractor.Initialize()
  renderWin.Render()
  renderInteractor.Start()

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

    def __init__(self,delta,wing):

        self.arrowColor = vtk.vtkColorTransferFunction()

        self.update_look_up_table()

        self.print_counter=0
        self.ren = vtk.vtkRenderer()

        self.geo_reader = vtk.vtkUnstructuredGridReader()
        self.geo_reader.SetFileName(wing)


        self.vec_reader = vtk.vtkStructuredPointsReader()
        self.vec_reader.SetFileName(delta)

        """ This is for drawing the wing,"""
        geo_Mapper=vtk.vtkDataSetMapper()
        geo_Mapper.SetInputConnection(self.geo_reader.GetOutputPort())

        geo_actor = vtk.vtkActor()
        geo_actor.SetMapper(geo_Mapper)

        self.ren.AddActor(geo_actor)

        """End of adding the wing """


        self.ren.AddActor(self.create_stream_line(25,150,0,0.5))


        #Add renderer to renderwindow and render
        self.renWin = vtk.vtkRenderWindow()
        self.renWin.AddRenderer(self.ren)
        self.renWin.SetSize(1920, 1080)

        iren = vtk.vtkRenderWindowInteractor()
        iren.SetRenderWindow(self.renWin)

        iren.AddObserver('RightButtonPressEvent', self.capture_image, 1.0)

        # Scalar Bar actor
        scalar_bar = vtk.vtkScalarBarActor()
        scalar_bar.SetOrientationToHorizontal()
        scalar_bar.SetLookupTable(self.arrowColor)
        scalar_bar.SetTitle("Color map")
        scalar_bar.SetLabelFormat("%5.2f")
        scalar_bar.SetMaximumHeightInPixels(300)
        scalar_bar.SetMaximumWidthInPixels(60)

        # Scalar Bar Widget
        scalar_bar_widget = vtk.vtkScalarBarWidget()
        scalar_bar_widget.SetInteractor(iren)
        scalar_bar_widget.SetScalarBarActor(scalar_bar)
        scalar_bar_widget.On()

        self.ren.SetBackground(0,0,0)
        self.renWin.Render()
        iren.Start()

def load_data_file(filename):
    reader = vtk.vtkStructuredPointsReader()
    reader.SetFileName(filename)
    reader.Update()
    Image = reader.GetOutput()
    array = VTKtoNumpy(Image)
    array = array.astype(float)
    return array

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

    def readSimulationData(self,_i):

        # self.counter+=1
        # if self.counter>3:
            # return
        # print "self.drawMutex=",self.drawMutex
        self.drawMutex.lock()
        # print "self.readFileSem.available()=",self.readFileSem.available()
        self.readFileSem.acquire()
        
        self.newFileBeingLoaded=True # this flag is used to prevent calling  draw function when new data is read from hard drive
        # print "LOCKED self.newFileBeingLoaded=",self.newFileBeingLoaded
        # self.drawMutex.lock()
        if _i >= len(self.ldsFileList):
            return
            
        fileName = self.ldsFileList[_i]
        
        self.simulationDataReader = vtk.vtkStructuredPointsReader()
        self.currentFileName = os.path.join(self.ldsDir,fileName)
        self.simulationDataReader.SetFileName(self.currentFileName)
        # print "path= ", os.path.join(self.ldsDir,fileName)
        self.simulationDataReader.Update()
        self.simulationData = self.simulationDataReader.GetOutput()
        # print "self.simulationData",self.simulationData
        # print 'dimensions=',self.simulationData.GetDimensions()
        
        #updating fieldDim each time we read data
        
        import CompuCell
        self.fieldDimPrevious=self.fieldDim
        # self.fieldDimPrevious=CompuCell.Dim3D()
        
        dimFromVTK=self.simulationData.GetDimensions()
        
        self.fieldDim=CompuCell.Dim3D(dimFromVTK[0],dimFromVTK[1],dimFromVTK[2])
        
        
        # 
        # self.fieldDim.x=dimFromVTK[0]
        # self.fieldDim.y=dimFromVTK[1]
        # self.fieldDim.z=dimFromVTK[2]
        
        
        
        # print 'self.fieldDim=',self.fieldDim
        
        
        # print "vtkStructuredPointsReader ERROR CODE=",self.simulationDataReader.GetErrorCode()
        
        
        self.currentStep = self.frequency * _i
        self.setCurrentStep(self.currentStep)
        # print "self.frequency=",self.frequency," _i=",self.frequency
        
        # print "\n\n\n\n FINISHED RUNNING readSimulationData step ",self.currentStep," \n\n\n"        
        self.drawMutex.unlock()
        self.readFileSem.release()

def readVTKSP(fn):
	'''
	reads a vtk structured points object from a file
	'''
	sp_reader = vtk.vtkStructuredPointsReader()
	sp_reader.SetFileName(fn)
	sp_reader.Update()
	sp = sp_reader.GetOutput()
	return sp

	def createReader(self,fileName):
		#I dont care about different types of readers for now
		#If i get a different file then i can mod reader here or something
		reader = vtk.vtkStructuredPointsReader()
		reader.SetFileName(fileName)
		reader.ReadAllVectorsOn()
		reader.ReadAllScalarsOn()
		reader.Update()
		self.determineCenterForCuttingPlane(reader)
		return reader 

def VTKSPtoNumpyFromFile(fn):
	'''
	reads a .vts file into a numpy array


	args:
		@a fn - string, filename of .sp file to read
	'''
	reader = vtk.vtkStructuredPointsReader()
	reader.SetFileName(fn)
	reader.Update()
	sp = reader.GetOutput()
	return VTKSPtoNumpy(sp)

def calculate_vtk_max_pointwise_difference(file1,file2,tol=1e-6):
    arrays = [0]*2

    reader = vtk.vtkStructuredPointsReader()
    for i, fname in enumerate([file1, file2]):
        reader.SetFileName(fname)
        reader.Update()
        data = reader.GetOutput().GetPointData()
        arrays[i] = np.array([vtk_to_numpy(data.GetArray(n)) for n in range(data.GetNumberOfArrays())])

    try:
        return np.allclose(arrays[0],arrays[1],rtol=0,atol=tol), np.max(np.abs(arrays[0]-arrays[1]))
    except:
        return False, np.inf

def readVTKVolume(name):
    try:
        reader = vtk.vtkStructuredPointsReader()
        reader.SetFileName( name )
        reader.Update()
        print "Input volume:", name
        vol = reader.GetOutput()
        reader = None
        return vol
    except:
        print "VTK volume reader 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 parse_vtk(vtk_file, sim_dict, data_fields):

    dim = sim_dict['dim']
    reader = vtk.vtkStructuredPointsReader()
    # reader = vtkUnstructuredGridReader()
    reader.SetFileName(vtk_file)
    reader.Update()
    field_data = reader.GetOutput()
    out_data = {}
    for data_field in data_fields:
        out_data[data_field] = vtk_to_numpy(
            field_data.GetPointData().GetArray(data_field)).reshape((dim['x'],
                                                                     dim['y']))
    reader.CloseVTKFile()
    return out_data

def main(argv):
  if len(argv) < 2:
    print "usage:",argv[0]," data.vtk"
    exit(1)
  data_fn = argv[1]
  reader = vtk.vtkStructuredPointsReader()
  reader.SetFileName(data_fn)
  reader.Update()
  data = reader.GetOutput()
  skin = vtk.vtkMarchingCubes()
  skin.ComputeNormalsOn()
  skin.ComputeGradientsOn()
  skin.SetValue(0, rgb[0][0])
  skin.SetInput(data)
  skin_mapper = vtk.vtkPolyDataMapper()
  skin_mapper.SetInputConnection(skin.GetOutputPort())
  skin_mapper.ScalarVisibilityOff()
  skin_actor = vtk.vtkActor()
  skin_property = vtk.vtkProperty()
  skin_property.SetColor(rgb[0][1], rgb[0][2], rgb[0][3])
  skin_property.SetOpacity(opacity[0][1])
  skin_actor.SetProperty(skin_property)
  skin_actor.SetMapper(skin_mapper)
  bone = vtk.vtkMarchingCubes()
  bone.ComputeNormalsOn()
  bone.ComputeGradientsOn()
  bone.SetValue(0, rgb[1][0])
  bone.SetInput(data)
  bone_mapper = vtk.vtkPolyDataMapper()
  bone_mapper.SetInputConnection(bone.GetOutputPort())
  bone_mapper.ScalarVisibilityOff()
  bone_actor = vtk.vtkActor()
  bone_actor.GetProperty().SetColor(rgb[1][1], rgb[1][2], rgb[1][3])
  bone_actor.GetProperty().SetOpacity(opacity[1][1])
  bone_actor.SetMapper(bone_mapper)
  renderer = vtk.vtkRenderer()
  renderWin = vtk.vtkRenderWindow()
  renderWin.AddRenderer(renderer)
  renderInteractor = vtk.vtkRenderWindowInteractor()
  renderInteractor.SetRenderWindow(renderWin)
  renderer.AddActor(skin_actor)
  #renderer.AddActor(bone_actor)
  renderer.SetBackground(0,0,0)
  renderWin.SetSize(400, 400)
  renderInteractor.Initialize()
  renderWin.Render()
  renderInteractor.Start()

    def create_reader (self): 
        "Create the corresponding reader."
        debug ("In VtkDataReader::create_reader ()")
        # set up the reader     
        if self.file_name == "":
            raise IOError, "No filename specifed for the data handler!"

        self.grid_type = find_data_type (self.file_name)
        if self.grid_type == "STRUCTURED_GRID":
            self.reader = vtk.vtkStructuredGridReader ()
        elif self.grid_type == "STRUCTURED_POINTS":
            self.reader = vtk.vtkStructuredPointsReader ()
        elif self.grid_type == "RECTILINEAR_GRID":
            self.reader = vtk.vtkRectilinearGridReader ()
        elif self.grid_type == "POLYDATA":
            self.reader = vtk.vtkPolyDataReader ()
        elif self.grid_type == "UNSTRUCTURED_GRID":
            self.reader = vtk.vtkUnstructuredGridReader () 

def _create_structured_points_pyvtk(x, y, z):
    
    """Creates a vtkStructuredPoints object given input data in the
    form of arrays.  This uses pyvtk to do the job and generates a
    temporary file in the process.

    Input Arguments:
       x -- Array of x-coordinates.  These should be regularly spaced.

       y -- Array of y-coordinates.  These should be regularly spaced.

       z -- Array of z values for the x, y values given.  The values
       should be computed such that the z values are computed as x
       varies fastest and y next.

    """

    import pyvtk
    import tempfile, os
    
    nx = len(x)
    ny = len(y)
    nz = len(z)
    assert nx*ny == nz, "len(x)*len(y) != len(z).  "\
           "You passed nx=%d, ny=%d,  nz=%d"%(nx, ny, nz)

    xmin, ymin = x[0], y[0]
    dx, dy= (x[1] - x[0]), (y[1] - y[0])

    # create a vtk data file
    sp = pyvtk.StructuredPoints ((nx, ny, 1), (xmin, ymin, 0), (dx, dy, 1))
    pd = pyvtk.PointData(pyvtk.Scalars(z, name='Scalars',
                                       lookup_table="default"))
    d = pyvtk.VtkData(sp, pd, "Surf data")
    file_name = tempfile.mktemp(suffix='.vtk')
    d.tofile(file_name, format='ascii')

    # read the created file - yes this is circuitous but works for now.
    reader = vtk.vtkStructuredPointsReader()
    reader.SetFileName(file_name)
    reader.Update()
    # cleanup.
    os.remove(file_name)
    return reader.GetOutput()

def read_data(file_in,dens_inc=False):
    """
    Reads in data from file_in.vtk into x, y, z, b, u
    """
    import vtk as vtk
    from vtk.util.numpy_support import vtk_to_numpy
    from numpy import rollaxis, reshape, array

    #Some vtk bookkeeping and object stuff
    reader = vtk.vtkStructuredPointsReader()
    reader.SetFileName(file_in)
    reader.ReadAllVectorsOn()
    reader.ReadAllScalarsOn()
    reader.Update()
    data = reader.GetOutput()
    dim = data.GetDimensions()
    nx = dim[0] - 1
    ny = dim[1] - 1
    nz = dim[2] - 1
    D = data.GetSpacing()
    N = (nx, ny, nz)

    #Read the raw data from the vtk file
    #u = vtk_to_numpy(data.GetCellData().GetArray('velocity')).reshape(nx, ny, nz, 3,order='F')
    #b = vtk_to_numpy(data.GetCellData().GetArray('cell_centered_B')).reshape(nx, ny, nz, 3,order='F') 

    ## #Convert to C-order and make a the 1st array index the vector component
    #b = array([b[...,0].ravel(order='F').reshape(nx,ny,nz), b[...,1].ravel(order='F').reshape(nx,ny,nz), b[...,2].ravel(order='F').reshape(nx,ny,nz)])
    #u = array([u[...,0].ravel(order='F').reshape(nx,ny,nz), u[...,1].ravel(order='F').reshape(nx,ny,nz), u[...,2].ravel(order='F').reshape(nx,ny,nz)])

    ## #Read the raw data from the vtk file
    u = vtk_to_numpy(data.GetCellData().GetArray('velocity')).reshape(nz,ny,nx,3)
    b = vtk_to_numpy(data.GetCellData().GetArray('cell_centered_B')).reshape(nz, ny, nx,3) 

    if(dens_inc):
        rho = vtk_to_numpy(data.GetCellData().GetArray('density')).reshape(nz,ny,nx)
        return (N,D,b.T,u.T,rho.T)
    #u = ascontiguousarray(array([u[...,0], u[...,1], u[...,2]]))
    #b = ascontiguousarray(array([b[...,0], b[...,1], b[...,2]]))
    #u = u.transpose(3,2,1,0)
    #b = b.transpose(3,2,1,0)
                        
    return (N, D, b.T, u.T)

def readvtk(path, datatype='polydata'):
    """Read VTK-file.

    Args:
        path: Path to file.
        type: 'imagedata', 'polydata', 'unstructeredgrid'

    """
    if datatype=='imagedata':
        reader = vtk.vtkStructuredPointsReader()
    elif datatype=='polydata':
        reader = vtk.vtkPolyDataReader()
    elif datatype=='unstructeredgrid':
        reader = vtk.vtkUnstructuredGridReader()
    else:
        print 'Invalid datatype'
    reader.SetFileName(path)
    reader.Update()
    return reader.GetOutput()

def load_maps(mol,rendmod,gfx,atomcol):
	mol.reader = vtk.vtkStructuredPointsReader()
	mol.reader.SetFileName(mol.mod.dfn)
	mol.reader.Update() #by calling Update() we read the file
	mol.iso = vtk.vtkMarchingContourFilter()
	mol.iso.UseScalarTreeOn()
	mol.iso.ComputeNormalsOn()
	mol.iso.SetInputConnection(mol.reader.GetOutputPort())
	mol.iso.SetValue(0,mol.mod.isov*mol.mod.sigavg[0]+mol.mod.sigavg[1])
	clean = vtk.vtkCleanPolyData()
  	clean.SetInputConnection(mol.iso.GetOutputPort())
 	clean.ConvertStripsToPolysOn()
	smooth = vtk.vtkWindowedSincPolyDataFilter()
 	smooth.SetInputConnection(clean.GetOutputPort())
 	smooth.BoundarySmoothingOn()
 	smooth.GenerateErrorVectorsOn()
  	smooth.GenerateErrorScalarsOn()
  	smooth.NormalizeCoordinatesOn()
  	smooth.NonManifoldSmoothingOn()
  	smooth.FeatureEdgeSmoothingOn()
  	smooth.SetEdgeAngle(90)
	smooth.SetFeatureAngle(90)
	smooth.Update()
	mol.mapper = vtk.vtkOpenGLPolyDataMapper()
	mol.mapper.SetInputConnection(smooth.GetOutputPort()) ### <- connection here
	mol.mapper.ScalarVisibilityOff()
	mol.mapper.Update()
	mol.acteur= vtk.vtkOpenGLActor()
	mol.acteur.SetMapper(mol.mapper)
	mol.acteur.GetProperty().SetColor(mol.col)
	if rendmod==5:
		mol.acteur.GetProperty().SetRepresentationToSurface()
	elif rendmod==6:
		mol.acteur.GetProperty().SetRepresentationToWireframe()
	elif rendmod==7:
		mol.acteur.GetProperty().SetRepresentationToPoints()
	else :
		mol.acteur.GetProperty().SetRepresentationToSurface()
	mol.acteur.GetProperty().SetInterpolationToGouraud()
	mol.acteur.GetProperty().SetSpecular(.4)
	mol.acteur.GetProperty().SetSpecularPower(10)
	gfx.renderer.AddActor(mol.acteur)
	gfx.renwin.Render()

def readDisplacement(filename):
    reader = vtk.vtkStructuredPointsReader();reader.SetFileName(filename);
    reader.Update();
    grid = reader.GetOutput();
    Dimensions=grid.GetExtent();
    NumPixX=Dimensions[1]+1;
    NumPixY=Dimensions[3]+1;
    centre = grid.GetCenter();
    bounds = grid.GetBounds();
    reader.SetScalarsName(reader.GetScalarsNameInFile(0));reader.Update();structuredPoints = reader.GetOutput();
    pd=structuredPoints.GetPointData();
    scalars=pd.GetScalars(reader.GetScalarsNameInFile(0));
    numPoints = scalars.GetSize();
    U=numpy.zeros([numPoints/2,1]);
    V=numpy.zeros([numPoints/2,1]);
    for i in range(0,(numPoints-1)/2):
      U[i]=scalars.GetValue(2*i);
      V[i]=scalars.GetValue(2*i+1);
    return {'U':U.reshape(NumPixY,NumPixX), 'V':V.reshape(NumPixY,NumPixX) }

def read_vtk(filename, skipGrid = False, reshape = False):

    reader = vtkStructuredPointsReader()
    reader.SetFileName(filename)
    # reader.SetDataByteOrderToBigEndian()
    # reader.ReadAllVectorsOn()
    # reader.ReadAllScalarsOn()
    reader.Update()

    data = reader.GetOutput()

    dim = data.GetDimensions()
    vec = list(dim)
    vec = [i-1 for i in dim]
    vec.append(3)
    # print dim
    # print vec
    # print data.GetPointData()

    try: 
        u = VN.vtk_to_numpy(data.GetPointData().GetArray('Re'))
    except AttributeError:
        #older file format
        u = VN.vtk_to_numpy(data.GetPointData().GetArray('scalar_field'))
        

    x = np.zeros(data.GetNumberOfPoints())
    y = np.zeros(data.GetNumberOfPoints())
    z = np.zeros(data.GetNumberOfPoints())

    if not skipGrid:
        # Looping over all the point takes a while... perhaps there is a faster way
        for i in range(data.GetNumberOfPoints()):
                x[i],y[i],z[i] = data.GetPoint(i)
    
    if reshape:
        u = u.reshape(dim,order='F')
        x = x.reshape(dim,order='F')
        y = y.reshape(dim,order='F')
        z = z.reshape(dim,order='F')
    
    
    return (u, np.array([x,y,z]))