14 replies to this topic

[Grumple]

I'm currently implementing a system to fire many rays against a triangle mesh, calculating intersections. Everything works great in general, except for the occasional case of rays firing between adjoining triangles without detecting a hit on either. I'm only using floats for all the math so I know my precision isn't as high as it could be, but the triangles are all using shared vertices, etc so they should be nicely locked together without 'real' gaps.

Can anyone suggest a nice way to handle this and detect this case where my rays would be striking the 'edge' of the triangle? Performance is a a big consideration. Cheers!

[Waterlimon]

I guess you should post the math for testing the intersection so someone can suggest a better algorithm or add something to it to make it work properly.

[Grumple]

Good call, Waterlimon, I was hoping for a general 'this algorithm handles edges nicely', but if anyone wants to suggest an improvement to my existing algorithm, thar she blows:

//return values:
// -1 triangle is degenerate, no processing done
// 0 no intersection
// 1 valid intersection
// 2 ray lies on triangle plane
//rayOrigin		 = x,y,z floats
//rayDir			= normalized x,y,z floats
//srcTriangleVert*  = x,y,z floats
//triNormal		 = x,y,z triangle normal
//destIntersectData = x,y,z intersect + distance to intersect. (< 0 if ray was headed in wrong direction)
int intersect_RayTriangleCPU_Native( float *rayOrigin, float *rayDir, float *srcTriangleVertA, float *srcTriangleVertB,
float *srcTriangleVertC, float *triNormal, float *destIntersectData )
{
	float	u[3], v[3];  // triangle vectors
	float	w0[3], w[3]; // ray vectors
	float	r, a, b;	 // params to calc ray-plane intersect
	float uu, uv, vv, wu, wv, D;
	float s, t;

	// get triangle edge vectors and plane normal
	u[0] = srcTriangleVertB[0] - srcTriangleVertA[0];
	u[1] = srcTriangleVertB[1] - srcTriangleVertA[1];
	u[2] = srcTriangleVertB[2] - srcTriangleVertA[2];

	v[0] = srcTriangleVertC[0] - srcTriangleVertA[0];
	v[1] = srcTriangleVertC[1] - srcTriangleVertA[1];
	v[2] = srcTriangleVertC[2] - srcTriangleVertA[2];

	w0[0] = rayOrigin[0] - srcTriangleVertA[0];
	w0[1] = rayOrigin[1] - srcTriangleVertA[1];
	w0[2] = rayOrigin[2] - srcTriangleVertA[2];

	//manually compute dot product and negate result
	a = -((triNormal[0] * w0[0]) + (triNormal[1] * w0[1]) + (triNormal[2] * w0[2]));
	//manually compute dot product
	b = ((triNormal[0] * rayDir[0]) + (triNormal[1] * rayDir[1]) + (triNormal[2] * rayDir[2]));

	if (fabs(b) < SMALL_NUM)	  // ray is parallel to triangle plane
	{
		if (a == 0.0)			 // ray lies in triangle plane
		{
			return 2;
		}
		else
		{
			return 0;			 // ray disjoint from plane
		}
	}

	// get intersect point of ray with triangle plane
	r = a / b;
	if (r < 0.0f)				   // ray goes away from triangle
	{
		return 0;				  // => no intersect
	}

	// for a segment, also test if (r > 1.0) => no intersect (we are firing raw rays, no endpoint)
	destIntersectData[0] = rayOrigin[0] + (rayDir[0] * r);  // intersect point of ray and plane
	destIntersectData[1] = rayOrigin[1] + (rayDir[1] * r);  // intersect point of ray and plane
	destIntersectData[2] = rayOrigin[2] + (rayDir[2] * r);  // intersect point of ray and plane

	//hijack the 4th element to return distance to intersection
	destIntersectData[3] = r;

	// is I inside srcTriangle?
	//manually compute dot products	  
	uu = (u[0] * u[0]) + (u[1] * u[1]) + (u[2] * u[2]);
	uv = (u[0] * v[0]) + (u[1] * v[1]) + (u[2] * v[2]);
	vv = (v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]);

	w[0] = destIntersectData[0] - srcTriangleVertA[0];
	w[1] = destIntersectData[1] - srcTriangleVertA[1];
	w[2] = destIntersectData[2] - srcTriangleVertA[2];

	wu = (w[0] * u[0]) + (w[1] * u[1]) + (w[2] * u[2]);
	wv = (w[0] * v[0]) + (w[1] * v[1]) + (w[2] * v[2]);
	D = uv * uv - uu * vv;

	// get and test parametric coords
	s = (uv * wv - vv * wu) / D;
	if (s < 0.0 || s > 1.0)		// I is outside T
	{
		return 0;
	}
	t = (uv * wu - uu * wv) / D;
	if ((t < 0.0) || ((s + t) > 1.0))  // I is outside T
	{
	   return 0;
	}

	return 1;					  // I is in T
}

