目录
1. 概述
射线和三角形的相交检测是游戏程序设计中一个常见的问题,最典型的应用就是拾取(Picking),本文介绍一个最常见的方法,这个方法也是DirectX中采用的方法,该方法速度快,而且存储空间少。先讲述理论,然后给出对应的代码实现。
2. 理论部分
2.1. 一个直观的方法
我想大多数人在看到这个问题时,可能都会想到一个简单而直观的方法:首先判断射线是否与三角形所在的平面相交,如果相交,再判断交点是否在三角形内。
但是,上面的方法效率并不很高,因为需要一个额外的计算,那就是计算出三角形所在的平面,而下面要介绍的方法则可以省去这个计算。
2.2. 本文的方法
接下来会涉及到一些数学知识,不过没关系,我会详细解释每一个步骤,不至于太晦涩,只要您不觉得烦就行了,好了开始!
射线的参数方程如下,其中O是射线的起点,D是射线的方向。
我们可以这样理解射线,一个点从起点O开始,沿着方向D移动任意长度,得到终点R,根据t值的不同,得到的R值也不同,所有这些不同的R值便构成了整条射线,比如下面的射线,起点是P0,方向是u,p0 + tu也就构成了整条射线。
三角形的参数方程如下,其中V0,V1和V2是三角形的三个点,u, v是V1和V2的权重,1-u-v是V0的权重,并且满足u>=0, v >= 0,u+v<=1。
确切的说,上面的方程是三角形及其内部所有点的方程,因为三角形内任意一点都可以理解为从顶点V0开始,沿着边V0V1移动一段距离,然后再沿着边V0V2移动一段距离,然后求他们的和向量。至于移动多大距离,就是由参数u和v控制的。
于是,求射线与三角形的交点也就变成了解下面这个方程-其中t,u,v是未知数,其他都是已知的
移项并整理,将t,u,v提取出来作为未知数,得到下面的线性方程组
现在开始解这个方程组,这里要用到两个知识点,一是克莱姆法则,二是向量的混合积。
令E1 = V1 - V0,E2 = V2 - V0,T = O - V0上式可以改写成
根据克莱姆法则,可得到t,u,v的解分别是
将这三个解联合起来写就是
根据混合积公式
上式可以改写成
令
得到最终的公式,这便是下面代码中用到的最终公式了,之所以提炼出P和Q是为了避免重复计算
3. 代码部分
理论部分阐述完毕,开始上代码,这份代码来自DirectX SDK中的Demo,名字叫做Picking(拾取),该函数位于文件Pick.cpp的最末尾。这个函数有一个特点,就是判断语句特别多,因为对于一个频繁被调用的函数来说,效率是最重要的,这么多判断就是为了在某个条件不满足时,及时返回,避免后续不必要的计算。
// Determine whether a ray intersect with a triangle
// Parameters
// orig: origin of the ray
// dir: direction of the ray
// v0, v1, v2: vertices of triangle
// t(out): weight of the intersection for the ray
// u(out), v(out): barycentric coordinate of intersection
bool IntersectTriangle(const Vector3& orig, const Vector3& dir,
Vector3& v0, Vector3& v1, Vector3& v2,
float* t, float* u, float* v)
{
// E1
Vector3 E1 = v1 - v0;
// E2
Vector3 E2 = v2 - v0;
// P
Vector3 P = dir.Cross(E2);
// determinant
float det = E1.Dot(P);
// keep det > 0, modify T accordingly
Vector3 T;
if( det >0 )
{
T = orig - v0;
}
else
{
T = v0 - orig;
det = -det;
}
// If determinant is near zero, ray lies in plane of triangle
if( det < 0.0001f )
return false;
// Calculate u and make sure u <= 1
*u = T.Dot(P);
if( *u < 0.0f || *u > det )
return false;
// Q
Vector3 Q = T.Cross(E1);
// Calculate v and make sure u + v <= 1
*v = dir.Dot(Q);
if( *v < 0.0f || *u + *v > det )
return false;
// Calculate t, scale parameters, ray intersects triangle
*t = E2.Dot(Q);
float fInvDet = 1.0f / det;
*t *= fInvDet;
*u *= fInvDet;
*v *= fInvDet;
return true;
}
3.1. 参数说明
输入参数:前两个参数orig和dir是射线的起点和方向,中间三个参数v0,v1和v2是三角形的三个顶点。
输出参数:t是交点对应的射线方程中的t值,u,v则是交点的纹理坐标值
3.2. 代码说明
变量的命名方式:为了方便阅读,代码中的变量命名与上面公式中的变量保持一致,如E1,E2,T等。
变量det表示矩阵的行列式值
27-35行用来确保det>0,如果det<0则令det = -det,并对T做相应的调整,这样做是为了方便后续计算,否则的话需要分别处理det>0和det<0两种情况。
第38行,注意浮点数和0的比较,一般不用 == 0的方式,而是给定一个Epsilon值,并与这个值比较。
第43行,这里实际上u还没有计算完毕,此时的值是Dot(P,T),如果Dot(P,T) > det, 那么u > 1,无交点。
第51行,要确保u + v <= 1,也即 [Dot(P,T) + Dot(Q, D)] / det 必须不能大于1,否则无交点。
第57-60行,走到这里时,表明前面的条件都已经满足,开始计算t, u, v的最终值。
3.3. 交点坐标
根据上面代码求出的t,u,v的值,交点的最终坐标可以用下面两种方法计算
O + Dt
(1 - u - v)V0 + uV1 + vV2
3.4. osg中的应用
在osg的LineSegmentIntersector.cpp文件的如下函数也用到上述的算法,参见3.2节解释去理解即可:
void intersect(const osg::Vec3& v0, const osg::Vec3& v1, const osg::Vec3& v2)
{
if (_settings->_limitOneIntersection && _hit) return;
// const StartEnd startend = _startEndStack.back();
// const osg::Vec3& ls = startend.first;
// const osg::Vec3& le = startend.second;
Vec3 T = _start - v0;
Vec3 E2 = v2 - v0;
Vec3 E1 = v1 - v0;
Vec3 P = _d ^ E2;
value_type det = P * E1;
value_type r,r0,r1,r2;
const value_type epsilon = 1e-10;
if (det>epsilon)
{
value_type u = (P*T);
if (u<0.0 || u>det)
{
return;
}
osg::Vec3 Q = T ^ E1;
value_type v = (Q*_d);
if (v<0.0 || v>det)
{
return;
}
if ((u + v) > det)
{
return;
}
value_type inv_det = 1.0/det;
value_type t = (Q*E2)*inv_det;
if (t<0.0 || t>_length) return;
u *= inv_det;
v *= inv_det;
r0 = 1.0-u-v;
r1 = u;
r2 = v;
r = t * _inverse_length;
}
else if (det<-epsilon)
{
value_type u = (P*T);
if (u>0.0 || u<det) return;
Vec3 Q = T ^ E1;
value_type v = (Q*_d);
if (v>0.0 || v<det) return;
if ((u+v) < det) return;
value_type inv_det = 1.0/det;
value_type t = (Q*E2)*inv_det;
if (t<0.0 || t>_length) return;
u *= inv_det;
v *= inv_det;
r0 = 1.0-u-v;
r1 = u;
r2 = v;
r = t * _inverse_length;
}
else
{
return;
}
// Remap ratio into the range of LineSegment
const osg::Vec3d& lsStart = _settings->_lineSegIntersector->getStart();
const osg::Vec3d& lsEnd = _settings->_lineSegIntersector->getEnd();
double remap_ratio = ((_start - lsStart).length() + r*_length)/(lsEnd - lsStart).length();
Vec3 in = lsStart*(1.0 - remap_ratio) + lsEnd*remap_ratio; // == v0*r0 + v1*r1 + v2*r2;
Vec3 normal = E1^E2;
normal.normalize();
LineSegmentIntersector::Intersection hit;
hit.ratio = remap_ratio;
hit.matrix = _settings->_iv->getModelMatrix();
hit.nodePath = _settings->_iv->getNodePath();
hit.drawable = _settings->_drawable;
hit.primitiveIndex = _primitiveIndex;
hit.localIntersectionPoint = in;
hit.localIntersectionNormal = normal;
if (_settings->_vertices.valid())
{
const osg::Vec3* first = &(_settings->_vertices->front());
hit.indexList.reserve(3);
hit.ratioList.reserve(3);
if (r0!=0.0f)
{
hit.indexList.push_back(&v0-first);
hit.ratioList.push_back(r0);
}
if (r1!=0.0f)
{
hit.indexList.push_back(&v1-first);
hit.ratioList.push_back(r1);
}
if (r2!=0.0f)
{
hit.indexList.push_back(&v2-first);
hit.ratioList.push_back(r2);
}
}
_settings->_lineSegIntersector->insertIntersection(hit);
_hit = true;
}