3D线面相交

Question

如果给定一条直线（用向量或两个点表示），如何找到该直线与平面相交的点？我找到了很多相关资源，但是我无法理解那里的方程式（它们似乎不是标准的代数方程）。我需要一个方程式（不管有多长），可以用标准编程语言（我正在使用Java）解释。

- jt78

---

可能是 http://stackoverflow.com/questions/4382591/line-plane-intersection 的重复问题。 - Cobra_Fast

这不是重复的，4382591并没有询问一般情况。 - ideasman42

Answer 1

以下是一个Python示例，用于查找线和平面的交点。

其中平面可以是点和法向量，也可以是4D向量（法向式）， 下面提供了两种情况的代码。

还要注意，此函数计算表示点在直线上的位置的值（在下面的代码中称为fac）。 您可能也想返回这个值，因为从0到1的值相交线段-这对调用者可能有用。

代码注释中注意到的其他细节。

---

注意：此示例使用纯函数，没有任何依赖项-以便于移植到其他语言。 使用Vector数据类型和运算符重载，可以更加简洁（包含在下面的示例中）。

# 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,
    )

如果平面被定义为一个四维向量（标准形式），我们需要找到平面上的一个点，然后像之前一样计算交点（参见p_co的分配）。

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

为了进一步参考，这是从Blender中提取并适应Python的内容。 math_geom.c中的isect_line_plane_v3()

---

为了更加清晰，这里提供使用mathutils API的版本（可以根据需要修改为其他支持运算符重载的数学库）。

# 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

- ideasman42

---

最终，真正有效的交点计算方法出现了！天啊，我花费了很长时间努力修复接触点代码，感谢你先生。 - Dan Bechard

从实际意义上讲，这种方法没有为线条设定起点/终点，因此在理论上是无限的。但在实践中，它受到浮点精度的限制。 - ideasman42

Answer 2

这里有一个Java方法，用于找到一条直线和一个平面的交点。虽然还有一些向量方法没有包含在内，但它们的功能相当容易理解。

/**
 * 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.normalize()) == 0) {
        return null;
    }

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

- ZGorlock

---

注意：LineDirection 必须 被规范化。 - smirkingman

我搜索了很长时间才找到这个答案！很棒的代码示例，易于理解。 - m_power

@smirkingman 你说得对，谢谢！在我的代码中，它在传递到方法之前已经被规范化了，我已经编辑了帖子以包括这一点。 - ZGorlock

只是更新一下，因为这种方法并不适用于所有情况。拿一个点 {0,0,0} 和一条射线 {-1,0,0}，在 x=10 的平面 YZ 上使用所以平面上的点为 {10,0,0}，法向量可以是{-1,0,0}或者{1,0,0}，即使从平面投影的射线被抛弃，它也能检测到撞击。 - Franck

@Franck 很好的观察！这段代码只适用于线，而不是射线。 - ZGorlock

Answer 3

您需要考虑三种情况：

• 平面与直线平行，且直线不在平面上（无交点）
• 平面与直线不平行（有一个交点）
• 平面包含直线（直线在平面上的每一点上都有交点）

您可以将直线表示为参数形式，例如此处：

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

这个讲座的前几页也是如此处理平面：

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

如果平面的法向量垂直于沿着直线的方向，则存在一个边缘情况，需要查看其是否有交点或位于平面内。

否则，您将得到一个交点，可以解决它。

我知道这不是代码，但要获得强大的解决方案，您可能需要将其置于应用程序的上下文中。

编辑：以下是仅有一个交点的示例。假设您从第一个链接中的参数化方程开始：

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

参数t可以是任何值。满足这些方程的所有(x, y, z)（无限集）构成该直线。然后，如果您有一个平面的方程式，例如：

x + 2y + 2z = 5

（摘自这里）您可以将上述x，y和z的方程代入平面方程中，现在平面方程仅包含参数t。解出t，这是该线在平面内所在的特定t值。然后，您可以通过回到线方程并重新代入t来解出x，y和z。

- John

---

谢谢，这很有帮助，但我还是有点困惑。我知道如何得到直线的方程，那很简单。但从那里开始，我就卡住了。用通俗易懂的语言来说，一个平面的“方程”究竟意味着什么？一旦解决了两个方程，我该怎么计算交点呢？ - jt78

你对你的平面有什么了解？比如平面上的点？ - John

几乎任何东西，我想。我至少知道平面上的4个点，并且它是独立定义的，所以我应该能够找到几乎任何东西。 - jt78

Answer 4

使用numpy和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)

- TimSC

Answer 5

如果你有两个点p和q来定义一条直线，以及一个平面的笛卡尔通式表示ax+by+cz+d=0，你可以使用参数方法。

如果你需要这个用于编码目的，下面是javascript代码片段：

/**
* 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)));

- DNax

---

“var t = - ( apx + bpy + cpz + d ) / tDenom;” 这句话中，实际上应该减去 d 吗？因此... var t = - ( apx + bpy + cpz - d ) / tDenom; - Russell Okamoto

d只是独立项，没有其他特殊之处。在这个例子中，你可以看到平面y=10被表示为d=-10（第一个控制台日志）。如果还不清楚，我可以写下公式的推导/解决方案。 - DNax

如果直线包含在平面内呢？这个函数似乎只返回一个点或null，在这种情况下，我期望它返回直线本身（因为那才是与平面相交的）。那么当返回null时，可以得出结论直线从未与平面相交吗？ - Mark Miller

这个公式中d被取反的原因是：给定平面的法向量n和平面上的两个点p和q，我们知道：dot(n, (p-q)) = 0。这可以扩展为（使用p和q作为点向量）dot(n,p)-dot(n,q)=0，其中-dot(n,q)变成了d，平面方程变为：dot(n,p)+d=0，这缩短为标准的[n|d]隐式平面形式。这里的重要结论是d=-dot(n,q)。 - pmw1234

Answer 6

基于this Matlab代码（不包括交叉点检查），使用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

结果

[-3.90909091  1.18181818 -0.27272727]

我进行了图形化检查，似乎可以工作。

我认为这是对before分享链接更加健壮的实现。

- BBSysDyn

Answer 7

这个问题很老，但由于有更方便的解决方案，我认为它可能会对某些人有帮助。

在齐次坐标中，平面和直线的交点非常优雅，但假设您只想要解决方案：

有一个描述平面的4x1向量p，使得对于平面上的任何齐次点，p^T*x=0。接下来计算线L=ab^T-ba^T的普拉克坐标，其中a={point_1; 1}，b={point_2; 1}，都是在线上的4x1向量。

计算：x=L*p={x0，x1，x2，x3}

x_intersect=({x0，x1，x2}/x3)，如果x3为零，则欧几里得意义下不存在交点。

- midjji

Answer 8

这是mathutils API的修复版本。

# 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 = -(p_no * w) / dot
        return p0 + (u * fac)

    return None

MathUtils使用向量运算来进行浮点数舍入，因此该解决方案不够精确。

- RLhoste

Answer 9

仅仅是为了扩展ZGorlock的答案，我已经完成了三维向量的点积、加法和缩放。这些计算的参考资料分别是点积, 两个三维向量相加和缩放。注意：Vec3D只是一个自定义类，它有三个点：x、y和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);
}

- Muhammad Qumail

---

这显然是不正确的：3D向量的点积应该是 (v0.x * v1.x + v0.y * v1.y + v0.z * v1.z)。 - Depau