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.
