ray casting, cast multiple rays at once

[ttsesm]

Hi guys,

I am trying to use the ray casting functionality for a PolyData, and I was looking at the corresponding example https://docs.pyvista.org/examples/01-filter/poly-ray-trace.html. However, it seems that at the moment you cover only the case of one ray each time, where you provide one origin and one direction. Thus I was wondering if it is possible to apply the casting procedure for multiple rays at once and ideally without a for loop, for example:

# Create source to ray trace
sphere = pv.Sphere(radius=0.85)

# Define line segments
start = [0, 0, 0]
stop = [[0.25, 1, 0.5],[0.5, 1, 0.25], [0, 0, 1]]

# Perform ray trace
points, ind = sphere.ray_trace(start, stop)

Thanks.

[akaszynski]

Sadly, at the moment this isn't something that vtk seems to support out of the box. Constructing the obbTree first takes a bit of time, but every ray trace afterwards will be a bit faster. I'd love to write an additional interface that vectorizes these functions, but I've had no luck getting vtk to work with cython.

[banesullivan]

FYI: #80 and pyvista/pyvista#478

[banesullivan]

You may be able to vectorize the sphere.ray_trace across an array of inputs with NumPy

[ttsesm]

Thanks for the feedback guys. @banesullivan do you mean something like with numpy.apply_along_axis. Actually I just found it, it seems interesting. I will give it a try.