[The_Neverending_Loop]

Just so I can understand correctly, Basically you have an object lets say its a model, and for the sake of argument lets say its a door.

This door is made up of a bunch of triangles, and your issue is that your code works correctly when it comes to picking seperate triangles, but sometimes when moving your mouse over the door, there are moments when the picking is not detected, and what you want is that whenever your mouse is over any part of that door, the door infact is recognized as being "highlighted" over the mouse.

Is this correct?

[Grumple]

Just so I can understand correctly, Basically you have an object lets say its a model, and for the sake of argument lets say its a door.

This door is made up of a bunch of triangles, and your issue is that your code works correctly when it comes to picking seperate triangles, but sometimes when moving your mouse over the door, there are moments when the picking is not detected, and what you want is that whenever your mouse is over any part of that door, the door infact is recognized as being "highlighted" over the mouse.

Is this correct?

Yeah exactly. In my case I have a triangle mesh such as a door, and I'm doing something very similar to picking in order to check where a ray intersects the mesh. Each triangle is tested individually. Normally it works perfectly, but as my picking ray crosses over two triangles of the door, it sometimes manages to pass between the two triangles, generating no hit with either. The triangles share vertices, so they don't really have gaps beyond the margin of error from float precision.

I've noticed that barycentric coord based solutions can technically detect edge cases such as this, but I wonder if in reality they have the same issues with float precision? I'm really hoping to avoid 'fudging' the numbers to make my triangles bigger than they are as a solution.

[The_Neverending_Loop]

I honestly haven't taken a thorough look at your code, but when you get the intersection point, how are you calculating if the point falls inside the triangle?

For example the method I use is that i calculat the angles that the point makes with respect to the 3 points of the triangle, if the sum of the angles it makes equals 2*pi or 360 degrees, its inside the triangle.

Is this what your using? or another method?

[Cornstalks]

I did something very similar in my ray tracer for ray-triangle intersection. For what it's worth, I got the ray-plane intersection from Wikipedia and the triangle intersection from SoftSurfer. While I used pretty much the same method that Grumple posted, I'll post my code just to give another perspective/implementation that will hopefully help it be a little more clear:

struct Triangle
{
    Vector3f a, b, c;
    Vector3f normal; // normalized


    bool intersects(const Ray& ray, float& t, Vector3f& n) const
    {
        // from http://softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm (and Plane.intersect())
        n = normal;

        float num = normal.dot(a - ray.start);
        float den = normal.dot(ray.direction);

        if (num != 0 && den == 0)
        {
            return false; // parallel, but never intersect
        }
        else if (num == 0 && den == 0)
        {
            t = 0; // parallel and in same plane
            return true;
        }
        else
        {
            t = num / den;
        }

        if (t < 0)
        {
            return false; // intersection point is behind the ray's start
        }

        Vector3f u = b - a;
        Vector3f v = c - a;
        Vector3f w = ray.start + ray.direction * t - a;

        float uv = u.dot(v);
        float wv = w.dot(v);
        float vv = v.squaredLength();
        float wu = w.dot(u);
        float uu = u.squaredLength();

        float s1 = (uv * wv - vv * wu) / (uv * uv - uu * vv);
        float t1 = (uv * wu - uu * wv) / (uv * uv - uu * vv);

        return s1 >= 0 && t1 >= 0 && (s1 + t1) <= 1;
    }
};

It's important to be aware that this degenerates when the ray lies in the triangle's plane, because I simply set t = 0, when really t = 0 could be outside the triangle. You could optionally return false if t = 0 and since a triangle is infinitely thin as well. It just depends on how you want to handle this situation.

