空间中直线段和三角形的相交算法

最新推荐文章于 2022-04-06 16:20:00 发布

惜双若雅

最新推荐文章于 2022-04-06 16:20:00 发布

阅读量1.9k

点赞数 2

分类专栏：综合

综合专栏收录该内容

3 篇文章 0 订阅

订阅专栏

原博客地址空间中直线段和三角形的相交算法

最近在看recast&detour源码的时候有遇到许多数学上的算法问题，特此记录，以便以后查看。

分析
问题：已经线段pq (p起点 q终点) ，三角形abc（a b c逆时针存储），判断pq与abc有无交点。

第一步：
判断三角形所在的面和线段所在的直线是否平行或者是否是从三角形面的背面进入。如果是就提前退出，否则进入第二步。

（使用面的法向量和线段的向量的点积）

第二步：
判断线段和三角形所在的面的交点是否在线段上。

平面的方程可以表示为：

X*normal = d （1）

其中 X 为平面上的点，d为原点到平面的距离。

以p为起点的，沿着pq方向的射线上的点的坐标可以表示为：

M = p + t*pq = p - t*qp （2）

假设M点即为平面与pq所在的直线的交点：(M点和a点均为在平面上的点，带入公式（1）中

(op→−t⋅qp→)⋅normal−→−−−−=oa−→⋅normal−→−−−−
(op→−t⋅qp→)⋅normal→=oa→⋅normal→

解得t为

t=ap−→⋅normal−→−−−−qp→⋅normal−→−−−−
t=ap→⋅normal→qp→⋅normal→

判断t是否（0,1）之间，如果是继续第三步，如果否返回false。

第三步：
判断是否交点在三角形内部。

三角形abc内的点M可以用三角形的重心坐标系表示：

M=(1−λ2−λ3)a+λ2b+λ3c
M=(1−λ2−λ3)a+λ2b+λ3c

其中如果lamda3 和 lamda2 都在（0,1）之间，那么表示M在abc的内部

改写成向量形式，且把M用公式（2）带入，得到：

op→−tqp→=(1−λ2−λ3)oa−→+λ2ob→+λ3oc→
op→−tqp→=(1−λ2−λ3)oa→+λ2ob→+λ3oc→

整理得：

λ2ab→+λ3ac→+tqp→=ap−→
λ2ab→+λ3ac→+tqp→=ap→

写成向量的形式：

⎛⎝⎜abxabyabzacxacyaczqpxqpyqpz⎞⎠⎟⎛⎝⎜λ2λ3t⎞⎠⎟=⎛⎝⎜apxapyapz⎞⎠⎟
(abxacxqpxabyacyqpyabzaczqpz)(λ2λ3t)=(apxapyapz)

求解lamda2 lamda3：（应用tips中的克莱姆法则以及混合积的运算法则）

λ2=(ap,ac,qp)(ab,ac,qp)=ac⋅(qp×ap)qp⋅(ab×ac)λ3=(ab,ap,qp)(ab,ac,qp)=ab⋅(ap×qp)qp⋅(ab×ac)=−ab⋅(qp×ap)qp⋅(ab×ac)
λ2=(ap,ac,qp)(ab,ac,qp)=ac⋅(qp×ap)qp⋅(ab×ac)λ3=(ab,ap,qp)(ab,ac,qp)=ab⋅(ap×qp)qp⋅(ab×ac)=−ab⋅(qp×ap)qp⋅(ab×ac)

源码
// 空间点 sp 起点 sq终点
// 三角形空间点 a b c
// 输出参数 t
static bool intersectSegmentTriangle(const float* sp, const float* sq,
                                   const float* a, const float* b, const float* c,
                                   float &t)
{
   float v, w;
   float ab[3], ac[3], qp[3], ap[3], norm[3], e[3];
   rcVsub(ab, b, a);
   rcVsub(ac, c, a);
   // 终点->起点的向量, 不是起点->终点
   rcVsub(qp, sp, sq);

   // Compute triangle normal. Can be precalculated or cached if
   // intersecting multiple segments against the same triangle
   // 法线方向为三角形正面方向
   rcVcross(norm, ab, ac);

   // Compute denominator d. If d <= 0, segment is parallel to or points
   // away from triangle, so exit early
   // d = 0.0 说明 qp 和 norm 垂直，说明三角形和 qp 平行。
   // d < 0.0 说明 qp 和 norm 是钝角说明是从三角形的背面进入和三角形相交的
   float d = rcVdot(qp, norm);
   if (d <= 0.0f) return false;

   // Compute intersection t value of pq with plane of triangle. A ray
   // intersects if 0 <= t. Segment intersects if 0 <= t <= 1. Delay
   // dividing by d until intersection has been found to pierce triangle
   rcVsub(ap, sp, a);
   t = rcVdot(ap, norm);
   if (t < 0.0f) return false;
   if (t > d) return false; // For segment; exclude this code line for a ray test

   // Compute barycentric coordinate components and test if within bounds
   rcVcross(e, qp, ap);
   v = rcVdot(ac, e);
   if (v < 0.0f || v > d) return false;
   w = -rcVdot(ab, e);
   if (w < 0.0f || v + w > d) return false;

   // Segment/ray intersects triangle. Perform delayed division
   t /= d;

   return true;
}
tips
克莱姆法则
一个线性方程组可以用矩阵（A为n*n的方阵）与向量（x，c均长度为n的列向量）的方程来表示：
Ax=c
Ax=c