Sign in to follow this  

Ray/Triangle intersection

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

Hi, I've spent the last few hours trying to get the seemingly-simple ray-triangle barycentric intersection routine going in my raytracer with absolutely no luck - and getting some very whacky results. Stepping through the code seems to yield some strange values for 'distance', which should contain the distance along the ray from the origin that the plane intersects at. Given that i am testing against a cube defined from z=24 to 30 and some distances are > 1000.0f im not sure whats going on - especially seeing as more or less the same distance routine is used elsewhere fine to calculate distance to planes. Would really appreciate some fresh eyes on this, its driving me insane :)
float OBJTriangleProxy::testIntersection (Ray &r, float maxDist)
{
	Vector3d O = r.getOrigin();
	Vector3d D = r.getDirection();
	
	Vector3d A = getVert (0);
	Vector3d B = getVert (1);
	Vector3d C = getVert (2);

	//Calculate average normal
	Vector3d N = getNormal(0) + getNormal(1) + getNormal(2);
	N.normalise();

	//Calculate distance from ray to the triangle plane
	float distance = ((O - A).dot (N)) / (D.dot (N));

	//Ray is behind the plane the triangle's plane
	if (distance < 0)
		return -1;
	
	//Work out dominant axis for projection of triangle points into 2d
	int dominant = 0;
	if (abs(N.getX()) > abs(N.getY()))
	{
		if (abs(N.getX() > abs(N.getZ())))
			dominant = 0; //x dominant
		else
			dominant = 2; //z dominant
	}
	else
	{
		if (abs(N.getY() > abs(N.getZ())))
			dominant = 1; //y dominant
		else
			dominant = 2; // z dominant
	}
	
	//Find point where the ray intersects the plane
	Vector3d P = O + D * distance;
	
	//Project points onto appropriate plane
	B = C - A;
	C = B - A;
	switch (dominant)
	{
	case 0:	//Project onto x
		B = Vector3d (B.getY(), A.getZ(), 0);
		C = Vector3d (C.getY(), C.getZ(), 0);
		P = Vector3d (P.getY(), P.getZ(), 0);
		break;
	case 1: //Project onto Y
		B = Vector3d (B.getX(), A.getZ(), 0);
		C = Vector3d (C.getX(), C.getZ(), 0);
		P = Vector3d (P.getX(), P.getZ(), 0);
		break;
	case 2: //Project onto Z
		B = Vector3d (B.getX(), A.getY(), 0);
		C = Vector3d (C.getX(), C.getY(), 0);
		P = Vector3d (P.getX(), P.getY(), 0);
		break;
	}
	
	//Calculate 2d barycentric components
	float u = ((P.getY() * C.getX()) - (P.getX()*B.getY())) / ((B.getY() * C.getX()) - (B.getX() - C.getY()));
	float v = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY()));	

	//Point in triangle?
	if (u >= 0 && v >= 0 && u + v <= 1)
		return distance;

	return -1;
}



Share this post


Link to post
Share on other sites
Don't average the vertex normals. Use the true geometric normal of the triangle.

N = normalize(cross(B-A, C-A));

Also make sure that these lines really do what you expect them to
B = C - A;
C = B - A;

Share this post


Link to post
Share on other sites
Presumbly not directly related to the problem, but: Check that the absolute value of D.dot(N) is greater than some epsilon before computing the distance. If you don't do so and the ray is nearly parallel to the plane, then a division by zero will happen.

Share this post


Link to post
Share on other sites
Oops, herein are 2 inconsistencies I've not seen during my first visit:

float u = ((P.getY() * C.getX()) - (P.getX()*B.getY())) / ((B.getY() * C.getX()) - (B.getX() - C.getY()));
float v = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY()));

Shouldn't the last part of u be like B.getX() * C.getY() or so? Furthurmore, P.getX()*B.getY() appearing in both formulas seems me curious, too.

Share this post


Link to post
Share on other sites
Thanks for the input!
I've fixed the normal calculation to use cross product (although this didn't matter as the test mesh im using uses the same normals for each vertex on a face anyway). Also, I forgot to multiply the distance by negative 1, and the calculation of the second triangle edge just before projection onto an axis-plane was incorrect (thanks pragma!).
Still getting wierd results, although starting to look a bit more normal.

The plane/ray intersection is working correctly now so im sure my problem is in the projection/baycentric part.

//Project points onto appropriate plane
Vector3d B2 = B;
B = C - A;
C = B2 - A;

switch (dominant)
{
case 0: //Project onto x
B = Vector3d (B.getY(), B.getZ(), 0);
C = Vector3d (C.getY(), C.getZ(), 0);
P = Vector3d (P.getY(), P.getZ(), 0);
break;
case 1: //Project onto Y
B = Vector3d (B.getX(), B.getZ(), 0);
C = Vector3d (C.getX(), C.getZ(), 0);
P = Vector3d (P.getX(), P.getZ(), 0);
break;
case 2: //Project onto Z
B = Vector3d (B.getX(), B.getY(), 0);
C = Vector3d (C.getX(), C.getY(), 0);
P = Vector3d (P.getX(), P.getY(), 0);
break;
}

//Calculate 2d barycentric components
float u = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY()));
float v = ((C.getY() * P.getX()) - (C.getX()*P.getY())) / ((B.getX() * C.getY()) - (B.getY() * C.getX()));
//Point in triangle?
if (u >= 0 && v >= 0 && u + v <= 1)
return distance;



Will have more of a play with it tonight :)

[Edited by - skotinator on March 27, 2008 12:48:46 AM]

