def testUnstructured(self):
        rt = vtk.vtkRTAnalyticSource()
        rt.SetWholeExtent(-5, 5, -5, 5, -5, 5)

        t = vtk.vtkThreshold()
        t.SetInputConnection(rt.GetOutputPort())
        t.ThresholdByUpper(-10)

        s = vtk.vtkSphere()
        s.SetRadius(2)
        s.SetCenter(0,0,0)

        c = vtk.vtkTableBasedClipDataSet()
        c.SetInputConnection(t.GetOutputPort())
        c.SetClipFunction(s)
        c.SetInsideOut(1)

        c.Update()

        self.assertEqual(c.GetOutput().GetNumberOfCells(), 64)

        eg = vtk.vtkEnSightGoldReader()
        eg.SetCaseFileName(VTK_DATA_ROOT + "/Data/EnSight/elements.case")
        eg.Update()

        pl = vtk.vtkPlane()
        pl.SetOrigin(3.5, 3.5, 0.5)
        pl.SetNormal(0, 0, 1)

        c.SetInputConnection(eg.GetOutputPort())
        c.SetClipFunction(pl)
        c.SetInsideOut(1)

        c.Update()
        data = c.GetOutputDataObject(0).GetBlock(0)
        self.assertEqual(data.GetNumberOfCells(), 75)

        rw = vtk.vtkRenderWindow()
        ren = vtk.vtkRenderer()
        rw.AddRenderer(ren)
        mapper = vtk.vtkDataSetMapper()
        mapper.SetInputData(data)
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
        ren.AddActor(actor)
        ac = ren.GetActiveCamera()
        ac.SetPosition(-7.9, 9.7, 14.6)
        ac.SetFocalPoint(3.5, 3.5, 0.5)
        ac.SetViewUp(0.08, 0.93, -0.34)
        rw.Render()
        ren.ResetCameraClippingRange()

        rtTester = vtk.vtkTesting()
        for arg in sys.argv[1:]:
            rtTester.AddArgument(arg)
        rtTester.AddArgument("-V")
        rtTester.AddArgument("tableBasedClip.png")
        rtTester.SetRenderWindow(rw)
        rw.Render()
        rtResult = rtTester.RegressionTest(10)

    def SphereDrill(self, m_holelist, m_holeRadius, m_quiet=False):
        """
        Drill sphere at locations specified by m_holelist.

        :param m_holelist:      [list]  A list of coordinates where holes are to be drilled
        :param m_holeRadius:    [float] The radius of the hole to drill
        :param m_quiet:         [bool]
        :return:
        """
        m_totalNumOfHoles = len(m_holelist)
        if not m_quiet:
            t = time.time()
            print "Drilling"

        for i in xrange(m_totalNumOfHoles):
            m_sphere = vtk.vtkSphere()
            m_sphere.SetCenter(m_holelist[i])
            m_sphere.SetRadius(m_holeRadius)

            clipper = vtk.vtkClipPolyData()
            clipper.SetInputData(self._data)
            clipper.SetClipFunction(m_sphere)
            clipper.Update()

            clipped = clipper.GetOutput()
            self._data.DeepCopy(clipped)

            if not m_quiet:
                print "\t%s/%s -- %.2f %%" % (i + 1, m_totalNumOfHoles, (i + 1) * 100 / float(m_totalNumOfHoles))
        if not m_quiet:
            print "Finished: Totaltime used = %.2f s" % (time.time() - t)
        pass

    def _Clip(self, pd):
        # The plane implicit function will be >0 for all the points in the positive side
        # of the plane (i.e. x s.t. n.(x-o)>0, where n is the plane normal and o is the 
        # plane origin).
        plane = vtkPlane()
        plane.SetOrigin(self.Iolet.Centre.x, self.Iolet.Centre.y, self.Iolet.Centre.z)
        plane.SetNormal(self.Iolet.Normal.x, self.Iolet.Normal.y, self.Iolet.Normal.z)
        
        # The sphere implicit function will be >0 for all the points outside the sphere.
        sphere = vtkSphere()
        sphere.SetCenter(self.Iolet.Centre.x, self.Iolet.Centre.y, self.Iolet.Centre.z)
        sphere.SetRadius(self.Iolet.Radius)
        
        # The VTK_INTERSECTION operator takes the maximum value of all the registered 
        # implicit functions. This will result in the function evaluating to >0 for all 
        # the points outside the sphere plus those inside the sphere in the positive 
        # side of the plane.
        clippingFunction = vtkImplicitBoolean()
        clippingFunction.AddFunction(plane)
        clippingFunction.AddFunction(sphere)
        clippingFunction.SetOperationTypeToIntersection() 
        
        clipper = vtkClipPolyData()
        clipper.SetInput(pd)
        clipper.SetClipFunction(clippingFunction)

        # Filter to get part closest to seed point
        connectedRegionGetter = vtkPolyDataConnectivityFilter()
        connectedRegionGetter.SetExtractionModeToClosestPointRegion()
        connectedRegionGetter.SetClosestPoint(*self.SeedPoint)
        connectedRegionGetter.SetInputConnection(clipper.GetOutputPort())
        connectedRegionGetter.Update()
        return connectedRegionGetter.GetOutput()

    def _add_ugrid_nodes_to_grid(self, name, diff_node_ids, nodes):
        """
        based on:
          _add_nastran_nodes_to_grid
        """
        nnodes = nodes.shape[0]
        assert nnodes > 0, nnodes
        # if nnodes == 0:
            # return
        nnodes = len(diff_node_ids)
        points = vtk.vtkPoints()
        points.SetNumberOfPoints(nnodes)

        for nid in diff_node_ids:
            node = nodes[nid, :]
            print('nid=%s node=%s' % (nid, node))
            points.InsertPoint(nid, *node)

            if 1:
                elem = vtk.vtkVertex()
                elem.GetPointIds().SetId(0, nid)
            else:
                elem = vtk.vtkSphere()
                sphere_size = self._get_sphere_size(dim_max)
                elem.SetRadius(sphere_size)
                elem.SetCenter(points.GetPoint(nid))

            self.alt_grids[name].InsertNextCell(elem.GetCellType(), elem.GetPointIds())
        self.alt_grids[name].SetPoints(points)

    def testStructured(self):
        rt = vtk.vtkRTAnalyticSource()
        rt.SetWholeExtent(-5, 5, -5, 5, -5, 5)
        rt.Update()
        i = rt.GetOutput()

        st = vtk.vtkStructuredGrid()
        st.SetDimensions(i.GetDimensions())

        nps = i.GetNumberOfPoints()
        ps = vtk.vtkPoints()
        ps.SetNumberOfPoints(nps)
        for idx in xrange(nps):
            ps.SetPoint(idx, i.GetPoint(idx))

        st.SetPoints(ps)

        s = vtk.vtkSphere()
        s.SetRadius(2)
        s.SetCenter(0,0,0)

        c = vtk.vtkTableBasedClipDataSet()
        c.SetInputData(st)
        c.SetClipFunction(s)
        c.SetInsideOut(1)

        c.Update()

        self.assertEqual(c.GetOutput().GetNumberOfCells(), 64)

    def openSurfaceAtPoint(self, polyData, seed):
        '''
        Returns a new surface with an opening at the given seed.
        '''

        someradius = 1.0

        pointLocator = vtk.vtkPointLocator()
        pointLocator.SetDataSet(polyData)
        pointLocator.BuildLocator()

        # find the closest point next to the seed on the surface
        # id = pointLocator.FindClosestPoint(int(seed[0]),int(seed[1]),int(seed[2]))
        id = pointLocator.FindClosestPoint(seed)

        # the seed is now guaranteed on the surface
        seed = polyData.GetPoint(id)

        sphere = vtk.vtkSphere()
        sphere.SetCenter(seed[0], seed[1], seed[2])
        sphere.SetRadius(someradius)

        clip = vtk.vtkClipPolyData()
        clip.SetInputData(polyData)
        clip.SetClipFunction(sphere)
        clip.Update()

        outPolyData = vtk.vtkPolyData()
        outPolyData.DeepCopy(clip.GetOutput())

        return outPolyData