Also is it possible to draw all the rays at once, with pvista.Line() again you need to go one by one. Then I was checking on the pvista.plot_arrows() (https://github.com/pyvista/pyvista/blob/a20729ae508e582fce2fc417081a88a2823f6199/tests/test_plotting.py#L407, https://github.com/pyvista/pyvista/blob/a20729ae508e582fce2fc417081a88a2823f6199/tests/test_plotting.py#L414) but I couldn't figure out how to include it in my plotter. Is it easy to provide an example. For instance:

import pyvista as pv

# Create source to ray trace
sphere = pv.Sphere(radius=0.85)

# Define line segment
start = [0, 0, 0]
stop = [0.25, 1, 0.5]

# Perform ray trace
points, ind = sphere.ray_trace(start, stop)

# Create geometry to represent ray trace
ray = pv.Line(start, stop) <---------- How to replace this line with plot_arrows()
intersection = pv.PolyData(points)

# Render the result
p = pv.Plotter()
p.add_mesh(sphere,
           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") <---------- How to replace this line with plot_arrows()
p.add_mesh(intersection, color="maroon",
           point_size=25, label="Intersection Points")
p.add_legend()
p.show()

[ttsesm]

I've tried to vectorize the sphere.ray_trace() by using both numpy.vectorize() and numpy.apply_along_axis(), e.g.

v_ray_trace = np.vectorize(sphere.ray_trace)
points, ind = v_ray_trace([0,0,1], [[0,0,1], [0,0,1]])

but it doesn't seem to work, I am getting TypeError: IntersectWithLine argument %Id: %V, similarly in numpy,apply_along_axis().

Any idea how we might be able to make it work.

[banesullivan]

plot_arrows plots vector data not necessary rays or lines... its meant for plotting a full 3D vector field where you have a point, direction, and magnitude for each arrow rather than start/stop coordinates.

pv.Line() creates a line mesh and it cannot be converted to an arrow. You coul manually create arrow geometry with: ray = pv.Arrow(start, stop). But that isn't really what I think you are looking for...

---

As for vectorizing the ray tracing...

do you mean something like with numpy.apply_along_axis. Actually I just found it, it seems interesting. I will give it a try.

Yep! That's what I was getting after. However, numpy.vectorize provides no performance benefit and can only handle scalar arguments so it won't really work in this case (my bad). actually numpy.vectorize is just a Python for-loop under the hood. Are you deadset on avoiding a for-loop? If so, you'll have to do some manual parrallelizing of this routine launching a bunch of threads each computing their own ray trace.

Also, the ray_trace method can have variable results depending on how many cells of the mesh are intersected if any, so it's not guaranteed that your output array would have a valid shape in NumPy

Here is a solution to perform batch ray tracing (a for-loop) that collects the first intersected cell IDs for the rays:

import numpy as np
import pyvista as pv

sphere = pv.Sphere(radius=0.85)

def ray_trace(start, stop):
    """Pass two same sized arrays of the start and stop coordinates for all rays"""
    assert start.shape == stop.shape
    assert start.ndim == 2
    # Launch this for loop in parallel if needed
    zeroth_cellids = []
    for i in range(len(start)):
        _, ids = sphere.ray_trace(start[i], stop[i])
        if len(ids) < 1:
            v = None
        else:
            v = ids[0]
        zeroth_cellids.append(v)
    return np.array(zeroth_cellids)

start = np.array([[0,0,0], [0,0,0]])
stop = np.array([[0.25, 1, 0.5], [0.5, 1, 0.25]])
cell_ids = ray_trace(start, stop)

cells = sphere.extract_cells(cell_ids)

p = pv.Plotter(notebook=0)
p.add_mesh(sphere,
           show_edges=True, opacity=0.5, color="w",
           lighting=False, label="Test Mesh")
p.add_arrows(start, stop, mag=1, color=True, opacity=0.5, )
p.add_mesh(cells, color="maroon",
           label="Intersection Cells")
p.add_legend()
p.show()

[ttsesm]

Thanks for the feedback. Indeed numpy.vectorize is not gonna help much, as you said it is a for loop under the hood.

The solution with the multiple threads is an option, and then I was checking on the the trimesh lib which also provides a ray casting method accelerated with embree, which seems a feasible solution as well.

add_arrows() seems to do what I want, but can you also show how to use the plot_arrows() on an existing plotter. It might get in hand somehow, if I have only the direction vectors instead of the stop points.

[banesullivan]

but can you also show how to use the plot_arrows() on an existing plotter. It might get in hand somehow, if I have only the direction vectors instead of the stop points.

plot_arrows() is not intended to be used on an existing plotter... it is a helper around BasePlotter.add_arrows, so you'd want to use that. Literally all plot_arrows() helper does is the following routine:

plotter = Plotter()
plotter.add_arrows(cent, direction)
plotter.show()

there is nothing special about it.

[Keou0007]

Has there been any further progress on this?

I'm trying to create a 2D contour plot from the Elevation data of a 3D surface. The issue is that any [x, y] coordinate could have multiple points of the surface along the z-axis and I want to plot the contour of only the maximum elevation at any point on the plane.

My solution so far, which is based on the distance between two surfaces example, is to create a 2D uniform grid above the surface then ray trace along -Z from every point and store the height of the first point found. The problem with this is that the python for loop can only get through around 200 ray traces per second. If I want my uniform grid to contain, say, 1 million points, then it's going to take hours to process, then times by a fairly large number of surfaces and I end up in months.

The best I've been able to find on google is the suggestion that this sort of thing would be much much faster done in C++, but that doesn't really help me much when the rest of my project is in Python. I have tried threading this operation in the past but ran into weird errors that seemed to suggest the VTKOBBTree wasn't thread safe.

[banesullivan]

There hasn't been any progress on this... 😞

Perhaps the note in #173 (comment) helps?

The solution with the multiple threads is an option, and then I was checking on the the trimesh lib which also provides a ray casting method accelerated with embree, which seems a feasible solution as well.

Otherwise, it may be a long time before any pyvista developers find the time to look into a way of doing this in batch with the VTK python bindings (I suspect it should be possible though)

[ttsesm]

Has there been any further progress on this?

I'm trying to create a 2D contour plot from the Elevation data of a 3D surface. The issue is that any [x, y] coordinate could have multiple points of the surface along the z-axis and I want to plot the contour of only the maximum elevation at any point on the plane.

My solution so far, which is based on the distance between two surfaces example, is to create a 2D uniform grid above the surface then ray trace along -Z from every point and store the height of the first point found. The problem with this is that the python for loop can only get through around 200 ray traces per second. If I want my uniform grid to contain, say, 1 million points, then it's going to take hours to process, then times by a fairly large number of surfaces and I end up in months.

The best I've been able to find on google is the suggestion that this sort of thing would be much much faster done in C++, but that doesn't really help me much when the rest of my project is in Python. I have tried threading this operation in the past but ran into weird errors that seemed to suggest the VTKOBBTree wasn't thread safe.

I personally went with trimesh and embree (it supports only until 2.17) were I was able to cast 18 million rays in less than >20 sec. Check here the discussion mikedh/trimesh#875 you might be able to use embree v.3 with this wrapper though https://github.com/sampotter/python-embree which might improve casting time even further (unfortunately I did not have the time to test it myself).

[Keou0007]

I got this working with the trimesh library. Unfortunately trimesh only supports meshio formats and meshio doesn't support VTK polydata, so I had to convert my surface into a .ply file first, but the results are impressive nonetheless.

Code

from pathlib import Path
from time import time

import trimesh
import pyvista as pv
import numpy as np

vtk_file = Path("test_file.vtp")

vista_mesh = pv.read(vtk_file)
vista_mesh.clear_arrays()

xmin, xmax, ymin, ymax, zmin, zmax = vista_mesh.bounds
points_per_metre = 20
dims = [
    int(round(xmax - xmin, 2) * points_per_metre) + 1,
    int(round(ymax - ymin, 2) * points_per_metre) + 1,
    1,
]
origin = [xmin, ymin, zmax + 1]
spacing = [1 / points_per_metre] * 3

ug = pv.UniformGrid(dims, spacing, origin)
ug["heights"] = np.empty(ug.n_points)
ug["heights"][:] = np.nan

vec = [0, 0, zmin - zmax - 1]
start = time()
for i in range(ug.n_points):
    p = ug.points[i]
    ip, ic = vista_mesh.ray_trace(p, p + vec, first_point=True)
    if ip.any():
        ug["heights"][i] = ip[2]

print(f"for-loop ray tracing with pyvista did {ug.n_points} points in {time()-start:.3f} s")

ug2 = pv.UniformGrid(dims, spacing, origin)
ug2["heights"] = np.empty(ug2.n_points)
ug2["heights"][:] = np.nan

vista_mesh.save("test_file.ply")
trimesh_mesh = trimesh.load("test_file.ply")
plane_points = ug2.points
vectors = np.array([[0, 0, -1]] * len(plane_points))

start = time()
loc, idr, idt = trimesh_mesh.ray.intersects_location(plane_points, vectors, multiple_hits=False)
print(f"\nvectorised ray tracing with trimesh did {ug2.n_points} points in {time()-start:.3f} s")

start = time()
for i, location in zip(idr, loc):
    ug2["heights"][i] = location[2]
print(f"plus {time()-start:.3f} s to unpack the results")

points_per_metre = 200
dims = [
    int(round(xmax - xmin, 2) * points_per_metre) + 1,
    int(round(ymax - ymin, 2) * points_per_metre) + 1,
    1,
]
origin = [xmin, ymin, zmax + 1]
spacing = [1 / points_per_metre] * 3

ug3 = pv.UniformGrid(dims, spacing, origin)
ug3["heights"] = np.empty(ug3.n_points)
ug3["heights"][:] = np.nan

plane_points = ug3.points
vectors = np.array([[0, 0, -1]] * len(plane_points))
start = time()
loc, idr, idt = trimesh_mesh.ray.intersects_location(plane_points, vectors, multiple_hits=False)
print(f"\nvectorised ray tracing with trimesh did {ug3.n_points} points in {time()-start:.3f} s")

start = time()
for i, location in zip(idr, loc):
    ug3["heights"][i] = location[2]
print(f"plus {time()-start:.3f} s to unpack the results")

and the output

for-loop ray tracing with pyvista did 11431 points in 127.226 s

vectorised ray tracing with trimesh did 11431 points in 6.147 s
plus 0.301 s to unpack the results

vectorised ray tracing with trimesh did 1122301 points in 1.232 s
plus 28.634 s to unpack the results

interestingly the majority of the initial 7 s is loading the file, which is then cached for the second run with many more traces.
The only downside seems to be a small number of traces that fail to find an intersection resulting in nan's in elevation data in the above code, but for my purpose of generating a contour plot, they can easily be ignored.

[banesullivan]

This is cool, thanks for posting a working solution, @Keou0007