def __init__(self, vtk_filename=None, vtk_data=None):
        """
        Initiate Viwer

        Parameters
        ----------
        vtk_filename : str
            Input VTK filename
        """

        QDialog.__init__(self)
        self.initUI()

        ren = vtk.vtkRenderer()
        self.vtkWidget.GetRenderWindow().AddRenderer(ren)
        iren = self.vtkWidget.GetRenderWindow().GetInteractor()

        if vtk_filename is not None:
            # VTK file
            reader = vtk.vtkUnstructuredGridReader()
            reader.SetFileName(vtk_filename)
            reader.Update()
            vtkdata = reader.GetOutput()

        if vtk_data is not None:
            vtkdata = vtk_data

        # VTK surface
        surface = vtk.vtkDataSetSurfaceFilter()
        surface.SetInput(vtkdata)
        surface.Update()

        mapper = vtk.vtkDataSetMapper()
        mapper.SetInput(surface.GetOutput())

        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
        actor.GetProperty().EdgeVisibilityOff()
        # actor.GetProperty().SetEdgeColor(1,1,1)
        # actor.GetProperty().SetLineWidth(0.1)
        ren.AddActor(actor)

        # annot. cube
        axesActor = vtk.vtkAnnotatedCubeActor()
        axesActor.SetXPlusFaceText('R')
        axesActor.SetXMinusFaceText('L')
        axesActor.SetYMinusFaceText('H')
        axesActor.SetYPlusFaceText('F')
        axesActor.SetZMinusFaceText('A')
        axesActor.SetZPlusFaceText('P')
        axesActor.GetTextEdgesProperty().SetColor(1, 1, 0)
        axesActor.GetCubeProperty().SetColor(0, 0, 1)
        self.axes = vtk.vtkOrientationMarkerWidget()
        self.axes.SetOrientationMarker(axesActor)
        self.axes.SetInteractor(iren)
        self.axes.EnabledOn()
        self.axes.InteractiveOn()

        ren.ResetCamera()
        iren.Initialize()

    def ReadOBJFile(self, filename):
        '''
        @param filename: str 
        @rtype: None
        '''
        reader = vtk.vtkOBJReader()
        reader.SetFileName(filename)
        
        try:
            reader.Update()
            cleanFilter = vtk.vtkCleanPolyData()
            cleanFilter.SetInput(reader.GetOutput())
            cleanFilter.ConvertLinesToPointsOn()
            cleanFilter.ConvertPolysToLinesOn()
            cleanFilter.ConvertStripsToPolysOn()
            cleanFilter.PointMergingOn()
            cleanFilter.Update()
            surfaceFilter = vtk.vtkDataSetSurfaceFilter()
            surfaceFilter.SetInput(cleanFilter.GetOutput())
            surfaceFilter.Update()
            self.setDataSet(surfaceFilter.GetOutput())
            del cleanFilter
            del surfaceFilter
        except Exception, e:
#            del reader
            print e
            raise IOError, "Could not read the OBJ file! "

def get_number_of_points_and_boundary_faces(pod_mode_dir,calibrate_coefficients):
	filename = pod_mode_dir+"/spatial_meanfield.vtk"
	reader = vtk.vtkUnstructuredGridReader()
	reader.ReadAllScalarsOn()
	reader.ReadAllVectorsOn()
	reader.ReadAllTensorsOn()
	reader.SetFileName(filename)
	reader.Update()
	num_points = reader.GetOutput().GetNumberOfPoints()

	if( (calibrate_coefficients=="none") or (calibrate_coefficients=="constant") ):
		geom_filter = vtk.vtkDataSetSurfaceFilter()
		geom_filter.SetInput(reader.GetOutput())
		geom_filter.Update()
		boundary_faces = vtk.vtkPolyData()
		boundary_faces = geom_filter.GetOutput()

		point_to_cell_data = vtk.vtkPointDataToCellData()
		point_to_cell_data.SetInput(geom_filter.GetOutput())
		point_to_cell_data.Update()
		num_boundary_faces = point_to_cell_data.GetOutput().GetNumberOfCells()
	else:
		num_boundary_faces = 1

	return num_points, num_boundary_faces

