Check if 3D Point is inside square based pyramid

Question

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 Check if a point is inside a rectangular shaped area (3D)? and here How do I find if a point exists in a3D solid?. I tried to change the math a bit but didn't work out for me (maybe because I don't understand some parts correctly).

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.

visualization

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

Python "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])

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.

I do not need the python code, I just added it for clarification purposes what I need. The math is enough.

Edit after I read the answers:
I posted this question on stackoverflow too. Below are more deep level explanations and as far as I could tell a different approach to some extent. Definitely interesting. On stackoverflow are more lightweight answer for those looking. https://stackoverflow.com/questions/68641598/check-if-a-3d-point-is-in-a-square-based-pyramid-or-not

Answer 1

All convex polyhedra can be described as the intersection of four or more half-spaces. We describe a half-space using a plane.

If $\vec{n} = (n_x, n_y, n_z)$ is the normal of the plane, and $d$ is the signed distance from the plane to origin, then point $\vec{p} = (x, y, z)$ is within the half-space if and only if $$\vec{p} \cdot \vec{n} \ge d \tag{1a}\label{BtV1a}$$ which is equivalent to $$x n_x + y n_y + z n_z \ge d \tag{1b}\label{BtV1b}$$ in Cartesian coordinates.

To define a plane or a half-space ($\vec{n}$ and $d$), you can use three points (three consecutive vertices of a polyhedron side polygon): $\vec{v}_1 = (x_1, y_1, z_1)$, $\vec{v}_2 = (x_2, y_2, z_2)$, $\vec{v}_3 = (x_3, y_3, z_3)$: $$\left\lbrace \, \begin{aligned} \vec{n} & = (v_1 - v_2) \times (v_3 - v_2) \\ d & = \vec{n} \cdot \vec{v}_1 = \vec{n} \cdot \vec{v}_2 = \vec{n} \cdot \vec{v}_3 \\ \end{aligned} \right. \tag{2a}\label{BtV2a}$$ where $\times$ is the vector cross product, and you can choose any of the three right sides for $d$ (because they are equivalent). This is equivalent to $$\left\lbrace \, \begin{aligned} n_x &= (y_1 - y_2)(z_3 - z_2) - (z_1 - z_2)(y_3 - y_2) \\ n_y &= (z_1 - z_2)(x_3 - x_2) - (x_1 - x_2)(z_3 - z_2) \\ n_z &= (x_1 - x_2)(y_3 - y_2) - (y_1 - y_2)(x_3 - x_2) \\ d &= n_x x_1 + n_y y_1 + n_z z_1 = n_x x_2 + n_y y_2 + n_z z_2 = n_x x_3 + n_y y_3 + n_z z_3 \\ \end{aligned} \right . \tag{2b}\label{BtV2b}$$ in Cartesian coordinates.

