Inaccuracies in mesh.contains

[LMescheder]

I experiencing some issues with the Trimesh.contains method. I tried to debug it myself, but I couldn't find any issues in the trimesh code. I suspect it's inaccuracies in PyEmbree.

Here is the code I'm using:

import numpy as np
import trimesh
from plyfile import PlyElement, PlyData

# Only for exporting the pointcloud
def export_pointcloud(vertices, out_file):
    assert(vertices.shape[1] == 3)
    vertices = vertices.astype(np.float32)
    vector_dtype = [('x', 'f4'), ('y', 'f4'), ('z', 'f4')]
    vertices = vertices.view(dtype=vector_dtype).flatten()
    plyel = PlyElement.describe(vertices, 'vertex')
    plydata = PlyData([plyel], text=True)
    plydata.write(out_file)

# Load mesh
mesh = trimesh.load('voxels.off', process=False)
points = np.random.rand(100000, 3) - 0.5
contained = mesh.contains(points)
export_pointcloud(points[contained], 'points.ply')

The result looks like this:

Here are the input / output files:
https://www.dropbox.com/sh/6i0222tk9lpkqgr/AAC1mZ0yblFsXxg9cEMmfneJa?dl=0

Of course, in this special case I can do it in a simpler way by exploiting the regular grid structure of the mesh. However, I also experience similar problems for other meshes (usually at higher resolutions).

I would appreciate any help!

[LMescheder]

I ended up writing my own cython code for checking if a point is inside a mesh, based on this answer stack overflow:
https://stackoverflow.com/a/6576840
As far as I can tell, it works quite robustly and yields more accurate results.
I would be happy to contribute the code (not now, but in a few weeks) if you are interested.

[mikedh]

Thanks for the report and reproducible example!

Yeah this looks like it's related to embree. Trimesh has two ray implementations with the same API, one with embree (which is hilariously faster, like 600K rays/sec) and one that's mostly straight numpy (12K RPS). The embree one only returns the first hit by default, and then trimesh has some wangling code to offset (to avoid a self-hit) then re-query. Avoiding the self hit means that in your example mesh with internal faces (from the non- surface voxelization) don't get counted correctly. The easiest thing to do is probably use the native ray- tracer, which gets the correct answer but is kinda slow. Once I push the changes that allow no- face PLY files to load, this works pretty well:

m = trimesh.load('voxels.off',
                    use_embree=False)
p = trimesh.load('points.ply')

c = m.contains(p.vertices)
trimesh.Scene([m, trimesh.PointCloud(p.vertices[c])]).show()

One other option for the future- you could marching- cubes convert your voxels to a mesh which wouldn't have the internal faces.

The nice thing about embree is that it is fast, although it doesn't handle this case all that well. A robust cython ray implementation could be great!

Trimesh currently doesn't have any built components, so adding any cython would be a substantial architectural change. The way we've been doing it is breaking up compiled libraries into small separate packages and setting up cibuildwheel (ala python-fcl) but it is a decent amount of work per- package. If there's a better way I'm also definitely open to it!

I think the contains logic would probably be mostly the same (count ray hits) but with more accurate rays.

Some other things we might do to help with this:

1. speed up current numpy/rtree ray tracer, particularly by making broad phase less stupid (the 2D algorithm you linked looks good).
2. We know which point queries are broken due to bidirectional queries, so fall back to native rays but only for those points.

[LMescheder]

Quick follow up: we released our cython code where we perform the inside/outside check. You can find it here. More specifically, we have implemented the function check_mesh_contains that performs the inside/outside check using the MeshIntersector class. Most code is in vanilla python, although there is a small cython kernel (the TriangleHash).

In our tests, this function gave more accurate results than the trimesh .contains method.
Maybe you find it useful given that this issue is still open.

[oliverstefanov]

I am also getting strange results sometimes with "mesh.contains".

I want to try your check_mesh_contains routine to compare but I'm getting the following error:

from .triangle_hash import TriangleHash as _TriangleHash
ImportError: attempted relative import with no known parent package

I included triangle_hash.pyx in my directory, but how can i import it properly in my main script ?