def save_lesion(output_file, input_file, field, limits):
    reader = v.vtkXMLUnstructuredGridReader()
    reader.SetFileName(input_file)
    reader.Update()

    threshold = v.vtkThreshold()
    threshold.SetInput(reader.GetOutput())

    if limits[0] is None:
        threshold.ThresholdByLower(limits[1])
    elif limits[1] is None:
        threshold.ThresholdByUpper(limits[0])
    else:
        threshold.ThresholdBetween(*limits)

    threshold.SetInputArrayToProcess(0, 0, 0, v.vtkDataObject.FIELD_ASSOCIATION_CELLS, field)
    threshold.Update()

    extract_surface = v.vtkDataSetSurfaceFilter()
    extract_surface.SetInput(threshold.GetOutput())
    extract_surface.Update()

    writer = v.vtkXMLPolyDataWriter()
    writer.SetFileName(output_file)
    writer.SetInput(extract_surface.GetOutput())
    writer.Write()

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 testLinear(self):
        pts = vtk.vtkPoints()
        pts.SetNumberOfPoints(4)
        pts.InsertPoint(0, (0, 0, 0))
        pts.InsertPoint(1, (1, 0, 0))
        pts.InsertPoint(2, (0.5, 1, 0))
        pts.InsertPoint(3, (0.5, 0.5, 1))

        te = vtk.vtkTetra()
        ptIds = te.GetPointIds()
        for i in range(4):
            ptIds.SetId(i, i)

        ghosts = vtk.vtkUnsignedCharArray()
        ghosts.SetName("vtkGhostLevels")
        ghosts.SetNumberOfTuples(4)
        ghosts.SetValue(0, 1)
        ghosts.SetValue(1, 1)
        ghosts.SetValue(2, 1)
        ghosts.SetValue(3, 0)

        grid = vtk.vtkUnstructuredGrid()
        grid.Allocate(1, 1)
        grid.InsertNextCell(te.GetCellType(), te.GetPointIds())
        grid.SetPoints(pts)
        grid.GetPointData().AddArray(ghosts)

        dss = vtk.vtkDataSetSurfaceFilter()
        dss.SetInputData(grid)
        dss.Update()
        self.assertEqual(dss.GetOutput().GetNumberOfCells(), 3)

