1. 概述

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

2. 常规算法

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

2.1. 理论推导

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

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

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

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

再将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组成的空间三角形，对于三角形内的任一点，有如下参数方程：

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

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

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

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

令\(E1 = V2 – V1，E2 = V3 – V1，T = O – V1\)，则上式可以改写成：

根据克莱姆法则，有：

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

上式可改写成：

令\(P=D \times E2, Q = T \times E1\)，进一步简化可得：

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] 三矢量的混合积