Accurately cutting loops out of 3D PolyData

I am trying to accurately cut loops out of a 3D PolyData with a loop which lies on the surface (or very close to it). Similar topics have been discussed a number of times before:

  • Punching a polygonal hole into a triangulated surface demonstrates how to do this for 2D using vtkImprintFilter
  • Splitting cells on a PolyData Extrudes the loop into a skirt and uses vtkIntersectionPolyDataFilter scribe the loop onto the surface. However, with that approach I have not been able to “cut out” the surface patch which is inside the loop - I don’t know how to select only the points/cells which are inside. Also, I have encountered simple cases where the intersection fails entirely.

Is there any updated guidance on this? It feels like the new vtkImprintFilter has the capacity to do this, but I have been unsuccessful getting it to work on 3D geometries. CC @will.schroeder who implemented this filter.

An example is shown below, using the surface.vtp and boundaryline.vtp files provided by @cobo on this post. This example produces the correct result when the target mesh and loop are co-planar, but not generally.

#!/usr/bin/env python
import vtk
import pyvista as pv
import numpy as np

target = pv.read("surface.vtp")
loop = pv.read("boundaryline.vtp")

# # make the data planar - the imprint works fine if these lines are un-commented!
# target.points[:, 2] = 0.
# loop.points[:, 2] = 0.

# ray-cast the loop points onto the target
origins = loop.points
direction = np.array([0, 0, 1])
length = 100.
intersections = []
for p in loop.points:
    # uses `vtkOBBTree.IntersectWithLine`
    intersection, _ = target.extract_surface().ray_trace(p + direction * length,
                                                         p - direction * length,
                                                         first_point=True)
    intersections.append(intersection)
intersections = np.array(intersections)

# update point data
loop_on_target = loop.copy(deep=True)
loop_on_target.points = intersections

# extract contour loop from loop polydata
loopExtract = vtk.vtkContourLoopExtraction()
loopExtract.SetInputDataObject(loop_on_target)
loopExtract.Update()

# clean
cleanLoop = vtk.vtkCleanPolyData()
cleanLoop.SetInputConnection(loopExtract.GetOutputPort())

# imprint
imp = vtk.vtkImprintFilter()
imp.SetTargetData(target)
imp.SetImprintConnection(cleanLoop.GetOutputPort())
imp.SetTolerance(0.01)
imp.BoundaryEdgeInsertionOn()
imp.TriangulateOutputOn()
imp.SetOutputTypeToImprintedRegion()
imp.Update()

# get result of filter pipeline
res = pv.filters._get_output(imp)

# plot
p = pv.Plotter()
p.add_mesh(target, style='wireframe')
p.add_mesh(res, color='green')
p.add_mesh(res.copy(), color='w', style='wireframe')
p.add_mesh(loop_on_target, color='red', style='wireframe', line_width=5)
p.view_xy()
p.show()

Correct result when target mesh and loop are co-planar. The green shaded region is the imprinted region.

Wrong result when the target mesh is 3D and the loop is projected onto the mesh.