Can someone explain to me in easier terms how to calculate if a 3d point in space is intersecting a plane.

Question

Given a point $P$ and a triangular Plane $T$ with three points $P_0$, $P_1$, and $P_2$, I want to find out if the $P$ intersects the $T$. I realize that there is a very similar question here, however I did not understand any of the answers on that page since I am fairly new to this type of math. The other resources that I read were illustrated on a 2D plane instead of 3D. Is there an easy method to calculate the intersection of a point on a 3d triangular Plane that can be translated to code(I am trying to use this in for a coding project I am working on)?

Answer 1

Given your points $\{P_1, P_2, P_3\}$, where $P_j=(x_j,y_j,z_j)$, the equation of the plane is found by solving a system of equations:

$$\begin{aligned} a x_1 + b y_1 + c z_1 = \rho \\ a x_2 + b y_2 + c z_2 = \rho \\ a x_3 + b y_3 + c z_3 = \rho \\ \end{aligned}$$

for $(a,b,c,\rho)$. To be a plane, at least one of $(a,b,c)$ must be non-zero. Once we have an equation for the plane, we can multiply each term by any nonzero constant, so the coefficients are only unique up to this constant.

The equation of a plane is then $$a x + b y + c z = \rho.$$ Given a point $\{x_0,y_0,z_0\}$, then it lies on that plane if and only if $$a x_0 + b y_0 + c z_0 = \rho.$$

Answer 2

There are two separate questions here:

1. How to calculate the distance between a point and a plane

2. How to find out if a point on a plane is within a triangle in that plane

Let's call the three triangle vertices $\vec{v}_0 = (x_0, y_0, z_0)$, $\vec{v}_1 = (x_1, y_1, z_1)$, and $\vec{v}_2 = (x_2, y_2, z_2)$, and the point $\vec{p} = (x_p, y_p, z_p)$.

The first question is answered by using the equation of a plane, $$\vec{n} \cdot \vec{p} = d$$ where $\vec{n}$ is the plane normal (perpendicular to the plane), and $d$ is the signed distance from the plane to origin in units of the normal vector length. (If $d = 0$, the plane goes through origin, $(0, 0, 0)$. If $d \gt 0$, the plane is in the same direction as the normal vector from origin, and if $d \lt 0$, the plane is in the opposite direction.)

