Mesh.intersect_with_line returns duplicated points

[baba-yaga]

Hi, Marco. I discovered the following unexpected behaviour of the mentioned procedure.

Draw a cube and identify its visible vertices from the cam:

from vedo import *

K = Cube().clean().tomesh().c('yellow').lighting('ambient').lw(3)

plt = Plotter()
plt += K
plt.show(interactive=False)
cam = plt.camera
pos = cam.GetPosition()
visible_ids = []
for i in range(8):
    pi = K.intersect_with_line(K.vertices[i], pos)
    print(f'i={i}, points of intersection: {pi}')
    if len(pi)==1:
        visible_ids.append(i)
plt.close()

Reproduce the plot with visible vertices marked in blue:

plt = Plotter()
plt += K.alpha(0.5)
plt += Points(K.vertices[visible_ids], r=20, c='blue')
plt.show(camera=cam).close()

No vertex is detected as visible! The code fails because the intersections pi for the front vertices contain duplicated points - the vertices themselves.

i=0, points of intersection: [[-0.5        -0.5        -0.5       ]
 [-0.36999699 -0.36999699  0.5       ]]
i=1, points of intersection: [[-0.5 -0.5  0.5]
 [-0.5 -0.5  0.5]
 [-0.5 -0.5  0.5]]
i=2, points of intersection: [[-0.5  0.5  0.5]
 [-0.5  0.5  0.5]
 [-0.5  0.5  0.5]]
....

The behaviour seems to be related to the angle between the view line and the gradient at the vertex. If you redefine
K = K.rotate_x(45).rotate_y(45).rotate_z(45)
the closest vertex is detected, but the others still not

[marcomusy]

Hi @baba-yaga this may be a precision problem (?)
Not sure about it.. maybe cleaning the points is a good enough solution?

from vedo import *

K = Cube().clean().c('yellow',0.2).lighting('off').lw(3)

plt = Plotter(interactive=False)
plt += K
plt.show()

pos = plt.camera.GetPosition()

hits = []
lines = []
for i in range(8):
    pi = K.intersect_with_line(K.vertices[i], pos)
    lines.append(Line(K.vertices[i], pos))
    hits.extend(pi)
    print(f'i={i}, points of intersection:\n {pi}\n')

hits = Points(hits, r=20, c='blue').clean()
plt += lines, hits

print("K.coordinates   ", K.coordinates)
print("hits.coordinates", hits.coordinates)

plt.interactive().close()

PS: no need to call .tomesh()

[baba-yaga]

Thanks, Marco!
Your suggestion to use clean() removes duplicates, but I think the problem is with the algorithm itself. In the following example the vertex 6 of the rotated cube has NO intersection with the cube!

K = Cube().clean().c('yellow').lighting('ambient').lw(3).rotate_x(45).rotate_y(45).rotate_z(45)

plt = Plotter()
plt += K
plt.show(interactive=False)
cam = plt.camera
pos = cam.GetPosition()
visible_ids = []
for i in range(8):
    # pi = K.intersect_with_line(K.vertices[i], pos)
    pip = Points(K.intersect_with_line(K.vertices[i], pos)).clean()
    pi = pip.vertices
    print(f'i={i}, points of intersection: {pi}')
    if len(pi)==1:
        visible_ids.append(i)
plt.close()

plt = Plotter()
plt += K.alpha(0.5)
plt += Points(K.vertices[visible_ids], r=20, c='blue')
plt.show(camera=cam).close()

Output:

i=0, points of intersection: [[-0.60355335 -0.60355335 -0.1464466 ]]
i=1, points of intersection: [[ 0.25       -0.75        0.35355338]]
i=2, points of intersection: [[0.10355338 0.10355338 0.85355335]]
i=3, points of intersection: [[-0.75        0.25        0.35355338]]
i=4, points of intersection: [[-0.10355338 -0.10355338 -0.85355335]
 [-0.0790121  -0.0790121   0.5953668 ]]
i=5, points of intersection: [[-0.25        0.75       -0.35355338]]
i=6, points of intersection: []
i=7, points of intersection: [[ 0.75       -0.25       -0.35355338]]

[marcomusy]

You are right I think this is an instability/bug in the vtk, so not much I can do to fix it... Eg

from vedo import *

K = Cube().clean()
K.rotate_x(45).rotate_y(45).rotate_z(45)
K.c('yellow', 0.25).lighting('off').lw(3)

plt = Plotter()
plt += K
plt.show(interactive=False)

# pos = np.array(plt.camera.GetPosition())
pos = np.array([1,2,-3]) # this works
# pos = np.array([1,2,3]) # this doesnt work
print(f'camera position: {pos}')

for i in range(8):
    pip = Points(K.intersect_with_line(K.vertices[i], pos)).clean()
    plt.add(pip.ps(10).c("blue5"))
    plt.add(Line(K.vertices[i], pos).c("red5"))
    print(f'i={i}, points of intersection: {pip.vertices}')
plt.show(axes=1, interactive=True).close()

[baba-yaga]

Thank you, Marco! Would it be too much to ask you to file a bug report to VTK project? Every time I take a deep breath to start learning VTK and every time it takes my breath away too quickly :)

[sergei9838]

A dirty hack (using your data above with a problematic pos = [1,2,3]:

from vedo import *


def intersection_points(msh: vedo.Mesh, ln: vedo.Line):
    """
        Add to the intersection the line's endpoints close enough to the mesh and clean the duplicates
    """
    int_pts = Points(msh.intersect_with_line(ln))
    close_endpts = [ln.vertices[i] for i in range(2) if ln.distance_to(msh)[i]<1e-6]
    if len(close_endpts) > 0:
        int_pts = merge(int_pts, Points(close_endpts))
    return int_pts.clean().vertices


plt = Plotter()
plt += K
plt += Point(pos, r=20, c="green")
for i in range(K.nvertices):
    L = Line(pos, K.vertices[i])
    plt += L.c("gray").lw(2).lighting("off")
    plt += Points(intersection_points(K, L), r=20, c="red")
plt.show().close()

No idea, though, why I get rainbow colours of the lines instead of grey. A modern trend...? ;)

[marcomusy]

That is because you pass L to intersection_points() and the "Distance" array is added:

print(L)

vedo.shapes.Line at (0x5582df03dd10)                                       
name          : Line
elements      : vertices=2 polygons=0 lines=1
position      : (0, 0, 0)
scaling       : (1.00000, 1.00000, 1.00000)
size          : average=2.02307, diagonal=4.04615
center of mass: (0.875000, 0.875000, 1.32322)
bounds        : x=(0.750, 1.00), y=(-0.250, 2.00), z=(-0.354, 3.00)
pointdata     : "Texture Coordinates" (float32), dim=2
pointdata *   : "Distance" (float64), dim=1, range=(0, 3.00)

just invert lines:

    plt += Points(intersection_points(K, L), r=20, c="red")
    plt += L.c("gray").lw(2).lighting("off")

[sergei9838]

Thank you, Marco!
Then I should use ln.clone() in the function as it should not mutate the arguments.
Have a nice weekend!
PS. Any news from vtk on this?

[marcomusy]

Sorry for the late reply... no news..