def visQuadFunc():
    """ vtk sample scene with iso contours """

    # VTK supports implicit functions of the form f(x,y,z)=constant. These
    # functions can represent things spheres, cones, etc. Here we use a
    # general form for a quadric to create an elliptical data field.
    quadric = vtk.vtkQuadric()
    quadric.SetCoefficients(.5, 1, .2, 0, .1, 0, 0, .2, 0, 0)

    # vtkSampleFunction samples an implicit function over the x-y-z range
    # specified (here it defaults to -1,1 in the x,y,z directions).
    sample = vtk.vtkSampleFunction()
    sample.SetSampleDimensions(30, 30, 30)
    sample.SetImplicitFunction(quadric)

    # Create five surfaces F(x,y,z) = constant between range specified. The
    # GenerateValues() method creates n isocontour values between the range
    # specified.
    contours = vtk.vtkContourFilter()
    contours.SetInputConnection(sample.GetOutputPort())
    contours.GenerateValues(8, 0.0, 1.2)

    contMapper = vtk.vtkPolyDataMapper()
    contMapper.SetInputConnection(contours.GetOutputPort())
    contMapper.SetScalarRange(0.0, 1.2)

    contActor = vtk.vtkActor()
    contActor.SetMapper(contMapper)

    # We'll put a simple outline around the data.
    outline = vtk.vtkOutlineFilter()
    outline.SetInputConnection(sample.GetOutputPort())

    outlineMapper = vtk.vtkPolyDataMapper()
    outlineMapper.SetInputConnection(outline.GetOutputPort())

    outlineActor = vtk.vtkActor()
    outlineActor.SetMapper(outlineMapper)
    outlineActor.GetProperty().SetColor(1, 0.5, 0)

    # extract data from the volume
    extract = vtk.vtkExtractVOI()
    extract.SetInputConnection(sample.GetOutputPort())
    extract.SetVOI(0, 29, 0, 29, 15, 15)
    extract.SetSampleRate(1, 2, 3)

    contours2 = vtk.vtkContourFilter()
    contours2.SetInputConnection(extract.GetOutputPort())
    contours2.GenerateValues(8, 0.0, 1.2)

    contMapper2 = vtk.vtkPolyDataMapper()
    contMapper2.SetInputConnection(contours2.GetOutputPort())
    contMapper2.SetScalarRange(0.0, 1.2)

    contActor2 = vtk.vtkActor()
    contActor2.SetMapper(contMapper2)

    return contActor, contActor2, outlineActor, contours, contours2

    def get_isosurface(self, iso_value=500):

        if not self.flag_read:
            sys.stderr.write('No Image Loaded!\n')
            return

        contour = vtk.vtkContourFilter()
        normals = vtk.vtkPolyDataNormals()
        stripper = vtk.vtkStripper()
        mapper = vtk.vtkPolyDataMapper()

        contour.SetInputData(self.reader)
        contour.SetValue(0, iso_value)

        normals.SetInputConnection(contour.GetOutputPort())
        normals.SetFeatureAngle(60.0)
        normals.ReleaseDataFlagOn()

        stripper.SetInputConnection(normals.GetOutputPort())
        stripper.ReleaseDataFlagOn()

        mapper.SetInputConnection(stripper.GetOutputPort())
        mapper.SetScalarVisibility(False)

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

        # Default colour, should be changed.
        actor.GetProperty().SetDiffuseColor(
            [247.0 / 255.0, 150.0 / 255.0, 155.0 / 255.0])  # Look like red
        actor.GetProperty().SetSpecular(0.3)
        actor.GetProperty().SetSpecularPower(20)

        return actor

    def __init__(self, module_manager, contourFilterText):

        # call parent constructor
        ModuleBase.__init__(self, module_manager)

        self._contourFilterText = contourFilterText
        if contourFilterText == 'marchingCubes':
            self._contourFilter = vtk.vtkMarchingCubes()
        else: # contourFilter == 'contourFilter'
            self._contourFilter = vtk.vtkContourFilter()

        module_utils.setup_vtk_object_progress(self, self._contourFilter,
                                           'Extracting iso-surface')

        # now setup some defaults before our sync
        self._config.isoValue = 128;

        self._viewFrame = None
        self._createViewFrame()

        # transfer these defaults to the logic
        self.config_to_logic()

        # then make sure they come all the way back up via self._config
        self.logic_to_config()
        self.config_to_view()

    def __init__ (self, mod_m): 
        debug ("In CustomGridPlane::__init__ ()")
        Common.state.busy ()
        Base.Objects.Module.__init__ (self, mod_m)
        self.act = None
        out = self.mod_m.GetOutput ()
        if out.IsA('vtkStructuredGrid'):
            self.plane = vtk.vtkStructuredGridGeometryFilter ()
        elif out.IsA ('vtkStructuredPoints') or out.IsA('vtkImageData'):
            if hasattr (vtk, 'vtkImageDataGeometryFilter'):
                self.plane = vtk.vtkImageDataGeometryFilter ()
            else:
                self.plane = vtk.vtkStructuredPointsGeometryFilter ()
        elif out.IsA ('vtkRectilinearGrid'):
            self.plane = vtk.vtkRectilinearGridGeometryFilter ()
        else:
            msg = "This module does not support the %s dataset."%(out.GetClassName())
            raise Base.Objects.ModuleException, msg

        self.cont_fil = vtk.vtkContourFilter ()
        self.mapper = self.map = vtk.vtkPolyDataMapper ()
        self.actor = self.act = vtk.vtkActor ()
        self._initialize ()
        self._gui_init ()
        self.renwin.Render ()
        Common.state.idle ()

 def getContour(self,value=0.5,arrName='',largestRegion=False):
     if arrName == '':
         arrName = self.waterArray[self.solver]
     contour=vtk.vtkContourFilter()
     contour.SetValue(0,value)
     self.pointData.SetActiveScalars(arrName)
     #contour.SetInputArrayToProcess(1, 0,0, vtkDataObject::FIELD_ASSOCIATION_POINTS , arrName);  #something like this may also work
     contour.SetInputConnection(self.outputPort)
     contour.Update()
     if largestRegion:
         conn=vtk.vtkConnectivityFilter()
         conn.SetInputConnection(contour.GetOutputPort())
         conn.SetExtractionModeToLargestRegion()
         conn.Update()
         connOutput = conn.GetOutput()
         IDs = []
         for iC in range(connOutput.GetNumberOfCells()):  #it has to be so complicated, no polylines are given
             cell = connOutput.GetCell(iC)
             for i in [0,1]:
                 ID = cell.GetPointId(i)
                 if not ID in IDs: IDs.append(ID)
         points=[connOutput.GetPoint(i) for i in IDs]
     else:
         cOutput = contour.GetOutput()
         points=[cOutput.GetPoint(i) for i in range(cOutput.GetNumberOfPoints())]
     
     if not points: sys.exit('vtkextract.py: No contours found - wrong initialization of water fraction?')
     
     points = np.array(points)
     return points

