Python code to check whether a point is inside a given mesh

Question

I am using the trimesh library, to check this. It has mesh.contains for performing these checks. I am trying to truncate an array of lines using the mesh.contains, such that the new lines only contain the points from inside the mesh.

Although this works fine for in some cases, it fails in multiple other cases. I am attaching a few pictures, the blue is the truncated line and the red is the original line. The gray shaded region is the mesh. In one case, it is able to detect points outside and truncate properly. However, there are issues in case of outside regions that are sandwiched between the domains. The points are outside, but this function is unable to eliminate them.

The code chunk is given below:

# create an empty array to store truncated lines
truncated_arr = []

# loop over each line in the array
for line in arr:
    new_line = []
    # loop over each point in the line and check if it's inside the mesh
    for point in line:
        if ((mesh1.contains([point])) or (mesh2.contains([point])) or (mesh3.contains([point]))):
            new_line.append(point)
    # append the truncated line to the new array
    if len(new_line) > 1:
        truncated_arr.append(new_line)

I am looking for ideas to overcome this. Is there a way to make this algorithm detect points that are outside the domain but is sandwinched between two wings of the domain, as seen in the picture?

Answer 1

A possible solution would be to use the mesh.ray.intersects_location() method to perform a ray casting for each point, and count the number of intersections to determine whether the point is inside or outside the mesh.

Here's a modified version of your code that implements this approach:

import numpy as np
import trimesh

def is_inside(point, meshes, direction=np.array([1.0, 0.0, 0.0]), tol=1e-5):
    count = 0

    for mesh in meshes:
        intersections, _ = mesh.ray.intersects_location(
            ray_origins=[point - direction * tol],
            ray_directions=[direction],
        )

        count += len(intersections) % 2

    return count % 2 == 1

# create an empty array to store truncated lines
truncated_arr = []

# the list of meshes
meshes = [mesh1, mesh2, mesh3]

# loop over each line in the array
for line in arr:
    new_line = []
    # loop over each point in the line and check if it's inside the meshes
    for point in line:
        if is_inside(point, meshes):
            new_line.append(point)
    # append the truncated line to the new array
    if len(new_line) > 1:
        truncated_arr.append(new_line)

Given a point, the is_inside() function ray-casts the point along an axis in the positive X direction, then it counts the number of times the ray intersects the mesh. If the number of intersections is odd, the point is assumed to be inside the mesh. You can change the direction of the ray by changing the direction parameter. The tol parameter is a small tolerance value to avoid starting the ray intersection exactly on the surface of the mesh.

Note that this approach might be slower than using the mesh.contains() method due to the additional ray intersection calculations. However, it should better handle points outside the domain but sandwiched between two wings of the domain.

Answer 2

The ray tracing proposed in the other answer seems overly resource intensive you already have an (presumably well-tuned) .contains function. I think instead you need to change your order of operations. I can't run this code, but it seems to me like if instead of chaining those three mesh.contain(.) or... you did all three separate then used numpy's set intersection on the three resulting lines you might have better luck.

Just a quick try, hope it works out better than trying to reimplement ray tracing. Not that it wouldnt work, idk but i bet it would, but it would be a ton of work to do and to be honest it's already been done for you. This is a logic error, not a deficiency in your tools.

Answer 3

Do you happen to have a volume mesh with that surface mesh as the boundary? I'm asking in case that surface mesh had been extracted. If you had a volume mesh, it could be easier to detect points outside the mesh. Moreover, it would be very cheap to test successive points.

Do you have the neighbours table?

Is it possible for you to choose the rays to intersect? In that case, it would be cheaper, more accurate and more informative to intersect the segment [point[i],point[i+1]] with the surface mesh rather than sending random rays from point[i]. You know the initial point point[0] is outside (presumably, though that's an edge case you can manage). Then, you only need to know if the following segments p[i] + tp[i+1] intersect the mesh for 0 < t <= 1. If the segment intersects, you also know exactly where, which avoids having gaps between the truncated lines and the mesh, as you would have by excluding segment ends.

If the segment intersects the mesh, you can keep in memory the triangle it intersected, put it in triangle_guesses. If you have neighbours, you can then use that as a starting point for a fast intersection search. Your segments are all packed together, it's very likely many are intersecting the same triangle, or triangles all close to each other.

From what I can tell, the RayMeshIntersector class may be modified to call contains_points with a check_direction argument other than the default. You could just code triangle-segment intersection yourself, too. Having neighbours, you'd even be faster than that algorithm because you could use previous guesses to restart the search.

If the algorithm is failing ray intersection, perhaps you have flat elements (vertex against opposite edge) or the mesh is not watertight.