Fast ray tri (with precomputation)
November 9, 2002
by Charles Bloom, cbloom@cbloom.com

Moller/Haines and others have good ray-tri intersection code for
the case of "no precomputation", but sometimes you want to do
precomputation.  For example, if you have only a few tris and
you want to shoot lots of rays at them, you should precompute
some things on your triangle to speed up the test.  Here's how.
I actually want segment-tri intersection, that is, my ray has
two end points.

We want to write a routine like

bool RayTriIntersect(const Triangle & tri,const Vec3 & from,
	const Vec3 & to, );

Compute the plane of the triangle.  The plane is a normal and
offset (4 floats).  Now you can do your first tests on your
segment :

dFm = tri.plane.Distance(from);
dTo = tri.plane.Distance(to);

Where plane.Distance is just (plane.normal * vec - plane.offset)

Now you need to check the signs of dFm and dTo to see if the
segment passes through the plane.  You have some options here
depending on where you want to treat your triangle as one-sided
or two-sided.  I'll assume that you want a one-sided triangle,
so only rays going "at" the triangle should count as hits, 
so :

if ( dTo <= 0 )
{
	// segment doesn't reach back side of triangle
	return false;
}
if ( dFm > 0 )
{
	// segment starts already on the back side
	return false;
}
assert( dTo > 0 && dFm <= 0 )

