# Splitting cells on a PolyData

I’d like to “split” cells on a PolyData surface, adding vertices and cells. For example, in the screenshot below, I am starting with the black mesh, and I’d like to add all the points where the pink lines intersect the mesh, and adjust the cells connectivity appropriately.

Is there an existing filter to do something like this, or a method that anybody could suggest?

If you just want to add know the point ids inside the square, you could use a `vtkAbstractPointLocator`. You can use any concrete subclass, but if you are in 2D, I’d recommend using `vtkStaticPointLocator2D`.

You can use the method `vtkAbstractPointLocator::FindPointsWithinRadius`, querying with the center of your rectangle with a radius equal to half its diagonal. This method outputs a list of point ids that are within the input radius. Then you can filter them out by checking if they fit inside the rectangle (threshold x and y coordinates separately).

I’m not super familiar with cell locators, but they exist, and you can take a look if you want to extract the cells as well. Look for `vtkAbstractCellLocator`.

The short answer is no: there is no such filter in VTK. So you will have to do the “hard work” manually by creating a new mesh (polydata), adding all the existing points and cells that are not intersected by the square (using a cell locator), then perform the intersection, adding the intersection points and the connectivity of the new cells.
VTK (in particular the vtkPolyData data structure) is clearly not the best tool to perform mesh manipulation but this is doable with a bit of work.

Thanks @Yohann_Bearzi and @Joachim_P!

I am interested in the general 3D case; the original figure I provided was not very clear.

More formally, the problem I need to solve is this - given a PolyData and a closed polyline which lies on that surface, build a new PolyData containing:

• All of the cells which are fully “outside” the polyline
• All of the cells which are fully “inside” the polyline
• New cells corresponding to those which are intersected by the polyline.

I am trying to avoid doing as much of the hard work as possible . I could be way off, but I feel like a lot of the machinery already exists in either vtkClipPolyData or vtkBooleanOperationPolyDataFilter.

I am also open to suggestions if if you are aware of other libraries which might be able to accomplish this.

My first idea was to use `vtkCookieCutter` but it doesn’t to work as expected - the other possibility is to cut using a loop of points:

``````import vtk
from vedo import *

mesh = Mesh('https://www.dropbox.com/s/bu9otn8mcbtbtll/dolfin_fine.vtk').lw(0.1)

pts = [
[ 8.53826183e-01,  1.90910461e-01, 0],
[ 8.55859667e-01,  8.72127538e-01, 0],
[ 7.50018833e-02,  8.68060570e-01, 0],
[ 1.01437173e-01,  6.07674966e-02, 0],
]

# ..this doesn't quite work:
# ccut.SetInputData(mesh.polydata())
# ccut.Update()
# cmesh = Mesh(ccut).lw(0.1).c('tomato')

# ..this uses vtkSelectPolyData with SetLoop() + vtkClipPolyData

``````

the `vtkCookieCutter` filter (commented out part above) gives me:

maybe I’m not using it properly…

In 3D the point loop method should work, see an example here (search “loop”).

…sorry I just realized I skipped your last message… maybe this fits better what you are looking for:

``````import vtk
from vedo import *

mesh = Mesh('https://www.dropbox.com/s/bu9otn8mcbtbtll/dolfin_fine.vtk').lw(0.1)

pts = [
[ 8.53826183e-01,  1.90910461e-01, 0],
[ 8.55859667e-01,  8.72127538e-01, 0],
[ 7.50018833e-02,  8.68060570e-01, 0],
[ 1.01437173e-01,  6.07674966e-02, 0],
]

line = Line(pts, closed=True).lw(3).c('red')

invert=1
onBoundary=0

ippd = vtk.vtkImplicitSelectionLoop()
ippd.SetLoop(line.polydata().GetPoints())
ippd.AutomaticNormalGenerationOn()

