raytri.cpp
/* Ray-Triangle Intersection Test Routines */
/* Different optimizations of my and Ben Trumbore's */
/* code from journals of graphics tools (JGT) */
/* http://www.acm.org/jgt/ */
/* by Tomas Moller, May 2000 */
#include <math.h>
#define EPSILON 1e-6
#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];
/* the original jgt code */
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, ray lies in plane of triangle */
det = DOT(edge1, pvec);
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;
/* prepare to test V parameter */
CROSS(qvec, tvec, edge1);
/* 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;
return 1;
}
/* code rewritten to do tests on the sign of the determinant */
/* the division is at the end in the code */
int intersect_triangle1(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, ray lies in plane of triangle */
det = DOT(edge1, pvec);
if (det > EPSILON)
{
/* 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;
}
else if(det < -EPSILON)
{
/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);
/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
/* printf("*u=%f\n",(float)*u); */
/* printf("det=%f\n",det); */
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;
}
else return 0; /* ray is parallell to the plane of the triangle */
inv_det = 1.0 / det;
/* calculate t, ray intersects triangle */
*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;
return 1;
}
/* code rewritten to do tests on the sign of the determinant */
/* the division is before the test of the sign of the det */
int intersect_triangle2(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, ray lies in plane of triangle */
det = DOT(edge1, pvec);
/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);
inv_det = 1.0 / det;
if (det > EPSILON)
{
/* 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;
}
else if(det < -EPSILON)
{
/* 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;
}
else return 0; /* ray is parallell to the plane of the triangle */
/* calculate t, ray intersects triangle */
*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;
return 1;
}
/* code rewritten to do tests on the sign of the determinant */
/* the division is before the test of the sign of the det */
/* and one CROSS has been moved out from the if-else if-else */
int intersect_triangle3(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, ray lies in plane of triangle */
det = DOT(edge1, pvec);
/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);
inv_det = 1.0 / det;
CROSS(qvec, tvec, edge1);
if (det > EPSILON)
{
*u = DOT(tvec, pvec);
if (*u < 0.0 || *u > det)
return 0;
/* calculate V parameter and test bounds */
*v = DOT(dir, qvec);
if (*v < 0.0 || *u + *v > det)
return 0;
}
else if(det < -EPSILON)
{
/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
if (*u > 0.0 || *u < det)
return 0;
/* calculate V parameter and test bounds */
*v = DOT(dir, qvec) ;
if (*v > 0.0 || *u + *v < det)
return 0;
}
else return 0; /* ray is parallell to the plane of the triangle */
*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;
return 1;
}
/* Sueda: inclusive borders */
#define ZERO -EPSILON
int intersect_triangle3_inc(
const double *orig, const double *dir,
const double *vert0, const double *vert1, const double *vert2,
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, ray lies in plane of triangle */
det = DOT(edge1, pvec);
/* calculate distance from vert0 to ray origin */
SUB(tvec, orig, vert0);
inv_det = 1.0 / det;
CROSS(qvec, tvec, edge1);
if (det > EPSILON)
{
*u = DOT(tvec, pvec);
if (*u < -ZERO || *u > det + ZERO)
return 0;
/* calculate V parameter and test bounds */
*v = DOT(dir, qvec);
if (*v < -ZERO || *u + *v > det + ZERO)
return 0;
}
else if(det < -EPSILON)
{
/* calculate U parameter and test bounds */
*u = DOT(tvec, pvec);
if (*u > ZERO || *u < det - ZERO)
return 0;
/* calculate V parameter and test bounds */
*v = DOT(dir, qvec) ;
if (*v > ZERO || *u + *v < det - ZERO)
return 0;
}
else return 0; /* ray is parallell to the plane of the triangle */
*t = DOT(edge2, qvec) * inv_det;
(*u) *= inv_det;
(*v) *= inv_det;
return 1;
}
