Skip to main content
GameDev.net gamedev.net
🔒 Locked

Ray-triangle intersection

Started by _Sauce_ Feb 21, 2010 at 5:00 AM 31 replies 13.9k views
Original Post
_Sauce_
_Sauce_
I'm working on my raytracer and in an effort to speed it up a bit, I'm ditching the algebraic method of ray-triangle intersection, and am replacing it with a parametric method, as described here. I'm gonna start talking some math stuff here so I made this diagram to help explain the method I'm using. Ro = origin of the ray. Rd = direction of the ray. P = point of intersection on the plane constructed from the three vertices of the triangle v1, v2, v3 = vertices of the triangle, in clockwise winding order. now according to the slides in the Princeton lecture, in order for P to be inside the triangle, a and b must both be in the range 0..1. My 3D math isn't the best but I recognised that a is simply the scalar projection of Z onto X, and b is the scalar projection of Z onto Y, where Z = P - v3, X = v2 - v3, and Y = v1 - v3. Correct? So here's my implementation... 
float RTTriangle::intersect(const Rayf &ray)
{
	float distance = mPlane.intersect(ray);
	if(distance == std::numeric_limits<float>::infinity())
	{
		return distance;
	}

	Vector3f P = ray.mOrigin + ray.mDir * res.distance;

	Vector3f Z(P - mVerts[2]); //P - v3
	Vector3f X(mVerts[1] - mVerts[2]); //v2 - v3

	X.normaliseSelf();
	float a = Z.dot(X);
	if(a < 0 || a > 1)
	{
		return std::numeric_limits<float>::infinity();
	}

	Vector3f Y(mVerts[0] - mVerts[2]); //v1 - v3

	Y.normaliseSelf();
	float b = Z.dot(Y);
	if(b < 0 || b > 1)
	{
		return std::numeric_limits<float>::infinity();
	}

	return res;
}
And the output... Doesn't quite look like the famous Utah Teapot, now does it? Where exactly did I go wrong here?
click
szecs
szecs
in case of CCW winding (I suggest you to use CCW)

if [(v3-p) cross (v1-p)] dot normal > 0
and [(v1-p) cross (v2-p)] dot normal > 0
and [(v2-p) cross (v3-p)] dot normal > 0

then we have a hit.

I use that for convex polygons (with any number of sides). Notice that you don't need to normalize that way at all, even the normal can be non-normalized.

There is a much simpler way too (google for methods, I can't recall the name), that will use some u,v parameters for checking.
It's much faster in case of triangles, and u,v coordinates can be used to calculate texture lookup coordinates with a very simple formula.
So I suggest you to use that.

[Edited by - szecs on February 21, 2010 6:44:49 AM]
szecs
szecs
By the way, your logic is wrong too: You would have to project with the direction of the other side of the triangle (and that way you only checked that the ray is inside a parallelogram), but you are projecting in a perpendicular direction to the side, so you check two sides and a some curve, I can't decide what kind at the moment (maybe ellipsis-arc, but more likely not).
haegarr
haegarr
Oh well, either I totally misunderstand something, or the slide that is cited by the OP is totally wrong!?

One can express p as
p = v3 + a * y + b * x
so that the condition for a triangle would be
a,b >=0 and a+b <= 1
The less restrict condition
0 <= a,b <= 1
would make up a parallelogram instead.

Now, looking at the slide, the image shows a and b being vectors but the text denotes them being scalars. Well ... okay.

Then, with a and b being scalars, the shown formula
P = a ( T2-T1 ) + b ( T3-T1 )
would indicate the same procedure as I've shown above: Adding 2 vectors that are along T2T1 and T3T1, resp., but then the difference vectors to P would not neccessarily be perpendicular to the sides.

Further, assume a triangle in 2D where
T1 = (0,0) ; T2 = (1,0) ; T3 = (0,1)
and a point of interest
P = ( 1,1 )
Doing the calculations gives me a=b=1, so that according to the slide, P is inside the triangle. But obviously, it isn't! Here the condition is wrong (see above).


Can anybody tell me what the slide is about?
haegarr
haegarr
Quote:
Original post by szecs
The slide is the same as the OP's drawing. But that seems plain wrong. Do I miss something?
Well, I came to the same conclusion. Still waiting for somebody to enlighten me ;)
_Sauce_
_Sauce_
Quote:
Original post by szecs
in case of CCW winding (I suggest you to use CCW)

if [(v3-p) cross (v1-p)] dot normal > 0
and [(v1-p) cross (v2-p)] dot normal > 0
and [(v2-p) cross (v3-p)] dot normal > 0

then we have a hit.

I use that for convex polygons (with any number of sides). Notice that you don't need to normalize that way at all, even the normal can be non-normalized.

