Sign in to follow this  

Ray-triangle intersection

This topic is 2852 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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?

Share this post


Link to post
Share on other sites
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]

Share this post


Link to post
Share on other sites
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).

Share this post


Link to post
Share on other sites
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?

Share this post


Link to post
Share on other sites
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 ;)

Share this post


Link to post
Share on other sites
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]

Share this post


Link to post
Share on other sites
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).

Share this post


Link to post
Share on other sites
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;
}

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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)

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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! :)

Share this post


Link to post
Share on other sites
Okay, so I'm still getting the weird issue with tris that simply do not line up where they should, despite having tried multiple ray-triangle algorithms. I'm not transforming anything with matrices yet, so where should I be looking for the culprit?

I'm guessing the best place to start would be my ray-plane intersection.

Share this post


Link to post
Share on other sites
Not to be offending, _Sauce_, but: Have you realized that the parametric solution was already mentioned in the 3rd answer, and detailed in the 10th answer? I have the impression that nobody besides szecs has realized that... Also the faults in the cited slide were already annotated in the 3rd answer...
Enough whining.

The thread has the title "Ray-triangle intersection", and that is what the solution proposed by me solves. The method shown by szecs at first (also he mentioned early to use the other one), as well as the method linked to by iMalc, as well as the flash applicated site you've linked to above, all solve the problem of "point in triangle" what is not exactly what you are looking for, because you have no such point. Instead, you have an infinite number of points, namely the ray, for which you want to find the one (if any) that is shared with a triangle. So please go back to the 10th answer. Compare that with the flash thingy and the paper iMalc has linked to, and see both the similarity but also the difference. Thanks.

Share this post


Link to post
Share on other sites
Quote:
Original post by _Sauce_
Okay, so I'm still getting the weird issue with tris that simply do not line up where they should, despite having tried multiple ray-triangle algorithms. I'm not transforming anything with matrices yet, ...

Does this mean that it is guaranteed that both the ray as well as the tris vertices are all given w.r.t. the same co-ordinate system when intersection is tested?

Share this post


Link to post
Share on other sites
Quote:
Original post by haegarr
Not to be offending, _Sauce_, but: Have you realized that the parametric solution was already mentioned in the 3rd answer, and detailed in the 10th answer? I have the impression that nobody besides szecs has realized that... Also the faults in the cited slide were already annotated in the 3rd answer...
Enough whining.

The thread has the title "Ray-triangle intersection", and that is what the solution proposed by me solves. The method shown by szecs at first (also he mentioned early to use the other one), as well as the method linked to by iMalc, as well as the flash applicated site you've linked to above, all solve the problem of "point in triangle" what is not exactly what you are looking for, because you have no such point. Instead, you have an infinite number of points, namely the ray, for which you want to find the one (if any) that is shared with a triangle. So please go back to the 10th answer. Compare that with the flash thingy and the paper iMalc has linked to, and see both the similarity but also the difference. Thanks.

I'm aware of this, however those solutions did not immediately seem to be correct, due to this strange issue I'm having regarding the odd transformations that are being applied to the triangles. I've looked through the mesh loading code, constructors for various classes, intersection tests, and I just cannot see anything that would cause what I'm seeing in my output. Each of the pentagons in that dodecahedron is 3 tris, and they all come out in the correct orientation and position with respect to any neighbouring tris that lie on the same plane. This thread is no longer about ray-triangle or point-in-triangle tests, it's about why the hell or how the hell these transformations could be occurring.


Quote:
Original post by haegarrDoes this mean that it is guaranteed that both the ray as well as the tris vertices are all given w.r.t. the same co-ordinate system when intersection is tested?

As far as I'm aware yes, they are all using the same coordinate system, but I can't be sure. I should be able to allocate a full day to sorting this out the day after tomorrow.

Share this post


Link to post
Share on other sites
I guess you should show us more of your code.

A wild guess without code: you are mixing x-y values somehow.

And by looking at your images: it's clear that you have some kind of transformation (the camera is pointing a bit from up to down).

Did you copy code snippets from somewhere?

Share this post


Link to post
Share on other sites

This topic is 2852 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this