3D Line-Plane Intersection

问题:

If given a line (represented by either a vector or two points on the line) how do I find the point at which the line intersects a plane? I've found loads of resources on this but I can't understand the equations there (they don't seem to be standard algebraic). I would like an equation (no matter how long) that can be interpreted by a standard programming language (I'm using Java).

回答1:

Here is a method in Java that finds the intersection between a line and a plane. There are vector methods that aren't included but their functions are pretty self explanatory.

/**
 * Determines the point of intersection between a plane defined by a point and a normal vector and a line defined by a point and a direction vector.
 *
 * @param planePoint    A point on the plane.
 * @param planeNormal   The normal vector of the plane.
 * @param linePoint     A point on the line.
 * @param lineDirection The direction vector of the line.
 * @return The point of intersection between the line and the plane, null if the line is parallel to the plane.
 */
public static Vector lineIntersection(Vector planePoint, Vector planeNormal, Vector linePoint, Vector lineDirection) {
    if (planeNormal.dot(lineDirection) == 0) {
        return null;
    }

    double t = (planeNormal.dot(planePoint) - planeNormal.dot(linePoint)) / planeNormal.dot(lineDirection);
    return linePoint.plus(lineDirection.scale(t));
}

回答2:

Heres 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)
    else:
        # 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)
    else:
        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)
    else:
        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)
    else:
        return None

回答3:

You'll need to consider three cases:

• Plane is parallel to line, and line does not lie in plane (no intersection)
• Plane is not parallel to line (one point of intersection)
• Plane contains the line (line intersects at every point on it)

You can express the line in paramaterized form, like here:

http://answers.yahoo.com/question/index?qid=20080830195656AA3aEBr

The first few pages of this lecture do the same for the plane:

http://math.mit.edu/classes/18.02/notes/lecture5compl-09.pdf

If the normal to the plane is perpendicular to the direction along the line, then you have an edge case and need to see whether it intersects at all, or lies within the plane.

Otherwise, you have one point of intersection, and can solve for it.

I know this isn't code but to get a robust solution you'll probably want to put this in the context of your application.

EDIT: Here's an example for which there's exactly one point of intersection. Say you start with the parameterized equations in the first link:

x = 5 - 13t
y = 5 - 11t
z = 5 - 8t

The parameter t can be anything. The (infinite) set of all (x, y, z) that satisfy these equations comprise the line. Then, if you have the equation for a plane, say:

x + 2y + 2z = 5

(taken from here) you can substitute the equations for x, y, and z above into the equation for the plane, which is now in only the parameter t. Solve for t. This is the particular value of t for that line that lies in the plane. Then you can solve for x, y, and z by going back up to the line equations and substituting t back in.

回答4:

Using numpy and python:

#Based on http://geomalgorithms.com/a05-_intersect-1.html
from __future__ import print_function
import numpy as np

epsilon=1e-6

#Define plane
planeNormal = np.array([0, 0, 1])
planePoint = np.array([0, 0, 5]) #Any point on the plane

#Define ray
rayDirection = np.array([0, -1, -1])
rayPoint = np.array([0, 0, 10]) #Any point along the ray

ndotu = planeNormal.dot(rayDirection) 

if abs(ndotu) < epsilon:
    print ("no intersection or line is within plane")

w = rayPoint - planePoint
si = -planeNormal.dot(w) / ndotu
Psi = w + si * rayDirection + planePoint

print ("intersection at", Psi)

回答5:

Based on this Matlab code (minus the checks for intersection), in Python

# n: normal vector of the Plane 
# V0: any point that belongs to the Plane 
# P0: end point 1 of the segment P0P1
# P1:  end point 2 of the segment P0P1
n = np.array([1., 1., 1.])
V0 = np.array([1., 1., -5.])
P0 = np.array([-5., 1., -1.])
P1 = np.array([1., 2., 3.])

w = P0 - V0;
u = P1-P0;
N = -np.dot(n,w);
D = np.dot(n,u)
sI = N / D
I = P0+ sI*u
print I

