3D 中线和三角形的交点
Intersection between line and triangle in 3D
我在 3D 中的某处有一条线和一个三角形 space。换句话说，三角形有 3 个点（每个 [x,y,z]），直线有两个点（也有 [x,y,z]）。
我需要找出一种方法，希望使用 C++，来确定直线是否穿过三角形。一条平行于三角形的线，并且有一个以上的公共点，应该算作 "does not intersect".
我已经编写了一些代码，但它不起作用，即使视觉表示清楚地显示了交叉点，我也总是出错。
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;
}
要在 3D 中找到直线和三角形之间的交点，请遵循以下方法：
- 计算支撑三角形的平面，
直线与支撑三角形的平面相交:
- 如果没有交集，则与三角形没有交集。
如果有交点，验证交点确实在三角形内：
- 三角形的每条边连同支撑三角形的平面的法线决定了三角形内部的半个space边界（对应的边界平面可以由法线和边顶点导出),
- 验证交点是否位于所有边半spaces 的内侧。
下面是一些示例代码，其中包含详细的计算，应该可行：
// 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 轴对齐边界框内。虽然这不足以保证相交，但它可以用来剔除明显不相交的点，避免进一步的复杂计算。
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)有不同的符号AND
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，« 快速、最小存储光线-三角形相交 »，
图形工具杂志，卷。 2, 1997, p. 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-A and E2 = C-A, u>=0, v>=0 and (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);
}
当函数returnstrue时，交点由R.Origin + t * R.Dir给出。三角形中交点的重心坐标为 u、v、1-u-v（对 Gouraud 着色或纹理映射很有用）。好消息是您可以免费获得它们！
请注意，代码是无分支的。
它被我在 ShaderToy
上的一些着色器使用
@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)
我在 3D 中的某处有一条线和一个三角形 space。换句话说，三角形有 3 个点（每个 [x,y,z]），直线有两个点（也有 [x,y,z]）。
我需要找出一种方法，希望使用 C++，来确定直线是否穿过三角形。一条平行于三角形的线，并且有一个以上的公共点，应该算作 "does not intersect".
我已经编写了一些代码，但它不起作用，即使视觉表示清楚地显示了交叉点，我也总是出错。
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;
}
要在 3D 中找到直线和三角形之间的交点，请遵循以下方法：
- 计算支撑三角形的平面，
直线与支撑三角形的平面相交:
- 如果没有交集，则与三角形没有交集。
如果有交点，验证交点确实在三角形内：
- 三角形的每条边连同支撑三角形的平面的法线决定了三角形内部的半个space边界（对应的边界平面可以由法线和边顶点导出),
- 验证交点是否位于所有边半spaces 的内侧。
下面是一些示例代码，其中包含详细的计算，应该可行：
// 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 轴对齐边界框内。虽然这不足以保证相交，但它可以用来剔除明显不相交的点，避免进一步的复杂计算。
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)有不同的符号AND
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，« 快速、最小存储光线-三角形相交 »， 图形工具杂志，卷。 2, 1997, p. 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-A and E2 = C-A, u>=0, v>=0 and (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);
}
当函数returnstrue时，交点由R.Origin + t * R.Dir给出。三角形中交点的重心坐标为 u、v、1-u-v（对 Gouraud 着色或纹理映射很有用）。好消息是您可以免费获得它们！
请注意，代码是无分支的。 它被我在 ShaderToy
上的一些着色器使用@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)