Bug: Extra faces created in cut_with_plane with incorrect pointdata

[JeffreyWardman]

I have a mesh that I slice via: mesh.cut_with_plane(origin=(0,0,0), normal=(1,0,0))

This mesh already had pointdata called Labels.

The mesh has labels 0-9. 0 is the background colour, 9 is closest to the plane cutting, then 8 is the object next to it, then 7 and so on. Only 0 and 9 should be on the plane intersection with the mesh.

Cutting with plane creates new vertices with Labels that are even numbered where labels 0 and 9 meet. That is, there are points for each time the boundary of the 0 and 9 labels change. The example is U shaped, where 9 is at the trough of the U and 0 is at the tips. Therefore there are at least 2 points where this issue occurs.

Now, to the issue: labels 2,4,6,8 are all generated in the intersection between the 0 and 9 labels for these newly created faces (note: I don't actually care about these new faces). Oddly enough, no odd numbers are added.

It looks like the pointdata is estimated incorrectly. Labels 1, 2, 3, 4, etc are all different instances of dogs for example. So if I say Dog 2 and 4 are next to each other so the point belongs to dog 3 (as the mean of 2 and 4 is 3) then that causes problems.

I'll try reproduce this with data I can share.

[JeffreyWardman]

An 8 gets created for the LHS but not RHS in this example even though points are only assigned labels 0 or 9.

import numpy as np
import vedo

sphere = vedo.Sphere()
sphere.pointdata["Labels"] = np.zeros(sphere.npoints).astype(np.uint16)
sphere.copy().cut_with_plane(origin=(0, 0, 0), normal=(0, 0, -1))

xmin, xmax, ymin, ymax, zmin, zmax = sphere.bounds()

sphere.pointdata["Labels"][np.nonzero(sphere.vertices[:, 1] > (ymin + ymax) / 2)[0]] = 0
sphere.pointdata["Labels"][np.nonzero(sphere.vertices[:, 1] <= (ymin + ymax) / 2)[0]] = 9

sphere.cmap("jet", sphere.pointdata["Labels"])

vedo.show(sphere, axes=1).close()

left = sphere.copy().cut_with_plane(
    origin=((xmin + xmax) / 2 + 0.001, (ymin + ymax) / 2 + 0.001, (zmin + zmax) / 2 + 0.001), normal=(1, 0, 0)
)
l_labels = left.pointdata["Labels"]
print(np.unique(l_labels))

left.show().close()

right = sphere.copy().cut_with_plane(
    origin=sphere.vertices[0],
    normal=(1, 0, 0),
    invert=True,
)
r_labels = right.pointdata["Labels"]
print(np.unique(r_labels))
right.show().close()

The reason seems to be if cut_with_plane has some value that the vertices don't directly touch, therefore the algorithm generates new faces for some reason.

[marcomusy]

Hi - yes this is how VTK works, it interpolates the values to the new created vertices (which must be created to match the cutting plane surface):

import numpy as np
import vedo

vedo.settings.interpolate_scalars_before_mapping = False

sphere = vedo.Sphere(res=12).lw(1)
sphere.pointdata["Labels"] = np.zeros(sphere.npoints).astype(np.uint16)

xmin, xmax, ymin, ymax, zmin, zmax = sphere.bounds()

ids0 = np.nonzero(sphere.vertices[:, 1]  > 0)[0]
ids9 = np.nonzero(sphere.vertices[:, 1] <= 0)[0]
sphere.pointdata["Labels"][ids0] = 0
sphere.pointdata["Labels"][ids9] = 9
print(np.unique(sphere.pointdata["Labels"]))

sphere.cmap("jet", "Labels", vmin=0, vmax=9)
# vedo.show(sphere, axes=1).close()

left = sphere.copy().cut_with_plane(origin=[0.1,0,0])
l_labels = left.pointdata["Labels"]
print(np.unique(l_labels))
left.show(axes=1).close()

[0 9]
[0 3 4 5 6 9]

this is an expected behavior..

[JeffreyWardman]

Is there a way to cut_with_plane without creating the new points and faces? Or to customise the interpolation method/get the newly generated point IDs?

[marcomusy]

you can grab the resulting array and threshold it manually in python/numpy. Or check out examples: examples/basic/mesh_threshold.py

[JeffreyWardman]

Is it possible to output the ids of the new vertices to make it a bit simpler? Or to assign them a pointdata?

[marcomusy]

the only thing that comes up to my mind is again the .add_ids() ...