def putMaskOnVTKGrid(data,grid,actorColor=None,deep=True):
  #Ok now looking 
  msk = data.mask
  imsk =  VN.numpy_to_vtk(msk.astype(numpy.int).flat,deep=deep)
  mapper = None
  if msk is not numpy.ma.nomask:
      msk =  VN.numpy_to_vtk(numpy.logical_not(msk).astype(numpy.uint8).flat,deep=deep)
      if actorColor is not None:
          grid2 = vtk.vtkStructuredGrid()
          grid2.CopyStructure(grid)
          geoFilter = vtk.vtkDataSetSurfaceFilter()
          if grid.IsA("vtkStructuredGrid"):
              grid2.GetPointData().SetScalars(imsk)
              #grid2.SetCellVisibilityArray(imsk)
              p2c = vtk.vtkPointDataToCellData()
              p2c.SetInputData(grid2)
              geoFilter.SetInputConnection(p2c.GetOutputPort())
          else:
              grid2.GetCellData().SetScalars(imsk)
              geoFilter.SetInputData(grid)
          geoFilter.Update()
          mapper = vtk.vtkPolyDataMapper()
          mapper.SetInputData(geoFilter.GetOutput())
          lut = vtk.vtkLookupTable()
          lut.SetNumberOfTableValues(1)
          r,g,b = actorColor
          lut.SetNumberOfTableValues(2)
          lut.SetTableValue(0,r/100.,g/100.,b/100.)
          lut.SetTableValue(1,r/100.,g/100.,b/100.)
          mapper.SetLookupTable(lut)
          mapper.SetScalarRange(1,1)
      if grid.IsA("vtkStructuredGrid"):
          grid.SetPointVisibilityArray(msk)
          #grid.SetCellVisibilityArray(msk)
  return mapper

    def doNonLinear(self, ghosts, ncells):
        pts = vtk.vtkPoints()
        pts.SetNumberOfPoints(10)
        pts.InsertPoint(0, (0, 0, 0))
        pts.InsertPoint(1, (1, 0, 0))
        pts.InsertPoint(2, (0.5, 1, 0))
        pts.InsertPoint(3, (0.5, 0.5, 1))
        pts.InsertPoint(4, (0.5, 0, 0))
        pts.InsertPoint(5, (1.25, 0.5, 0))
        pts.InsertPoint(6, (0.25, 0.5, 0))
        pts.InsertPoint(7, (0.25, 0.25, 0.5))
        pts.InsertPoint(8, (0.75, 0.25, 0.5))
        pts.InsertPoint(9, (0.5, 0.75, 0.5))

        te = vtk.vtkQuadraticTetra()
        ptIds = te.GetPointIds()
        for i in range(10):
            ptIds.SetId(i, i)

        grid = vtk.vtkUnstructuredGrid()
        grid.Allocate(1, 1)
        grid.InsertNextCell(te.GetCellType(), te.GetPointIds())
        grid.SetPoints(pts)
        grid.GetPointData().AddArray(ghosts)

        ugg = vtk.vtkUnstructuredGridGeometryFilter()
        ugg.SetInputData(grid)

        dss = vtk.vtkDataSetSurfaceFilter()
        dss.SetNonlinearSubdivisionLevel(2)
        dss.SetInputConnection(ugg.GetOutputPort())
        dss.Update()
        self.assertEqual(dss.GetOutput().GetNumberOfCells(), ncells)

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

 def ConvertDataSetToSurface(algorithmOutputPort):
     dataSetSurfaceFilter = vtk.vtkDataSetSurfaceFilter()
     dataSetSurfaceFilter.SetInputConnection(algorithmOutputPort)
     dataSetSurfaceFilter.Update()
     polyData = vtk.vtkPolyData()
     polyData.ShallowCopy(dataSetSurfaceFilter.GetOutput())
     return polyData

    def __init__(self, renderer, data, type, index):
        PeacockActor.__init__(self, renderer)
        self.data = data
        self.type = type
        self.index = index

        self.solid_visible = False
        self.edges_visible = False

        self.mesh = data.GetBlock(type).GetBlock(index)

        self.geom = vtk.vtkDataSetSurfaceFilter()
        if vtk.VTK_MAJOR_VERSION <= 5:
            self.geom.SetInput(self.mesh)
        else:
            self.geom.SetInputData(self.mesh)
        self.geom.Update()

        self.mapper = vtk.vtkDataSetMapper()
        if vtk.VTK_MAJOR_VERSION <= 5:
            self.mapper.SetInput(self.mesh)
        else:
            self.mapper.SetInputData(self.mesh)

        self.actor = vtk.vtkActor()
        self.actor.SetMapper(self.mapper)
        self.actor.GetProperty().SetPointSize(5)
        self.actor.GetProperty().SetEdgeColor(0,0,0)
        self.actor.GetProperty().SetAmbient(0.3);