We can easily calculate $\vec{n}$ and $d$ from the three triangle vertices: $$\begin{aligned} \vec{n} &= (\vec{v}_1 - \vec{v}_0) \times (\vec{v}_2 - \vec{v}_0) \\ d &= \vec{n} \cdot \vec{v}_0 = \vec{n} \cdot \vec{v}_1 = \vec{n} \cdot \vec{v}_2 \\ \end{aligned}$$ where $\times$ represents vector cross product, and $\cdot$ vector dot product. Because we can scale $\vec{n}$ freely to any nonzero length (it'll just change $d$, and the pair will still represent the same plane), it is very useful to scale $\vec{n}$ to unit length, $\lVert\vec{n}\rVert = 1$. (I like to denote such unit length vectors as $\hat{n}$.)

A function that calculates the unit normal vector $\hat{n}$ and the signed distance $d$ from three points can be written in e.g. simplified Python as follows:

from math import sqrt

def Plane(x0,y0,z0, x1,y1,z1, x2,y2,z2):
    # g = v1 - v0
    gx = x1 - x0
    gy = y1 - y0
    gz = z1 - z0

    # h = v2 - v0
    hx = x2 - x0
    hy = y2 - y0
    hz = z2 - z0

    # Cross product n = g × h
    nx = gy*hz - gz*hy
    ny = gz*hx - gx*hz
    nz = gx*hy - gy*hx

    # nn = ||n||, Euclidean length 
    nn = sqrt(nx*nx + ny*ny + nz*nz)
    if nn == 0:
        raise ValueError("Degenerate triangle: line or point!")

    # Scale n to unit length
    nx = nx / nn
    ny = ny / nn
    nz = nz / nn

    # Find d - just use any of the three vertices; say v0.
    d = nx*x0 + ny*y0 + nz*z0

    # Return nx, ny, nz, and d.
    return (nx, ny, nz, d)

Note that this identifies the plane, and only needs to be calculated once.

Essentially, given a point $\vec{p} = (x, y, z)$, the unit normal vector to the plane $\hat{n}$, and the signed distance $d$, $$h = \lvert \hat{n} \cdot \vec{p} - d \rvert$$ tells you the distance between the point and the plane.

When working with floating-point numbers in a computer program, we cannot test for exact equality, because the numbers themselves are not exact. We must include a margin of rounding error. This is typically done by using a machine epsilon, a value that represents the largest number (say, in a point coordinate) that we consider still "zero", when those rounding errors are taken into account.

So, to test for point-in-plane, you can use something like the following simplified Python code:

def point_in_plane(x,y,z, nx,ny,nz, d, epsilon=0.00001):
    return (abs(nx*x + ny*y + nz*z - d) <= epsilon)

which returns True if the point is within epsilon distance of the plane (in units of the unit normal length, so we do expect the plane normal $\hat{n}$ to be an unit vector here, as given by the previous function), and False otherwise.

That takes care of the first part.

---

If we assume the point has already been determined to be close enough to the plane, we can use the vector cross product again to check if the point is within any convex polygon (and triangles are definitely convex), given just the vertices.

You see, if we have two successive vertices on the convex polygon, $\vec{v}_i$ and $\vec{v}_{i+1}$, the cross product between the vector between the two and the vector to the point in question, $(\vec{v}_{i+1} - \vec{v}_i)\times(\vec{p} - \vec{v}_i)$, will point in the same direction as the plane normal when the three ($\vec{v}_i$, $\vec{v}_{i+1}$, and $\vec{p}$) are in counterclockwise order, and in the opposite direction when they are in clockwise order. (It will be a zero vector when the three are collinear.) If we use dot product on the result and the plane normal, we'll get a positive number or negative number depending on whether the point is counterclockwise or clockwise, with respect to the edge. So, ignoring any zero cases, if the cases are either all positive or all negative, but not a mix, we know the point is surrounded by the edges, and therefore inside (or on) an edge. If there is a mix, then we know the point must be outside the convex polygon.

This function can be written in simplified Python as follows:

def point_in_convex_polygon(p, v, n, epsilon=0.00001):

    # Number of counterclockwise and clockwise results
    nccw = 0
    ncw  = 0

    vnext = v[len(v)-1 ]  # Last vertex
    for i in range(0, len(v)):  # i = 0, ..., len(v)-1
        vprev = vnext
        vnext = v[i]

        gx = vnext[0] - vprev[0]
        gy = vnext[1] - vprev[1]
        gz = vnext[2] - vprev[2]

        hx = p[0] - vprev[0]
        hy = p[1] - vprev[1]
        hz = p[2] - vprev[2]

        # s = n.(g×h)
        s = n[0]*(gy*hz-gz*hy) + n[1]*(gz*hx-gx*hz) + n[2]*(gx*hy-gy*hx)

        if s > epsilon:
            nccw = nccw + 1
            if ncw > 0:
                return False
        elif s < -epsilon:
            ncw = ncw + 1
            if nccw > 0:
                return False

    return True

Technically, we compare a product of coordinates to epsilon here, so mathematically we should be using epsilon squared instead; but the operations we're performing are the ones that are producing the rounding error with floating-point numbers, so using the same epsilon actually makes more sense.

In a real world program, you would precompute $\hat{n}$ and $d$ whenever a triangle is moved or added, for efficiency reasons, but if speed is not an issue, you could implement the is-point-in-3D-triangle test in simplified Python as

def point_in_3d_triangle(p, v0, v1, v2, epsilon=0.0001):

    # Plane unit normal and signed distance from origin
    nx, ny, nz, d = Plane(v0[0], v0[1], v0[2],
                          v1[0], v1[1], v1[2],
                          v2[0], v2[1], v2[2])

    # Point on the plane?
    if not point_in_plane(p[0], p[1], p[2], nx, ny, nz, d, epsilon):
        return False

    # Point within the triangle?
    return point_in_convex_polygon(p, [v0, v1, v2], n, epsilon)

---

For triangles (and parallelograms), there is an even better test, via barycentric coordinates $(u, v)$, which correspond to point $$\vec{p} = \vec{v}_0 + u \bigr( \vec{v}_1 - \vec{v}_0 \bigr) + v \bigr( \vec{v}_2 - \vec{v}_0 \bigr) \tag{1a}\label{None1a}$$ i.e. $$\left\lbrace ~ \begin{aligned} x &= x_0 + u ( x_1 - x_0 ) + v ( x_2 - x_0 ) \\ y &= y_0 + u ( y_1 - y_0 ) + v ( y_2 - y_0 ) \\ z &= z_0 + u ( z_1 - z_0 ) + v ( z_2 - z_0 ) \\ \end{aligned} \right . \tag{1b}\label{None1b}$$ Now, $(u, v)$ are within the triangle if and only if $$\begin{array}{c} 0 \le u \le 1 \\ 0 \le v \le 1 \\ 0 \le u + v \le 1 \\ \end{array} \tag{2}\label{None2}$$ Thing is, we have three equations, but only two unknowns. Essentially, we need to exclude the axis with the smallest angle to the plane as irrelevant.

For each triangle, or after moving a triangle, calculate $$\begin{aligned} c_{xy} &= x_0 ( y_2 - y_1 ) + x_1 ( y_0 - y_2 ) + x_2 ( y_1 - y_0 ) \\ c_{xz} &= x_0 ( z_2 - z_1 ) + x_1 ( z_0 - z_2 ) + x_2 ( z_1 - z_0 ) \\ c_{yz} &= y_0 ( z_2 - z_1 ) + y_1 ( z_0 - z_2 ) + y_2 ( z_1 - z_0 ) \\ \end{aligned} \tag{3}\label{None3}$$ and pick the one with the largest magnitude. Then,

$$\begin{aligned} u &= \displaystyle \frac{(y_0 - y_2) x + (x_2 - x_0) y + (x_0 y_2 - x_2 y_0)}{c_{xy}} \\ ~ &= \displaystyle \frac{(z_0 - z_2) x + (x_2 - x_0) z + (x_0 z_2 - x_2 z_0)}{c_{xz}} \\ ~ &= \displaystyle \frac{(z_0 - z_2) y + (y_2 - y_0) z + (y_0 z_2 - z_2 y_0)}{c_{yz}} \\ \end{aligned} \tag{4a}\label{None4a}$$ and $$\begin{aligned} v &= \displaystyle \frac{(y_1 - y_0) x + (x_0 - x_1) y + (x_1 y_0 - x_0 y_1)}{c_{xy}} \\ ~ &= \displaystyle \frac{(z_1 - z_0) x + (x_0 - x_1) z + (x_1 z_0 - x_0 z_1)}{c_{xz}} \\ ~ &= \displaystyle \frac{(z_1 - z_0) y + (y_0 - y_1) z + (y_1 z_0 - y_0 z_1)}{c_{yz}} \\ \end{aligned} \tag{4b}\label{None4b}$$

However, since we only need to use one of the three above, and they're linear functions of the coordinates of the point $\vec{p} = (x, y, z)$, we can actually just use $$\left\lbrace ~ \begin{aligned} u &= \vec{p} \cdot \vec{u} + u_0 \\ v &= \vec{p} \cdot \vec{v} + v_0 \\ \end{aligned} \right . \quad \iff \quad \left\lbrace ~ \begin{aligned} u &= x u_x + y u_y + z u_z + u_0 \\ v &= x v_x + y v_y + z v_z + v_0 \\ \end{aligned} \right . \tag{5}\label{None5}$$ Ah-ha! This will be seriously efficient. (But hopefully at this point you can see why I had to write such a wall of text to explain: just seeing this last one would have kept at least myself just scratching my head with a Wha?, unless I understood how this came to be, and that this is the better method.)

What we need to do in a computer program, is whenever the triangle is changed, to calculate $c_{xy}$, $c_{xz}$, and $c_{yz}$.

If $\lvert c_{xy} \rvert \ge \lvert c_{xz} \rvert$ and $\lvert c_{xy} \rvert \ge \lvert c_{yz} \rvert$, define $$\begin{aligned} u_x &= \displaystyle \frac{y_0 - y_2}{c_{xy}} \\ u_y &= \displaystyle \frac{x_2 - x_0}{c_{xy}} \\ u_z &= 0 \\ u_0 &= \displaystyle \frac{x_0 y_2 - x_2 y_0}{c_{xy}} \\ \end{aligned} \quad \text{and} \quad \begin{aligned} v_x &= \displaystyle \frac{y_1 - y_0}{c_{xy}} \\ v_y &= \displaystyle \frac{x_0 - x_1}{c_{xy}} \\ v_z &= 0 \\ v_0 &= \displaystyle \frac{x_1 y_0 - x_0 y_1}{c_{xy}} \\ \end{aligned} \tag{6xy}\label{None6xy}$$ Otherwise, if $\lvert c_{xz} \rvert \ge \lvert c_{xy} \rvert$ and $\lvert c_{xz} \rvert \ge \lvert c_{yz} \rvert$, define $$\begin{aligned} u_x &= \displaystyle \frac{z_0 - z_2}{c_{xz}} \\ u_y &= 0 \\ u_z &= \displaystyle \frac{x_2 - x_0}{c_{xz}} \\ u_0 &= \displaystyle \frac{x_0 z_2 - x_2 z_0}{c_{xz}} \\ \end{aligned} \quad \text{and} \quad \begin{aligned} v_x &= \displaystyle \frac{z_1 - z_0}{c_{xz}} \\ v_y &= 0 \\ v_z &= \displaystyle \frac{x_0 - x_1}{c_{xz}} \\ v_0 &= \displaystyle \frac{x_1 z_0 - x_0 z_1}{c_{xz}} \\ \end{aligned} \tag{6xz}\label{None6xz}$$ Otherwise, $\lvert c_{yz} \rvert \ge \lvert c_{xy} \rvert$ and $\lvert c_{yz} \rvert \ge \lvert c_{xz} \rvert$, and you can define $$\begin{aligned} u_x &= 0 \\ u_y &= \displaystyle \frac{z_0 - z_2}{c_{yz}} \\ u_z &= \displaystyle \frac{y_2 - y_0}{c_{yz}} \\ u_0 &= \displaystyle \frac{y_0 z_2 - y_2 z_0}{c_{yz}} \\ \end{aligned} \quad \text{and} \quad \begin{aligned} v_x &= 0 \\ v_y &= \displaystyle \frac{z_1 - z_0}{c_{xz}} \\ v_z &= \displaystyle \frac{y_0 - y_1}{c_{xz}} \\ v_0 &= \displaystyle \frac{y_1 z_0 - y_0 z_1}{c_{xz}} \\ \end{aligned} \tag{6yz}\label{None6yz}$$

Here is a Public Domain (Creative Commons Zero) -licensed proper example in Python3 that implements a Triangle class (and a helper Vector class implementing basic 3D vector algebra) optimized for testing containmen:

# SPDX-License-Identifier: CC0-1.0
from math import sqrt

class Triangle(tuple):
    """3D Triangle class optimized for point-in-triangle testing"""

    def __new__(cls, a, b, c):
        a = Vector(a[0], a[1], a[2])
        b = Vector(b[0], b[1], b[2])
        c = Vector(c[0], c[1], c[2])

        # First, calculate the unit normal vector (cross product).
        n = ((b - a) ^ (c - a)).unit
        if not n:
            raise ValueError("Degenerate triangle")

        # Calculate the signed distance from origin (dot product).
        d = n | a

        # Calculate the three possible divisors.
        div_xy = a.x*(c.y-b.y) + b.x*(a.y-c.y) + c.x*(b.y-a.y)
        div_xz = a.x*(c.z-b.z) + b.x*(a.z-c.z) + c.x*(b.z-a.z)
        div_yz = a.y*(c.z-b.z) + b.y*(a.z-c.z) + c.y*(b.z-a.z)
        abs_xy = abs(div_xy)
        abs_xz = abs(div_xz)
        abs_yz = abs(div_yz)
        if abs_xy >= abs_xz and abs_xy >= abs_yz:
            # d_xy has the largest absolute value; using xy plane
            u_axis = Vector((a.y-c.y)/div_xy, (c.x-a.x)/div_xy, 0)
            v_axis = Vector((b.y-a.y)/div_xy, (a.x-b.x)/div_xy, 0)
            u_offset = (a.x*c.y - a.y*c.x)/div_xy
            v_offset = (a.y*b.x - a.x*b.y)/div_xy
        elif abs_xz >= abs_xy and abs_xz >= abs_yz:
            # d_xz has the largest absolute value; using xz plane
            u_axis = Vector((a.z-c.z)/div_xz, 0, (c.x-a.x)/div_xz)
            v_axis = Vector((b.z-a.z)/div_xz, 0, (a.x-b.x)/div_xz)
            u_offset = (a.x*c.z - a.z*c.x)/div_xz
            v_offset = (a.z*b.x - a.x*b.z)/div_xz
        else:
            # d_yz has the largest absolute value; using yz plane
            u_axis = Vector(0, (a.z-c.z)/div_yz, (c.y-a.y)/div_yz)
            v_axis = Vector(0, (b.z-a.z)/div_yz, (a.y-b.y)/div_yz)
            u_offset = (a.y*c.z - a.z*c.y)/div_yz
            v_offset = (a.z*b.y - a.y*b.z)/div_yz

        return tuple.__new__(cls, (a, b, c, n, d, u_axis, u_offset, v_axis, v_offset))

    @property
    def vertex(self):
        """Tuple of all three vertices"""
        return (self[0], self[1], self[2])

    @property
    def normal(self):
        """Triangle plane unit normal vector"""
        return self[3]

    @property
    def d(self):
        """Triangle plane signed distance from origin"""
        return self[4]

    @property
    def plane(self):
        """Triangle plane (nx, ny, nz, d)"""
        return self[3][0], self[3][1], self[3][2], self[4]

    @property
    def u_axis(self):
        """U axis vector"""
        return self[5]

    @property
    def u_offset(self):
        """U axis offset"""
        return self[6]

    @property
    def v_axis(self):
        """V axis vector"""
        return self[7]

    @property
    def v_offset(self):
        """V axis offset"""
        return self[8]

    def distance_to_point(self, p):
        """Distance between point p and the triangle plane"""
        return abs((Vector(p[0], p[1], p[2]) | self.normal) - self.d)

    def contains(self, p, epsilon=0.0001):
        """Returns (u,v) coordinates of point p if it is inside this triangle, None otherwise."""
        # Make sure p is a Vector.
        p = Vector(p[0], p[1], p[2])

        # Outside the plane of the triangle?
        if abs((p | self.normal) - self.d) > epsilon:
            return None

        # U coordinate. We can use exact zero for u,v coordinates.
        u = (p | self.u_axis) + self.u_offset
        if u < 0 or u > 1:
            return None

        # V coordinate.
        v = (p | self.v_axis) + self.v_offset
        if v < 0 or v > 1:
            return None

        # This is a triangle, not a parallelogram; u+v <= 1.  (We know u+v >= 0.)
        if u + v > 1:
            return None

        # Yes, p is within this triangle.
        return (u, v)


class Vector(tuple):
    """3D Vector class with basic vector algebra support"""

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

    def __str__(self):
        return "(%s,%s,%s)" % (str(self[0]), str(self[1]), str(self[2]))

    def __bool__(self):
        """Nonzero vectors are True, zero vectors False."""
        return (self[0]*self[0] + self[1]*self[1] + self[2]*self[2] > 0.0)

    def __abs__(self):
        """abs(v): Absolute value is Euclidean length"""
        return sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])

    def __neg__(self):
        """-v: Vector negation"""
        return tuple.__new__(self.__class__, (-self[0], -self[1], -self[2]))

    def __pos__(self):
        """+v: Vector itself"""
        return self

    def __add__(self, other):
        """v+w: Vector addition (left side term)"""
        return tuple.__new__(self.__class__, (self[0]+other[0], self[1]+other[1], self[2]+other[2]))

    def __radd__(self, other):
        """w+v: Vector addition (right side term)"""
        return tuple.__new__(self.__class__, (other[0]+self[0], other[1]+self[1], other[2]+self[2]))

    def __sub__(self, other):
        """v-w: Vector subtraction (left side term)"""
        return tuple.__new__(self.__class__, (self[0]-other[0], self[1]-other[1], self[2]-other[2]))

    def __rsub__(self, other):
        """w-v: Vector subtraction (right side term)"""
        return tuple.__new__(self.__class__, (other[0]-self[0], other[1]-self[1], other[2]-self[2]))

    def __mul__(self, scalar):
        """v*n: Vector-scalar multiplication (left side term)"""
        if isinstance(scalar, (int, float)):
            return tuple.__new__(self.__class__, (self[0]*scalar, self[1]*scalar, self[2]*scalar))
        else:
            return NotImplemented

    def __rmul__(self, scalar):
        """n*v: Vector-scalar multiplication (right side term)"""
        if isinstance(scalar, (int, float)):
            return tuple.__new__(self.__class__, (scalar*self[0], scalar*self[1], scalar*self[2]))
        else:
            return NotImplemented

    def __truediv__(self, scalar):
        """v/n: Vector-scalar division (left side term)"""
        if isinstance(scalar, (int, float)):
            return tuple.__new__(self.__class__, (self[0]/scalar, self[1]/scalar, self[2]/scalar))
        else:
            return NotImplemented

    def __or__(self, other):
        """v|w: Vector dot product (left side term)"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return self[0]*other[0] + self[1]*other[1] + self[2]*other[2]
        else:
            return NotImplemented

    def __ror__(self, other):
        """w|v: Vector dot product (right side term)"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return other[0]*self[0] + other[1]*self[1] + other[2]*self[2]
        else:
            return NotImplemented

    def __xor__(self, other):
        """v^w: Vector cross product (left side term)"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return tuple.__new__(self.__class__, (self[1]*other[2] - self[2]*other[1],
                                                  self[2]*other[0] - self[0]*other[2],
                                                  self[0]*other[1] - self[1]*other[0]))
        else:
            return NotImplemented

    def __rxor__(self, other):
        """w^v: Vector dot product (right side term)"""
        if isinstance(other, (tuple, list)) and len(other) >= 3:
            return tuple.__new__(self.__class__, (other[1]*self[2] - other[2]*self[1],
                                                  other[2]*self[0] - other[0]*self[2],
                                                  other[0]*self[1] - other[1]*self[0]))
        else:
            return NotImplemented

    def __and__(self, notimplemented): return NotImplemented

    def __rand__(self, notimplemented): return NotImplemented

    def __rtruediv__(self, notimplemented): return NotImplemented

    @property
    def x(self):
        return self[0]

    @property
    def y(self):
        return self[1]

    @property
    def z(self):
        return self[2]

    @property
    def norm(self):
        """Absolute value or Euclidean length"""
        return sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])

    @property
    def normsqr(self):
        """v|v: Euclidean length squared"""
        return self[0]*self[0] + self[1]*self[1] + self[2]*self[2]

    @property
    def unit(self):
        """Vector scaled to unit length"""
        n = sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])
        return tuple.__new__(self.__class__, (self[0]/n, self[1]/n, self[2]/n))


if __name__ == '__main__':
    from random import Random
    from sys import exit

    rng = Random()

    # Coordinate range for the random points.
    Cmin = -9
    Cmax = +9

    N = 100
    P = 500
    epsilon = 0.0001
    for i in range(0, N):

        # Generate three random vertices, at least 0.01 distance away from each other.
        while True:
            v0 = Vector(rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax))
            v1 = Vector(rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax))
            v2 = Vector(rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax), rng.uniform(Cmin,Cmax))
            if (v1-v0).norm >= 0.01 and (v2-v0).norm >= 0.01 and (v2-v1).norm >= 0.01:
                break
        tri = Triangle(v0, v1, v2)

        # Unit normal vector to the triangle plane
        n = ((v1 - v0) ^ (v2 - v0)).unit

        for j in range(0, P):
            # Generate a random point (u, v) within a triangle
            while True:
                u = rng.uniform(0, 1)
                v = rng.uniform(0, 1)
                if u + v <= 1:
                    break
            # Calculate its position.
            p = v0 + u*(v1 - v0) + v*(v2 - v0)

            # Ask the Triangle class for the (u,v) coordinates.
            c = tri.contains(p)
            if c is None:
                print("Error: %s is at (%.3f,%.3f) inside the triangle %s, %s, %S, but .contains() failed." % (str(p), u, v, str(v0), str(v1), str(v2)))
                exit(1)

            # Find q definitely outside the triangle
            if Cmin < 0.1 or Cmax > 0.1:
                while True:
                    r = rng.uniform(Cmin, Cmax)
                    if r < -0.1 or r > 0.1:
                        break
                q = p + r*n

                c = tri.contains(q)
                if c is not None:
                    print("Error: %s is outside the triangle %s, %s, %s, but .contains returned %s." % (str(q), str(v0), str(v1), str(v2), str(c)))
                    exit(1)

    print("Tested %d random triangles and %d points inside and outside successfully." % (N, P))

If you save it as triangle.py, you can see its interface documentation using pydoc3 triangle in the same directory. You can run it to do self-checks with random triangles and points (both inside and outside). You can even include it in your own projects by putting it in the same directory, and adding from triangle import Vector, Triangle to your script.

After creating an instance of the Triangle class, you can use its .contains(point) member function to test if point is within the triangle or not. If not, it will return None, otherwise the $(u,v)$ Barycentric coordinates to point within the triangle.