There is a much simpler way too (google for methods, I can't recall the name), that will use some u,v parameters for checking.
It's much faster in case of triangles, and u,v coordinates can be used to calculate texture lookup coordinates with a very simple formula.
So I suggest you to use that.


A quick implementation of your method, doesn't seem to work. It looks less chaotic than my other attempts but it's still far from correct.

float RTTriangle::intersect(const Rayf &ray){	float distance = mPlane.intersect(ray);		if(distance != std::numeric_limits<float>::infinity())	{				Vector3f intersection(ray.mOrigin + ray.mDir * res.distance),			a(mVerts[0] - intersection),			b(mVerts[1] - intersection),			c(mVerts[2] - intersection),			normal = mPlane.getNormal();		if((c.cross(a)).dot(normal) > 0)		{			if((a.cross(b)).dot(normal) > 0)			{				if((b.cross(c)).dot(normal) > 0)				{					return distance;				}			}		}	}		return std::numeric_limits<float>::infinity();}




[Edited by - _Sauce_ on February 22, 2010 7:26:14 AM]
click
szecs
szecs
Quote:
Original post by haegarr
Quote:
Original post by szecs
The slide is the same as the OP's drawing. But that seems plain wrong. Do I miss something?
Well, I came to the same conclusion. Still waiting for somebody to enlighten me ;)
Sorry, I totally misunderstood your first sentence (haven't slept too much).
szecs
szecs 
bool Poin_On_Poly(int sides,float *point,float *poly_pts, float *norm){	float v1[3], v2[3], n[3];	float skal;	for( int i = 0; i < sides; i++ )	{	v1[0] = poly_pts[  i          *3 + 0] - point[0];		v1[1] = poly_pts[  i          *3 + 1] - point[1];		v1[2] = poly_pts[  i          *3 + 2] - point[2];		v2[0] = poly_pts[((i+1)%sides)*3 + 0] - point[0];		v2[1] = poly_pts[((i+1)%sides)*3 + 1] - point[1];		v2[2] = poly_pts[((i+1)%sides)*3 + 2] - point[2];		vekt_szorz(v1,v2,n);		skal = skal_szorz(n,norm);		if( skal > 0 ) //I have no idea why not skal < 0			return false;	}	return true;}
haegarr
haegarr
If you (_Sauce_) are still looking for a parametric method, then I can contribute this basic method:

Assume a ray in parameter form
r(k) := r0 + k * rt
with its origin and track vectors both given in the same co-ordinate system as the triangle.

The triangle itself with its 3 vertices
{ p0, p1, p2 }
can be expressed in parameter form as
t(u,v) := p0 + u * ( p1 - p0 ) + v * ( p2 - p0 )
where
u,v >= 0 and u+v <= 1

If both have a common point then
r0 + k * rt == p0 + u * ( p1 - p0 ) + v * ( p2 - p0 )
and this equation system of 3 equations needs to be solved for { k, u, v } once. Care must be taken for k because it would become infinity if the ray is parallel to the plane where the triangle lies in. In such case you have to abort further calculations with a "miss", of course.
szecs
szecs
If you use haegarr's method, getting the intersection point's texture coordinates becomes very easy:
t1 = p1_text-p0_text;t2 = p2_text-p0_text;intersection_text = p0_text + t1*u + t2*v;

You have to multiply intersection_text with the sizes of the texture image to obtain the texel.

So I suggest you to use haegarr's method.
My method is good for convex polygons and other kinds of testing: for example picking.
Melekor
Melekor
Quote:
Original post by _Sauce_
now according to the slides in the Princeton lecture, in order for P to be inside the triangle, a and b must both be in the range 0..1. My 3D math isn't the best but I recognised that a is simply the scalar projection of Z onto X, and b is the scalar projection of Z onto Y, where Z = P - v3, X = v2 - v3, and Y = v1 - v3. Correct?


This is where you are going wrong. The a, b in the slide aren't the scalar projections, they are fractions of the triangle side length.

If you interpret a and b as scalar projections instead, then your test needs to be changed to:

0 < a < |Y|
0 < b < |X|

or equivalently,

0 < a/|Y| < 1
0 < b/|X| < 1

also, the slide fails to mention the additional constraint

0 < a/|Y| + b/|X| < 1
_Sauce_
_Sauce_
Okay, so it appears there's some really funky stuff going on inside my raytracer. I've done a bit of investigating and tested on a less complex model (a dodecahedron). I used the surface normals for shading in order to get a better idea of what's going on.



As you can see, the edges of the tris don't line up at all. In fact, all of those pentagons are inverted around the origin. The blue one on top there belongs on the bottom of the model, towards the back. left is right, up is down, and in is out.

Kinda sucky. Now I gotta figure out just why that is. The model data is fine (I've viewed the mesh in meshlab just fine) so somewhere along the line things are getting messed up.
click
_Sauce_
_Sauce_
Quote:
Original post by iMalc
I found this using google, and it looks pretty good, doing just what you appear to want.


As much as I appreciate all the suggestions, I'd prefer to get at least one of these algorithms working correctly before attempting to use a more efficient one. It's like changing multiple variables in a scientific experiment and then trying to figure out what effect each one had.

I'm currently using Szec's method, which isn't quite working (see failed dodecahedron) but after I've done some uni work, I'll draw up a diagram of what the algorithm is doing and work out exactly why it's broken (if it's even the algorithm that's the problem)
click
c1984
c1984
RayTriangle

The problem is that the red angles aren't 90 degrees. To achieve the real u and v coordinates the lines should be parallel to a and b like the green ones.
_Sauce_
_Sauce_
You're right c1984.

I've finally found some material on the barycentric method I was attempting to use before, including a nifty flash demo illustrating it a little more clearly.

http://www.blackpawn.com/texts/pointinpoly/default.html

After all this time googling "ray-triangle intersection", I tried "point in triangle" instead. The number of results increased significantly! :)
click

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.