def splitDomainBetweenEndoAndEpi(
        pdata_domain,
        r=0.99,
        verbose=1):

    myVTK.myPrint(verbose, "*** splitDomainBetweenEndoAndEpi ***")

    bounds = pdata_domain.GetBounds()
    assert (r > 0.)
    assert (r < 1.)
    origin = [(1./2)*bounds[0]+(1./2)*bounds[1],
              (1./2)*bounds[2]+(1./2)*bounds[3],
              (1.-r)*bounds[4]+(  r )*bounds[5]]
    #myVTK.myPrint(verbose-1, "bounds = "+str(bounds))
    #myVTK.myPrint(verbose-1, "origin = "+str(origin))

    (pdata_domain,
     cap) = myVTK.clipPDataUsingPlane(
         pdata_mesh=pdata_domain,
         plane_O=origin,
         plane_N=[0,0,1],
         verbose=verbose-1)

    connectivity0 = vtk.vtkPolyDataConnectivityFilter()
    connectivity0.SetExtractionModeToSpecifiedRegions()
    connectivity0.AddSpecifiedRegion(0)
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        connectivity0.SetInputData(pdata_domain)
    else:
        connectivity0.SetInput(pdata_domain)
    connectivity0.Update()
    pdata0 = connectivity0.GetOutput()
    assert (pdata0.GetNumberOfPoints())

    connectivity1 = vtk.vtkPolyDataConnectivityFilter()
    connectivity1.SetExtractionModeToSpecifiedRegions()
    connectivity1.AddSpecifiedRegion(1)
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        connectivity1.SetInputData(pdata_domain)
    else:
        connectivity1.SetInput(pdata_domain)
    connectivity1.Update()
    pdata1 = connectivity1.GetOutput()
    assert (pdata1.GetNumberOfPoints())

    if (myVTK.getPDataSurfaceArea(pdata0,0) < myVTK.getPDataSurfaceArea(pdata1,0)):
        return pdata0, pdata1
    else:
        return pdata1, pdata0

    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 cleanMesh(mesh, connectivityFilter=False):
    try:
        t = time.clock()
        connect = vtk.vtkPolyDataConnectivityFilter()
        clean = vtk.vtkCleanPolyData()

        if (connectivityFilter):
            if vtk.vtkVersion.GetVTKMajorVersion() >= 6:
                connect.SetInputData( mesh )
            else:
                connect.SetInput( mesh )
            connect.SetExtractionModeToLargestRegion()
            clean.SetInputConnection( connect.GetOutputPort() )
        else:
            if vtk.vtkVersion.GetVTKMajorVersion() >= 6:
                clean.SetInputData( mesh )
            else:
                clean.SetInput( mesh )

        clean.Update()
        print ("Surface cleaned")
        m2 = clean.GetOutput()
        print "    ", m2.GetNumberOfPolys(), "polygons"
        elapsedTime(t)
        clean = None
        connect = None
        return m2
    except:
        print "Surface cleaning 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 extractlargestregion(polydata):
    """Extract largest of several disconnected regions."""
    connect = vtk.vtkPolyDataConnectivityFilter()
    connect.SetInput(polydata)
    connect.SetExtractionModeToLargestRegion()
    connect.Update()
    return connect.GetOutput()

 def __init__(self):
     """Setup the pipeline.
     """
     self.reader = vtkSTLReader()
     
     self.cutter = vtkCutter()
     self.cutter.SetInputConnection(self.reader.GetOutputPort())
     
     self.regionPicker = vtkPolyDataConnectivityFilter()
     self.regionPicker.SetInputConnection(self.cutter.GetOutputPort())
     self.regionPicker.SetExtractionModeToClosestPointRegion()
     
     self.scaler = vtkTransformPolyDataFilter()
     self.scaler.SetInputConnection(self.regionPicker.GetOutputPort())
     
     self.GetFileName = self.reader.GetFileName
     self.SetFileName = self.reader.SetFileName
     
     self.GetOutputPort = self.scaler.GetOutputPort
     self.GetOutput = self.scaler.GetOutput
     self.Update = self.scaler.Update
     
     self.iolet = None
     self.fileUnitLength = None
     
     return

    def extractSurfaces(self):
        self.PrintLog("Extracting surfaces.")

        def bnorm(data):
            bounds = data.GetBounds()
            return math.sqrt(sum((bounds[2*i+1]-bounds[2*i])**2 for i in range(3)))

        # Get bounds of entire surface
        bounds_norm = bnorm(self.Surface)

        # Currently assuming that this is the region numbers that were produced by coloring!
        # TODO: If necessary, get from array to be more robust.
        region_ids = (0, 1)

        # Extract each surface in turn to find the smallest one
        subsurfaces = {}
        for k in region_ids:
            connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
            connectivityFilter.SetInputData(self.ColoredSurface)
            connectivityFilter.SetExtractionModeToSpecifiedRegions()
            connectivityFilter.AddSpecifiedRegion(k)
            connectivityFilter.ColorRegionsOff()
            connectivityFilter.SetScalarConnectivity(0)
            connectivityFilter.Update()
            subsurfaces[k] = connectivityFilter.GetOutput()

        # The inner surface has smaller bounds
        if bnorm(subsurfaces[region_ids[0]]) < bnorm(subsurfaces[region_ids[1]]):
            self.InnerRegionId = region_ids[0]
            self.OuterRegionId = region_ids[1]
        else:
            self.InnerRegionId = region_ids[1]
            self.OuterRegionId = region_ids[0]
        self.InnerSurface = subsurfaces[self.InnerRegionId]
        self.OuterSurface = subsurfaces[self.OuterRegionId]

    def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: No input surface.')

        if (self.GroupId != -1) & (self.GroupIdsArrayName!=''):
            self.Surface.GetPointData().SetActiveScalars(self.GroupIdsArrayName)

        connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
        connectivityFilter.SetInput(self.Surface)
        connectivityFilter.ColorRegionsOff()       
        connectivityFilter.SetExtractionModeToLargestRegion()
        if self.GroupId != -1:
            connectivityFilter.ScalarConnectivityOn()
            scalarRange = [self.GroupId,self .GroupId]
            connectivityFilter.SetScalarRange(scalarRange)
        connectivityFilter.Update()

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

            self.Surface = cleaner.GetOutput()

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

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