def threshold(polydata, arrayname, valuerange=[0, 1], iscelldata=True,
              allscalars=True):
    """Extract those cells from polydata whose pointdata/celldata values are
    within a specified range.

    For pointdata, cells are included if scalar values of all cell points are
    within the range (allscalars=True) or if the scalar value of at least one of
    the cell points is within the range (allscalar=False).

    """
    thresholdfilter = vtk.vtkThreshold()
    thresholdfilter.SetInput(polydata)
    thresholdfilter.ThresholdBetween(valuerange[0], valuerange[1])
    if iscelldata:
        thresholdfilter.SetInputArrayToProcess(0, 0, 0,
                                           vtk.vtkDataObject.
                                           FIELD_ASSOCIATION_CELLS, arrayname)
    else:
        thresholdfilter.SetInputArrayToProcess(0, 0, 0,
                                           vtk.vtkDataObject.
                                           FIELD_ASSOCIATION_POINTS, arrayname)
    if allscalars:
        thresholdfilter.AllScalarsOn()
    else:
        thresholdfilter.AllScalarsOff()
    thresholdfilter.Update()
    surfacefilter = vtk.vtkDataSetSurfaceFilter()
    surfacefilter.SetInput(thresholdfilter.GetOutput())
    surfacefilter.Update()
    return surfacefilter.GetOutput()

    def _plotInternalCustomBoxfill(self):
        """Implements the logic to render a custom boxfill."""
        self._mappers = []

        self._customBoxfillArgs = self._prepContours()
        tmpLevels = self._customBoxfillArgs["tmpLevels"]
        tmpColors = self._customBoxfillArgs["tmpColors"]
        tmpOpacities = self._customBoxfillArgs["tmpOpacities"]

        style = self._gm.fillareastyle

        luts = []
        geos = []
        wholeDataMin, wholeDataMax = vcs.minmax(self._originalData1)
        _colorMap = self.getColorMap()
        assert(style != 'solid' or len(tmpLevels) == 1)
        for i, l in enumerate(tmpLevels):
            # Ok here we are trying to group together levels can be, a join
            # will happen if: next set of levels continues where one left off
            # AND pattern is identical

            # TODO this should really just be a single polydata/mapper/actor:
            for j, color in enumerate(tmpColors[i]):
                mapper = vtk.vtkPolyDataMapper()
                lut = vtk.vtkLookupTable()
                th = vtk.vtkThreshold()
                th.ThresholdBetween(l[j], l[j + 1])
                th.SetInputConnection(self._vtkPolyDataFilter.GetOutputPort())
                geoFilter2 = vtk.vtkDataSetSurfaceFilter()
                geoFilter2.SetInputConnection(th.GetOutputPort())
                # Make the polydata output available here for patterning later
                geoFilter2.Update()
                geos.append(geoFilter2)
                mapper.SetInputConnection(geoFilter2.GetOutputPort())
                lut.SetNumberOfTableValues(1)
                r, g, b, a = self.getColorIndexOrRGBA(_colorMap, color)
                if style == 'solid':
                    tmpOpacity = tmpOpacities[j]
                    if tmpOpacity is None:
                        tmpOpacity = a / 100.
                    else:
                        tmpOpacity = tmpOpacities[j] / 100.
                    lut.SetTableValue(0, r / 100., g / 100., b / 100., tmpOpacity)
                else:
                    lut.SetTableValue(0, 1., 1., 1., 0.)
                mapper.SetLookupTable(lut)
                mapper.SetScalarRange(l[j], l[j + 1])
                luts.append([lut, [l[j], l[j + 1], False]])
                # Store the mapper only if it's worth it?
                # Need to do it with the whole slab min/max for animation
                # purposes
                if not (l[j + 1] < wholeDataMin or l[j] > wholeDataMax):
                    self._mappers.append(mapper)

        self._resultDict["vtk_backend_luts"] = luts
        if len(geos) > 0:
            self._resultDict["vtk_backend_geofilters"] = geos

def delaunay3d(points):
    """Construct a 3D Delaunay triangulation from a set of points."""
    delaunay = vtk.vtkDelaunay3D()
    delaunay.SetInput(points)
    delaunay.Update()
    surface = vtk.vtkDataSetSurfaceFilter()
    surface.SetInput(delaunay.GetOutput())
    surface.Update()
    return surface.GetOutput()

def write_stl(ugrid, filename):
    surface_filter = vtk.vtkDataSetSurfaceFilter()
    surface_filter.SetInputData(ugrid)

    triangle_filter = vtk.vtkTriangleFilter()
    triangle_filter.SetInputConnection(surface_filter.GetOutputPort())

    writer = vtk.vtkSTLWriter()
    writer.SetFileName(filename)
    writer.SetInputConnection(triangle_filter.GetOutputPort())
    writer.Write()

def cellthreshold(polydata, arrayname, start=0, end=1):
    threshold = vtk.vtkThreshold()
    threshold.SetInputData(polydata)
    threshold.SetInputArrayToProcess(0,0,0,vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS,arrayname)
    threshold.ThresholdBetween(start,end)
    threshold.Update()

    surfer = vtk.vtkDataSetSurfaceFilter()
    surfer.SetInputConnection(threshold.GetOutputPort())
    surfer.Update()
    return surfer.GetOutput()

def pointthreshold(polydata, arrayname, start=0, end=1,alloff=0):
    threshold = vtk.vtkThreshold()
    threshold.SetInputData(polydata)
    threshold.SetInputArrayToProcess(0,0,0,vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS,arrayname)
    threshold.ThresholdBetween(start,end)
    if (alloff):
        threshold.AllScalarsOff()
    threshold.Update()

    surfer = vtk.vtkDataSetSurfaceFilter()
    surfer.SetInputConnection(threshold.GetOutputPort())
    surfer.Update()
    return surfer.GetOutput()

 def _createPolyDataFilter(self):
     """Overrides baseclass implementation."""
     self._vtkPolyDataFilter = vtk.vtkDataSetSurfaceFilter()
     if self._useCellScalars:
         # Sets data to point instead of just cells
         c2p = vtk.vtkCellDataToPointData()
         c2p.SetInputData(self._vtkDataSet)
         c2p.Update()
         # For contouring duplicate points seem to confuse it
         self._vtkPolyDataFilter.SetInputConnection(c2p.GetOutputPort())
     else:
         self._vtkPolyDataFilter.SetInputData(self._vtkDataSet)
     self._resultDict["vtk_backend_filter"] = self._vtkPolyDataFilter

