Re-close the surface of combined meshs

Hi guys,

I’m trying to do a clip on any closed mesh using a sphere, and get a closed mesh after clip.
And here is my try:

  1. clip the mesh using vtk.vtkClipPolyData() with a sphere
  2. get the sphere polys inside the to-be-cut mesh using vtk.vtkSelectEnclosedPoints() and vtk.vtkMultiThreshold()
  3. Combine the two polydata using vtk.vtkAppendPolyData()

But after that, I get a unclosed polygon data (as figure shows, cut a sphere source using a sphere clip).

PS: I have tried vtk.vtkBooleanOperationPolyDataFilter() to do this job, but it’s too slow…
Anyone could give a clue? How can I do the sphere clip job? Or how can I make the polydata closure?

Here is my code:

import vtk 

def GetSpherePD( radius, center ):
    """
    :return: The PolyData representation of a sphere.
    """
    reduction=0.5
    source = vtk.vtkSphereSource()
    source.SetRadius(radius)
    source.SetThetaResolution(100)
    source.SetPhiResolution(100)
    source.SetCenter(center)
    source.Update()
    return source.GetOutput()

center=[0,0,0]
polyData = GetSpherePD(0.5,center)
cut_center=[0.5,0.,0.]
cut_radius=0.5
cutpolyData = GetSpherePD(cut_radius,cut_center)

boundaryMapper = vtk.vtkPolyDataMapper()
ori_transform=vtk.vtkTransform()
OriTransformFileter=vtk.vtkTransformPolyDataFilter()
OriTransformFileter.SetTransform(ori_transform)
OriTransformFileter.SetInputData(polyData)
OriTransformFileter.Update()

sphere= vtk.vtkSphere()
sphere.SetRadius(cut_radius)
sphere.SetCenter(cut_center)

clipper = vtk.vtkClipPolyData()
clipper.SetInputData(OriTransformFileter.GetOutput())
clipper.SetClipFunction(sphere)
clipper.SetValue(0.)
clipper.Update()

select = vtk.vtkSelectEnclosedPoints()
select.SetInputData(cutpolyData)
select.SetSurfaceData(polyData)
select.Update()
threshold = vtk.vtkMultiThreshold()
# Outside points have a 0 value in ALL points of a cell
outsideId = threshold.AddBandpassIntervalSet(
    0, 0,
    vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "SelectedPoints",
    0, 1)
# Inside points have a 1 value in ALL points of a cell
insideId = threshold.AddBandpassIntervalSet(
    1, 1,
    vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "SelectedPoints",
    0, 1)
# Border points have a 0 or a 1 in at least one point of a cell
borderId = threshold.AddIntervalSet(
    0, 1,
    vtk.vtkMultiThreshold.OPEN, vtk.vtkMultiThreshold.OPEN,
    vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "SelectedPoints",
    0, 0)
threshold.SetInputData(select.GetOutput())
threshold.OutputSet(outsideId)
threshold.OutputSet(insideId)
threshold.OutputSet(borderId)
threshold.Update()

surfaceFilter=vtk.vtkDataSetSurfaceFilter()
surfaceFilter.SetInputData(threshold.GetOutput().GetBlock(insideId).GetBlock(0))
# surfaceFilter.SetInputConnection(threshold.GetOutputPort())
surfaceFilter.Update()

normFilter= vtk.vtkPolyDataNormals()
normFilter.SetInputData(surfaceFilter.GetOutput())
normFilter.SetComputeCellNormals(1)
normFilter.SetComputePointNormals(1)
normFilter.SetFlipNormals(1)
normFilter.Update()

combined_polydata = vtk.vtkAppendPolyData()
combined_polydata.AddInputData(normFilter.GetOutput())
combined_polydata.AddInputData(clipper.GetOutput())
combined_polydata.Update()

combined_mapper=vtk.vtkDataSetMapper()
combined_mapper.SetInputData(combined_polydata.GetOutput())
combinedActor=vtk.vtkActor()
combinedActor.SetMapper(combined_mapper)

renderer = vtk.vtkRenderer()
renderer.AddActor(combinedActor)

renderWindow = vtk.vtkRenderWindow()
renderWindow.AddRenderer(renderer)
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(renderWindow)

interactor.Start()```