[LMescheder]

@oliverstefanov
You should probably follow the instructions from the Occupancy Networks Repo to build the whole package and then import check_mesh_contains using

from im2mesh.utils.libmesh import check_mesh_contains

If you experience any issues, please open a separate issue there.

[oliverstefanov]

@LMescheder Ok thanks!

[benoitrosa]

Hi,

For easy testing, you can download the inside_mesh.py and triangle_hash.pyx files from the Occupancy Network repo, and compile the pyx file. You just have to create a setup.py file with the following lines :

`from distutils.core import setup
from Cython.Build import cythonize

setup(
ext_modules = cythonize("triangle_hash.pyx")
)`

And then compile it using cython : python setup.py build_ext --inplace

The compilation creates a triangle_hash.so which can be imported in python normally, making the inside_mesh module work as a standalone one for using the check_mesh_contains method

[mikedh]

So if a mesh is loaded with trimesh.load(filename, use_embree=False), does anyone have a reproducible example of incorrect results? Or is this isolated to embree?

[benoitrosa]

After looking a bit more in depth, I found that the mesh I was using didn't have any normals computed. I didn't look into the details of the computations but I thought that trimesh was using ray tracing, which should not require face normals ?

Anyway, here are two dummy STL files which I used. They are simple STLs generated using PyMesh : pymesh.generate_icosphere(radius, center). I don't remember the radius I used, but the center is (0,0,0). Dummy.stl was directly generated and saved as is, and dummy2 was postprocessed using CloudCompare to compute normals.

I used the following code :

`import trimesh
import numpy as np

mesh = trimesh.load('dummy.stl')
mesh2 = trimesh.load('dummy.stl', use_embree=False)

print mesh.contains(np.array([[0,0,0],[200,200,200]]))
print mesh2.contains(np.array([[0,0,0],[200,200,200]]))

mesh3 = trimesh.load('dummy2.stl')
mesh4 = trimesh.load('dummy2.stl', use_embree=False)

print mesh3.contains(np.array([[0,0,0],[200,200,200]]))
print mesh4.contains(np.array([[0,0,0],[200,200,200]]))`

The output is :

[False False] [False False] [ True False] [ True False]

One can see, that the mesh without normals behaves incorrectly irregardless of the pyembree setting. The mesh with normals behaves good also irregardless of the pyembree setting. Perhaps having some check for mesh normals (if they are required) could be a good idea ? For the same meshes, the is_watertight property is OK (True for all)

dummy_stls.zip

[oliverstefanov]

@benoitrosa Bonjour Benoit, do you confirm that precomputing normals on the mesh resolves the issues encountered with mesh.contains ?
Thanks for your insight

[benoitrosa]

I just re-tested with another (more complex) mesh and indeed, the mesh.contains function behaves correctly when the mesh normals are pre-computed.

