0

If I have a square based pyramid of which I know every point (all 5), what is the algorithm to determine if a given point is inside that pyramid or not? I saw a similar question here Algorithm for checking if 3D point inside convex polyhedron (square pyramid) and here How to check whether the point is in the tetrahedron or not?. The answers given there did not really help and the math did not work out for me. I have seen approaches with creating a normal to each plane (4 triangles + 1 square) and checking on which side of the plane the point lies... however I do not know the correct math for that.

If someone could help me with the math/algorithm for that (preferably with explanation), it would really help.

visualization

Example Input:

# Pyarmid Points
A = np.array([2, 2, 4]) # pyramid top
B = np.array([1, 1, 1]) 
C = np.array([1, 3, 3])
D = np.array([3, 3, 1])
E = np.array([3, 1, 1])

# Example points
point_inside = np.array([2, 2, 2])
point_outside = np.array([3, 3, 3])

Example Expectation:

def check_point_in_pyramid(A, B, C, D, E, Point):

    # math...
    
    if point_inside_pyramid:

        return true


    return false

So for point_inside the function return would be true and for the point_outside it would be false.

Steve Piercy
  • 13,693
  • 1
  • 44
  • 57
Nils
  • 3
  • 5

2 Answers2

1

I am taking the right as positive and left as negative here:

As per your pyramid, it is defined by 5 faces and therefore can also be defined by 5 PLANES. A plane is simply a flat (without a volume; only area) surface. , defined by the equation ax + by + cz + d = 0. However, to answer your question we only need the plane's normal, which is given by the three coeff's of the equation above N = (a, b, c).

Assuming we used plane ACD from your image. We can also find the plane's normal by taking the cross product like this: (D - A) x (C - A), which gives us the normal, N.

Now we also need a direction vector, that is pointing from the 3D point, which we are testing, towards any point on the plane. To simplify, we can just use either one of the points: A, C or D (let's use A) as they are still considered to part of the plane.

Like this: directionVector = point_outside - A

Now we have something on the means like this:

Here

Now if we take the dot product of the the direction vector and normal (float value = directionVector . N) we will only get a numeric value.

However the sign of the numeric value tells us whether the point is inside or outside (Which side of the plane the 3D point is on).

So if the value > 0, the point is inside the pyramid (to the right side of the plane) and if value < 0, the point is outside the pyramid (to the left side of the plane).

Now if repeat this process with all the other 4 planes of the pyramid. If all of them return a positive value, then the point is indeed inside the pyramid. However, if even one of them returns a negative value, then the point is on the outside the pyramid.

Note 1: Make sure the planes' normals are facing out of the pyramid and not inside.

Note 2: If you reverse the positive and negative sides (if left as positive and right as negative now), make sure to reverse the inequalities too.

Dstarred
  • 321
  • 2
  • 10
0

Yes, you can solve this problem using normals.

To get normal, calculate cross product of two edges for every plane. Be careful wtih normal direction - use vertices in strict order to get either all external, or all internal normals.

Then check signs of dot product for all pairs of plane normal with vector (vertex in plane - point to check).

In case of internal point all dot products will have the same sign (or zero when point belongs to the plane).

Example for outside point - note distinct signs:

nabc = np.cross(B - A, C - A)
nacd = np.cross(C - A, D - A)
nade = np.cross(D - A, E - A)
naeb = np.cross(E - A, B - A)
nbcd = np.cross(D - B, C - B)

r = np.dot(point_outside - A, nabc)
print(r)
r = np.dot(point_outside - A, nacd)
print(r)
r = np.dot(point_outside - A, nade)
print(r)
r = np.dot(point_outside - A, naeb)
print(r)
r = np.dot(point_outside - B, nbcd)
print(r)

8
-4
-4
8
8
MBo
  • 77,366
  • 5
  • 53
  • 86