探索射线与三角形求交的方法

Tomas Moller的Fast,Minimum Storage Ray/Triangle Intersection ,好像是99年的吧,这篇在讲到射线与三角形所在平面相交后,如何判断点在三角形内,运用了三角形重心坐标系的方法(详情看paper),同时运用行列式的克莱姆法则去求u,v( 重心坐标运用此两个来表示) ,不过又对其行列式进行了优化:

　　　|A B C| = -(A叉乘C)*B ＝-(C叉乘B)*A　，这样来减少计算量．

其后又看了另一个paper，讲到：　让射线与三角形所在平面相交后产生的点，在此平面一个方向上发射线，如果交点为奇数，则为在三角形里，否则在三角形外面．其它的继续探讨．．．．．

先把Moller的程序搞过来：

#include

<iostream>

using

namespace std
;

#define

EPSILON
0.000001

#define

CROSS(dest, v1,v2
) /

dest[0] = v1[1]*v2[2]-v1[2]*v2
[1]; /

dest[1] = v1[2]*v2[0]-v1[0]*v2
[2]; /

dest[2] = v1[0]*v2[1]-v1[1]*v2
[0];

#define

DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2
[2])

#define

SUB(dest,v1,v2
) /

dest[0] = v1[0]-v2
[0]; /

dest[1] = v1[1]-v2
[1]; /

dest[2] = v1[2]-v2
[2];

int

intersect_triangle

(double orig[3], double dir
[3],

double vert0[3],double vert1[3], double vert2
[3],

double *t, double *u, double *v
)
{

double edge1[3],edge2[3],tvec[3],pvec[3],qvec
[3];

double det,inv_det
;

// find vectors for two edges sharing vert0

SUB(edge1,vert1,vert0
);

SUB(edge2,vert2,vert0
);

// begin calculating determinant - also used to calculate U parameter

CROSS(pvec,dir,edge2
);

// if determinant is near zero , the ray lies in the plane or parallel to the plane

det = DOT(edge1,pvec
);

#ifdef

TEST_CULL
// define TEST_CULL if culling is desired

if(det<EPSILON)

return 0;

// calculate distance from vert0 to ray origin

SUB(tvec,orig,vert0);

// calculate U parameter and test bounds

*u = DOT(tvec,pvec);

if(*u<0.0 || *u>det)

return 0;

// prepare to test V parameter

CROSS(qvec,tvec,edge1);

// calculate V parameter and test bounds

*v = DOT(dir,qvec);

if(*v<0.0 || *u + *v >det)

return 0;

// calculate t, scale parameters, ray intersects triangle

*t = DOT(edge2,qvec);

inv_det = 1.0/det;

*t *= inv_det;

*u *= inv_det;

*v *= inv_det;

#else

if( det > -EPSILON && det < EPSILON
)

return
0;

inv_det = 1.0 / det
;

// calculate distance from vert0 to ray origin

SUB(tvec,orig,vert0
);

// calculate U parameter and test bounds

*

u = DOT(tvec,pvec)*inv_det
;

if(*u < 0.0 || *u
> 1.0)

return
0;

// calculate V parameter and test bounds

*

v = DOT(dir,qvec)*inv_det
;

if(*v<0.0 || *u + *v
> 1.0)

return
0;

// calculate t ,ray intersects triangle

*

t = DOT(edge2,qvec)*inv_det
;

#endif

return
1;
}

int

main
()
{

double origin
[3] = { 0,0,0};

double dir
[3] = {1.0,1.0,1.0};

double vert1
[3] ={4,5,6};

double vert2
[3] = {5,6,8};

double vert3
[3] = {7,8,9};

double t,u,v
;

if(intersect_triangle(origin,dir,vert1,vert2,vert3,&t,&u,&v
))

cout<<" The ray has intersected the triangle"<<endl
;

else

cout<<" The ray has not intersected the triangle"<<endl
;

cin.get
();

return
1;