def create_actors_for_skin_and_bone(reader):
    actors_list = []
    for contour_val, color, opacity in SKIN_BONE_LIST:
	contour = vtk.vtkContourFilter()
	contour.SetInput(reader.GetOutput())
	contour.SetNumberOfContours(1)
	contour.SetValue(contour_val[0], contour_val[1])

	normals = vtk.vtkPolyDataNormals()
	normals.SetInput(contour.GetOutput())
	normals.SetFeatureAngle(60)
	normals.ConsistencyOff()
	normals.SplittingOff()

	mapper = vtk.vtkPolyDataMapper()
	mapper.SetInput(normals.GetOutput())
	mapper.ScalarVisibilityOff()

	actor = vtk.vtkActor()
	actor.SetMapper(mapper)
	actor.GetProperty().SetColor(color)
        actor.GetProperty().SetOpacity(opacity)
	actor.RotateX(-90)
        actors_list.append(actor)
    return actors_list

def render_image(vtk_reader):
    contour = vtk.vtkContourFilter()
    contour.SetInput(vtk_reader.GetOutput())
    contour.GenerateValues(20,20,200)

    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInput(contour.GetOutput())
    mapper.ScalarVisibilityOn()
    mapper.SetScalarRange(0,255)

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

    renderer = vtk.vtkRenderer()
    renderWindow = vtk.vtkRenderWindow()
    renderWindow.AddRenderer(renderer)
    renderWindowInteractor = vtk.vtkRenderWindowInteractor()
    renderWindowInteractor.SetRenderWindow(renderWindow)

    renderer.AddActor(actor)
    renderer.SetBackground(.5, .5, .5)

    renderWindow.SetSize(600, 600)
    renderWindow.Render()
    renderWindowInteractor.Start()

    def setUp(self):
        self.vtk_iso = vtkContourFilter()
        # self.vtk_iso.SetInput(...)

        self.vtk_dnorm = vtkPolyDataNormals()
        self.vtk_subdiv = vtkLinearSubdivisionFilter()
        self.vtk_dmap = vtkPolyDataMapper()

 def __init__(self, module_manager):
     SimpleVTKClassModuleBase.__init__(
         self, module_manager,
         vtk.vtkContourFilter(), 'Processing.',
         ('vtkDataSet',), ('vtkPolyData',),
         replaceDoc=True,
         inputFunctions=None, outputFunctions=None)

 def BuildPolyBallSurface(self):
     #Build a surface for displaying the polyball
     if self.PolyBall == None:
         return
     
     #Sample the polyball
     sampler = vtk.vtkSampleFunction()
     sampler.SetImplicitFunction(self.PolyBall)
     
     #Set the bounds to be slightly larger than those of the mesh
     meshBounds = self.Mesh.GetBounds()
     meshCenter = self.Mesh.GetCenter()
     polyBallBounds = [0, 0, 0, 0, 0, 0]
     for i in range(0,3):
         length = 1.2*(meshBounds[2*i+1] - meshCenter[i])
         polyBallBounds[2*i] = meshCenter[i] - length
         polyBallBounds[2*i+1] = meshCenter[i] + length
     
     sampler.SetModelBounds(polyBallBounds)
     sampler.SetSampleDimensions(self.PolyBallResolution)
     sampler.ComputeNormalsOff()
     sampler.Update()
     
     #Extract the isosurface at 0
     contour = vtk.vtkContourFilter()
     contour.SetInput(sampler.GetOutput())
     contour.SetValue(0,0.)
     contour.Update()
     
     #Set the new model as the mapper input
     self.PolyBallActor.GetMapper().SetInput(contour.GetOutput())