def MaskWithPatch(id,t,c,r,maskArray,centerlines,voronoi):
   patch = ExtractPatch(id,centerlines)

   tubeFunction = vtkvmtk.vtkvmtkPolyBallLine()
   tubeFunction.SetInput(patch)
   tubeFunction.SetPolyBallRadiusArrayName(radiusArrayName)

   lastSphere = vtk.vtkSphere()
   lastSphere.SetRadius(r*1.5)
   lastSphere.SetCenter(c)

   for i in range(voronoi.GetNumberOfPoints()):
      point = [0.0,0.0,0.0]
      voronoiVector = [0.0,0.0,0.0]

      voronoi.GetPoint(i,point)
      voronoiVector[0] = point[0]-c[0]
      voronoiVector[1] = point[1]-c[1]
      voronoiVector[2] = point[2]-c[2]
      voronoiVectorDot = vtk.vtkMath.Dot(voronoiVector,t)   

      tubevalue = tubeFunction.EvaluateFunction(point)
      spherevalue = lastSphere.EvaluateFunction(point)
      if (spherevalue<0.0) & (voronoiVectorDot<0.0): continue
      elif (tubevalue<=0.0):
         maskArray.SetTuple1(i,1)

    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.ren.SetBackground(0.1, 0.2, 0.4)
        self.vtkWidget.GetRenderWindow().AddRenderer(self.ren)
        self.iren = self.vtkWidget.GetRenderWindow().GetInteractor()
 
        sphere1 = vtk.vtkSphere()
        sphere1.SetCenter(0.9, 0, 0)
        sphere2 = vtk.vtkSphere()
        sphere2.SetCenter(-0.9, 0, 0)

        implicitBoolean = vtk.vtkImplicitBoolean()
        implicitBoolean.AddFunction(sphere1)
        implicitBoolean.AddFunction(sphere2)
        implicitBoolean.SetOperationTypeToUnion()
        #implicitBoolean.SetOperationTypeToIntersection()

        # Sample the function
        sample = vtk.vtkSampleFunction()
        sample.SetSampleDimensions(50, 50, 50)
        sample.SetImplicitFunction(implicitBoolean)
        sample.SetModelBounds(-3, 3, -3, 3, -3, 3)

        # Create the 0 isosurface
        contours = vtk.vtkContourFilter()
        contours.SetInputConnection(sample.GetOutputPort())
        contours.GenerateValues(1, 1, 1)
 
        # Create a mapper
        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInputConnection(contours.GetOutputPort())
 
        # Create an actor
        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
 
        self.ren.AddActor(actor)
        self.ren.ResetCamera()

        self._initialized = False

    def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: no Surface.')

        if self.WidgetType == "box":
            self.ClipFunction = vtk.vtkPlanes()
        elif self.WidgetType == "sphere":
            self.ClipFunction = vtk.vtkSphere()

        self.Clipper = vtk.vtkClipPolyData()
        self.Clipper.SetInput(self.Surface)
        self.Clipper.SetClipFunction(self.ClipFunction)
        self.Clipper.GenerateClippedOutputOn()
        self.Clipper.InsideOutOn()
        
        if not self.vmtkRenderer:
            self.vmtkRenderer = vmtkrenderer.vmtkRenderer()
            self.vmtkRenderer.Initialize()
            self.OwnRenderer = 1

        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInput(self.Surface)
        mapper.ScalarVisibilityOff()
        self.Actor = vtk.vtkActor()
        self.Actor.SetMapper(mapper)
        self.vmtkRenderer.Renderer.AddActor(self.Actor)

        if self.WidgetType == "box":
            self.ClipWidget = vtk.vtkBoxWidget()
            self.ClipWidget.GetFaceProperty().SetColor(0.6,0.6,0.2)
            self.ClipWidget.GetFaceProperty().SetOpacity(0.25)
        elif self.WidgetType == "sphere":
            self.ClipWidget = vtk.vtkSphereWidget()
            self.ClipWidget.GetSphereProperty().SetColor(0.6,0.6,0.2)
            self.ClipWidget.GetSphereProperty().SetOpacity(0.25)
            self.ClipWidget.GetSelectedSphereProperty().SetColor(0.6,0.0,0.0)
            self.ClipWidget.GetSelectedSphereProperty().SetOpacity(0.75)
            self.ClipWidget.SetRepresentationToSurface()
            self.ClipWidget.SetPhiResolution(20)
            self.ClipWidget.SetThetaResolution(20)

        self.ClipWidget.SetInteractor(self.vmtkRenderer.RenderWindowInteractor)
        self.Display()

        self.Transform = vtk.vtkTransform()
        self.ClipWidget.GetTransform(self.Transform)

        if self.OwnRenderer:
            self.vmtkRenderer.Deallocate()

        if self.CleanOutput == 1:
            cleaner = vtk.vtkCleanPolyData()
            cleaner.SetInput(self.Surface)
            cleaner.Update()
            self.Surface = cleaner.GetOutput()

        if self.Surface.GetSource():
            self.Surface.GetSource().UnRegisterAllOutputs()

    def getImplicitSphere(self):
        impSphere = vtk.vtkSphere()  # implicit sphere
        centerPoint = self.getCentralPoint()
        radius = self.getRadius()

        impSphere.SetRadius(radius)
        impSphere.SetCenter(centerPoint)
        return impSphere

    def testImage(self):
        r = vtk.vtkRTAnalyticSource()
        r.SetWholeExtent(-5, 5, -5, 5, -5, 5)
        r.Update()

        s = vtk.vtkSphere()
        s.SetRadius(2)
        s.SetCenter(0,0,0)

        c = vtk.vtkTableBasedClipDataSet()
        c.SetInputConnection(r.GetOutputPort())
        c.SetClipFunction(s)
        c.SetInsideOut(1)

        c.Update()

        self.assertEqual(c.GetOutput().GetNumberOfCells(), 64)