Share this post


Link to post
Share on other sites
I can't believe i missed that!
Well it is printing out correctly proportioned triangles now, but they are not correctly positioned.
for instance

v 30.0 0.0 20.0
v 15.0 10.0 20.0
v 10.0 0.0 20.0

Is the right shape, but has the right most edge (which should be at x=30) at x=0. If i change to

v 50.0 0.0 20.0
v 15.0 10.0 20.0
v 10.0 0.0 20.0

the triangle is deformed as you'd expect, but the extra length shifts the left hand edge, not the right - which remains fixed at x=0.
Generally, for the XY plane the triangle is rendered into, an extension of X or Y of a certain amount results in the triangle stretching in the opposite direction to accomodate - while keeping the most centered point on the origin (ie x=0 y=0).
I think my barycentric equations are to blame

Share this post


Link to post
Share on other sites
Okay, lets have a deeper look.

Your plane equation looks like
( P - A ) . n = 0
and your ray equation looks like
P(k) := O + k * d
(where O is not a Null 0 but an Ooh ;) ).

That yields in the
( O + k * d - A ) . n = 0
<=>
k = ( ( A - O ) . n ) / ( d . n )
If I made no mistake, then there is already a difference to your implementation: Your distance is negated!

Your calculation of the dominant and intersection point seems me correct.

Then you compute the difference vectors. For any reason you chosed to "exchange" the variables. I would not do so for clarity, but I've followed your decision here:
b := C - A
c := B - A
However, you've forgotten to relate P to A, too. IMHO here a
p := P - A
is missed, isn't it?

Your projection stuff seems me correct.

Now you have to compute u,v from linear combinations. From the 2 possibilities, you've chosen
u * c + v * b = p
or, as the equal equation system
u * cx + v * bx = px
u * cy + v * by = py
Eliminating u,v, resp., and re-ordering terms, let me assume that
u = ( px * by - py * bx ) / ( cx * by - cy * bx )
v = ( px * cy - py * cx ) / ( bx * cy - by * cx )
are the results. Your implementation isn't identical, but it's just an expansion of the fractions by -1/-1, so it doesn't really matter.

Share this post


Link to post
Share on other sites
Haegar!
Thankyou so much for your time explaining that through.
I didn't entirely understand the purpose of calculating the line segments of the triangle edges relative to A, so I didn't think to do the same to the intersection point :)
I was negating distance because I was using (O-A).N not (A-0).N. Guess thats how i saw the formula written to begin with :P
Anyway, here is complete & working, hopefully clearly commented code for anyone else that happens to stumble across this looking for a full implementation of the barycentric intersection technique.

/**
Function to return the distance distance along a ray an intersection with the triangle occurs at
Ray R The ray to compute intersection with
float maxDist The maximum acceptable distance along the ray an intersection can occur at
return float +ve the distance along the ray an intersection occurs at or -ve no intersection
**/

float OBJTriangleProxy::testIntersection (Ray &r, float maxDist)
{
//Get ray origin and direction to simplify equations later on
Vector3d O = r.getOrigin();
Vector3d D = r.getDirection();

//Get triangle vertices
Vector3d A = getVert (0);
Vector3d B = getVert (1);
Vector3d C = getVert (2);

//Calculate average normal
Vector3d N = (B - A).cross (C-A);
N.normalise();

//Calculate distance from ray to the plane the triangle lays in
float distance = ((A - O).dot (N)) / (D.dot (N));

//Ray is behind the plane the triangle's plane
if (distance < 0)
return -1;

//Work out dominant axis for projection of triangle points into 2d
int dominant = 0;
if (abs(N.getX()) > abs(N.getY()))
{
if (abs(N.getX() > abs(N.getZ())))
dominant = 0; //x dominant
else
dominant = 2; //z dominant
}
else
{
if (abs(N.getY() > abs(N.getZ())))
dominant = 1; //y dominant
else
dominant = 2; // z dominant
}

//Find point where the ray intersects the plane
Vector3d P = O + D * distance;

//Obtain points relative to A
B = B - A;
C = C - A;
P = P - A;

//Project vectors onto dominant axis plane
switch (dominant)
{
case 0: //Project onto x
B = Vector3d (B.getY(), B.getZ(), 0);
C = Vector3d (C.getY(), C.getZ(), 0);
P = Vector3d (P.getY(), P.getZ(), 0);
break;
case 1: //Project onto Y
B = Vector3d (B.getX(), B.getZ(), 0);
C = Vector3d (C.getX(), C.getZ(), 0);
P = Vector3d (P.getX(), P.getZ(), 0);
break;
case 2: //Project onto Z
B = Vector3d (B.getX(), B.getY(), 0);
C = Vector3d (C.getX(), C.getY(), 0);
P = Vector3d (P.getX(), P.getY(), 0);
break;
}

//Calculate 2d barycentric components
float u = ((P.getY() * B.getX()) - (P.getX()*B.getY())) / ((C.getY() * B.getX()) - (C.getX() * B.getY()));
float v = ((C.getY() * P.getX()) - (C.getX()*P.getY())) / ((B.getX() * C.getY()) - (B.getY() * C.getX()));

//Point in triangle?
if (u >= 0 && v >= 0 && u + v <= 1)
return distance;

//No intersection
return -1;
}





There is obviously massive potential for precalulation here that i've skipped over, but for the sake of gaining some understanding of the algorithm, it is better to keep it all in one routine to begin with.

[Edited by - skotinator on March 28, 2008 12:53:17 AM]

Share this post


Link to post
Share on other sites

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