Now we know that the plane is hit.  We can now make the
point on the plane (we don't actually want to do this, but
we'll write it out) :

float t = dFm / (dTo - dFm)
onPlane = lerp(from,to,t);

Now we want to check if this point is on the triangle.  We do
that using barycentric coordinates.  To compute the barycentric
coords, you can use scaled planes.  You do :

edge1 = tri.vert1 - tri.vert0;
edge2 = tri.vert2 - tri.vert1;
edgePerp1 = tri.normal CROSS edge1
edgePerp2 = tri.normal CROSS edge2
edgePlane1 = { edgePerp1, edgePerp1*tri.vert0 }
edgePlane2 = { edgePerp2, edgePerp2*tri.vert1 }

Now a plane is just a vector and a constant; we can scale this
normal to change the distances we want, so we make

baryPlane1 = edgePlane1 / (edgePlane1.Distance(tri.vert2))
baryPlane2 = edgePlane2 / (edgePlane2.Distance(tri.vert0))

The result are these two "baryPlanes".  These give you the
barycentric coords just by doing dots :

u = baryPlane1.Distance(onPlane);
v = baryPlane2.Distance(onPlane);

These baryPlanes are pretty neat; they are the planes that lie
along the triangle edges and normal; their vectors are scaled
so that the opposite point on the triangle returns a distance
of 1.  Once we have these barycentric coordinates, we can
check collision easily :

if ( u < 0 || v < 0 || (u+v) > 1 )
{
	// no collision
	return false;
}
// collision !

And we also have the time of hit and the bary-coords, which is
good info to have.  Now lets put it all together and delay the
divide :

struct Triangle
{
	Plane	plane;
	Plane	bary1;
	Plane	bary2;
};

bool RayTriIntersect(const Triangle & tri,const Vec3 & from,
	const Vec3 & to)
{
	float dTo = tri.plane.Distance(to);
	if ( dTo <= 0 )
	{
		// segment doesn't reach back side of triangle
		return false;
	}

	float dFm = tri.plane.Distance(from);
	if ( dFm > 0 )
	{
		// segment starts already on the back side
		return false;
	}
	assert( dTo > 0 && dFm <= 0 );

	float denom = dTo - dFm;
	assert( denom > 0 );
	Vec3 temp = dTo * from + dFm * to;

	float uTimesDenom = (tri.bary1.vector DOT temp) - tri.bary1.offset * denom;

	if ( uTimesDenom < 0 || uTimesDenom > denom )
	{
		// off triangle
		return false;
	}

	float vTimesDenom = (tri.bary2.vector DOT temp) - tri.bary2.offset * denom;

	if ( vTimesDenom < 0 || (uTimesDenom+vTimesDenom) > denom )
	{
		// off triangle
		return false;
	}

	// it's on the triangle, compute hit info :
	float inv = 1 / denom;
	*pT = dFm * inv;
	*pU = uTimesDenom * inv;
	*pV = vTimesDenom * inv;

	return true;
}

On my machine (1 GHz P3), I measure about 165 clocks for Moller's optimized ray-tri #1 
(late division).  See :

http://www.ce.chalmers.se/staff/tomasm/raytri/

The code here takes about 120 clocks.  The exact numbers depend on how much my extra
early-outs help.  The 120 clock number is with no help from the early-outs ; in
practice, the number is much lower, more like 50-60 clocks average in typical use.

Most of all, I like the mental construct of the "barycentric planes".

The size of the pre-computed "Triangle" structure is very small.  It's only 12 floats,
compared to 9 floats to store the triangle vertices.  Thus, the memory-fetch cost of
this routine is not significantly higher.  Note that I don't use the triangle's actual
vertices at all in the ray-intersection routine.

Finally, note that in the real world there should be higher-level acceleration 
structures used for ray intersection of all kinds.

====================================================================================
== Addendum : some notes on Moller-Haines ray-triangle test :

I've been thinking about the classic Moller-Haines ray-triangle test.
The basic ray-tri test looks like this :

bool intersect_triangle(const Vec3& orig, const Vec3& dir,
			const Vec3& vert0, const Vec3& vert1, const Vec3& vert2)
{
	const Vec3 edge1 = v1 - v0;
	const Vec3 edge2 = v2 - v0;

	/* begin calculating determinant - also used to calculate U parameter */
	// ^ = cross product
	Vec3 pvec = dir ^ edge2;

	/* if determinant is near zero, ray lies in plane of triangle */
	const float det = fabsf( edge1 * pvec ); // = Vec3U::TripleProduct(dir,edge1,edge2);

	/* calculate distance from vert0 to ray origin */
	const Vec3 tvec = orig - vert0;

	/* calculate U parameter and test bounds */
	const float u = tvec * pvec; // = Vec3U::TripleProduct(dir,edge2,tvec);
	if ( u < 0 || u > det )
		return false;

	/* prepare to test V parameter */

	/* calculate V parameter and test bounds */
	const float v = Vec3U::TripleProduct(dir,tvec,edge1); // = dir * (tvec ^ edge1)
	if ( v < 0 || u + v > det )
		return false;

	return true;
}


It's interesting what's going on here geometrically. First let's back up
and think about how you do a ray-triangle test. The obvious way is to
take your ray, find the point where it hits the plane of the triangle,
and make the point on that plane. Now you have to do a 2d
point-in-triangle on that plane. You do this with Barycentric
coordinates. They good way to do that is using areas. Let's say your
verts are (v0,v1,v2), and your point on the plane is (P). You compute
Barycentric coordinates like this :

u = Area(v0,v1,P) / Area(v0,v1,v2)
v = Area(v0,P,v2) / Area(v0,v1,v2)
w = Area(P,v1,v2) / Area(v0,v1,v2) = 1 - u - v

and the areas are signed and the winding matters. Now, computing these
Area's in 3d is actually a little bit of a pain, because you would have
to project down to 3d and such. In fact, the easier thing to do is to
compute the *volume* of the parallelopiped defined by the triangle
extruded along its normal direction. Let me switch to coordinates where
v0 = 0 by doing :

e1 = v1 - v0
e2 = v2 - v0
Q =  P  - v0

u = Area(0,e1,Q) / Area(0,e1,e2)
v = Area(0,Q,e2) / Area(0,e1,e2)

now if we have N = e1 ^ e2 , which is proportional to the triangle normal, we can do :

u = Volume(N,e1,Q) / Volume(N,e1,e2)
v = Volume(N,Q,e2) / Volume(N,e1,e2)

note that the sign and magnitude of N actually don't matter since we
divide through by them. Now the Volume of 3 vectors is just a "triple
product" :

Volume(a,b,c) = a * (b ^ c) = Epsilon_ijk * a_i b_j c_k

which is the determinant of the matrix with a,b,c as columns.

Now, you can notice that what direction we use for N actually doesn't
matter, it just introduces a scaling factor, since the ratio of Volumes
is all that matters. That is, we can shear our parallelopiped so it's
not at right angles. 

u = Volume(D,e1,Q) / Volume(D,e1,e2)
v = Volume(D,Q,e2) / Volume(D,e1,e2)

This works for any D. In particular, we can use the original ray
direction. The only time this breaks down is if D is in the plane of the
triangle, but in that case you ray-tri is falling apart anyway.

Next, note that if we add anything to Q in the direction of D, it
doesn't affect the answer. That's because we cross product it out - 

Volume(D,e1,Q + tD) = Volume(e1,Q + tD,D) = e1 * (Q+tD) X D = e1 * Q X D = Volume(e1,Q,D) = Volume(D,e1,Q)

That means we can use the original ray source point, we don't have to
use the point Q that's on the triangle surface. This gives us our final
equations :

u = Volume(dir,e1,orig-v0) / Volume(dir,e1,e2)
v = Volume(dir,orig-v0,e2) / Volume(dir,e1,e2)

We can understand this geometrically as the volumes of the
parallelopides from the triangle verts, in the direction of the ray
normal, through the ray origin. If you look along the ray direction, you
see the triangle in 2d; the ray origin must lie inside that triangle to
have an intersection, which is exactly this test.

This trivially leads us to the Moller-Haines tests :

bool intersect()
{
	det = Volume(dir,e1,e2);

	tvec = orig - v0;

	u = Volume(dir,e1,tvec);

	if ( u < 0 || u > det )
		return false;

	v = Volume(dir,tvec,e2);

	if ( v < 0 || v > det )
		return false;

	return true;
}