• Advertisement

• ### Popular Now

• 15
• 15
• 11
• 9
• 10
• Advertisement
• Advertisement
• Advertisement

# Ray-triangle intersection

This topic is 2950 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

##### Share on other sites
Advertisement
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

##### 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

##### 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

##### Share on other sites
But that seems plain wrong. Do I miss something?

See: the point successes the test, but is sooo outside.

All points succeed that are inside the blue quad:

#### Share this post

##### Share on other sites
Quote:
 Original post by szecsThe 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

##### Share on other sites
I guess you have to do it in all three directions, but in that case, my method seems a bit faster.

#### Share this post

##### Share on other sites
Quote:
 Original post by szecsin case of CCW winding (I suggest you to use CCW)if [(v3-p) cross (v1-p)] dot normal > 0and [(v1-p) cross (v2-p)] dot normal > 0and [(v2-p) cross (v3-p)] dot normal > 0then 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

##### Share on other sites
Quote:
Original post by haegarr
Quote:
 Original post by szecsThe 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

##### 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

##### Share on other sites

• Advertisement