def test_connectivity(pd):
    conn = vtk.vtkPolyDataConnectivityFilter()
    conn.SetInput(pd)
    conn.SetExtractionModeToAllRegions()
    conn.Update()

    return conn.GetNumberOfExtractedRegions()

    def Execute(self):

        if self.Image == None:
            self.PrintError('Error: No Image.')

        extent = self.Image.GetExtent()
        translateExtent = vtk.vtkImageTranslateExtent()
        translateExtent.SetInputData(self.Image)
        translateExtent.SetTranslation(-extent[0],-extent[2],-extent[4])
        translateExtent.Update()

        if (self.ArrayName != ''):
            translateExtent.GetOutput().GetPointData().SetActiveScalars(self.ArrayName)

        marchingCubes = vtk.vtkMarchingCubes()
        marchingCubes.SetInputConnection(translateExtent.GetOutputPort())
        marchingCubes.SetValue(0,self.Level)
        marchingCubes.Update()

        self.Surface = marchingCubes.GetOutput()

        if self.Connectivity == 1:
            connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
            connectivityFilter.SetInputData(self.Surface)
            connectivityFilter.SetExtractionModeToLargestRegion()
            connectivityFilter.Update()
            self.Surface = connectivityFilter.GetOutput()

def labelregions(polydata):
    """Label connected regions by assigning a 'RegionId' as pointdata."""
    connect = vtk.vtkPolyDataConnectivityFilter()
    connect.SetInput(polydata)
    connect.SetExtractionModeToAllRegions()
    connect.ColorRegionsOn()
    connect.Update()
    return connect.GetOutput()

def extractclosestpointregion(polydata, point=[0, 0, 0]):
    """Extract region closest to specified point."""
    connect = vtk.vtkPolyDataConnectivityFilter()
    connect.SetInput(polydata)
    connect.SetExtractionModeToClosestPointRegion()
    connect.SetClosestPoint(point)
    connect.Update()
    return connect.GetOutput()

def computeComponents(vertices, faces):
    pts, tris = arrayToPolydata(vertices, faces)
    polys = vtk.vtkPolyData()
    polys.SetPoints(pts)
    polys.SetPolys(tris)
    polys.Update()
    
    cleanFilter = vtk.vtkCleanPolyData()
    #cleanFilter.PointMergingOff()
    cleanFilter.ConvertStripsToPolysOff()
    cleanFilter.ConvertPolysToLinesOff()
    cleanFilter.ConvertLinesToPointsOff()
    cleanFilter.SetInput(polys)
    cleanFilter.Update()
    
    connFilter = vtk.vtkPolyDataConnectivityFilter()
    connFilter.ScalarConnectivityOff()
    connFilter.SetExtractionModeToAllRegions()
    connFilter.SetInput(cleanFilter.GetOutput())
    connFilter.Update()
    
    nregions = connFilter.GetNumberOfExtractedRegions()
    print "Regions extracted: %d" % (nregions)
    
    components = []
    for i in range(nregions):
        connFilter = vtk.vtkPolyDataConnectivityFilter()
        connFilter.ScalarConnectivityOff()
        connFilter.SetExtractionModeToSpecifiedRegions()
        connFilter.AddSpecifiedRegion(i)
        connFilter.SetInput(cleanFilter.GetOutput())
        connFilter.Update()
    
        cFilter = vtk.vtkCleanPolyData()
