1. 概述

任何复杂的三维模型都可以视作空间三角面片的集合，很容易碰到的一个问题就是空间射线与三角形相交的问题，例如拾取、遮蔽检测等。这里就总结下该问题的两种算法实现。

2. 常规算法

一种很常规的思路就是先计算射线与三角面片的交点，再看该交点是否再三角形内部。

2.1. 理论推导

对于空间一条射线，令起点为O，其方向为D，根据射线的参数公式，其上任意一点P（也就是要求的交点）为：

P=O+tD(1)(1)P=O+tD

其中t>0，根据t的取值不同，可得射线上不同的点，也就是关键在于求未知量t的值。

已知空间三角面片三个顶点为v1，v2，v3，那么很容易可以求得三角面片的法向量n。显然面上的向量(v1-P)与n是垂直的，则它们的点积为0：

(v1−P)⋅n=0(2)(2)(v1−P)⋅n=0

将式(1)代入式(2)，求得未知量t为：

t=(v1−O)⋅nD⋅nt=(v1−O)⋅nD⋅n

再将t代入到(1)式中，即可得到射线与该三点组成的平面了。

接下来就是判断这个交点是否在三角形面之内了，由于是空间三角形，所以比较好的算法是文献[2]中提到的同向法，摘录如下：

2.2. 具体实现

具体的C/C++实现代码如下：

