How to find if a line segment/ray intersects a plane in the 3-D space using the vertices of the plane and line

Question

I have the four corner vertices of a plane and vertices of the start and end point of a line segment. I want to be able to check if the line intersects the plane in the 3-D Space.

I was initially using this https://www.youtube.com/watch?v=_W3aVWsMp14 video for my math by converting it into parametrics and then solving.

However, i'm a little lost on how to actually solve and format the equation.

import numpy as np
from sympy import Point3D, Plane, Line3D
import sympy as sp


def find_intersection_point(line_vertices, plane_vertices):

    t, tt = sp.symbols('t tt')

    line_start = np.array(line_vertices[0])
    line_end = np.array(line_vertices[1])
    line_direction = line_end - line_start

    line_x = line_start[0]*t + line_direction[0]
    line_y = line_start[1]*t + line_direction[1]
    line_z = line_start[2]*t + line_direction[2]

    print(line_x)
    print(line_y)
    print(line_z)

    # Find two vectors in the plane
    vector1 = np.subtract(plane_vertices[1], plane_vertices[0])
    vector2 = np.subtract(plane_vertices[2], plane_vertices[0])

    # Find the normal vector
    normal_vector = np.cross(vector1, vector2)
    print("normal vector: ", normal_vector)

    # normal_vector[0]*(line_x-plane_vertices)



def get_line_vetices(line):
    line_vertices = [[0 for i in range(3)] for j in range(2)]
    line_vertices = np.empty((2, 3))
    for i in range(2):
        line_vertices[i][0] = line.Shape.Vertexes[i].X
        line_vertices[i][1] = line.Shape.Vertexes[i].Y
        line_vertices[i][2] = line.Shape.Vertexes[i].Z
    return line_vertices


def get_plane_vertices(plane):
    plane_vertices = [[0 for i in range(3)] for j in range(4)]
    plane_vertices = np.empty((4, 3))
    for i in range(4):
        plane_vertices[i][0] = plane.Shape.Vertexes[i].X
        plane_vertices[i][1] = plane.Shape.Vertexes[i].Y
        plane_vertices[i][2] = plane.Shape.Vertexes[i].Z
    return plane_vertices


# Define the vertices of the plane as [x, y, z]
plane_vertices = [
    [-26.551074, 29.194627000000004, 0.0],
    [45.091893, 29.194627000000004, 0.0],
    [45.091893, -34.722119000000006, 0.0],
    [-26.551074, -34.722119000000006, 0.0]
]

plane_vertices = np.array(plane_vertices)

# Define the start and end points of the line segment
line_vertices = [
    [-112.944015, -52.98405, 0.0],
    [136.40151699999996, 39.730347999999985, 0.0]
]

line_vertices = np.array(line_vertices)

find_intersection_point(line_vertices, plane_vertices)

This code down here was another method that i was experimenting with but didn't end up working out.

# plane Points
a1 = Point3D(-26.551074, 29.194627000000004, 0.0)
a2 = Point3D(45.091893, 29.194627000000004, 0.0)
a3 = Point3D(45.091893, -34.722119000000006, 0.0)
# line Points
p0 = Point3D(-112.944015, -52.98405, 0.0)  # point in line
line_direction = np.subtract(line_vertices[0], line_vertices[1])

print(line_direction)

# create plane and line
plane = Plane(a1, a2, a3)

line = Line3D(p0, direction_ratio=line_direction)


print(f"plane equation: {plane.equation()}")
print(f"line equation: {line.equation()}")

# find intersection:
# Find intersection
intr = plane.intersection(line)
if intr:
    intersection = intr[0]
    print("It intersects")
    print(f"Intersection: {intersection}")
else:
    print("No intersection found.")

Answer 1

Having three plane points a1,a2,a3 and line points p0, p1, you can make parametric equation, using vectors

 a12 = a2 - a1
 a13 = a3 - a1
 p01 = p1 - p0