def FindClippingPointOnParentArtery(centerlines,parentCenterlines,toll):
   divergingPointID = -1
   divergingPoint = [0.0,0.0,0.0]
   divergingPointMISR = -1

   clippingPointID = -1
   clippingPoint = [0.0,0.0,0.0]
   
   cell0PointIds = vtk.vtkIdList()
   cell1PointIds = vtk.vtkIdList()

   centerlines.GetCellPoints(0,cell0PointIds) #this is the cl that goes through the aneurysm
   centerlines.GetCellPoints(1,cell1PointIds)
 
   for i in range(0,min(cell0PointIds.GetNumberOfIds(),cell1PointIds.GetNumberOfIds())):
      cell0Point = centerlines.GetPoint(cell0PointIds.GetId(i))
      cell1Point = centerlines.GetPoint(cell1PointIds.GetId(i))

      distanceBetweenPoints = math.sqrt(vtk.vtkMath.Distance2BetweenPoints(cell0Point,cell1Point))
      if (distanceBetweenPoints>toll):
        divergingPointID = cell1PointIds.GetId(i)
        divergingPoint = centerlines.GetPoint(cell1PointIds.GetId(i))
        divergingPointMISR = centerlines.GetPointData().GetArray(radiusArrayName).GetTuple1(cell1PointIds.GetId(i))
        break
  
   MISphere = vtk.vtkSphere()
   MISphere.SetCenter(divergingPoint)
   MISphere.SetRadius(divergingPointMISR)
   tempPoint = [0.0,0.0,0.0]

   for i in range(divergingPointID,0,-1):
      value = MISphere.EvaluateFunction(centerlines.GetPoint(i))
      if (value>=0.0):
         tempPoint = centerlines.GetPoint(i) 
         break
       
   locator = vtk.vtkPointLocator()
   locator.SetDataSet(parentCenterlines)
   locator.BuildLocator()
 
   clippingPointID = locator.FindClosestPoint(tempPoint)
   clippingPoint = parentCenterlines.GetPoint(clippingPointID)

   return clippingPoint,divergingPoint

    def testRectilinear(self):
        rt = vtk.vtkRTAnalyticSource()
        rt.SetWholeExtent(-5, 5, -5, 5, -5, 5)
        rt.Update()
        i = rt.GetOutput()

        r = vtk.vtkRectilinearGrid()
        dims = i.GetDimensions()
        r.SetDimensions(dims)
        exts = i.GetExtent()
        orgs = i.GetOrigin()

        xs = vtk.vtkFloatArray()
        xs.SetNumberOfTuples(dims[0])
        for d in range(dims[0]):
            xs.SetTuple1(d, orgs[0] + exts[0] + d)
        r.SetXCoordinates(xs)

        ys = vtk.vtkFloatArray()
        ys.SetNumberOfTuples(dims[1])
        for d in range(dims[1]):
            ys.SetTuple1(d, orgs[1] + exts[2] + d)
        r.SetYCoordinates(ys)

        zs = vtk.vtkFloatArray()
        zs.SetNumberOfTuples(dims[2])
        for d in range(dims[2]):
            zs.SetTuple1(d, orgs[2] + exts[4] + d)
        r.SetZCoordinates(zs)

        s = vtk.vtkSphere()
        s.SetRadius(2)
        s.SetCenter(0,0,0)

        c = vtk.vtkTableBasedClipDataSet()
        c.SetInputData(r)
        c.SetClipFunction(s)
        c.SetInsideOut(1)

        c.Update()

        self.assertEqual(c.GetOutput().GetNumberOfCells(), 64)