[Grumple]

Thanks for posting that, Cornstalks. Your implementation is essentially identical (although admittedly more readable =P) than mine. Referring to your code, in situations where I should be detecting intersection with the triangle, I'm getting values of s1 along the lines of -0.0008, or (s1 + t1) of 1.0002.

I guess I could accept fudging my boundaries a bit if I could easily relate s1 and t1 to the units of measurement the overall world adheres to. For example if I can guarantee I am adding a set number of centimeters to my triangle edges as opposed to a fraction of the unknown triangle size.

Have you run into this issue with your implementation?

[jefferytitan]

That's a strange one, and I'm not sure how this edge case is happening. Can you characterise the nature of the bug, e.g. find an angle between triangles or camera angle where this happens repeatably? You could perhaps create a test harness to try a bunch of options with known outcomes (e.g. two triangles meeting at a known point, loop through different angles and test points) until you find a situation where it occurs.

One slightly unrelated question, do you actually *need* to ray test against a raw mesh? Often games use simpler models for collision etc for speed reasons. A simpler object may yield faster code, e.g. ray testing against transformed cubes, spheres, cylinders, etc is very easy.

[Cornstalks]

With really small errors like that, it could very easily be floating point rounding errors. Try replacing the intersection code with double instead of float and see if the error is smaller (this is just to make sure it is indeed floating point accuracy errors).

I'm going to go ahead and assume it is indeed floating point rounding issues, and to be honest there isn't a fool-proof way around this. One commonly used solution would be to use some epsilon to allow for a controllable degree of error (which is similar to what you mentioned). For example:

float epsilon = 0.001f;
return (s1 >= 0.0f - epsilon) && (t1 >= 0.0f - epsilon) && ((s1 + t1) <= 1.0f + epsilon);

Another thing would be to look at your compiler's optimization settings. There are usually some flags that can be set in the compiler that will trade accuracy for speed (usually I've found it to be accuracy is favored over speed by default, and you have to explicitly say to favor speed over accuracy). Double check them and make sure you aren't sacrificing accuracy for speed.

I've never had precision errors this small cause any problems for me though, so I can't give any more solutions. It also depends on your data. If you're mixing large numbers and small numbers in calculations, the errors will be much more noticeable. Try to make sure your floats don't differ by too many orders of magnitude if you want more precision.

[Grumple]

With really small errors like that, it could very easily be floating point rounding errors. Try replacing the intersection code with double instead of float and see if the error is smaller (this is just to make sure it is indeed floating point accuracy errors).

I'm going to go ahead and assume it is indeed floating point rounding issues, and to be honest there isn't a fool-proof way around this. One commonly used solution would be to use some epsilon to allow for a controllable degree of error (which is similar to what you mentioned). For example:

float epsilon = 0.001f;
return (s1 >= 0.0f - epsilon) && (t1 >= 0.0f - epsilon) && ((s1 + t1) <= 1.0f + epsilon);

Another thing would be to look at your compiler's optimization settings. There are usually some flags that can be set in the compiler that will trade accuracy for speed (usually I've found it to be accuracy is favored over speed by default, and you have to explicitly say to favor speed over accuracy). Double check them and make sure you aren't sacrificing accuracy for speed.

I've never had precision errors this small cause any problems for me though, so I can't give any more solutions. It also depends on your data. If you're mixing large numbers and small numbers in calculations, the errors will be much more noticeable. Try to make sure your floats don't differ by too many orders of magnitude if you want more precision.

Thanks a lot for the tips! I actually went off and did pretty much exactly what you are suggesting yesterday, so I had a good chuckle when I signed in and read your post this morning. The issue is definitely related to floating point precision. For anyone else who runs into this issue, I've found that going to double precision data helps immensely.

I am mixing somewhat large and small numbers. My values tend to be in the thousands of whole units with accuracy mattering down to 3 decimal places.

To really get a handle on things I implemented two versions of the number crunching algorithm....one using all floats (floats in, floats for intermediate calc results, and floats out). My input/ouput data cant be doubles due to storage/framework constraints, and I didn't want to add the overhead of transferring the input float data into doubles before running calcs. So for the 'high precision' version I just used floats in, and floats out, but with doubles for all intermediate calculation results. The results are cast down to floats for return only after all calcs are completed.

