3D中直线和三角形的交点

Question

我有一条线和一个三角形，在3D空间的某个地方。换句话说，这个三角形有三个点([x,y,z]每个点)，这条线有两个点(也是[x,y,z]每个点)。

我需要找到一种方法，最好使用C++，来确定这条线是否穿过了这个三角形。与三角形平行且具有多个共同点的直线应被计算为“不相交”。

我已经编写了一些代码，但它并不起作用，而且即使可视化表示清楚地显示了相交，我仍然得到false。

ofVec3f P1, P2;
P1 = ray.s;
P2 = ray.s + ray.t;

ofVec3f p1, p2, p3;
p1 = face.getVertex(0);
p2 = face.getVertex(1);
p3 = face.getVertex(2);

ofVec3f v1 = p1 - p2;
ofVec3f v2 = p3 - p2;

float a, b, c, d;

a = v1.y * v2.z - v1.z * v2.y;
b = -(v1.x * v2.z - v1.z * v2.x);
c = v1.x * v2.y - v1.y * v2.x;
d = -(a * p1.x + b * p1.y + c * p1.z);

ofVec3f O = P1;
ofVec3f V = P2 - P1;

float t;

t = -(a * O.x + b * O.y + c * O.z + d) / (a * V.x + b * V.y + c * V.z);

ofVec3f p = O + V * t;

float xmin = std::min(P1.x, P2.x);
float ymin = std::min(P1.y, P2.y);
float zmin = std::min(P1.z, P2.z);

float xmax = std::max(P1.x, P2.x);
float ymax = std::max(P1.y, P2.y);
float zmax = std::max(P1.z, P2.z);


if (inside(p, xmin, xmax, ymin, ymax, zmin, zmax)) {
    *result = p.length();
    return true;
}
return false;

这里是inside()的定义：

bool primitive3d::inside(ofVec3f p, float xmin, float xmax, float ymin, float ymax, float zmin, float zmax) const {
    if (p.x >= xmin && p.x <= xmax && p.y >= ymin && p.y <= ymax && p.z >= zmin && p.z <= zmax)
        return true;

    return false;
}

- Kaito Kid

---

inside() 函数是做什么用的？ - Fureeish

抱歉，我编辑了我的问题以添加它。 - Kaito Kid

你是在寻找线段/三角形、直线/三角形还是射线/三角形的交点呢？从代码来看，我会说是射线/三角形，但当前的描述是直线/三角形。 - user3146587

Answer 1

1）如果你只是想知道直线是否与三角形相交（而不需要实际的交点）：

设p1，p2，p3为你的三角形。

在直线上选取两个点q1，q2，这两个点要在两个方向上离得很远。

令SignedVolume(a,b,c,d)表示四面体a,b,c,d的带符号体积。

如果SignedVolume(q1,p1,p2,p3)和SignedVolume(q2,p1,p2,p3)有不同的符号，并且SignedVolume(q1,q2,p1,p2)、SignedVolume(q1,q2,p2,p3)和SignedVolume(q1,q2,p3,p1)都有相同的符号，则存在交点。

SignedVolume(a,b,c,d) = (1.0/6.0)*dot(cross(b-a,c-a),d-a)

2) 如果您想要交集，当第1步测试通过时

以参数形式编写线的方程：p(t) = q1 + t*(q2-q1)

编写平面方程：dot(p-p1,N) = 0 其中

N = cross(p2-p1, p3-p1)

将p(t)注入平面方程中：dot(q1 + t*(q2-q1) - p1, N) = 0

展开得：dot(q1-p1,N) + t dot(q2-q1,N) = 0

推导得：t = -dot(q1-p1,N)/dot(q2-q1,N)

交点为：q1 + t*(q2-q1)

3）一种更加高效的算法

我们现在研究的是以下算法：

Möller和Trumbore，“快速，最小存储光线-三角形相交”，Journal of Graphics Tools，vol. 2，‎1997，p. 21–28 （也可参见：https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm）

该算法最终比1)和2)中的算法更简单（指令较少），但稍微复杂一些。让我们一步步推导。

符号：

• O = 光线起点,

• D = 光线方向向量,

• A,B,C = 三角形的顶点

光线上的任意点P可以表示为：P = O + tD

三角形上的任意点P可以表示为：P = A + uE1 + vE2 其中 E1 = B-A，E2 = C-A，u≥0，v≥0 且 (u+v)≤1

将P的两个表达式代入得：

O + tD = A + uE1 + vE2 

或者：

uE1 + vE2 -tD = O-A

矩阵形式:

            [u]
 [E1|E2|-D] [v] = O-A
            [t]

(其中[E1 | E2 | -D]是以E1，E2和-D作为其列的3x3矩阵)

使用Cramer公式解决：

   [a11 a12 a13][x1]   [b1]
   [a12 a22 a23][x2] = [b2]
   [a31 a32 a33][x3]   [b3]