def main():
    pointSource = vtk.vtkPointSource()
    pointSource.SetNumberOfPoints(20)
    pointSource.Update()

    idFilter = vtk.vtkIdFilter()
    idFilter.SetInputConnection(pointSource.GetOutputPort())
    idFilter.SetIdsArrayName("OriginalIds")
    idFilter.Update()

    surfaceFilter = vtk.vtkDataSetSurfaceFilter()
    surfaceFilter.SetInputConnection(idFilter.GetOutputPort())
    surfaceFilter.Update()

    poly_input = surfaceFilter.GetOutput()

    # Create a mapper and actor
    mapper = vtk.vtkPolyDataMapper()

    if vtk.VTK_MAJOR_VERSION <= 5:
        mapper.SetInputConnection(poly_input.GetProducerPort())
    else:
        mapper.SetInputData(poly_input)
    mapper.ScalarVisibilityOff()

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

    # Visualize
    renderer = vtk.vtkRenderer()
    renderWindow = vtk.vtkRenderWindow()
    renderWindow.AddRenderer(renderer)

    areaPicker = vtk.vtkAreaPicker()
    renderWindowInteractor = vtk.vtkRenderWindowInteractor()
    renderWindowInteractor.SetPicker(areaPicker)
    renderWindowInteractor.SetRenderWindow(renderWindow)

    renderer.AddActor(actor)
    #renderer.SetBackground(1,1,1) # Background color white

    renderWindow.Render()

    style = vtk.vtkRenderWindowInteractor()
    #style = myInteractorStyle()
    #style = InteractorStyle()
    #style = QVTKRenderWindowInteractor()
    #style.SetPoints(poly_input)
    renderWindowInteractor.SetInteractorStyle(style)

    renderWindowInteractor.Start()

    def _plotInternalBoxfill(self):
        """Implements the logic to render a non-custom boxfill."""
        # Prep mapper
        mapper = vtk.vtkPolyDataMapper()
        self._mappers = [mapper]

        if self._gm.ext_1 and self._gm.ext_2:
            mapper.SetInputConnection(self._vtkPolyDataFilter.GetOutputPort())
            self._resultDict["vtk_backend_geofilters"] = \
                [self._vtkPolyDataFilter]
        else:
            thr = vtk.vtkThreshold()
            thr.SetInputConnection(self._vtkPolyDataFilter.GetOutputPort())
            if not self._gm.ext_1 and not self._gm.ext_2:
                thr.ThresholdBetween(self._contourLevels[0],
                                     self._contourLevels[-1])
            elif self._gm.ext_1 and not self._gm.ext_2:
                thr.ThresholdByLower(self._contourLevels[-1])
            elif not self._gm.ext_1 and self._gm.ext_2:
                thr.ThresholdByUpper(self._contourLevels[0])

            geoFilter2 = vtk.vtkDataSetSurfaceFilter()
            geoFilter2.SetInputConnection(thr.GetOutputPort())
            mapper.SetInputConnection(geoFilter2.GetOutputPort())
            self._resultDict["vtk_backend_geofilters"] = [geoFilter2]

        # Colortable bit
        # make sure length match
        numLevels = len(self._contourLevels)
        while len(self._contourColors) < numLevels:
            self._contourColors.append(self._contourColors[-1])

        lut = vtk.vtkLookupTable()
        lut.SetNumberOfTableValues(numLevels)
        _colorMap = self.getColorMap()
        for i in range(numLevels):
            r, g, b, a = self.getColorIndexOrRGBA(_colorMap, self._contourColors[i])
            lut.SetTableValue(i, r / 100., g / 100., b / 100., a / 100.)

        mapper.SetLookupTable(lut)
        if numpy.allclose(self._contourLevels[0], -1.e20):
            lmn = self._min - 1.
        else:
            lmn = self._contourLevels[0]
        if numpy.allclose(self._contourLevels[-1], 1.e20):
            lmx = self._mx + 1.
        else:
            lmx = self._contourLevels[-1]
        mapper.SetScalarRange(lmn, lmx)
        self._resultDict["vtk_backend_luts"] = [[lut, [lmn, lmx, True]]]