Check if a point is inside an arbitrary hexahedron

Question

I am working on a 3D finite element code, where i face the following problem:

If I take an arbitrary point (say x), how do I figure out, which element it is in?

This can be simplified to: How do I check if an arbitrary point (x) lies inside or outside of an (hexahedral) element?

What I already found:

• Limited to cubes: How to determine a point is inside or outside a cube?
• Limited to rectangular shapes: https://math.stackexchange.com/questions/1472049/check-if-a-point-is-inside-a-rectangular-shaped-area-3d

Contrary to the two approaches above, my problem does not assume right angles nor parallel faces.

Problem sketch:

Notation: (again: though the sketch shows a regular shape, our hexahedron is assumed to be of general shape)

• 8-node hexahedron topology, nodes: 0,..,7

• axis: r,s,t

                    t
                    |
           4--------|-------------7
          /|        |            /|
         / |        |           / |
        /  |        |          /  |
       /   |        |         /   |
      /    |        |        /    |
     /     |        |       /     |
    5----------------------6      |
    |      |        |      |      |
    |      |        o------|---------s
    |      |       /       |      |
    |      0------/--------|------3
    |     /      /         |     /
    |    /      /          |    /
    |   /      /           |   /
    |  /      /            |  /
    | /      r             | /
    |/                     |/
    1----------------------2
  

Data that we have available:

• coordinates of the nodes (vectors P0 to P7)
• coordinates of the point we want to check (lets say Px)

Additionaly we assume the nodes are ordered as sketched above.

My approach/solution so far:

1. calculate the surface (outward) normal vectors

   Use cross products, eg. for the r_pos_normal_vec (pointing out of the plane)

   r_pos_normvec = (P2-P1) x (P5-P1)

   and for the r_neg_normal_vec

   r_neg_normvec = (P4-P0) x (P3-P0)

   similarly for the s and t directions

2. check two opposite corner nodes (I chose node0 and node 6)

   • For node0

     • calculate vector from P0 to Px:

       P0x = Px - P0

     • calculate inner prodcut of P0x and surfaces adjacent to node 0

       <P0x, r_neg_normal_vec>

       <P0x, s_neg_normal_vec>

       <P0x, t_neg_normal_vec>

   • For node1

     • same scheme as for node 0, whereas P1 instead of P0 and the positive counterparts of the normal vectors are used

3. Iff all 6 (3 from node0 and 3 from node1) inner products result in negative value -> the point is inside the hexahedron.

Question:

I implemented the functionality described above in my code and ran some tests. It seems to work, from the math side I am quite confident.

Please discuss my approach, I am happy for any hints/clues/recommendations/bug fixes ... Is there some way to make this faster? Alternative solutions?

Note:

• To speed up the algorithm a box check can be done first:
  • Construct a rectangular box around the hexahedron:
  • Get the min and max values of the node coordinates in each direction.
  • If the point to check (x) is outside this (larger) box, it cannot be inside the hexahedron.

Answer 1

For any convex polyhedron, establish the implicit equations of the faces (f.i. plane by three points), of the form ax+by+cz+d=0.

When you plug the coordinates of a known point inside the volume (such as the center) in the expression ax+by+cz+d, you will get a set of signs. An arbitrary point is inside if it yields the same signs.

---

Update:

For maximum performance, you can consider also using an axis-aligned bounding box for quick rejection. This only makes sense if many of the points are outside. Make sur to use a shortcut evaluation so that early rejection can happen.

Note that a rejection test such as X<Xmin is nothing but the above sign test against the plane of equation X-Xmin=0.

Answer 2

I personnally prefer your method, however there also is a way to approach the problem if the hexahedral restricted to parallelepiped. So you can transfer the coordinate of P in the frame $(0; e_1; e_2; e_3)$ to $(P_0, P_0P_1,P_0P_3,P_0P_4)$. We call it $(a,b,c)$, then the point is in that parallelepiped if $a,b,c > 0 a \in [0,1], a+b+c \in [0,1]$.

Answer 3

Because you mentioned that you want to be able to handle arbitrary hexahedrons, I think that your process might be improved if you adjust your s, r, and t normals to account for having slightly warped faces. I would do this by making the following change to r normals (and similar for s and t):

r_pos_normvec = (P6-P1) x (P5-P2)

r_neg_normvec = (P7-P0) x (P4-P3)

This would be important for a case where you shifted node 6 towards node 7 (say 0.9xP6) and had a point at 0.95xP6. Without the warping correction, I believe you would erroneously determine the point as inside the hexahedron.

Answer 4

Here is a python example :

def point_is_in_hexa(point,centers,normals):
    vect=[]
    prod=[]
    for i in range(6):
        vect.append(point-centers[i])
    vect= np.array(vect)
    for i in range(6):
        prod.append(np.dot(vect[i],normals[i]))
    prod=np.array(prod)
    if all(prod <= 0):
        is_in_hexa = 1
    else:
       is_in_hexa = -1
    return is_in_hexa

https://github.com/fgomez03/hexacheck