提供：

       |b1 a12 a13|   |a11 a12 a13|
  x1 = |b2 a22 a23| / |a21 a22 a23|
       |b3 a32 a33|   |a31 a32 a33|

       |a11 b1 a13|   |a11 a12 a13|
  x2 = |a21 b2 a23| / |a21 a22 a23|
       |a31 b3 a33|   |a31 a32 a33|

       |a11 a12 b1|   |a11 a12 a13|
  x3 = |a21 a22 b2| / |a21 a22 a23|
       |a31 a32 b3|   |a31 a32 a33|

现在我们得到：

  u = (O-A,E2,-D) / (E1,E2,-D)
  v = (E1,O-A,-D) / (E1,E2,-D)
  t = (E1,E2,O-A) / (E1,E2,-D)

其中(A,B,C)表示以A,B,C为列向量的3x3矩阵的行列式。

现在我们使用以下恒等式：

  (A,B,C) = dot(A,cross(B,C))  (develop the determinant w.r.t. first column)

  (B,A,C) = -(A,B,C)           (swapping two vectors changes the sign)

  (B,C,A) =  (A,B,C)           (circular permutation does not change the sign)

现在我们得到：

u = -(E2,O-A,D)  / (D,E1,E2)
v =  (E1,O-A,D)  / (D,E1,E2)
t = -(O-A,E1,E2) / (D,E1,E2)  

使用:

N=cross(E1,E2);

AO = O-A; 

DAO = cross(D,AO)

我们最终得到以下代码（这里使用GLSL编写，易于转换为其他语言）：

bool intersect_triangle(
    in Ray R, in vec3 A, in vec3 B, in vec3 C, out float t, 
    out float u, out float v, out vec3 N
) { 
   vec3 E1 = B-A;
   vec3 E2 = C-A;
         N = cross(E1,E2);
   float det = -dot(R.Dir, N);
   float invdet = 1.0/det;
   vec3 AO  = R.Origin - A;
   vec3 DAO = cross(AO, R.Dir);
   u =  dot(E2,DAO) * invdet;
   v = -dot(E1,DAO) * invdet;
   t =  dot(AO,N)  * invdet; 
   return (det >= 1e-6 && t >= 0.0 && u >= 0.0 && v >= 0.0 && (u+v) <= 1.0);
}
 

当函数返回true时，交点由R.Origin + t * R.Dir给出。在三角形中，交点的重心坐标为u、v和1-u-v（用于Gouraud着色或纹理映射）。好消息是，你可以免费获得它们！

请注意，该代码是无分支的。它被一些ShaderToy上的我的着色器使用。

链接如下：

• https://www.shadertoy.com/view/tl3XRN
• https://www.shadertoy.com/view/3ltSzM

- BrunoLevy

---

平面的方程不应该是：dot(p,N) - dot(p1,N) = 0 吗？ - Garrett

最终算法中，希望更清楚地看到交点在u、v和t方面的位置。 - Chet

晚来一步，但是其他评论中还没有提到：在返回语句中，det >= 1e-6 这个条件应该改为 fabs(det) >= 1e-6 吧？我想这是检查矩阵列之间的依赖性，所以符号应该不重要。我检查了你的算法（用C语言实现），这是唯一让我无法通过测试的问题。 - Aldo

如果你想要分段版本而不是射线版本，只需在结尾添加一个检查，即t是否小于线段长度，并将线段转换为Dir为单位向量的射线。如果你想要直线版本，完全摆脱t >= 0的检查。 - Bryce Wagner

代码并不是无分支的。原因是 && 也使用了分支代码。 - Pirulax

Answer 2

@BrunoLevi: 你的算法似乎无法工作，请看下面的Python实现：

def intersect_line_triangle(q1,q2,p1,p2,p3):
    def signed_tetra_volume(a,b,c,d):
        return np.sign(np.dot(np.cross(b-a,c-a),d-a)/6.0)

    s1 = signed_tetra_volume(q1,p1,p2,p3)
    s2 = signed_tetra_volume(q2,p1,p2,p3)

    if s1 != s2:
        s3 = signed_tetra_volume(q1,q2,p1,p2)
        s4 = signed_tetra_volume(q1,q2,p2,p3)
        s5 = signed_tetra_volume(q1,q2,p3,p1)
        if s3 == s4 and s4 == s5:
            n = np.cross(p2-p1,p3-p1)
            t = -np.dot(q1,n-p1) / np.dot(q1,q2-q1)
            return q1 + t * (q2-q1)
    return None

我的测试代码如下：

q0 = np.array([0.0,0.0,1.0])
q1 = np.array([0.0,0.0,-1.0])
p0 = np.array([-1.0,-1.0,0.0])
p1 = np.array([1.0,-1.0,0.0])
p2 = np.array([0.0,1.0,0.0])

print(intersect_line_triangle(q0,q1,p0,p1,p2))

提供：

[ 0.  0. -3.] 

与期望的不同

[ 0.  0. 0.]

查看该行

t = np.dot(q1,n-p1) / np.dot(q1,q2-q1)