#include <iostream>
using namespace std;
#define EPSILON 0.000001
// 3D vector
class Vector3d
{
public:
  Vector3d()
  {
  }
  ~Vector3d()
  {
  }
  Vector3d(double dx, double dy, double dz)
  {
    x = dx;
    y = dy;
    z = dz;
  }
  // 矢量赋值
  void set(double dx, double dy, double dz)
  {
    x = dx;
    y = dy;
    z = dz;
  }
  // 矢量相加
  Vector3d operator + (const Vector3d& v) const
  {
    return Vector3d(x + v.x, y + v.y, z + v.z);
  }
  // 矢量相减
  Vector3d operator - (const Vector3d& v) const
  {
    return Vector3d(x - v.x, y - v.y, z - v.z);
  }
  //矢量数乘
  Vector3d Scalar(double c) const
  {
    return Vector3d(c*x, c*y, c*z);
  }
  // 矢量点积
  double Dot(const Vector3d& v) const
  {
    return x * v.x + y * v.y + z * v.z;
  }
  // 矢量叉积
  Vector3d Cross(const Vector3d& v) const
  {
    return Vector3d(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
  }
  double _x()
  {
    return x;
  }
  double _y()
  {
    return y;
  }
  double _z()
  {
    return z;
  }
private:
  double x, y, z;
};
// v1 = Cross(AB, AC)
// v2 = Cross(AB, AP)
// 判断矢量v1和v2是否同向
bool SameSide(Vector3d A, Vector3d B, Vector3d C, Vector3d P)
{
  Vector3d AB = B - A;
  Vector3d AC = C - A;
  Vector3d AP = P - A;
  
  Vector3d v1 = AB.Cross(AC);
  Vector3d v2 = AB.Cross(AP);
  // v1 and v2 should point to the same direction
  //return v1.Dot(v2) >= 0 ;
  return v1.Dot(v2) > 0;
}
// 判断点P是否在三角形ABC内(同向法)
bool PointinTriangle1(Vector3d A, Vector3d B, Vector3d C, Vector3d P)
{
  return SameSide(A, B, C, P) && SameSide(B, C, A, P) && SameSide(C, A, B, P);
}
//ray-triangle intersection algorithm  (通过平面方程计算)
//参数说明：V1,V2,V3,三角形三点；O,射线原点；D,射线方向
bool ray_triangle_intersection1(Vector3d V1, Vector3d V2, Vector3d V3, Vector3d O, Vector3d D, Vector3d *I)
{
  bool rv = false;
  //v1(n1,n2,n3);
  //平面方程: na * (x – n1) + nb * (y – n2) + nc * (z – n3) = 0 ;
  double na = (V2._y() - V1._y())*(V3._z() - V1._z()) - (V2._z() - V1._z())*(V3._y() - V1._y());
  double nb = (V2._z() - V1._z())*(V3._x() - V1._x()) - (V2._x() - V1._x())*(V3._z() - V1._z());
  double nc = (V2._x() - V1._x())*(V3._y() - V1._y()) - (V2._y() - V1._y())*(V3._x() - V1._x());
  //平面法向量
  Vector3d nv(na, nb, nc);
  //平面法向量与射线方向向量差积
  double vpt = D.Dot(nv);
  if (vpt == 0)
  {
    rv = false;  //此时直线与平面平行
  }
  else
  {
    Vector3d P = V1 - O;
    double t = P.Dot(nv) / vpt;
    *I = O + D.Scalar(t);
    if (PointinTriangle1(V1, V2, V3, *I))
    {
      rv = true;
    }
    else
    {
      rv = false;
    }
  }
  return rv;
}
int main()
{
  Vector3d V1(0, 0, 0);
  Vector3d V2(50, 0, 0);
  Vector3d V3(0, 50, 0);
  Vector3d O(5, 10, -10);
  Vector3d P(10, 10, 10);
  Vector3d D = P - O;
  Vector3d I;
  if (ray_triangle_intersection1(V1, V2, V3, O, D, &I)) {
    cout << I._x() << '\t' << I._y() << '\t' << I._z() << endl;
  } 
}

3. 优化算法

仔细思考常规算法的思路，在计算射线与平面的交点的时候，实际是将射线的参数方程与平面的参数方程联立求值即可。那么如果知道空间三角形的参数方程，将其与射线的参数方程联立，不就可以直接求得交点了吗？Tomas Moller的论文《Fast, Minimum Storage Ray Triangle Intersection》提出了一种优化算法，正是基于这个思路，并且给出了合理的解法。

3.1. 理论推导

对于三个顶点为V1，V2，V3组成的空间三角形，对于三角形内的任一点，有如下参数方程：

P=(1−u−v)V1+uV2+vV3(3)(3)P=(1−u−v)V1+uV2+vV3

u, v是V2和V3的权重，1-u-v是V1的权重，并且满足u>=0, v >= 0,u+v<=1。这个参数方程的具体解释可参考文献[5]，摘录如下：

将射线公式(1)与三角形公式(3)联立起来，有：

(1−u−v)V1+uV2+vV3=O+tD(1−u−v)V1+uV2+vV3=O+tD

很显然，u、v、t都是未知数，移项并整理，可得如下线性方程组：

[−DV2−V1V3−V1]⎡⎢⎣tuv⎤⎥⎦=O−V1[−DV2−V1V3−V1][tuv]=O−V1

可以使用克莱姆法则来求解这个线性方程组，大家可以复习下线性代数(文献[6])，我这里也将其摘录如下：

令E1=V2−V1，E2=V3−V1，T=O−V1E1=V2−V1，E2=V3−V1，T=O−V1，则上式可以改写成：

[−DE1E2]⎡⎢⎣tuv⎤⎥⎦=T[−DE1E2][tuv]=T

根据克莱姆法则，有：

⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩t=1|−DE1E2||TE1E2|u=1|−DE1E2||−DTE2|v=1|−DE1E2||−DE1T|{t=1|−DE1E2||TE1E2|u=1|−DE1E2||−DTE2|v=1|−DE1E2||−DE1T|

接下来就要用到向量的混合积公式（具体可参看文献[7]）了，对于三向量a，b，c，有：

|abc|=a×b⋅c=−a×c⋅b=−c×b⋅a|abc|=a×b⋅c=−a×c⋅b=−c×b⋅a

上式可改写成：

⎧⎪ ⎪⎨⎪ ⎪⎩t=1D×E2⋅E1(T×E1⋅E2)u=1D×E2⋅E1(D×E2⋅T)v=1D×E2⋅E1(T×E1⋅D){t=1D×E2⋅E1(T×E1⋅E2)u=1D×E2⋅E1(D×E2⋅T)v=1D×E2⋅E1(T×E1⋅D)

令P=D×E2,Q=T×E1P=D×E2,Q=T×E1，进一步简化可得：

⎧⎪ ⎪⎨⎪ ⎪⎩t=1P⋅E1(Q⋅E2)u=1P⋅E1(P⋅T)v=1P⋅E1(Q⋅D){t=1P⋅E1(Q⋅E2)u=1P⋅E1(P⋅T)v=1P⋅E1(Q⋅D)

3.2. 具体实现

具体的C/C++实现代码如下：

#include <iostream>
using namespace std;
#define EPSILON 0.000001
// 3D vector
class Vector3d
{
public:
  Vector3d()
  {
  }
  ~Vector3d()
  {
  }
  Vector3d(double dx, double dy, double dz)
  {
    x = dx;
    y = dy;
    z = dz;
  }
  // 矢量赋值
  void set(double dx, double dy, double dz)
  {
    x = dx;
    y = dy;
    z = dz;
  }
  // 矢量相加
  Vector3d operator + (const Vector3d& v) const
  {
    return Vector3d(x + v.x, y + v.y, z + v.z);
  }
  // 矢量相减
  Vector3d operator - (const Vector3d& v) const
  {
    return Vector3d(x - v.x, y - v.y, z - v.z);
  }
  //矢量数乘
  Vector3d Scalar(double c) const
  {
    return Vector3d(c*x, c*y, c*z);
  }
  // 矢量点积
  double Dot(const Vector3d& v) const
  {
    return x * v.x + y * v.y + z * v.z;
  }
  // 矢量叉积
  Vector3d Cross(const Vector3d& v) const
  {
    return Vector3d(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);
  }
  double _x()
  {
    return x;
  }
  double _y()
  {
    return y;
  }
  double _z()
  {
    return z;
  }
private:
  double x, y, z;
};
//ray-triangle intersection algorithm
//参数说明：V1,V2,V3,三角形三点；O,射线原点；D,射线方向。
bool ray_triangle_intersection(Vector3d V1, Vector3d V2, Vector3d V3, Vector3d O, Vector3d D, Vector3d *I)
{
  //Find vectors for two edges sharing V1
  Vector3d e1 = V2 - V1;
  Vector3d e2 = V3 - V1;
  //Begin calculating determinant - also used to calculate u parameter
  Vector3d P = D.Cross(e2);
  //if determinant is near zero, ray lies in plane of triangle
  double det = e1.Dot(P);
  //NOT CULLING
  if (det > -EPSILON && det < EPSILON)
  {
    return false;
  }
  double inv_det = 1.f / det;
  //calculate distance from V1 to ray origin
  Vector3d T = O - V1;
  //Calculate u parameter and test bound
  double u = T.Dot(P) * inv_det;
  //The intersection lies outside of the triangle
  if (u < 0.f || u > 1.f)
  {
    return false;
  }
  //Prepare to test v parameter
  Vector3d Q = T.Cross(e1);
  //Calculate V parameter and test bound
  double v = D.Dot(Q) * inv_det;
  //The intersection lies outside of the triangle
  if (v < 0.f || u + v  > 1.f)
  {
    return false;
  }
  double t = e2.Dot(Q) * inv_det;
  //ray intersection
  if (t > EPSILON)
  {
    *I = O + D.Scalar(t);
    return true;
  }
  return false;
}
int main()
{
  Vector3d V1(0, 0, 0);
  Vector3d V2(50, 0, 0);
  Vector3d V3(0, 50, 0);
  Vector3d O(5, 10, -10);
  Vector3d P(10, 10, 10);
  Vector3d D = P - O;
  Vector3d I;
  if (ray_triangle_intersection(V1, V2, V3, O, D, &I)) {
    cout << I._x() << '\t' << I._y() << '\t' << I._z() << endl;
  }
}

可以看到这种优化算法无论是代码量还是时间、空间复杂度都由于原来的常规算法，最直观的体现就是判断语句多，能够即使返回避免后续运算。

4. 参考

[1] Möller–Trumbore intersection algorithm

[2] 判断点是否在三角形内

[3] 射线与平面的相交检测(Ray-Plane intersection test)

[4] 射线和三角形的相交检测（ray triangle intersection test）

[5] 三角形方程？ - 高崎汀步的回答 - 知乎

[6] 克莱姆法则

[7] 三矢量的混合积

分类: 计算几何

标签: 相交 , 射线 , 三角形 , C++ , 算法

空间射线与三角形相交算法的两种实现

空间射线与三角形相交算法的两种实现

1. 概述

2. 常规算法

2.1. 理论推导

2.2. 具体实现

3. 优化算法

3.1. 理论推导

3.2. 具体实现

4. 参考

热门文章

最新文章

相关课程

相关电子书

相关实验场景

热门

活动广场

任务中心

开发者评测

高校计划

乘风者计划

训练营

阿里云MVP

话题

直播

下载

镜像站

技术资料

插件

空间射线与三角形相交算法的两种实现

1. 概述

2. 常规算法

2.1. 理论推导

2.2. 具体实现

3. 优化算法

3.1. 理论推导

3.2. 具体实现

4. 参考

热门文章

最新文章

相关课程

相关电子书

相关实验场景