How to check whether a point is inside a mesh

[baranaydogan]

I wonder if there is a way in Geogram to query whether a point is inside a mesh or not.

I noticed that under third_party/PoissonRecon there is an OctNode class, which has an "isInside" member function. But I am not sure whether this is meant for this task. And if so, there is also no documentation/comments nor an example for how to use it.

Any suggestions?

[BrunoLevy]

Hi, it depends whether you need performance or not.
If you do not need performance, just launch a ray from your point, and count the number of intersections with the mesh, if the number is odd, then you are inside.
If the ray passes exactly through a vertex or an edge, then you need to take additional care ... or if you do not want to, just throw another ray in a random direction (until no degenerate configuration is encountered).
Something very similar to that is done here

Now if you need performance (if you have many points to test), you may use an Axis-Aligned Bbox tree to accelerate ray-triangle intersections.

[baranaydogan]

Yes, performance is actually crucial here. I will need to query points in the order of millions (or maybe even billions in some cases).

I went through the AABB class. Is it so that, your suggestion is to first initialize MeshFacetsAABB, then call, ray_all_intersections with a function that adds the number of intersections? If any of the intersections occur on a vertex/edge (checked using the barycentric coordinates(?), throw new rays until no degeneracy occurs. If number of intersections is odd then the point is inside.

If this is the most efficient way do this with Geogram, it will be useful to have a "bool isInside(const vec3& p)" member function under MeshFacetsAABB in the future. Forgive the shameless request 😄

[BrunoLevy]

There is an important question: do you need the result to be exact or not ? I mean, is it acceptable to have some kind of "z fighting" near the boundary ?

If you do not need exactness, then you can simply use ray_nearest_intersection() (and that is the end of the story, if it returns false, then you are outside). You may use as you suggested the barycentric coordinates returned in the Intersection to detect whether the intersection is very near a vertex or edge (keep a small epsilon).

If you need exact intersections, then it is a completely different story: you can no longer use ray_all_intersections() (because it computes the invert of the ray direction coordinates using floating point coordinates which is not exact), so you will need to re-implement an exact traversal on your own.

A key component is this function that tells whether there is an intersection between a ray and a triangle (and whether this intersection is degenerate). You can reuse it as is (it is a template, I use it with exact points, but it will also work with vec3. Then it uses the exact implementation of the used geometric predicates for vec3).

If you have millions / billions points (but a rather small number of triangles in your surface, like 100K), there is another option:

• tetrahedralize your mesh here
• use MeshCellsAABB, function containing_tet (if it returns -1 then you are outside the mesh)
  This one is exact and simple to use (I'd suggest to start from there, of course, provided that your surface mesh does not have zillions triangles !)

If your points are organized on a grid, then it is also possible to "rasterize" the mesh in the grid, which is much much faster !

[baranaydogan]

Thanks! Great suggestions!

I would say that rather than "exactness", "robustness" is more important in this case, e.g., if I know that errors within an epsilon near the boundary can happen, that is fine; but a few flipped triangles should not cause an issue.

I should have mentioned earlier that I am working with brain images and surface meshes, and many times they come with non-manifold vertices/triangles, and many times they can be hard to fix despite the mesh repair routines in Geogram. To give a concrete example for an application, consider the signed Euclidean distance transform on an image grid. (This is not the only place I need this isInside function but it tells about the importance of robustness I think.)

Until now, I have been using a combination of approaches for checking "isInside" (that includes a form a rasterization of the mesh into voxels). I have found that the approach in libigl's fast winding number is quite fast and robust to degeneracies on the mesh too. But this is the only place I am using libigl, so it would make my software lighter if I can find a Geogram replacement.

I think, overall, the tetrahedralization idea can be more suitable for my needs but there I hit a license issue with tetgen. I was not able to find mg-tetra, which is mentioned in your code though, maybe it has a more permissible license(?) Can you share a link or clarify how to compile Geogram with mg-tetra?

[BrunoLevy]

The "voxelization" is really the way to go, that is, outer loop is "for each triangle" instead of "for each voxel"). It is much much faster.

Why not re-implementing "fast winding numbers" ? I think it is not too difficult, and the method is fully explained here. I think it is not too much code.

Another interesting list of references here

About the tetrahedralization approach (but re-implementing fast winding number is probably better)

• regarding mg-tetra, it is part of a commecial product, MeshGems
• there is a new tet mesher with GPL / LGPL license here (did not test it)
• then I'd recommend to do it like that:

SignedDistance:
For each tetrahedron t
   Get t's bbox
      for x = xmin to xmax
         for y = ymin to ymax
             for z = zmin to zmax
                  if(x,y,z is inside t)
                      setvoxel(x,y,z,distanceToSurface(x,y,z))
   for x = 0 to NX-1
      for y = 0 to NY-1
         for z = 0 to NZ-1
              if voxel(x,y,z) was not set
                   setvoxel(x,y,z,-distanceToSurface(x,y,z))

To compute distanceToSurface(), you can use the AABB and it will work. To make it faster, you can compute the (unsigned) distance transform in the grid (there is a super fast algoritm to do that in 6 "sweepings" in the grid, I do not remember the reference but it will be easy to find with a few google queries).

To test whether a point p is inside a tet (p1,p2,p3,p4), compute the four signed volumes of (p,p2,p3,p4), (p1, p, p3, p4), (p1, p2, p, p4), (p1, p2, p3, p) and compare their signs. If they are all the same, then p is inside (p1,p2,p3,p4).
If performance is an issue, you can:

• rasterize the tetrahedra in parallel
• compute the signed volumes in an incremental manner (by factoring out x,y,z from the expressions of the determinants, and compute the deltas used to increment them when incrementing x,y and z)

[BrunoLevy]

Oh I did not remember that I had this code (I mean AABB for computing distance to surface + tet for inside/outside test) ready in Graphite.
See here
(it is not optimized, but it may be a starting point)

[baranaydogan]

Great! I think the function in Graphite is indeed doing all that is needed for me.

If more efficiency is needed, one can look into the voxelization based idea with a fast winding number implementation too (unfortunately, I don't have the time to implement it myself now, even though that would be fun!)

Thanks again for your comments.