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

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;