for plane:

 x = a1.x + u * a12.x  + v * a13.x
 y = a1.y + u * a12.y  + v * a13.y
 z = a1.z + u * a12.z  + v * a13.z

for line:

 x = p0.x + t * p01.x
 y = p0.y + t * p01.y
 z = p0.z + t * p01.z

Now put line equations into the plane ones:

 a1.x + u * a12.x  + v * a13.x - p0.x - t * p01.x = 0
 a1.y + u * a12.y  + v * a13.y - p0.y - t * p01.y = 0
 a1.z + u * a12.z  + v * a13.z - p0.z - t * p01.z = 0

and solve this linear equation system for unknowns u, v, t

If system has zero discriminant, then this line lies in the plane or it is parallel to the plane.

Otherwise you have parameters for intersection point and can calculate it.

If you have a ray, then check that parameter t >= 0
If you have line segment, then check that t parameter is in range 0..1
If you need some part of the plane - check u,v parameters. Say for rectangle (more general - parallelogram) defined by a1,a2,a3 points (and implicit four vertex) parameters should be in range 0..1

Example with sympy

import sympy as sp
u,v, t = sp.symbols('u, v, t')

a1x,a1y,a1z = 1,0,0
a12x,a12y,a12z = -1,1,0
a13x,a13y,a13z = -1,0,1
p0x,p0y,p0z = 0,0,0
p01x,p01y,p01z = 2,2,2

eq1 = sp.Eq(a1x + u * a12x  + v * a13x - p0x - t * p01x, 0)
eq2 = sp.Eq(a1y + u * a12y  + v * a13y - p0y - t * p01y, 0)
eq3 = sp.Eq(a1z + u * a12z  + v * a13z - p0y - t * p01y, 0)
ans = sp.solve((eq1, eq2, eq3), (u, v, t))
print(ans)

>>>
{u: 1/3, v: 1/3, t: 1/6}

P.S. Made formulas with sympy. At first get denominanor (should be the same everywhere) and check for zero value.

import sympy
from sympy.abc import u, v, t
a1x,a1y,a1z =sympy.symbols('a1x,a1y,a1z')
a12x,a12y,a12z = sympy.symbols('a12x,a12y,a12z')
a13x,a13y,a13z = sympy.symbols('a13x,a13y,a13z')
p0x,p0y,p0z = sympy.symbols('p0x,p0y,p0z')
p01x,p01y,p01z = sympy.symbols('p01x,p01y,p01z')

print(sympy.solve([a1x + u * a12x  + v * a13x - p0x - t * p01x,
             a1y + u * a12y  + v * a13y - p0y - t * p01y,
             a1z + u * a12z  + v * a13z - p0y - t * p01y],
             u, v, t, dict=True))

[{u: (-p01y*(a13y - a13z)*(a1x - p0x) + (a1y - p0y)*(a13x*p01y - a13z*p01x) - (a1z - p0y)*(a13x*p01y - a13y*p01x))/(a12x*a13y*p01y - a12x*a13z*p01y - a12y*a13x*p01y + a12y*a13z*p01x + a12z*a13x*p01y - a12z*a13y*p01x), 
v: (p01y*(a12y - a12z)*(a1x - p0x) - (a1y - p0y)*(a12x*p01y - a12z*p01x) + (a1z - p0y)*(a12x*p01y - a12y*p01x))/(a12x*a13y*p01y - a12x*a13z*p01y - a12y*a13x*p01y + a12y*a13z*p01x + a12z*a13x*p01y - a12z*a13y*p01x), 
t: ((a1x - p0x)*(a12y*a13z - a12z*a13y) - (a1y - p0y)*(a12x*a13z - a12z*a13x) + (a1z - p0y)*(a12x*a13y - a12y*a13x))/(a12x*a13y*p01y - a12x*a13z*p01y - a12y*a13x*p01y + a12y*a13z*p01x + a12z*a13x*p01y - a12z*a13y*p01x)}]