3D 中线和三角形的交点

新手上路,请多包涵

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

我需要找出一种方法,希望使用 C++ 来确定这条线是否穿过三角形。一条平行于三角形的线,并且有多个共同点,应该算作“不相交”。

我已经编写了一些代码,但它不起作用,即使视觉表示清楚地显示了交叉点,我也总是会出错。

 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 发布,翻译遵循 CC BY-SA 4.0 许可协议

阅读 1.1k
2 个回答

1)如果您只想知道直线 是否 与三角形相交(不需要实际的交点):

p1,p2,p3 表示你的三角形

选择两个点 q1,q2 在线上两个方向都很远。

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

If SignedVolume(q1,p1,p2,p3) and SignedVolume(q2,p1,p2,p3) have different signs AND SignedVolume(q1,q2,p1,p2) , SignedVolume(q1,q2,p2,p3) and SignedVolume(q1,q2,p3,p1) have the same sign, then there是一个路口。

 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,《快速、最小存储射线三角交点》,图形工具杂志,第一卷。 2,1997 年,第21–28

(也可以看看:)

https://en.wikipedia.org/wiki/M%C3%B6ller%E2%8%93Trumbore_intersection_algorithm

该算法最终更简单(比我们在 1)和 2)中所做的指令更少),但理解起来要复杂得多。让我们一步一步推导出来。

符号:

  • O = 射线的原点,

  • D = 光线的方向矢量,

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

射线上的任意一点 P 可以写成 P = O + tD

三角形上的任意点 P 可以写成 P = A + uE1 + vE2 其中 E1 = B-AE2 = 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 矩阵)

使用克莱默公式求解:

    [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 给出。三角形中交点的重心坐标为 uv1-u-v (用于 Gouraud 着色或纹理映射)。好消息是您可以免费获得它们!

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

原文由 BrunoLevy 发布,翻译遵循 CC BY-SA 4.0 许可协议

要在 3D 中找到直线和三角形之间的交点,请遵循以下方法:

  • 计算支撑三角形的平面,

  • 将直线与支撑三角形的平面相交:

    • 如果没有交点,则与三角形没有交点。

    • 如果存在交点,请验证交点确实位于三角形中:

      • 三角形的每条边与支撑三角形的平面的法线一起确定了一个包围三角形内部的半空间(对应的边界平面可以从法线和边的顶点导出),
      • 验证交点是否位于所有边缘半空间的内侧。

这是一些示例代码,其中包含应该可以工作的详细计算:

 // 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] 的 3D 轴对齐边界框中。这与可能的其他错误相结合可以解释为什么结果不正确。
  • 可以验证交点是否在(整个)三角形的 3D 轴对齐边界框中。虽然这不足以保证相交,但它可以用于剔除明显不相交的点并避免进一步的复杂计算。

原文由 user3146587 发布,翻译遵循 CC BY-SA 3.0 许可协议

撰写回答
你尚未登录,登录后可以
  • 和开发者交流问题的细节
  • 关注并接收问题和回答的更新提醒
  • 参与内容的编辑和改进,让解决方法与时俱进
推荐问题