#.........这里部分代码省略.........
        else:
            # slow, but more generic
            self._StencilGenerator = vtk.vtkImplicitFunctionToImageStencil()

        self._StencilGenerator.SetOutputOrigin(self.getImageOrigin())
        self._StencilGenerator.SetOutputSpacing(self.getImageSpacing())

        # set extent of stencil - taking into account transformation
        self._StencilGenerator.SetOutputWholeExtent(e_t)

        if is_identity:
            # use DG's fast routines
            if roi_type == "box":
                self._StencilGenerator.SetShapeToBox()
            elif roi_type == "cylinder":
                if roi_orientation == "X":
                    self._StencilGenerator.SetShapeToCylinderX()
                elif roi_orientation == "Y":
                    self._StencilGenerator.SetShapeToCylinderY()
                elif roi_orientation == "Z":
                    self._StencilGenerator.SetShapeToCylinderZ()
            elif roi_type == "ellipsoid":
                self._StencilGenerator.SetShapeToEllipsoid()
            self._StencilGenerator.SetBounds(b)
        else:
            # use JG's slow routines
            if roi_type == "box":
                obj = vtk.vtkBox()
                obj.SetTransform(t1)
                obj.SetBounds(b)
            elif roi_type == "cylinder":
                cyl = vtk.vtkCylinder()
                cyl.SetRadius(1.0)

                xc, yc, zc = (b[1] + b[0]) * 0.5, (b[3] + b[2]) * 0.5, (b[5] + b[4]) * 0.5
                diam_a, diam_b, diam_c = (b[1] - b[0]), (b[3] - b[2]), (b[5] - b[4])

                # The cylinder is infinite in extent, so needs to be cropped by using the intersection
                # of three implicit functions -- the cylinder, and two cropping
                # planes
                obj = vtk.vtkImplicitBoolean()
                obj.SetOperationTypeToIntersection()
                obj.AddFunction(cyl)

                clip1 = vtk.vtkPlane()
                clip1.SetNormal(0, 1, 0)
                obj.AddFunction(clip1)

                clip2 = vtk.vtkPlane()
                clip2.SetNormal(0, -1, 0)
                obj.AddFunction(clip2)

                t2 = vtk.vtkTransform()
                t2.Translate(xc, yc, zc)

                if roi_orientation == "X":
                    # cylinder is infinite in extent in the y-axis
                    t2.Scale(1, diam_b / 2.0, diam_c / 2.0)
                    t2.RotateZ(90)
                    r = diam_a / 2.0
                elif roi_orientation == "Y":
                    # cylinder is infinite in extent in the y-axis
                    t2.Scale(diam_a / 2.0, 1, diam_c / 2.0)
                    r = diam_b / 2.0
                elif roi_orientation == "Z":
                    # cylinder is infinite in extent in the y-axis
                    t2.Scale(diam_a / 2.0, diam_b / 2.0, 1)
                    t2.RotateX(90)
                    r = diam_c / 2.0

                clip1.SetOrigin(0, r, 0)
                clip2.SetOrigin(0, -r, 0)

                # combine transforms
                t2.SetInput(self.__Transform)

                obj.SetTransform(t2.GetInverse())

            elif roi_type == "ellipsoid":
                obj = vtk.vtkSphere()
                obj.SetRadius(1.0)

                xc, yc, zc = (b[1] + b[0]) * 0.5, (b[3] + b[2]) * 0.5, (b[5] + b[4]) * 0.5
                diam_a, diam_b, diam_c = (b[1] - b[0]), (b[3] - b[2]), (b[5] - b[4])

                t2 = vtk.vtkTransform()
                t2.Translate(xc, yc, zc)
                t2.Scale(diam_a / 2.0, diam_b / 2.0, diam_c / 2.0)

                # combine transforms
                t2.SetInput(self.__Transform)

                obj.SetTransform(t2.GetInverse())

            self._StencilGenerator.SetInput(obj)

        _t1 = time.time()
        self._StencilGenerator.Update()
        _t2 = time.time()
        return self._StencilGenerator.GetOutput()