#        cFilter.ConvertStripsToPolysOff()
#        cFilter.ConvertPolysToLinesOff()
#        cFilter.ConvertLinesToPointsOff()
        cFilter.SetInput(connFilter.GetOutput())
        cFilter.Update()
    
        v,f = polydataToArray(cFilter.GetOutput())
        components.append([v,f])
        
    return components

def connectivityFilter(polyData):
    connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
    connectivityFilter.SetInputConnection(polyData.GetOutputPort())
    connectivityFilter.SetExtractionModeToAllRegions()
    connectivityFilter.ColorRegionsOn()
    connectivityFilter.Update()
    global n_regions
    n_regions = connectivityFilter.GetNumberOfExtractedRegions()
    print ("number of regions extracted: " + str(n_regions) + "\n")
    return connectivityFilter

def concomp(opts, argv):
	p = nv.readVTK(argv[0])
	c = vtk.vtkPolyDataConnectivityFilter()
	c.SetInput(p)
	c.SetExtractionModeToLargestRegion()
	c.Update()
	d = vtk.vtkCleanPolyData()
	d.SetInput(c.GetOutput())
	d.Update()
	p = d.GetOutput()
	nv.writeVTK(argv[1],p)

  def applyFilters(self, state):

    surface = None
    surface = state.inputModelNode.GetPolyDataConnection()

    if state.decimation:
      triangle = vtk.vtkTriangleFilter()
      triangle.SetInputConnection(surface)
      decimation = vtk.vtkDecimatePro()
      decimation.SetTargetReduction(state.reduction)
      decimation.SetBoundaryVertexDeletion(state.boundaryDeletion)
      decimation.PreserveTopologyOn()
      decimation.SetInputConnection(triangle.GetOutputPort())
      surface = decimation.GetOutputPort()

    if state.smoothing:
      if state.smoothingMethod == "Laplace":
        smoothing = vtk.vtkSmoothPolyDataFilter()
        smoothing.SetBoundarySmoothing(state.boundarySmoothing)
        smoothing.SetNumberOfIterations(state.laplaceIterations)
        smoothing.SetRelaxationFactor(state.laplaceRelaxation)
        smoothing.SetInputConnection(surface)
        surface = smoothing.GetOutputPort()
      elif state.smoothingMethod == "Taubin":
        smoothing = vtk.vtkWindowedSincPolyDataFilter()
        smoothing.SetBoundarySmoothing(state.boundarySmoothing)
        smoothing.SetNumberOfIterations(state.taubinIterations)
        smoothing.SetPassBand(state.taubinPassBand)
        smoothing.SetInputConnection(surface)
        surface = smoothing.GetOutputPort()

    if state.normals:
      normals = vtk.vtkPolyDataNormals()
      normals.AutoOrientNormalsOn()
      normals.SetFlipNormals(state.flipNormals)
      normals.SetSplitting(state.splitting)
      normals.SetFeatureAngle(state.featureAngle)
      normals.ConsistencyOn()
      normals.SetInputConnection(surface)
      surface = normals.GetOutputPort()

    if state.cleaner:
      cleaner = vtk.vtkCleanPolyData()
      cleaner.SetInputConnection(surface)
      surface = cleaner.GetOutputPort()

    if state.connectivity:
      connectivity = vtk.vtkPolyDataConnectivityFilter()
      connectivity.SetExtractionModeToLargestRegion()
      connectivity.SetInputConnection(surface)
      surface = connectivity.GetOutputPort()

    state.outputModelNode.SetPolyDataConnection(surface)
    return True

    def colorSurfaceRegions(self):
        self.PrintLog("Coloring surface regions.")

        connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
        connectivityFilter.SetInputData(self.DoubleSurface)
        connectivityFilter.ColorRegionsOn()
        connectivityFilter.SetExtractionModeToAllRegions()
        connectivityFilter.Update()

        assert connectivityFilter.GetNumberOfExtractedRegions() == 2

        self.ColoredSurface = connectivityFilter.GetOutput()

    def PickCallback(self, obj):
        picker = vtk.vtkPointPicker()
        eventPosition = self.vmtkRenderer.RenderWindowInteractor.GetEventPosition()
        result = picker.Pick(float(eventPosition[0]),float(eventPosition[1]),0.0,self.vmtkRenderer.Renderer)
        # We treat result as a boolean value (which can also be nonetype if a selection outside the render windw is made)
        # If the value of result is False or None, then an invalid actor was selected in the scene. 
        if result is None:
            return
        elif result == 0:
            return

        # find the value of "RegionId" stored in the vtkPolyData PointData array. This is a unique number
        # which is generated for each connected surface point by the allRegionConnectivity 
        # vtkPolyDataConnectivityFilter function
        pickedPointId = picker.GetPointId()
        self.CurrentPickedActor = picker.GetActor()
        self.CurrentPickedPolyData = self.CurrentPickedActor.GetMapper().GetInputAsDataSet()
        regionId = self.CurrentPickedPolyData.GetPointData().GetArray('RegionId').GetValue(pickedPointId)
        
        # extract only the desired region Id from the dataset
        isolatedFilter = vtk.vtkPolyDataConnectivityFilter()
        isolatedFilter.SetInputData(self.CurrentPickedPolyData)
        isolatedFilter.SetExtractionModeToSpecifiedRegions()
        isolatedFilter.AddSpecifiedRegion(regionId)

        # remove all points not in the extracted cells
        cleanedFilter = vtk.vtkCleanPolyData()
        cleanedFilter.SetInputConnection(isolatedFilter.GetOutputPort())
        cleanedFilter.Update()

        outputSurface = cleanedFilter.GetOutput()

        # create a new mapper and actor for this extracted surface object
        mapper = vtk.vtkPolyDataMapper()
        mapper.SetInputData(outputSurface)
        mapper.ScalarVisibilityOff()
        mapper.Update()

        actor = vtk.vtkActor()
        actor.SetMapper(mapper)
        actor.GetProperty().SetColor(1.0, 0.0, 0.0) # red color
        actor.GetProperty().SetDiffuse(1.0)
        actor.GetProperty().SetSpecular(0.2)
        actor.GetProperty().SetRepresentationToSurface()
        actor.GetProperty().EdgeVisibilityOff()

        # place the extracted surface object's actor instance in the necessary render windows
        # and lists to log it's selection        
        self.PastPickedActorList.append(actor)
        self._OutputPolyDataList.append(outputSurface)
        self.vmtkRenderer.Renderer.AddActor(actor)
        self.vmtkRenderer.RenderWindow.Render()
        return

 def GetTargetSurface(self):
     """Get the connected part of the clipped surface that lies
     closest to the first clipping plane's origin, using
     vtkPolyDataConnectivityFilter.
     
     """
     filter = vtk.vtkPolyDataConnectivityFilter()
     filter.SetExtractionModeToClosestPointRegion()
     filter.SetClosestPoint(self.Planes[0].GetOrigin())
     filter.SetInput(self.Clipped)
     filter.Update()
     self.Output = filter.GetOutput()
     return

    def Execute(self):

        if self.Surface == None:
            self.PrintError('Error: No input surface.')

        if (self.GroupId != -1) and (self.GroupIdsArrayName!=''):
            self.Surface.GetPointData().SetActiveScalars(self.GroupIdsArrayName)

        barycenter = [0.0,0.0,0.0]
        if self.Method == 'closest' and self.ClosestPoint == None:
            n = self.ReferenceSurface.GetNumberOfPoints()
            for i in range(n):
                point = self.ReferenceSurface.GetPoint(i)
                barycenter[0] += point[0]
                barycenter[1] += point[1]
                barycenter[2] += point[2]
            barycenter[0] /= n
            barycenter[1] /= n
            barycenter[2] /= n

        connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
        connectivityFilter.SetInputData(self.Surface)
        connectivityFilter.ColorRegionsOff()       
        if self.Method == 'largest':
            connectivityFilter.SetExtractionModeToLargestRegion()
        elif self.Method == 'closest':
            connectivityFilter.SetExtractionModeToClosestPointRegion()
            if self.ClosestPoint:
                connectivityFilter.SetClosestPoint(self.ClosestPoint)
            else:
                connectivityFilter.SetClosestPoint(barycenter)
        elif self.Method == 'all':
            connectivityFilter.SetExtractionModeToAllRegions()
            connectivityFilter.ColorRegionsOn()

        if self.GroupId != -1:
            connectivityFilter.ScalarConnectivityOn()
            scalarRange = [self.GroupId,self.GroupId]
            connectivityFilter.SetScalarRange(scalarRange)
        connectivityFilter.Update()

        self.Surface = connectivityFilter.GetOutput()

        if self.CleanOutput == 1:
            cleaner = vtk.vtkCleanPolyData()
            cleaner.SetInputConnection(connectivityFilter.GetOutputPort())
            cleaner.Update()

            self.Surface = cleaner.GetOutput()