3D Line-Plane Intersection

Here is a Python example which finds the intersection of a line and a plane.

Where the plane can be either a point and a normal, or a 4d vector (normal form),
In the examples below (code for both is provided).

Also note that this function calculates a value representing where the point is on the line, (called fac in the code below).
You may want to return this too, because values from 0 to 1 intersect the line segment – which may be useful for the caller.

Other details noted in the code-comments.

---

Note: This example uses pure functions, without any dependencies – to make it easy to move to other languages.
With a Vector data type and operator overloading, it can be more concise (included in example below).

# intersection function
def isect_line_plane_v3(p0, p1, p_co, p_no, epsilon=1e-6):
    """
    p0, p1: Define the line.
    p_co, p_no: define the plane:
        p_co Is a point on the plane (plane coordinate).
        p_no Is a normal vector defining the plane direction;
             (does not need to be normalized).

    Return a Vector or None (when the intersection can't be found).
    """

    u = sub_v3v3(p1, p0)
    dot = dot_v3v3(p_no, u)

    if abs(dot) > epsilon:
        # The factor of the point between p0 -> p1 (0 - 1)
        # if 'fac' is between (0 - 1) the point intersects with the segment.
        # Otherwise:
        #  < 0.0: behind p0.
        #  > 1.0: infront of p1.
        w = sub_v3v3(p0, p_co)
        fac = -dot_v3v3(p_no, w) / dot
        u = mul_v3_fl(u, fac)
        return add_v3v3(p0, u)

    # The segment is parallel to plane.
    return None

# ----------------------
# generic math functions

def add_v3v3(v0, v1):
    return (
        v0[0] + v1[0],
        v0[1] + v1[1],
        v0[2] + v1[2],
    )


def sub_v3v3(v0, v1):
    return (
        v0[0] - v1[0],
        v0[1] - v1[1],
        v0[2] - v1[2],
    )


def dot_v3v3(v0, v1):
    return (
        (v0[0] * v1[0]) +
        (v0[1] * v1[1]) +
        (v0[2] * v1[2])
    )


def len_squared_v3(v0):
    return dot_v3v3(v0, v0)


def mul_v3_fl(v0, f):
    return (
        v0[0] * f,
        v0[1] * f,
        v0[2] * f,
    )

If the plane is defined as a 4d vector (normal form), we need to find a point on the plane, then calculate the intersection as before (see p_co assignment).

def isect_line_plane_v3_4d(p0, p1, plane, epsilon=1e-6):
    u = sub_v3v3(p1, p0)
    dot = dot_v3v3(plane, u)

    if abs(dot) > epsilon:
        # Calculate a point on the plane
        # (divide can be omitted for unit hessian-normal form).
        p_co = mul_v3_fl(plane, -plane[3] / len_squared_v3(plane))

        w = sub_v3v3(p0, p_co)
        fac = -dot_v3v3(plane, w) / dot
        u = mul_v3_fl(u, fac)
        return add_v3v3(p0, u)

    return None

For further reference, this was taken from Blender and adapted to Python.
isect_line_plane_v3() in math_geom.c

---

For clarity, here are versions using the mathutils API (which can be modified for other math libraries with operator overloading).

# point-normal plane
def isect_line_plane_v3(p0, p1, p_co, p_no, epsilon=1e-6):
    u = p1 - p0
    dot = p_no * u
    if abs(dot) > epsilon:
        w = p0 - p_co
        fac = -(plane * w) / dot
        return p0 + (u * fac)

    return None


# normal-form plane
def isect_line_plane_v3_4d(p0, p1, plane, epsilon=1e-6):
    u = p1 - p0
    dot = plane.xyz * u
    if abs(dot) > epsilon:
        p_co = plane.xyz * (-plane[3] / plane.xyz.length_squared)

        w = p0 - p_co
        fac = -(plane * w) / dot
        return p0 + (u * fac)

    return None