g3Mapper.SetScalarRange(output.GetScalarRange())
g3Actor = vtk.vtkActor()
g3Actor.SetMapper(g3Mapper)
g3Actor.AddPosition(0,0,15)
gf4 = vtk.vtkDataSetSurfaceFilter()
gf4.SetInputConnection(gf2.GetOutputPort())
gf4.UseStripsOn()
g4Mapper = vtk.vtkPolyDataMapper()
g4Mapper.SetInputConnection(gf4.GetOutputPort())
g4Mapper.SetScalarRange(output.GetScalarRange())
g4Actor = vtk.vtkActor()
g4Actor.SetMapper(g4Mapper)
g4Actor.AddPosition(0,15,15)
# create pipeline - unstructured grid
#
s = vtk.vtkSphere()
s.SetCenter(output.GetCenter())
s.SetRadius(100.0)
#everything
eg = vtk.vtkExtractGeometry()
eg.SetInputData(output)
eg.SetImplicitFunction(s)
gf5 = vtk.vtkDataSetSurfaceFilter()
gf5.SetInputConnection(eg.GetOutputPort())
g5Mapper = vtk.vtkPolyDataMapper()
g5Mapper.SetInputConnection(gf5.GetOutputPort())
g5Mapper.SetScalarRange(output.GetScalarRange())
g5Actor = vtk.vtkActor()
g5Actor.SetMapper(g5Mapper)
g5Actor.AddPosition(0,0,30)
gf6 = vtk.vtkDataSetSurfaceFilter()

    def makeModels(self):
        """
        make vtk model
        """

        # Here we create two ellipsoidal implicit functions and boolean them
        # together to form a "cross" shaped implicit function.
        quadric = vtk.vtkQuadric()
        quadric.SetCoefficients(.5, 1, .2, 0, .1, 0, 0, .2, 0, 0)

        sample = vtk.vtkSampleFunction()
        sample.SetSampleDimensions(50, 50, 50)
        sample.SetImplicitFunction(quadric)
        sample.ComputeNormalsOff()

        trans = vtk.vtkTransform()
        trans.Scale(1, .5, .333)
        sphere = vtk.vtkSphere()
        sphere.SetRadius(self.radius)
        sphere.SetTransform(trans)

        trans2 = vtk.vtkTransform()
        trans2.Scale(.25, .5, 1.0)
        sphere2 = vtk.vtkSphere()
        sphere2.SetRadius(self.radius)
        sphere2.SetTransform(trans2)

        self.sphere_geom_1 = sphere
        self.sphere_geom_2 = sphere2

        union = vtk.vtkImplicitBoolean()
        union.AddFunction(sphere)
        union.AddFunction(sphere2)
        union.SetOperationType(0)

        # Here is where it gets interesting. The implicit function is used to
        # extract those cells completely inside the function. They are then
        # shrunk to helpr show what was extracted.
        extract = vtk.vtkExtractGeometry()
        extract.SetInputConnection(sample.GetOutputPort())
        extract.SetImplicitFunction(union)
        shrink = vtk.vtkShrinkFilter()
        shrink.SetInputConnection(extract.GetOutputPort())
        shrink.SetShrinkFactor(self.shrink_factor)
        dataMapper = vtk.vtkDataSetMapper()
        dataMapper.SetInputConnection(shrink.GetOutputPort())

        self.shrink_geom = shrink

        # data actor
        self.data_actor = vtk.vtkActor()
        self.data_actor.SetMapper(dataMapper)

        # The outline gives context to the original data.
        outline = vtk.vtkOutlineFilter()
        outline.SetInputConnection(sample.GetOutputPort())
        outlineMapper = vtk.vtkPolyDataMapper()
        outlineMapper.SetInputConnection(outline.GetOutputPort())

        # outline actor
        self.outline_actor = vtk.vtkActor()
        self.outline_actor.SetMapper(outlineMapper)
        outlineProp = self.outline_actor.GetProperty()
        outlineProp.SetColor(0, 0, 0)

    def _createImplicit(self, implicitType, implicitName, bounds, primaryInput):
        if implicitType in self._implicitTypes and \
               implicitName not in self._implicitsDict:

            pi = primaryInput
            rwi = self.slice3dVWR.threedFrame.threedRWI
            implicitInfoBounds = None
            
            if implicitType == "Plane":
                implicitWidget = vtk.vtkImplicitPlaneWidget()
                implicitWidget.SetPlaceFactor(1.25)
                if pi != None:
                    implicitWidget.SetInput(pi)
                    implicitWidget.PlaceWidget()
                    b = pi.GetBounds()
                    implicitWidget.SetOrigin(b[0], b[2], b[4])
                    implicitInfoBounds = b
                elif bounds != None:
                    implicitWidget.PlaceWidget(bounds)
                    implicitWidget.SetOrigin(bounds[0], bounds[2], bounds[4])
                    implicitInfoBounds = bounds
                else:
                    # this can never happen
                    pass

                implicitWidget.SetInteractor(rwi)
                implicitWidget.On()

                # create the implicit function
                implicitFunction = vtk.vtkPlane()
                # sync it to the initial widget
                self._syncPlaneFunctionToWidget(implicitWidget)
                # add it to the output
                self.outputImplicitFunction.AddFunction(implicitFunction)

                # now add an observer to the widget
                def observerImplicitPlaneWidget(widget, eventName):
                    # sync it to the initial widget
                    ret = self._syncPlaneFunctionToWidget(widget)
                    # also select the correct grid row
                    if ret != None:
                        name, ii = ret
                        row = self.findGridRowByName(name)
                        if row >= 0:
                            self._grid.SelectRow(row)

                oId = implicitWidget.AddObserver('EndInteractionEvent',
                                                 observerImplicitPlaneWidget)
                    
            elif implicitType == "Sphere":
                implicitWidget = vtk.vtkSphereWidget()
                implicitWidget.SetPlaceFactor(1.25)
                implicitWidget.TranslationOn()
                implicitWidget.ScaleOn()
                #implicitWidget.HandleVisibilityOn()
                
                if pi != None:
                    implicitWidget.SetInput(pi)
                    implicitWidget.PlaceWidget()
                    b = pi.GetBounds()
                    implicitInfoBounds = b
                    #implicitWidget.SetOrigin(b[0], b[2], b[4])
                elif bounds != None:
                    implicitWidget.PlaceWidget(bounds)
                    implicitInfoBounds = bounds
                    #implicitWidget.SetOrigin(bounds[0], bounds[2], bounds[4])
                else:
                    # this can never happen
                    pass

                implicitWidget.SetInteractor(rwi)
                implicitWidget.On()

                # create the implicit function
                implicitFunction = vtk.vtkSphere()
                # sync it to the initial widget
                self._syncSphereFunctionToWidget(implicitWidget)
                # add it to the output
                self.outputImplicitFunction.AddFunction(implicitFunction)

                # now add an observer to the widget
                def observerImplicitSphereWidget(widget, eventName):
                    # sync it to the initial widget
                    ret = self._syncSphereFunctionToWidget(widget)
                    # also select the correct grid row
                    if ret != None:
                        name, ii = ret
                        row = self.findGridRowByName(name)
                        if row >= 0:
                            self._grid.SelectRow(row)

                oId = implicitWidget.AddObserver('EndInteractionEvent',
                                                 observerImplicitSphereWidget)



            if implicitWidget:
                # set the priority so it gets interaction before the
                # ImagePlaneWidget.  3D widgets have default priority 0.5,
                # so we assign our widgets just a tad higher. (voiwidget