从法线中减去p1对我来说没有意义，您想要从q1投影到三角形的平面上，因此您需要沿着法线进行投影，并且距离应该与从q1到平面的距离和q1-q2沿法线的距离的比例成正比，对吗？

以下代码修复了这个问题:

n = np.cross(p2-p1,p3-p1)
t = np.dot(p1-q1,n) / np.dot(q2-q1,n)
return q1 + t * (q2-q1)

- Jochemspek

---

赞同这个答案。这个答案中纠正过的数学内容应该被合并到被接受的答案中。 - TiberiumFusion

非常感谢您注意到了这一点，我已经修正了另一个答案。 - BrunoLevy

Answer 3

要在3D中找到一条直线和一个三角形的交点，请按照以下步骤进行：

1.计算支持三角形的平面；

2.将直线与支持三角形的平面相交： - 如果没有相交，则不存在与三角形的交点； - 如果有相交，验证交点确实位于三角形内部： - 三角形的每条边及其支持平面法向量确定一个半空间，该半空间包围三角形内部（相应的包围平面可以从法向量和边顶点派生）； - 验证交点是否位于所有边半空间的内部。

以下是一些带有详细计算的示例代码，应该能够正常工作：

// Compute the plane supporting the triangle (p1, p2, p3)
//     normal: n
//     offset: d
//
// A point P lies on the supporting plane iff n.dot(P) + d = 0
//
ofVec3f v21 = p2 - p1;
ofVec3f v31 = p3 - p1;

ofVec3f n = v21.getCrossed(v31);
float d = -n.dot(p1);

// A point P belongs to the line from P1 to P2 iff
//     P = P1 + t * (P2 - P1)
//
// Find the intersection point P(t) between the line and
// the plane supporting the triangle:
//     n.dot(P) + d = 0
//                  = n.dot(P1 + t (P2 - P1)) + d
//                  = n.dot(P1) + t n.dot(P2 - P1) + d
//
//     t = -(n.dot(P1) + d) / n.dot(P2 - P1)
//
ofVec3f P21 = P2 - P1;
float nDotP21 = n.dot(P21);

// Ignore line parallel to (or lying in) the plane
if (fabs(nDotP21) < Epsilon)
    return false;

float t = -(n.dot(P1) + d) / nDotP21;
ofVec3f P = P1 + t * P21;

// Plane bounding the inside half-space of edge (p1, p2): 
//     normal: n21 = n x (p2 - p1)
//     offset: d21 = -n21.dot(p1)
//
// A point P is in the inside half-space iff n21.dot(P) + d21 > 0
//

// Edge (p1, p2)
ofVec3f n21 = n.cross(v21);
float d21 = -n21.dot(p1);

if (n21.dot(P) + d21 <= 0)
    return false;

// Edge (p2, p3)
ofVec3f v32 = p3 - p2;
ofVec3f n32 = n.cross(v32);
float d32 = -n32.dot(p2);

if (n32.dot(P) + d32 <= 0)
    return false;

// Edge (p3, p1)
ofVec3f n13 = n.cross(-v31);
float d13 = -n13.dot(p3);

if (n13.dot(P) + d13 <= 0)
    return false;

return true;

关于问题发布的代码的一些评论：

• 当可用时，应优先使用ofVec3f的预定义操作（如几何乘积的.dot()和.cross()等），因为这样更易读、避免了实现错误等问题。
• 代码最初遵循上述方法，但后来仅检查交点是否在线段[P1，P2]的三维轴对齐边界框中。这与可能存在的其他错误结合起来可能会导致结果不正确。
• 可以验证交点是否在（整个）三角形的三维轴对齐边界框中。虽然这还不足以保证交点存在，但可以用于剔除明显没有相交的点并避免进一步的复杂计算。

- user3146587

Answer 4

我有一种不同的方法来做这件事，我发现在我的渲染器中比BrunoLevy给出的第一种方法要快得多。（我还没有实现第二种方法）

点A、B、C是三角形的顶点 O是射线的起点 D是射线的方向（不需要归一化，只需比三角形更靠近原点）

检查方向（D+O）是否在四面体ABC、O内。

bool SameSide(vec3 A, vec3 B, vec3 C, vec3 D, vec3 p)
{
    vec3 normal = cross(B - A, C - A);
    float dotD = dot(normal, D - A);
    float dotP = dot(normal, p - A);
    return signbit(dotD) == signbit(dotP);
}

bool LineIntersectTri(vec3 A, vec3 B, vec3 C, vec3 O, vec3 D)
{
    return SameSide(A, B, C, O, O+D) &&
           SameSide(B, C, O, A, O+D) &&
           SameSide(C, O, A, B, O+D) &&
           SameSide(O, A, B, C, O+D);               
}

如果D变化，其他所有东西保持不变（例如在光线投射渲染器中），那么法线和点积就不需要重新计算；这就是我发现它快得多的原因。

这段代码来自于这个答案 https://dev59.com/voLba4cB1Zd3GeqPcEJl#25180294。

- Ben