faster trimesh.proximity.closest_point for large query point sets

[martinResearch]

Thank you for this amazing library.

I would like to compute the distance of a large number of point to a relatively simple mesh.
I am using trimesh.proximity.closest_point, but unfortunately it is quit slow : it takes 30 second when querying 500000 points on my machine.

One major bottleneck seem to be the line candidates = [list(rtree.intersection(b)) for b in bounds] in proximity.py.
I believe this could run much faster if the loop was implemented in c++. Ideally rtree.intersection would accept a list of boxes and implement the loop in c++ and return the list of list of candidates (similarly to scipy.spatial.KDTree.query_ball_point), but that does not seem supported by rtree.

Here is the code I use to get some timings

import trimesh
import time
import numpy as np
import matplotlib.pyplot as plt

vertices = np.random.rand(20, 3)
faces = (np.random.rand(10, 3) * 20).astype(np.int)


mesh = trimesh.Trimesh(faces=faces, vertices=vertices)

nb_points_list = [1000, 10000, 100000, 500000]
durations_list = []

for nb_points in nb_points_list:
    points = np.random.rand(nb_points, 3)

    start = time.time()
    trimesh.proximity.closest_point(mesh, points)
    duration = time.time()-start
    durations_list.append(duration)
    print(f"trimesh.proximity.closest_point wit {nb_points} points took {duration} seconds")

plt.plot(nb_points_list, durations_list)
plt.xlabel("nb qurey points")
plt.xlabel("duration in seconds")
plt.show()

I get this timings
trimesh.proximity.closest_point wit 1000 points took 0.11966729164123535 seconds
trimesh.proximity.closest_point wit 10000 points took 0.6350014209747314 seconds
trimesh.proximity.closest_point wit 100000 points took 6.338748455047607 seconds
trimesh.proximity.closest_point wit 500000 points took 32.05847191810608 seconds

It would be great if we could list the potential ways we could accelerate the function. I looked at alternative libararies available from python to do the nearest points query but did not find a good alternative yet.

[martinResearch]

added an issue on the rtree repo Toblerity/rtree#178

[mikedh]

Agreed, vectorizing rtree.Index.intersection would be awesome here (and elsewhere).

[martinResearch]

tried https://pypi.org/project/pysdf/ that is amazingly fast.
with 500000 points I get

trimesh.proximity.closest_point took 31.534400939941406 seconds
kdtree.query   took 0.25 seconds
pysdf   took 0.058998 seconds

pysdf does not return the nearest face index for each point but could probably be modified to do that, in which case it could maybe be used a a dependency to accelerate trimesh's function

[mikedh]

Whoa that is awesome! Yeah possibly some packaging work to do to get that included as a dep but really impressive results.

[martinResearch]

I figured the distance computed by pysdf is actually approximate as, according to the readme it computes "only distance to faces adjacent to nearest neighbors". For each point it searches for the nearest vertex and then enumerate only the adjacent faces as face candidates for the nearest face.

[martinResearch]

The module python-prtree has been recently improved to implements a rtree on 3d data . The code is in C++ with pybind11 binding that provide a vectorized interface . It could be good alternative to Toblerity/rtree.

[kickd0wn]

@mikedh Do you plan to implement/add https://github.com/atksh/python_prtree or a similar solution to trimesh? It would be very nice to use this function to analyze huge 3D point clouds.

Unfortunately, my programming skills are not good enough to implement it myself.
I hope that you can help me and other users with this. Thanks!

[guatavita]

@kickd0wn If you are still looking for something faster you can use IGL

signed_dist, _, _ = igl.signed_distance(coords, mesh.vertices, mesh.faces)

[cosw0t]

One more suggestion, with open3d:

import open3d as o3d

scene = o3d.t.geometry.RaycastingScene()
scene.add_triangles(o3d.t.geometry.TriangleMesh.from_legacy(mesh.as_open3d))
ret_dict = scene.compute_closest_points(o3d.core.Tensor.from_numpy(coords.astype(np.float32)))

[mikedh]

Oh nice, also should we update to_open3d so it doesn't need that function labeled legacy?