Result

[-3.90909091  1.18181818 -0.27272727]

I checked it graphically it seems to work,

This I believe is a more robust implementation of the link shared before

回答6:

This question is old but since there is such a much more convenient solution I figured it might help someone.

Plane and line intersections are quite elegant when expressed in homogeneous coordinates but lets assume you just want the solution:

There is a vector 4x1 p which describes the plane such that p^T*x =0 for any homogeneous point on the plane. Next compute the plucker coordinates for the line L=ab^T - ba^T where a = {point_1; 1}, b={point_2;1}, both 4x1 on the line

compute: x=L*p = {x0,x1,x2,x3}

x_intersect=({x0,x1,x2}/x3) where if x3 is zero there is no intersection in the euclidean sense.

回答7:

If you have two points p and q that define a line, and a plane in the general cartesian form ax+by+cz+d = 0, you can use the parametric method.

If you needed this for coding purposes, here's a javascript snippet:

/**
* findLinePlaneIntersectionCoords (to avoid requiring unnecessary instantiation)
* Given points p with px py pz and q that define a line, and the plane
* of formula ax+by+cz+d = 0, returns the intersection point or null if none.
*/
function findLinePlaneIntersectionCoords(px, py, pz, qx, qy, qz, a, b, c, d) {
    var tDenom = a*(qx-px) + b*(qy-py) + c*(qz-pz);
    if (tDenom == 0) return null;

    var t = - ( a*px + b*py + c*pz + d ) / tDenom;

    return {
        x: (px+t*(qx-px)),
        y: (py+t*(qy-py)),
        z: (pz+t*(qz-pz))
    };
}

// Example (plane at y = 10  and perpendicular line from the origin)
console.log(JSON.stringify(findLinePlaneIntersectionCoords(0,0,0,0,1,0,0,1,0,-10)));

// Example (no intersection, plane and line are parallel)
console.log(JSON.stringify(findLinePlaneIntersectionCoords(0,0,0,0,0,1,0,1,0,-10)));

回答8:

Just to expand on ZGorlock's answer, I have done the dot product, plus and scalse of 3D Vectors. The references for these calculations are Dot Product, Add two 3D vectors and Scaling. Note: Vec3D is just a custom class which has points: x, y and z.

/**
 * Determines the point of intersection between a plane defined by a point and a normal vector and a line defined by a point and a direction vector.
 *
 * @param planePoint    A point on the plane.
 * @param planeNormal   The normal vector of the plane.
 * @param linePoint     A point on the line.
 * @param lineDirection The direction vector of the line.
 * @return The point of intersection between the line and the plane, null if the line is parallel to the plane.
 */
public static Vec3D lineIntersection(Vec3D planePoint, Vec3D planeNormal, Vec3D linePoint, Vec3D lineDirection) {
    //ax × bx + ay × by
    int dot = (int) (planeNormal.x * lineDirection.x + planeNormal.y * lineDirection.y);
    if (dot == 0) {
        return null;
    }

    // Ref for dot product calculation: https://www.mathsisfun.com/algebra/vectors-dot-product.html
    int dot2 = (int) (planeNormal.x * planePoint.x + planeNormal.y * planePoint.y);
    int dot3 = (int) (planeNormal.x * linePoint.x + planeNormal.y * linePoint.y);
    int dot4 = (int) (planeNormal.x * lineDirection.x + planeNormal.y * lineDirection.y);

    double t = (dot2 - dot3) / dot4;

    float xs = (float) (lineDirection.x * t);
    float ys = (float) (lineDirection.y * t);
    float zs = (float) (lineDirection.z * t);
    Vec3D lineDirectionScale = new Vec3D( xs, ys, zs);

    float xa = (linePoint.x + lineDirectionScale.x);
    float ya = (linePoint.y + lineDirectionScale.y);
    float za = (linePoint.z + lineDirectionScale.z);

    return new Vec3D(xa, ya, za);
}