def GetXandt(filelist):
  time = []
  X_ns = []
  X_fs = []
  for files in filelist:

      data = vtktools.vtu(files) 
      
      time.append(data.GetScalarField("Time")[0])
      
      # Get X
      data.ugrid.GetPointData().SetActiveScalars('Temperature')
      data = data.ugrid
      
      contour = vtk.vtkContourFilter()
      if vtk.vtkVersion.GetVTKMajorVersion() <= 5:
        contour.SetInput(data)
      else:
        contour.SetInputData(data)
      contour.SetValue(0, 0.0)
      contour.Update()
      polydata = contour.GetOutput()

      bounding_box = polydata.GetBounds()
   
      X_ns.append(bounding_box[1])
      X_fs.append(bounding_box[0])

  return time, X_ns, X_fs

    def __init__(self,data_reader,origin,normal,camera_normal):
        self.axial=self.get_matrix(data_reader,origin,normal,camera_normal)
        # Extract a slice in the desired orientation
        self.reslice = vtk.vtkImageReslice()
        self.reslice.SetInput(data_reader.get_data_set())
        self.reslice.SetOutputDimensionality(2)
        self.reslice.SetResliceAxes(self.axial)
        self.reslice.SetInterpolationModeToLinear()
        
        self.contour=vtk.vtkContourFilter()
        self.contour.SetInputConnection(self.reslice.GetOutputPort())

        self.contour.GenerateValues(25, data_reader.get_scalar_range())
        self.contour.ComputeScalarsOn()
        self.contour.ComputeGradientsOn()      
        self.cutmapper=vtk.vtkPolyDataMapper()
        self.cutmapper.SetInputConnection(self.contour.GetOutputPort())  
        self.actor=vtk.vtkActor()
        self.actor.SetMapper(self.cutmapper)
        self.actor.PokeMatrix(self.axial)
        c="c"
        if origin[0]==c or origin[1]==c or origin[2]==c:
			origin=self.center
        
        origin=(float(origin[0]),float(origin[1]),float(origin[2]))
        self.actor.SetOrigin(origin)		  

    def test_contours(self):
        cell = vtk.vtkUnstructuredGrid()
        cell.ShallowCopy(self.Cell)

        np = self.Cell.GetNumberOfPoints()
        ncomb = pow(2, np)

        scalar = vtk.vtkDoubleArray()
        scalar.SetName("scalar")
        scalar.SetNumberOfTuples(np)
        cell.GetPointData().SetScalars(scalar)

        incorrectCases = []
        for i in range(1,ncomb-1):
            c = Combination(np, i)
            for p in range(np):
                scalar.SetTuple1(p, c[p])

            gradientFilter = vtk.vtkGradientFilter()
            gradientFilter.SetInputData(cell)
            gradientFilter.SetInputArrayToProcess(0,0,0,0,'scalar')
            gradientFilter.SetResultArrayName('grad')
            gradientFilter.Update()

            contourFilter = vtk.vtkContourFilter()
            contourFilter.SetInputConnection(gradientFilter.GetOutputPort())
            contourFilter.SetNumberOfContours(1)
            contourFilter.SetValue(0, 0.5)
            contourFilter.Update()

            normalsFilter = vtk.vtkPolyDataNormals()
            normalsFilter.SetInputConnection(contourFilter.GetOutputPort())
            normalsFilter.SetConsistency(0)
            normalsFilter.SetFlipNormals(0)
            normalsFilter.SetSplitting(0)

            calcFilter = vtk.vtkArrayCalculator()
            calcFilter.SetInputConnection(normalsFilter.GetOutputPort())
            calcFilter.SetAttributeTypeToPointData()
            calcFilter.AddVectorArrayName('grad')
            calcFilter.AddVectorArrayName('Normals')
            calcFilter.SetResultArrayName('dir')
            calcFilter.SetFunction('grad.Normals')
            calcFilter.Update()

            out = vtk.vtkUnstructuredGrid()
            out.ShallowCopy(calcFilter.GetOutput())

            numPts = out.GetNumberOfPoints()
            if numPts > 0:
                dirArray = out.GetPointData().GetArray('dir')
                for p in range(numPts):
                    if(dirArray.GetTuple1(p) > 0.0): # all normals are reversed
                        incorrectCases.append(i)
                        break

        self.assertEquals(','.join([str(i) for i in incorrectCases]), '')

    def test_method_signature(self):
        """Check if VTK method signatures are parsed correctly."""
        p = self.p

        # Simple tests.
        o = vtk.vtkProperty()
        self.assertEqual([(['string'], None)],
                         p.get_method_signature(o.GetClassName))
        if hasattr(vtk, 'vtkArrayCoordinates'):
            self.assertEqual([([('float', 'float', 'float')], None),
                              ([None], (['float', 'float', 'float'],)),
                              ([None], ('float', 'float', 'float'))],
                             p.get_method_signature(o.GetColor))
        else:
            self.assertEqual([([('float', 'float', 'float')], None),
                              ([None], (('float', 'float', 'float'),))],
                             p.get_method_signature(o.GetColor))
        if hasattr(vtk, 'vtkArrayCoordinates'):
            self.assertEqual([([None], ('float', 'float', 'float')),
                ([None], (['float', 'float', 'float'],))],
                             p.get_method_signature(o.SetColor))

        else:
            self.assertEqual([([None], ('float', 'float', 'float')),
                              ([None], (('float', 'float', 'float'),))],
                             p.get_method_signature(o.SetColor))

        # Get VTK version to handle changed APIs.
        vtk_ver = vtk.vtkVersion().GetVTKVersion()

        # Test vtkObjects args.
        o = vtk.vtkContourFilter()
        sig = p.get_method_signature(o.SetInput)
        if len(sig) == 1:
            self.assertEqual([([None], ['vtkDataSet'])],
                             sig)
        elif vtk_ver[:3] in ['4.2', '4.4']:
            self.assertEqual([([None], ['vtkDataObject']),
                              ([None], ('int', 'vtkDataObject')),
                              ([None], ['vtkDataSet']),
                              ([None], ('int', 'vtkDataSet'))
                              ], sig)
        elif vtk_ver[:2] == '5.' or vtk_ver[:3] == '4.5':
            self.assertEqual([([None], ['vtkDataObject']),
                              ([None], ('int', 'vtkDataObject')),
                              ], sig)

        self.assertEqual([(['vtkPolyData'], None),
                          (['vtkPolyData'], ['int'])],
                         p.get_method_signature(o.GetOutput))

        # Test if function arguments work.
        self.assertEqual([(['int'], ('int', 'function'))],
                         p.get_method_signature(o.AddObserver))
        # This one's for completeness.
        self.assertEqual([([None], ['int'])],
                         p.get_method_signature(o.RemoveObserver))

 def test_to_tvtk_returns_tvtk_object(self):
     # Given
     v = vtk.vtkContourFilter()
     # When
     x = tvtk.to_tvtk(v)
     # Then
     self.assertEqual(x.class_name, 'vtkContourFilter')
     self.assertTrue(isinstance(x, tvtk_base.TVTKBase))
     self.assertTrue(isinstance(x, tvtk.ContourFilter))
     self.assertTrue(v is x._vtk_obj)

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 __init__(self, parent = None):
        super(VTKFrame, self).__init__(parent)

        self.vtkWidget = QVTKRenderWindowInteractor(self)
        vl = QtGui.QVBoxLayout(self)
        vl.addWidget(self.vtkWidget)
        vl.setContentsMargins(0, 0, 0, 0)
 
        self.ren = vtk.vtkRenderer()
        self.vtkWidget.GetRenderWindow().AddRenderer(self.ren)
        self.iren = self.vtkWidget.GetRenderWindow().GetInteractor()
 
        # Create cylinder
        cylinder = vtk.vtkCylinder()
        cylinder.SetCenter(0, 0, 0)
        cylinder.SetRadius(1.0)
        
        # Create plane
        plane = vtk.vtkPlane()
        plane.SetOrigin(0, 0, 0)
        plane.SetNormal(0, -1, 0)

        # Cut the cylinder
        cuted_cylinder = vtk.vtkImplicitBoolean()
        cuted_cylinder.SetOperationTypeToIntersection()
        #cuted_cylinder.SetOperationTypeToUnion()
        cuted_cylinder.AddFunction(cylinder)
        cuted_cylinder.AddFunction(plane)

        # Sample 
        sample = vtk.vtkSampleFunction()
        sample.SetImplicitFunction(cuted_cylinder)
        sample.SetModelBounds(-1.5 , 1.5 , -1.5 , 1.5 , -1.5 , 1.5)
        sample.SetSampleDimensions(60, 60, 60)
        sample.SetComputeNormals(0)

        #
        surface = vtk.vtkContourFilter()
        #surface.SetInput(sample.GetOutput())
        surface.SetInputConnection(sample.GetOutputPort())
 
        # Create a mapper
        mapper = vtk.vtkPolyDataMapper()
        #mapper.SetInput(surface.GetOutput())
        mapper.SetInputConnection(surface.GetOutputPort())
 
        # Create an actor
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
 
        self.ren.AddActor(actor)
        self.ren.ResetCamera()

        self._initialized = False