#.........这里部分代码省略.........

ren1 = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# create implicit function primitives
cone = vtk.vtkCone()
cone.SetAngle(20)
vertPlane = vtk.vtkPlane()
vertPlane.SetOrigin(.1,0,0)
vertPlane.SetNormal(-1,0,0)
basePlane = vtk.vtkPlane()
basePlane.SetOrigin(1.2,0,0)
basePlane.SetNormal(1,0,0)
iceCream = vtk.vtkSphere()
iceCream.SetCenter(1.333,0,0)
iceCream.SetRadius(0.5)
bite = vtk.vtkSphere()
bite.SetCenter(1.5,0,0.5)
bite.SetRadius(0.25)
# combine primitives to build ice-cream cone
theCone = vtk.vtkImplicitBoolean()
theCone.SetOperationTypeToIntersection()
theCone.AddFunction(cone)
theCone.AddFunction(vertPlane)
theCone.AddFunction(basePlane)
theCream = vtk.vtkImplicitBoolean()
theCream.SetOperationTypeToDifference()
theCream.AddFunction(iceCream)
theCream.AddFunction(bite)

sTube = positions[0]
sTubeX = sTube[:,0]
sTubeY = sTube[:,1]
sTubeZ = sTube[:,2]

tPairs = []

for i in range(1,np.shape(positions[0])[0]):
    tPairs.append((i-1,i))


