Algorithm for determining whether a point is inside a 3D mesh

Question

What is a fast algorithm for determining whether or not a point is inside a 3D mesh? For simplicity you can assume the mesh is all triangles and has no holes.

What I know so far is that one popular way of determining whether or not a ray has crossed a mesh is to count the number of ray/triangle intersections. It has to be fast because I am using it for a haptic medical simulation. So I cannot test all of the triangles for ray intersection. I need some kind of hashing or tree data structure to store the triangles in to help determine which triangle are relevant.

Also, I know that if I have any arbitrary 2D projection of the vertices, a simple point/triangle intersection test is all necessary. However, I'd still need to know which triangles are relevant and, in addition, which triangles lie in front of a the point and only test those triangles.

Solution 1

I solved my own problem. Basically, I take an arbitrary 2D projection (throw out one of the coordinates), and hash the AABBs (Axis Aligned Bounding Boxes) of the triangles to a 2D array. (A set of 3D cubes as mentioned by titus is overkill, as it only gives you a constant factor speedup.) Use the 2D array and the 2D projection of the point you are testing to get a small set of triangles, which you do a 3D ray/triangle intersection test on (see Intersections of Rays, Segments, Planes and Triangles in 3D) and count the number of triangles the ray intersection where the z-coordinate (the coordinate thrown out) is greater than the z-coordinate of the point. An even number of intersections means it is outside the mesh. An odd number of intersections means it is inside the mesh. This method is not only fast, but very easy to implement (which is exactly what I was looking for).

Solution 2

This is algorithm is efficient only if you have many queries to justify the time for constructing the data structure.

Divide the space into cubes of equal size (we'll figure out the size later). For each cube know which triangles has at least a point in it. Discard the cubes that don't contain anything. Do a ray casting algorithm as presented on wikipedia, but instead o testing if the line intersects each triangle, get all the cubes that intersect with the line, and then do ray casting only with the triangles in these cubes. Watch out not to test the same triangle more than one time because it is present in two cubes.
Finding the proper cube size is tricky, it shouldn't be neither to big or too small. It can only be found by trial and error. Let's say number of cubes is c and number of triangles is t.
The mean number of triangles in a cube is t/c
k is mean number of cubes that intersect the ray
line-cube intersections + line-triangle intersection in those cubes has to be minimal
c+k*t/c=minimal => c=sqrt(t*k)
You'll have to test out values for the size of the cubes until c=sqrt(t*k) is true
A good starting guess for the size of the cube would be sqrt(mesh width)
To have some perspective, for 1M triangles you'll test on the order of 1k intersections

Comments

[rwols]

"Do a ray casting algorithm", in which direction?

[rwols]

If i'm understanding correctly, the direction for the ray test is the direction of the chosen 2D projection?

[gnzlbg]

What do you do when the ray intersects a vertex or an edge of a triangle? Don't you get the wrong number of intersections?

[Franz]

If you precisely hit an edge or vertex, you would have to check the surface normals (iactually the sign of their value in the direction component) of all participating triangles. More mundane solutions could be choosing another coordinate direction (if you have the according hashtable available) or adding small distortions to the original point and doing a majority vote.

[dgrat]

I implemented something similar. There are three reference implementations of different rasterizers. github.com/3DStuff/voxelizer

[Willem Van Onsem]

@rwols: any direction, as long as the ray goes through an odd number of triangles, it is in the mesh.

[alpha027]

What is the complexity of the algorithm ?