两个矩阵是否相交的算法_直线-三角形相交算法总结

算法一： 直接与平面相交 + 交点是否处于三角形内部

直线与三角形相交，首先三角形是在平面上，直线与平面的交点比较容易求得。

平面公式：

直线可以写成：

代入：

当我们计算出来 t 之后，再代回原公式得到交点P：

这里有一点点危险的地方

可能为0， 这个时候就意味着 N 和 AB 垂直，那么直线与平面平行，它与平面没有交点。

当直线与平面相交之后，我们需要判断交点是否在三角形内， 之前在二维平面上，我们用过这样的办法： 判断

和三角形边

的法向量的 n 的点乘是否小于0来看P点是否位于三角形边的‘左侧’。

实际上这个思路也可以推广到三维：

n 是三角形的normal，也是平面的normal是

结果。如果P点在三角形内侧，那么

的向量 n' 会方向相同，也就是它们的点乘会大于0.

查看另外两边，也会有类似的效果，所以三角形和射线相交可以这样写：

bool ray_triangle_intersect(const Vec3f A, const Vec3f B, const Vec3f p0, const Vec3f p1, const Vec3f p2){

  // 1. find the p point
  Vec3f n = cross(p1 - p0, p2 - p0);
  float d = -n*p1;

  Vec3f AB = B - A;
  if (fabsf(n*AB) < epsilon) return false;

  float t = (-d - n*A )/(n * AB);
  // t < 0 means opposite direction of AB
  if (t < 0 ) return false;
  Vec3f p = A + AB*t;

  if (cross(p1-p0, p-p0) * n < 0 ) return false;
  if (cross(p2-p1, p-p1) * n < 0 ) return false;
  if (cross(p0-p2,p-p2) * n < 0 ) return false;
  return true;
}

当然我们交点P也是直接求出来了。

算法二： 重心坐标系

重心坐标系我之前也写过 → 三角形重心坐标

实际上重心坐标系也叫面积坐标，如果三角形面积为1的话，那么P点分割的三角形有以下性质：

如果想验证也很简单：

所以这给我们提供了一个比较简单的重心坐标系的算法，再结合算法一,我们可以很容易的算出 u, v.

// AB: ray, p0p1p2: triangle
bool bary_centric_coordinate(const Vec3f A, const Vec3f B, const Vec3f p0, const Vec3f p1, const Vec3f p2, float &u, float &v){

  // 1. find the p point, same as above
  Vec3f n = cross(p1 - p0, p2 - p0);
  float area = n.norm();
  float d = -n*p1;

  Vec3f AB = B - A;
  if (fabsf(n*AB) < epsilon) return false;

  float t = (-d - n*A )/(n * AB);
  if (t < 0 ) return false;
  Vec3f p = A + AB*t;

  // 2. find u,v
  Vec3f uvector, vvector;
  if (cross(p2-p1, p-p1) * n < 0 ) return false;
  if ((vvector = cross(p1-p0,p-p0)) * n < 0) return false;
  if ((uvector = cross(p0-p2,p-p2)) * n < 0) return false;

  u = uvector.norm()/area;
  v = vvector.norm()/area;

  return true;
}

实际上代码跟算法一很多部分都是一致的，只是这里我们加了计算 u, v 而已。

Möller-Trumbore 算法

Möller-Trumbore算法我感觉和重心坐标系差不多，只是加入了更多线性代数的优化？三角形依旧是 ABC， 假设光线源点为O，方向为D，有相交点满足：

矩阵形式：

利用Cramer's rule，可知：

其中 T = O - A, E1 = B - A, E2 = C - A.

又线性代数中行列式的性质：

继续：

其中

.

所以整个计算就是我们无需再去计算 三角形平面的一些性质，取而代之我们用以上式子就可以计算出 t, u, v.

// OD: ray, p0p1p2: triangle
bool ray_triangle_intersect_mt(const Vec3f O, const Vec3f D, const Vec3f p0, const Vec3f p1, const Vec3f p2, float &t, float &u, float &v){
  Vec3f e1 = p1-p0;
  Vec3f e2 = p2-p0;
  Vec3f pvec = cross(D, e2);
  float det = e1*pvec;
  if (det < epsilon) return false;

  Vec3f tvec = O - p0;
  u = tvec*pvec*(1./det);
  if (u < 0 || u > 1) return false;

  Vec3f qvec = cross(tvec, e1);
  v = D*qvec*(1./det);
  if (v < 0 || u + v > 1) return false;

  t = e2*qvec*(1./det);
  return t > epsilon;
}

结束

直线与三角形相交是如此重要，是因为有了这些算法，在‘光线追踪’中，我们可以放model来玩，同时在 ‘光栅化’中，利用三角形的重心坐标来做各种插值也算光栅化的基石之一。

代码：

KrisYu/miscellaneous​github.com

参考：

Fast, minimum storage ray-triangle intersection.

Ray Tracing: Rendering a Triangle