spheres = [None]*np.size(sTubeX)
sphereActs = [None]*np.size(sTubeX)
for i in range(len(spheres)):
    spheres[i] = vtk.vtkSphere()
    spheres[i].SetCenter(sTubeX[i], sTubeY[i], sTubeZ[i])
    spheres[i].SetRadius(scalars[0])


cyls = [None]*len(tPairs)
for c in range(len(cyls)):
    icyl = vtk.vtkCylinder()
    p1 = vtk.vtkPlane()
    p1.SetOrigin(sTubeX[tPairs[c][0]],sTubeY[tPairs[c][0]],sTubeZ[tPairs[c][0]])
    p1.SetNormal(sTubeX[tPairs[c][1]] - sTubeX[tPairs[c][0]],sTubeY[tPairs[c][1]] - sTubeY[tPairs[c][0]],sTubeZ[tPairs[c][1]] - sTubeZ[tPairs[c][0]])
    p2 = vtk.vtkPlane()
    p2.SetOrigin(sTubeX[tPairs[c][1]],sTubeY[tPairs[c][1]],sTubeZ[tPairs[c][1]])
    p2.SetNormal(sTubeX[tPairs[c][0]] - sTubeX[tPairs[c][1]],sTubeY[tPairs[c][0]] - sTubeY[tPairs[c][1]],sTubeZ[tPairs[c][0]] - sTubeZ[tPairs[c][1]])
    ccyl = vtk.vtkImplicitBoolean()
    ccyl.SetOperationTypeToIntersection()

