pyvista ray trace detects only one intersection

Question

I try to do detect the intersections of a ray with a polygon using pyvista. In this case the polygon is a cube. The cube is defined by points in 2d and then extruded. Why does this simple example where a ray is intersecting with the cube detect only one intersection?

import numpy as np
import pyvista as pv

points = [(0, 0), (0, 10), (10, 10), (10, 0)]

def points_2d_to_poly(points, z):
        """Convert a sequence of 2d coordinates to polydata with a polygon."""
        faces = [len(points), *range(len(points))]
        poly = pv.PolyData([p + (z,) for p in points], faces=faces)
        return poly

polygon = points_2d_to_poly(points, 0.0)
polygon_extruded = polygon.extrude((0, 0, 6.0), capping=True)
poly_surface = polygon_extruded.extract_surface()
poly_mesh = poly_surface.triangulate()

# Define line segment
start = [5, 2, -10]
stop = [5, 2, 10]

# Perform ray trace
points, ind = poly_mesh.ray_trace(start, stop, first_point=False)

# Create geometry to represent ray trace
ray = pv.Line(start, stop)
intersection = pv.PolyData(points)

# Render the result
p = pv.Plotter()
p.add_mesh(poly_mesh, show_edges=True, opacity=0.5, color="w", lighting=False, label="Test Mesh")
p.add_mesh(ray, color="blue", line_width=5, label="Ray Segment")
p.add_mesh(intersection, color="maroon", point_size=25, label="Intersection Points")
p.add_legend()
p.show()

Answer 1

poly_mesh normals are missing, and this is required per the documentation of the vtk method used. https://vtk.org/doc/nightly/html/classvtkAbstractCellLocator.html#a027818794762a37868a3dccc53ad6e81

This method assumes that the data set is a vtkPolyData that describes a closed surface, and the intersection points that are returned in 'points' alternate between entrance points and exit points.

So if the normals are computed first, it works as expected.

import numpy as np
import pyvista as pv

points = [(0, 0), (0, 10), (10, 10), (10, 0)]

def points_2d_to_poly(points, z):
        """Convert a sequence of 2d coordinates to polydata with a polygon."""
        faces = [len(points), *range(len(points))]
        poly = pv.PolyData([p + (z,) for p in points], faces=faces)
        return poly

polygon = points_2d_to_poly(points, 0.0)
polygon_extruded = polygon.extrude((0, 0, 6.0), capping=True)
poly_surface = polygon_extruded.extract_surface()
# Note this line and following line are the only changes
poly_mesh = poly_surface.triangulate().compute_normals()

# Define line segment
start = [5, 2, -10]
stop = [5, 2, 10]

# Perform ray trace
points, ind = poly_mesh.ray_trace(start, stop, first_point=False)

# Create geometry to represent ray trace
ray = pv.Line(start, stop)
intersection = pv.PolyData(points)

# Render the result
p = pv.Plotter()
p.add_mesh(poly_mesh, show_edges=True, opacity=0.5, color="w", lighting=False, label="Test Mesh")
p.add_mesh(ray, color="blue", line_width=5, label="Ray Segment")
p.add_mesh(intersection, color="maroon", point_size=25, label="Intersection Points")
p.add_legend()
p.show()