In my environment, I'm finding I need a high epsilon for floats, but not toooo bad. Allowing +/- of 0.001f on either side of the normal range seems to eliminate instances of rays shooting through triangles. However, I start to see occasional misses again if I try to go another decimal place down (ie +/- 0.0001f).

The double precision implementation is quite nice but with one obvious drawback. It does seem to add add an average of 10-15% to my execution time. However, using the same data set it allows my epsilon to go down to at least +/- 0.00001 without generating any misses. I didn't try any smaller values as this is well within my accuracy requirement.

Keep in mind the above results probably do vary depending on the data being used. As Cornstalks suggests, mixing very large and very small values would like require the double precision intermediate calculations to be useful.

[Grumple]

That's a strange one, and I'm not sure how this edge case is happening. Can you characterise the nature of the bug, e.g. find an angle between triangles or camera angle where this happens repeatably? You could perhaps create a test harness to try a bunch of options with known outcomes (e.g. two triangles meeting at a known point, loop through different angles and test points) until you find a situation where it occurs.

One slightly unrelated question, do you actually *need* to ray test against a raw mesh? Often games use simpler models for collision etc for speed reasons. A simpler object may yield faster code, e.g. ray testing against transformed cubes, spheres, cylinders, etc is very easy.

Sorry, I didnt notice this post earlier. To originally narrow the issue down I did something very similar to what you suggest. I basically created an automatic ray generation system that would sweep over a mesh repeatedly until I got a ray that didnt hit any triangles. Then it would just highlight that ray and focus the camera on where it crossed the mess. It quickly became clear that regardless of just about all other variables, the misses always happened for rays crossing the seams between triangles.

You are completely right about generalized collisions for games, but this is actually more of a science project....the exact ray/mesh collision points are critical. =)

Cheers!

[Cornstalks]

My values tend to be in the thousands of whole units with accuracy mattering down to 3 decimal places.

That could be an issue right there. 32-bit IEEE754 floats (i.e. normal floats) are only accurate to about 7 significant digits. So if you have 1234.567, the 7 there is is getting right on the border of the float's accuracy, and it may or may not be able to properly represent it. Start doing math with a few of these numbers as floats, and you're almost guaranteed that for the result xxxx.xxy, the y will be noticeably off.

[Grumple]

My values tend to be in the thousands of whole units with accuracy mattering down to 3 decimal places.

That could be an issue right there. 32-bit IEEE754 floats (i.e. normal floats) are only accurate to about 7 significant digits. So if you have 1234.567, the 7 there is is getting right on the border of the float's accuracy, and it may or may not be able to properly represent it. Start doing math with a few of these numbers as floats, and you're almost guaranteed that for the result xxxx.xxy, the y will be noticeably off.

Yeah that makes a lot of sense, Cornstalks...this is certainly pushing the limits of floating point accuracy. I guess the real surprise for me is how quickly said accuracy falls off past that 7th digit. o_O

[jefferytitan]

My values tend to be in the thousands of whole units with accuracy mattering down to 3 decimal places.

That could be an issue right there. 32-bit IEEE754 floats (i.e. normal floats) are only accurate to about 7 significant digits. So if you have 1234.567, the 7 there is is getting right on the border of the float's accuracy, and it may or may not be able to properly represent it. Start doing math with a few of these numbers as floats, and you're almost guaranteed that for the result xxxx.xxy, the y will be noticeably off.

Yeah that makes a lot of sense, Cornstalks...this is certainly pushing the limits of floating point accuracy. I guess the real surprise for me is how quickly said accuracy falls off past that 7th digit. o_O

Yeah, it does drop off very quickly, particularly with multiplication or division, or addition/subtraction if the exponents are very different. How about using floats for coordinates, but using doubles for calculations in the intersection tests? Using doubles for both would likely introduce more subtle bugs, e.g. loss of precision when combining two doubles that use every possible digit of the mantissa.

Also if possible work with coordinates close to (0,0,0). Both floats and doubles suffer from precision loss when dealing with numbers with a large absolute value.