def gen_slice(args):

    verbose = args.verbose
    origin = args.origin
    normal = args.normal
    output_file = args.output_file[0]
    input_file =  args.input_file[0]
    variable =  args.variable[0]
    on_sphere = args.sphere
    sphere_radius = args.sphere_radius
    depth = args.sphere_radius

    xi = []
    yi = []
    zi = []
    vtu_files = []
    path = os.path.dirname(input_file)
    
    if (not on_sphere):
        if (normal == (0,0,1)):
            index_1 = 0
            index_2 = 1
            x_axis = "X"
            y_axis = "Y"
        elif(normal == (0,1,0)):
            index_1 = 0
            index_2 = 2
            x_axis = "X"
            y_axis = "Depth"
        elif(normal == (1,0,0)):
            index_1 = 1
            index_2 = 2
            x_axis = "Y"
            y_axis = "Depth"
        else:
            print "Normal should be one of (1,0,0), (0,1,0), or (0,0,1)"
            sys.exit(-1)
        # check origin
        if (origin == None):
            print "Origin should be a set of three corrdinate, eg. (0,0,-10)"
            sys.exit(-1)
    else:
        import fluidity.spheretools as st
        index_1 = 0
        index_2 = 1
        index_3 = 2
        x_axis = "Longitude"
        y_axis = "Latitude"


    # get vtus from pvtu file. They are in lines such as:
    # <Piece Source="restratA-np64-adapt-C_99/restratA-np64-adapt-C_99_0.vtu" />
    # As PVTU is XML, parse as such and grab the Piece elements
    xml_root = etree.parse(input_file,xml_parser)
    find = etree.XPath("//Piece")
    peices = find(xml_root)
    for p in peices:
        name = p.attrib['Source']
        vtu_files.append(os.path.join(path,name))
      
    if (on_sphere):
        sphere = vtk.vtkSphere()
        sphere.SetCenter(0,0,0)
        sphere.SetRadius(sphere_radius-10)
    else:
        plane = vtk.vtkPlane()
        plane.SetOrigin(origin[0], origin[1], origin[2])
        plane.SetNormal(normal[0],normal[1],normal[2])


    getMagnitude = False
    component = -1
    # Special cases of variables
    if (variable == "VelocityMagnitude"):
        getMagnitude = True
        variable = "Velocity_projection"
    if (variable == "VelocityU"):
        getMagnitude = False
        variable = "ProjectedVelocity"
        component = 0
    if (variable == "VelocityV"):
        getMagnitude = False
        variable = "ProjectedVelocity"
        component = 1   
    if (variable == "BedShearStressMagnitude"):
        getMagnitude = True
        variable = "BedShearStress_projection"
    if (variable == "MaxBedShearStressMagnitude"):
        getMagnitude = True
        variable = "MaxBedShearStress_projection"
    if (variable == "AveBedShearStressMagnitude"):
        getMagnitude = True
        variable = "AveBedShearStress_projection"
    if (variable == "MaxVelocityMagnitude"):
        getMagnitude = True
        variable = "MaxVelocity_projection"
    if (variable == "AveVelocityMagnitude"):
        getMagnitude = True
        variable = "AveVelocity_projection"
    if (variable == "AveVelocityU"):
#.........这里部分代码省略.........