clipper = vtk.vtkExtractPolyDataGeometry()
clipper.SetInputData(mesh.polydata())
clipper.SetImplicitFunction(ippd)
clipper.SetExtractInside(invert)
clipper.SetExtractBoundaryCells(onBoundary)
clipper.Update()
cmesh = Mesh(clipper).lw(1).c('t')

show([[mesh, line], [cmesh, line]], N=2, axes=1)

``````

Thanks for the responses, and for introducing me to some new VTK filters and vedo!

I am a bit hesitant to use clipping with an implicit function due to the discussion in https://www.vtkjournal.org/browse/publication/797.

I did find a solution using vtkIntersectionPolyDataFilter with the following approach:

1. Create a polyline from the points
2. Extrude the polyline to form a skirt which intersects the surface (if the points originally lie exactly on the surface, then you only need to extrude by some small epsilon)
3. Intersect the original surface with vtkIntersectionPolyDataFilter. This returns the intersection polyline, and the two surfaces which are remeshed to form triangular cells from the split cells.

The example below (uses pyvista) adds the edges of a trapezoid projected onto a sphere.

``````import pyvista as pv
import vtk

surf = pv.Sphere(center=(0, 0, 0), radius=200)

polyline = pv.PolyLine([(10, 0, 0),
(10, 50, 0),
(10, 25, 100),
(10, 0, 100)], closed=True)
skirt = polyline.extrude((400, 0, 0)).triangulate().clean().extract_surface()

f = vtk.vtkIntersectionPolyDataFilter()
f.SetInputDataObject(0, surf)
f.SetInputDataObject(1, skirt)
f.SetComputeIntersectionPointArray(True)
f.Update()

intersection_polyline = pv.wrap(f.GetOutputDataObject(0))  # intersection of the two surfaces
surf_remeshed = pv.wrap(f.GetOutputDataObject(1))  # the surface which has been remeshed along the split
skirt_remeshed = pv.wrap(f.GetOutputDataObject(2))  # the skirt which has been remeshed along the split

p = pv.Plotter(shape=(1, 2))
p.subplot(0, 0)

p.subplot(0, 1)
p.show()
``````

1 Like

I had the same problem with vtkCookieCutter. There was a bug in an internal method dealing with line intersections which is now fixed, like 4 weeks ago. You could try using/building the latest vtk version from gitlab and see if that fixes things for you.

1 Like

@whophil, If you’d like to try another alternative, might I suggest TriSurfaceCutter?

All you need to do is turn on Embedding and turn off Removal.

``````    import tsc
...
cutter = tsc.tscTriSurfaceCutter()
cutter.SetInputConnection(triangulation.GetOutputPort())
cutter.SetLoopsConnection(1, loops.GetOutputPort())
cutter.RemoveOff()
cutter.EmbedOn()
...
``````

It is not yet integrated in VTK, so you have to build it separately and use the optional plugin with ParaView or in a Python/C++ application.

Cons:

1. It can only process triangulations. All other linear cells are passed through to output.
2. Always embeds the loops in the xy-plane. Imagine it to be a silhouette of the loop on the surface mesh as seen from +/- Z-axis.
3. It’s not yet a part of VTK.

Pros:

1. It uses custom-made almost robust line-line intersection functions on a 2D plane.
2. Precise geometric cuts, unlike filters that do scalar-based clipping.
3. It colours the embedded edges with a Constrained cell data array.
4. Similarly colours points that were acquired from loop polygons (named Acquired).
5. It’s faster than vtkCookieCutter (comes with a price, only triangles)
• Here are the options in ParaView.
1. Embed: On | Remove: Off | Invert: NA

2. Embed: On | Remove: On | Invert: On

3. Embed: On | Remove: On | Invert: Off

4. Embed: Off | Remove: On | Invert: On

5. Embed: Off | Remove: On | Invert: Off

4, 5 look that way since no cell-line intersection was computed (Embed: Off)

I think the first set of options is sufficient for you!