def TestDataType(dataType, reader, writer, ext, numTris, useSubdir=False):
    s = GetSource(dataType)

    filename = VTK_TEMP_DIR + "/%s.p%s" % (dataType, ext)

    writer.SetInputConnection(s.GetOutputPort())
    npieces = 16
    writer.SetNumberOfPieces(npieces)
    pperrank = npieces // nranks
    start = pperrank * rank
    end = start + pperrank - 1
    writer.SetStartPiece(start)
    writer.SetEndPiece(end)
    writer.SetFileName(filename)
    if useSubdir:
        writer.SetUseSubdirectory(True)
    #writer.SetDataModeToAscii()
    writer.Write()

    if contr:
        contr.Barrier()

    reader.SetFileName(filename)

    cf = vtk.vtkContourFilter()
    cf.SetValue(0, 130)
    cf.SetComputeNormals(0)
    cf.SetComputeGradients(0)
    cf.SetInputConnection(reader.GetOutputPort())
    cf.UpdateInformation()
    cf.GetOutputInformation(0).Set(vtk.vtkStreamingDemandDrivenPipeline.UPDATE_NUMBER_OF_PIECES(), nranks)
    cf.GetOutputInformation(0).Set(vtk.vtkStreamingDemandDrivenPipeline.UPDATE_PIECE_NUMBER(), rank)
    cf.Update()

    ntris = cf.GetOutput().GetNumberOfCells()
    da = vtk.vtkIntArray()
    da.InsertNextValue(ntris)

    da2 = vtk.vtkIntArray()
    da2.SetNumberOfTuples(1)
    contr.AllReduce(da, da2, vtk.vtkCommunicator.SUM_OP)

    if rank == 0:
        print(da2.GetValue(0))
        import os
        os.remove(filename)
        for i in range(npieces):
            if not useSubdir:
                os.remove(VTK_TEMP_DIR + "/%s_%d.%s" % (dataType, i, ext))
            else:
                os.remove(VTK_TEMP_DIR + "/%s/%s_%d.%s" %(dataType, dataType, i, ext))

    assert da2.GetValue(0) == numTris

def createGeometry(reader, threshold = 0):
    isoSurface = vtk.vtkContourFilter()
    isoSurface.SetInputConnection(reader.GetOutputPort())
    isoSurface.SetValue(0, threshold)

    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(isoSurface.GetOutputPort())

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

    return actor

def get_interface_depth(file):
  vtu=vtktools.vtu(file)
  data = vtu.ugrid
  data.GetPointData().SetActiveScalars("Dense::MaterialVolumeFraction")
  contour = vtk.vtkContourFilter ()
  contour.SetInput(data)
  contour.SetValue(0, 0.5)
  contour.Update()
  polydata = contour.GetOutput()
  bounding_box = polydata.GetBounds()
  interface_depth = bounding_box[2]
  return interface_depth