Note that I did compute the normals using Cloudcompare. The fix_normals() function of trimesh does not do anything in this case (but that's perhaps another bug due to the format of the input mesh, all normals are with [0,0,0] values)

[oliverstefanov]

It is still not working for me, I get wrong results for points located beyong the mesh in the y+ direction only...!? Computing or recomputing the normals makes no difference.

I tried installing the inside_mesh/triangle_hash as @benoitrosa describes but I am not getting a triangle_hash.so file, only .cpp file although the compiling seems ok (after including numpy in the setup.py). So I don't know how to import the check_mesh_contains method into my python code...

Any help or insight would be greatly appreciated because this unreliability in the mesh.contains function kinda ruins everything I've done so far in my program :-(
@mikedh any ideas?

[benoitrosa]

@oliverstefanov There was a mistake in the command needed. It should be python setup.py build_ext. When you compile it normally creates a build folder, with subdirectories containing the built files (the .o and the .so). If you copy the so file in the folder where the inside_mesh.py is, then you can import the inside_mesh module and use the check_mesh_contains file. Hope that helps !

[LMescheder]

@benoitrosa
For simplicity, I just created a public gist for the standalone version of check_mesh_contains.

Compiling it with python setup.py build_ext --inplace should do the trick

[oliverstefanov]

Thank you for your response. When compiling I am getting .pyd/.obj/.lib/.exp (as I am on Windows?) in the build directory. I moved everything to my main directory along with triangle_hash.pyx and inside_mesh.py. I included "import inside_mesh" in my code, which leads to the following error:

from .triangle_hash import TriangleHash as _TriangleHash ImportError: attempted relative import with no known parent package

[LMescheder]

@oliverstefanov Can you try to replace the line with

from triangle_hash import TriangleHash as _TriangleHash

as I did in my gist?

[oliverstefanov]

Oh ok thanks, but I think i've tried this before, I still get an error:

from triangle_hash import TriangleHash as _TriangleHash ImportError: cannot import name 'TriangleHash' from 'triangle_hash'

[mikedh]

Hey @LMescheder that self contained bit is awesome! What do you think about putting it on pypi, kind of like triangle or python-fcl?

I made a first pass at packaging your gist here: https://github.com/mikedh/contains

It compiles and works nicely, with cibuildwheel and automatic pypi releases it could be pretty slick. I'm happy to transfer that repo to you or add you to it, and apologies for the automatic license in there (delete away).

[LMescheder]

@mikedh Sure, I can look into that if you add me or transfer the repo.

[mikedh]

Sweet! Sorry this fell off my radar, transferred!

[AtinAngrish]

Hi @mikedh , is this bug fixed or should I be using the gist supplied by @LMescheder ?

Thanks,

A

[avlonder]

Is the Pyembree issue solved by temporally removing the face that was hit instead of translating the ray origin beyond the face?

[AymericFerreira]

Hello

I still have problems with contains. Normal contains is too slow in my case and pyembree return wrong results.

I tried check_mesh_contains but the computing time is very strange. I tested on a mesh generated by an algorithm
And on a mesh generated by trimesh (cylinders)

Computing time is crazy for the cylinder mesh despite he contains 2 times less vertices and faces. I tried to check the type of variable inside the mesh but they all seems identical.

Here is my script

import trimesh
from inside_mesh import check_mesh_contains
import numpy as np
import time

def timeme(method):
    def wrapper(*args, **kw):
        startTime = int(round(time.time() * 1000))
        result = method(*args, **kw)
        endTime = int(round(time.time() * 1000))

        print(f'{endTime - startTime} ms')
        return result

    return wrapper

points = np.random.uniform(size=(1000000, 3)).astype(dtype=np.int64)


mesh = trimesh.load_mesh('mesh.ply')
cyl = trimesh.load_mesh('cyl2.ply')

@timeme
def check(mesh, points):
    check_mesh_contains(mesh, points)

print(f'Mesh vertices : {mesh.vertices.shape}')
print(f'Mesh faces {mesh.faces.shape}')
print('Mesh : ')
check(mesh, points)

print(f'Cylinder vertices {cyl.vertices.shape}')
print(f'Cylinder faces : {cyl.faces.shape}')
print(f'Cylinder : ')
check(cyl, points)

And here the result :

Mesh vertices : (120651, 3)
Mesh faces (238646, 3)
Mesh : 
194 ms
Cylinder vertices (76032, 3)
Cylinder faces : (147456, 3)
Cylinder : 
10370 ms

Why this function in the cylinder mesh is so long ?
I added the 2 meshes in an archive.
I post it here because I suspect it can be bind to generated mesh from trimesh but I don't know how to correct this
Sadly I can't use non pyembree contains because it's too low and memory consumption goes crazy

meshes.zip

[aluo-x]

For issue #263 and #331 I think the following could be done while still utilizing embree, in pseudo code

Given Face array F of shape [F, 3]
Given Vertex array V of shape [V, 3]
Collect vertices for each face V_new = V[F]  of shape [F, 3, 3]
V_soft = []
For each face index F_i in range(F):
    N = normal(F_i)
    V_new_i = expand_face_verts(V_new[F_i], N, tol = 0.1)
    V_soft.append(V_new_i)
# F_new is a new face array of consecutive numbers
# V_soft is an array of shape [F, 3]
# ray trace as usual
# use the collected face index, index back into the original vertex array V, compute bary etc.

This doesn't solve the multi-hit problem, but does give some tolerance for each ray to hit. And I think all of this could be implemented in a vectorized way without needing any loops.

Also to note, embree3 now has an example to gather all hits along a single ray. I've taken a brief look, and it seems to be just tracking the most recent hit, then offsetting the ray. Unfortunately the current wrapper for embree that Trimesh uses is for embree2, and the next_hit example relies on features only introduced in embree 3.7.0,

[aluo-x]

More experimentation w/ embree3 using the python-embree wrapper. Even with the ROBUST mode enabled, the first intersection location output of embree3 is not reliable when computed with either UV or tfar (changes even between runs).

I think the current Trimesh implementation is probably as good as it gets:

• Use embree to check for intersection
• get the id of the triangle that a ray hits
• Perform location computation in python.
• Avoid second collisions

The second/third/non-first collision is totally unreliable.
@mikedh Not sure if you have any thoughts on this. Did similar reasoning cause you to implement collision location computation in python?

[mikedh]

Nice thanks for experimenting! Yeah the "second hit" logic with offsetting is not great. Does robust mode solve the "leaking through on-vertex or on-edge" issues with embree? Embree seems like it is pretty much only designed for first-hits as far as I can tell.

If they actually solved the leaks, I was wondering if something like this would work:

• query the first hit on the mesh
• for multiple hits, rather than offsetting remove the hit triangle through a mask, dynamic meshes, etc, before running the ray check again.

[aluo-x]

As far as I can tell, not when the triangles are large. Individual rays fluctuate between runs (the id it hits changes).

I've attempted the following to get second hit working:

• Changing the origin (hit into the bug where python-embree mutates the array each time)
• Changing tnear = tfar * ray_dist + ray_orig (hit the bug where ray_dist sometimes is on the order of 10^30 even for successful intersects)
• Currently attempting to use your code, compute the origin in python each time... Will check back in a few hours

[aluo-x]

I have it working now. For my use case the results appear to be correct with embree3 (when they were wrong with embree2).

[aluo-x]

Just to note, I'm not certain this solution gives points that are on the mesh.

I'm comparing this purely in my use case where I retrieve all hits along a ray. When using embree2 and the trimesh ray-offset, there were multiple holes in the ray-traced solution (very obvious when performing visualization as sudden gaps in the z-dist along the ray).

I no longer see this issue with embree3 + python-embree. The trick is to use scene.set_flags(4) which activates ROBUST mode. The only reliable information that embree returns is the triangle index, and the tfar and uv are both changing between runs (by orders of magnitude).

So the current pipeline is as follows:
Download embree 3.12.1 for linux.
You will need the embree3 wrapper by sampotter, I installed it on Linux with GCC 10 as my compiler. Set the correct include_dirs and library_dirs for the embree location in setup.py.

1. Create a TrimeshShapeModel object - see link here -- or my modified version here
2. Modify def _make_scene(self):, add a line scene.set_flags(4), this activates ROBUST mode
3. The TrimeshShapeModel can be set once and reused, np.copy is not necessary on the ray origin and direction (was another unrelated bug). I use the exact same ray origin and direction each time, but simply change the tnear value.
4. Override tnear with a list of offsets, I use (1 + 1e-5) as my tnear scaling along the ray direction. I observed duplicate intersections with smaller values. Depending on your use case, an even higher value (1+1e-4) or (1+1e-3) may work better.
5. I call intersect with the same exact origin & direction each time, I simply filter the output with a updated copy of "bad rays". "bad rays" are rays which have up until this point returned at least one -1 as the triangle intersection index. The number of "bad rays" can never decrease. For "bad rays", I set the tnear to a high value - 100 in my case
6. Use the Trimesh planes_lines function to compute the actual intersection. After using the bad rays array to filter the triangles, filter again with the valid output of plane_lines
7. I used np.float32 for tnear, ray origin, and ray direction

[sinAshish]

Is this issue solved? or the gist by @LMescheder is better? @mikedh