When working with floating-point numbers, you should use maximum of the three right sides for the signed distance $d$, $$d = \max\bigl( n_x x_1 + n_y y_1 + n_z z_1 ,\, n_x x_2 + n_y y_2 + n_z z_2 ,\, n_x x_3 + n_y y_3 + n_z z_3 \bigr)$$ to minimize rounding errors. (For the base of the pyramid, you'll want to use the fourth planar vertex here also, to ensure it too is considered "inside" the pyramid.)

Using the notation in $\eqref{BtV1a}$ and $\eqref{BtV1b}$, the plane normal points to "inside". Depending on the order of the three points above, you can get $\vec{n}$ that points "inside" or "outside". If you know you have point $\vec{p} = (x, y, z)$ "inside" the half-space, just calculate $$\vec{p} \cdot \vec{n} - d = x n_x + y n_y + z n_z - d \tag{3}\label{BtV3}$$ If this is negative, you need to negate both $\vec{n}$ and $d$ (i.e, negate all four scalars $n_x$, $n_y$, $n_z$, and $d$). Then, you can be sure $\vec{n}$ points "inside" the half-space.

In your program, you need to calculate the plane/half-space for each side of the rectangular pyramid; five in total. You can use either the centroid of the vertices (average of the five vertex coordinates) –– because the pyramid is convex, this is always inside it ––, or one of the vertices to verify your normals all point "inside".

Here is an example Python 3 implementation, pyramid.py:

# SPDX-License-Identifier: CC0-1.0
# -*- encoding: utf8 -*-

from math import sqrt

class Halfspace(tuple):
    """Simple 3D half-space/plane"""

    def __new__(cls, x, y, z, d):
        return tuple.__new__(cls, (float(x), float(y), float(z), float(d)))

    @property
    def normal(self):
        """Normal vector as a tuple"""
        return self[0:3]

    @property
    def x(self):
        """Normal vector x coordinate"""
        return self[0]

    @property
    def y(self):
        """Normal vector y coordinate"""
        return self[1]

    @property
    def z(self):
        """Normal vector z coordinate"""
        return self[2]

    @property
    def d(self):
        """Signed distance from origin"""
        return self[3]

    def contains(self, point):
        """Test if specified point is inside this halfspace"""
        return self[0]*point[0] + self[1]*point[1] + self[2]*point[2] >= self[3]

    def __contains__(self, point):
        """Returns True if specified point is inside this halfspace"""
        return self[0]*point[0] + self[1]*point[1] + self[2]*point[2] >= self[3]

    def distance(self, point):
        """Signed distance to plane, in units of normal vector length.
           This is positive inside the half-space, negative outside."""
        return self[0]*point[0] + self[1]*point[1] + self[2]*point[2] - self[3]

    @classmethod
    def fourpoint(cls, p, v1, v2, v3, *vmore):
        """Construct a halfspace from points v1, v2, v3, such that p is inside it"""

        ux, uy, uz = v1[0]-v2[0], v1[1]-v2[1], v1[2]-v2[2]
        vx, vy, vz = v3[0]-v2[0], v3[1]-v2[1], v3[2]-v2[2]
        nx, ny, nz = uy*vz-uz*vy, uz*vx-ux*vz, ux*vy-uy*vx

        # Use the max of the signed distance from origin, to minimize rounding errors
        d = max(nx*v1[0] + ny*v1[1] + nz*v1[2],
                nx*v2[0] + ny*v2[1] + nz*v2[2],
                nx*v3[0] + ny*v3[1] + nz*v3[2])

        # If additional plane vertices are specified, verify they are "inside" also.
        for v in vmore:
            d = max(d, nx*v[0] + ny*v[1] + nz*v[2])

        # Optional: we normalize (nx, ny, nz) to unit length here.
        #           This is only necessary if the .distance() method is used/needed.
        nn = sqrt(nx*nx + ny*ny + nz*nz)
        nx, ny, nz, d = nx/nn, ny/nn, nz/nn, d/nn

        # Check the normal vector direction, and correct (swap) it if necessary.
        if nx*p[0] + ny*p[1] + nz*p[2] < d:
            nx, ny, nz, d = -nx, -ny, -nz, -d

        return cls(nx, ny, nz, d)

class Pyramid(tuple):
    """Simple rectangular 3D pyramid"""

    def __new__(cls, p0, p1, p2, p3, p4):
        """Create a new pyramid with apex at p0, and base p1-p2-p3-p4"""

        # h0 is the base of the pyramid.  We assume it is planar.
        h0 = Halfspace.fourpoint(p0, p1,p2,p3,p4)

        # h1-h4 are the sides of the pyramid.
        h1 = Halfspace.fourpoint(p3, p1,p2,p0)
        h2 = Halfspace.fourpoint(p4, p2,p3,p0)
        h3 = Halfspace.fourpoint(p1, p3,p4,p0)
        h4 = Halfspace.fourpoint(p2, p4,p1,p0)

        return tuple.__new__(cls, (h0, h1, h2, h3, h4))

    def contains(self, point):
        """Check if point is inside the pyramid"""
        return (self[0].contains(point) and self[1].contains(point) and
                self[2].contains(point) and self[3].contains(point) and
                self[4].contains(point))

    def __contains__(self, point):
        """Returns True if point is inside this pyramid"""
        return (self[0].contains(point) and self[1].contains(point) and
                self[2].contains(point) and self[3].contains(point) and
                self[4].contains(point))

    def distance(self, point):
        """Return the distance between point and the surface of the pyramid,
           if the point is inside the pyramid.  Otherwise return -1 (outside)."""
        # Note: This requires Halfspace.fourpoint() normalizes the normal to unit length.
        d = min(self[0].distance(point), self[1].distance(point),
                self[2].distance(point), self[3].distance(point),
                self[4].distance(point))
        if d >= 0:
            return d
        else:
            return -1

Note that the vertices for the base must be listed in (cyclic) order, i.e. tracing the vertices in order as one travels the perimeter of the face. If you specify them in "butterfly" order, you'll get very odd results.

For example, to construct a pyramid P with base $(-1,0,0)$, $(0,-1,0)$, $(1,0,0)$, $(0,1,0)$, and apex at $(0,0,1)$, and to check if point $(0,0,0.5)$ is inside the pyramid, use (in the same directory as above pyramid.py)

import pyramid

P = Pyramid((0,0,1), (-1,0,0), (0,-1,0), (1,0,0), (0,1,0))
if (0,0,0.5) in P:
    print("Point (0, 0, 0.5) is inside pyramid P.")

In general, this half-space based test tends to be the most efficient test for convex polyhedra. (When there are many faces, it is often useful to define one or more spheres that are completely inside the polyhedron, and/or one (or more) that contain the polyhedron. This is because point-in-sphere test is even more trivial; the former allow shortcuts to "inside", and the latter to "outside".)

For this pyramid case, the test involves 15 multiplications, 10 additions, and 5 comparisons (between floating-point scalars), which is quite efficient.