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

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;