Sign in to follow this  
harmless

Yet another Ray-Triangle Intersection Test Question

Recommended Posts

Hi everyone, I already have a function that checks for intersection between a ray and triangle and the pseudo code works like this: 1. derive the plane where the triangle is on 2. check for the ray intersection with the above plane 3. if intersection exists, calculate the intersection points 4. check if the point is inside the triangle Now my problem is on item 2. If ray is parallel to plane AND lies on the plane, and assuming the the origin of ray is not inside the triangle, then the intersection point will be somewhere along the edges of the triangle. Currently, if this is the case, i just perform ray-segment intersection check with ray against each of the 3 edges ofthe triangle. It does work though. however, i am not too sure if this is the best way to approach this problem. I appreciate any advise on this problem. thanks! [Edited by - harmless on December 6, 2004 1:48:50 AM]

Share this post


Link to post
Share on other sites
Just assume that if ray is parallel to plane and lies on that plane, it does not intersect that plane. It's OK for any sane use on computers. Your ray and your plane is already non-perfect.

Even if you'll especially handle that case you'll have 'problems' with precision (if ray "should be" parallel, but only slightly inparallel due to inprecision, it will not work)

In 3D that issue corresponds to handling border of triangle especially.

Share this post


Link to post
Share on other sites
Quote:
Original post by harmless
Yes i am aware of the imprecision issue. but i am concerned about the rare possibility that ray is perfectly parallel to triangle, in which case it hits one of the edges of the triangle.

What you are coding? Raytracer? Collision detection? None of such applications requirs what you want to do. Even more, it is more valid assumption to assume tha parallel rays do not intersect plane at all no matter what. Unless you want to use your 3D code for 2D things, it is totally unnecessary to handle singular cases. You must only check not to divide by zero in intersection code.

If you turn your camera by any extremely small amount, image shouldn't change much, right? But if you somehow placed camera that there is some rays parallel to triangle, any small turn of camera will make that there will be no rays parallel to triangle.

Also, if your triangles is in mesh, parallel ray will intersect neighbouring triangles anyway.

[Edited by - Dmytry on December 6, 2004 6:21:35 AM]

Share this post


Link to post
Share on other sites
Quote:
Original post by harmless
Yes i am aware of the imprecision issue. but i am concerned about the rare possibility that ray is perfectly parallel to triangle, in which case it hits one of the edges of the triangle.


Parallel isn't enough. Must be colinear.

Share this post


Link to post
Share on other sites
just define an intersection to be a point on your ray where the ray goes from one side of a surface to the other one. voila, no conflict.

that way, a colinear ray and triangle arnt intersecting, and as other have mentioned, youre not going to see the slightest bit of difference in the endresult.

the reason why you want to define such cases as no intersection is because its a physicly impossible situation anyway, and hence might give rise to singularities/oddities in further calculations.

Share this post


Link to post
Share on other sites
If this is for a collision system as for a game, technically you'd never be parallel and lie on the plane. Assuming that all objects are closed, you would hit the poly joining the other poly that is parallel to your ray first.


|
A|
|
|_________
B
|
|
|
R


If you have ray R, it's parallel to plane A, but technically would hit the edge of plane B first. If that makes more sense.

Share this post


Link to post
Share on other sites
Guest Anonymous Poster
If I'm not totally wrong then there is a better way to calculate triangle-ray intersection. There is a method that uses barycentric-coordinates to calculate the exact intersection point without calculating the normal of the triangle.

As soon as planes are used in intersection detection then you are heading for trouble. Just look at all the issues with precision when using k-DOPs or all the instable triangle-triangle intersection implementations that are availible on the net. The problem is that a normal n can't span u x v (in theory it does).

/__fold

Share this post


Link to post
Share on other sites
Try the Moller-Trumbore method here, there are other algorithms that *may* give better performance but this one is pretty good.

(some code in OCaml ;)

let closest_ray_intersect tri ray tmin tmax =
let origin = Ray.org ray in
let direction = Ray.dir ray in
let epsilon = 0.000001 in
let edge1 = Vector3.sub tri.p1 tri.p0 in
let edge2 = Vector3.sub tri.p2 tri.p0 in
let pvec = Vector3.cross direction edge2 in
let det = Vector3.dot edge1 pvec in

if (det > ~-.epsilon) && (det < epsilon) then
Intersection.None
else
let inv_det = 1.0 /. det in
let tvec = Vector3.sub origin tri.p0 in
let u = (Vector3.dot tvec pvec) *. inv_det in

if u < 0.0 || u > 1.0 then
Intersection.None
else
let qvec = Vector3.cross tvec edge1 in
let v = (Vector3.dot direction qvec) *. inv_det in

if v < 0.0 || (v +. u) > 1.0 then
Intersection.None
else
let t = (Vector3.dot edge2 qvec) *. inv_det in
if not ((t < tmin) || (t > tmax)) then
Intersection.Hit (t, Vector2.make u v )
else
Intersection.None

Share this post


Link to post
Share on other sites
I posted the most stable method for determining ray-tri here (at the bottom)

note that, triangles which are coplanar with the ray will return no intersection, as there is no unique value t for the intersection (technically, it rejects because the last test is > 0.0f instead of >= 0.0f, which means that it will miss any intersection that falls directly on the t1,t2 edge... I view this as a benefit, however, because for 2/3 of cases in mesh where ray intersects exactly a triangle edge, the test will return true for exactly one triangle, and not both that share the edge).

You can easily change the code to get exactly the behavior you want, by noting that a determinant == 0.0f means it hits an edge. 2 determinants == 0.0f means it hits the vertex shared by the two corresponding edges. All 3 == 0.0f means coplanar, the 'bad news' case.

This is much better than moller trumbore, IMO. I wrote an SSE version of this algorithm that determines ray-tri on batches of indexed triangle lists in 67 cycles per triangle, with no branching except to loop (which is predicted correctly on every iteration except the last).

Share this post


Link to post
Share on other sites
That is interesting, I note though that theres nothing obvious in your method which gives u,v bilinear barycentrics (or some other parameterization), can this be added without breaking the performance? (obviously applications vary, but it's not a fair comparison since the moller-trumbore above gives those coordinates as a by-product).

[Edited by - JuNC on December 6, 2004 3:24:02 PM]

Share this post


Link to post
Share on other sites
I also listed code to calculate t for r0 + t(r1 - r0), which is used to easily calculate the exact intersection point. Do you need u and v to calculate anything other than the intersection point? Also, remember that (in general) the number of tests will be much larger than the number of intersection results, so it is more important, speed-wise, to make the test code fast.

Straying from mathematics and delving into processor design, the moller-trombore "trivial reject test" is not so trivial. By the time you compute (u < 0.0f || u > 1.0f), you've already done a cross product, 2 dots, 3 vector subs, and a scalar-vec multiply. THEN, the 'trivial reject' is taken approximately 2/3 of the time, which is really bad for branch-predictors. The great irony is that the case where moller-trombore might be faster, where the majority of triangles tested intersect successfully, the div lies on the branch-predict path, which means that if the 'trivial reject' does succeed, the pipeline has to wait for the div to finish before it can flush! That's going to be a huge penalty, probably ballpark 50-cycles... so it's really no win unless the % of successful intersection results is up around 80-90%.

I could actually do a benchmark, but it would be pretty meaningless to test an optimized version of one algorithm against the textbook version of a different one.

[edit] And although, I wish I were smart enough to discover something like this, it's not "my" method :) I originally saw it in an article by segura and feito , although I think they made errors in their paper that I describe in another post, and I suspect that this method turns out to be fundamentally equivelent to plucker space... the equations are suspiciously similar.[/edit]

[Edited by - ajas95 on December 6, 2004 3:11:40 PM]

Share this post


Link to post
Share on other sites
Yeah I noticed the t code, but obviously for a raytracer you'll need some parameterization of the triangle as well. As I said, it depends on your application whether you need coordinates on the triangle and comparing an algorithm which doesn't give you them 'for free' and one that does isn't very fair.

If I didn't need the u,v then 67 cycles without branches sounds very attractive, with that requirement then it'll take a bit more convincing - I haven't checked the math so it's possible that some parameterization falls out automatically from the dets and you can just store them.

EDIT: I'm not arguing that the method you gave isn't faster - I'm sure it is. Also your point about not needing u,v for a large number of the tests is taken. In fact, I think I'll actually incorporate the method for my large scale mesh intersections - is your SSE optimized code available or do I have to do it myself? :)

Share this post


Link to post
Share on other sites
Guest Anonymous Poster
ajas95: I think Möller uses his code mainly for ray-tracing. He probably does some sort of Octree first so he increases the possibility of an intersection and this limits his test to the triangles in a few voxels. I guess that he also had in mind the possibility of saving results from previous tests and update does with the same basic transformations as for the objects. Also, if many rays are parallel to each other and then checked against the same triangle then the code before the first rejection test needs to calculated only once.

/__fold

Share this post


Link to post
Share on other sites
JuNC,
what you said about intermediate determinants being re-used, that may very well be true. If you notice, my code to calculate t is just the top equation of moller's intermediate:

[t u v] = 1.0f / det(-d, e1, e2) * [ det(s,e1,e2), det(s,e2,-d), det(-d,e1,s) ]

Also, it's somewhat misleading to say that Moller's method calculates u,v as a "bi-product" of intersection. Mollers method IS to calculate u,v :) do you see what I mean? You could more accurately say that intersection is a biproduct of determining u and v (by checking 0 <= u,v <= 1). Also, just because u and v are calculated sequentially, (i.e. implying that you don't have to calculate v if u is outside a range), that is not how you would actually do it in an optimized version. Dots and crosses each require about 5 instructions, each of which is dependent on the previous one. So you would calculate these things interleaved. By the end, I can already see that you're not saving anything by trying to branch halfway through. You might save yourself a div, but if you understand my note above, that savings might be an illusion.

Anyway, segura's method can also use pre-computed info to speed it up. It's equivalent (I suspect) to storing the plucker coords of each edge. My implementation assumes that you don't have this info, as my mesh changes every frame. If you're doing many tests over the same mesh every frame, then it would be a benefit to calculate the cross product for each vertex (i.e. v0 x v1, v1 x v2 and v2 x v0). Then you would only need to compute 3 dot products to test intersection... hmmm, I should really try that. However your storage requirements skyrocket, it completely negates the benefit of indexed tri-lists. Although you really only need 2 of the cross-products per triangle... Or you really only need one cross-product per edge, but then you'd need to also store directionality of each edge for each triangle.

Alright, enough stream of conciousness :) It would be complicated to do it in a space-efficient manner. If you do it in a non-space-efficient manner, you risk cache-effects, although it becomes computationally very fast.

Hard to say without actually doing it.

Yeah, I can email you the code. The intersection code is commented, but unless you have a thorough understanding of SSE and the algorithm itself, it will be very difficult to follow, as the operations are heavily interwoven.

__fold : when I wrote it, it was with the assumption that there was no precomputed info. It seemed Moller was making the same assumption with his (no precomputed normal). He also seemed aware that branching might not be such a benefit ("However, on modern CPUs where conditional statements tend to be fairly expensive, this test may actually decrease performance" -- Moller, RTR 2nd Ed.)

Share this post


Link to post
Share on other sites
Quote:

Also, it's somewhat misleading to say that Moller's method calculates u,v as a "bi-product" of intersection. Mollers method IS to calculate u,v :) do you see what I mean? You could more accurately say that intersection is a biproduct of determining u and v (by checking 0 <= u,v <= 1). Also, just because u and v are calculated sequentially, (i.e. implying that you don't have to calculate v if u is outside a range), that is not how you would actually do it in an optimized version. Dots and crosses each require about 5 instructions, each of which is dependent on the previous one. So you would calculate these things interleaved. By the end, I can already see that you're not saving anything by trying to branch halfway through. You might save yourself a div, but if you understand my note above, that savings might be an illusion.


Yes, the method IS to calculate u,v - so therefore having them is 'free'. If I don't have them I'll have to do extra calculations to work them out (admittedly I've never bothered checking the performance of doing it this way since higher-level optimizations tend to make more difference in the field I'm interested in).

I'm totally in sync with what you're saying about branching - my point is more that can I get the exact same information out of your algorithm as I can out of the MT method with *less* work (in whatever way we want to define that - branching, FP ops or whatever)?

There was some optimization done by Moller based on feedback from various architectures - he reorganized which tests were done first to avoid pipeline stalls. Of course, those kind of optimizations change from CPU to CPU and between generations of the same CPU.

I can't really precompute much information for triangles because I'm aiming at ray tracing (in various guises like photon mapping and MLT) many-hundred-million face meshes and memory is a massive problem (and caches are a pain).

Share this post


Link to post
Share on other sites
okay, well we can quantify work based on % of tris tested which intersect the ray (this assumes that there is no aggregate benefit to the 'early' exit test 0 <= u <= 1.0f, which I'm fairly convinced of... on the other hand, if you allow his algorithm to run to completion it can be written with no branches also). Moller's method does 2 crosses, 4 dots, 3 scalar * vector muls, 3 subs, and a div. Segura's method does 2 crosses, 3 dots, 4 subs... a difference of 3 muls, a dot and a div (subs are pretty much free :)

So it really depends on how dependencies and registers work out... but the div is going to cost about 30 cycles, 12 or so for the dot, and another 12 for the three muls. This is totally a ball-park guess, but I would guess Mollers method to take an extra 50 cycles or so (remember the big assumption above, however).

So, if m is tests, and n is probability of intersection:


segura vs moller
67m + mn * 120 = 120m

67 + 120n = 120

n = (120 - 67)/120
n = .44


So if your probability of intersection is < 44%, it's probably better to run Segura's test, then on each intersecting result, run moller's test to calculate u,v,t. I'll accept that all of this (except the # of individual ops) is open for debate.

100 of millions? whew. I guess your concern is probably creating very cache-efficient index/vertex storage, and space partitioning.

Share this post


Link to post
Share on other sites
Actually, I just thought of a good optimization for MT to eliminate the div with a branch taken almost everytime... it boils down to calculating s.p and d.q (i.e a*u and a*v), and then checking vs a and 0... i.e.

xmm0 = a*u a*v -a*u -a*v
xmm1 = a a 0 0

cmpps xmm1, xmm0, 1 (less_than)
movmskps eax, xmm1
cmp eax, 1111b
jnz no_intersect

// test u+v < 1
// do div

so, say this might cut MT down to 95 cycles or so (on average, but now it's branching, so it will vary).

Using my math above... that would make the intersection expectation > 27% to make MT worth it. However, because of the branch it becomes MORE expensive when there actually is an intersection :
just a thought

Share this post


Link to post
Share on other sites
thanks to everyone for such an enormous response! seems like i got much more answers than i could ask for.

based on your replies, i come to realize that it does make sense rays collinear (grhodes_at_work you're right!) to triangle should not intersect. in which case answers my main question.

as for other better algorithms such as moller and trumbore and the ones using barycentric coordinates, i'm aware of their existence. i actually have the geometric tools book by eberly/schneider (bought it last week) and these two are explained. however, i still don't understand them yet. i am just taking it one at a time though, using the easier to understand (although less efficient) algorithm first.

with regards to possible issues i may encounter when using the algorithm i am using (planes used for intersection), at this moment i am happy i don't encounter any trouble with my intersection function. however, i'll take this into account as soon as i get into much more complicated applications.

again, many thanks and i appreciate your continued support!

[Edited by - harmless on December 7, 2004 12:48:34 AM]

Share this post


Link to post
Share on other sites
Guest Anonymous Poster
harmless: One thing you have to realize is that even if you insert the origin in that plane equation and find it to be outside the plane then the ray actually might in some situations be coplanar with the plane due to numerical errors. I had that problem with a tri-tri intersection test code.


ajas95: Have you done some robustnest check for your code? The vectors v0, v1 and v2 might become very small in some cases and in others very big depending on the triangles. Doesn't long thin triangles give your code problems? (cause it gives big differences between the length of the v-vectors and then you make a crossproduct between those vectors) The numerical errors in the planes that you use for checking the signed distance against also seems to grow lienarly with the distance between the ray origin and the triangle. But I guess it all depends on what you use the code for and your implementation is nice. :)

I think Möllers original code also had some numerical problems and I believe that others has extended his work to fix the problems with his code.

/__fold

Share this post


Link to post
Share on other sites
Quote:
Original post by Anonymous Poster
harmless: One thing you have to realize is that even if you insert the origin in that plane equation and find it to be outside the plane then the ray actually might in some situations be coplanar with the plane due to numerical errors. I had that problem with a tri-tri intersection test code.


/__fold


what do you mean by inserting the origin in the plane? would you care to explain further? thanks!

Share this post


Link to post
Share on other sites
Guest Anonymous Poster
First you derive the plane where the triangle is on and you wanted to now what to do if the ray is located inside that plane. Now that plane will have some numercial errors in it. I assume you insert the origin of the ray into that plane equation and check if it's inside that plane. If the distance is < epsilon then you assume that the origin of the ray is inside the plane. I was trying to point out that a test like this does only assures you that the origin of the ray is in that plane and it actually doesn't tell you anything about how it's realted to the triangle.

To sum it up, there is no fail-safe way to test if a ray is coplanar with a triangle and it's better as suggested before to not handle this case at all.

Consider this example:
A triangle is located in the xz-plane. The normal should be (0,1,0) but then due to numerical errors let say it's (0.0000001, 1, 0) (using floats).

The origin of a ray starts in (n, 0, 0) (where n is any real number) and it has a direction along the positive x-axis (1,0,0).

Since the plane goes throught the origin then d value is 0 so the plane equation is Pi = 0.0000001x + 1y + 0z. Now insert the origin of the ray. We have 0.0000001 * n + 1 * 0 + 0 * 0 = 0.0000001 * n.

Now in theory, we should have Pi = 0 for any value of n but reality we have an error that grows linearly with the distance between the ray and the triangle.

Finally insert for example the ray origin (-10000, 0, 0) in Pi. We get Pi = -0.001. The abs(-0.001) is surely greater than a predefined epsilon. However, we know that the ray is coplaner with the triangle we had from the start.

This example has only a small error along one axis. The situation is often much worse in reality.

What I was trying to say is that it's more intresting that your code can handle this situation without crashing then finding the intersection for a ray that is "coplanar" with the triangle.

/__fold

Share this post


Link to post
Share on other sites
Quote:
Original post by Anonymous Poster
First you derive the plane where the triangle is on and you wanted to now what to do if the ray is located inside that plane. Now that plane will have some numercial errors in it. I assume you insert the origin of the ray into that plane equation and check if it's inside that plane. If the distance is < epsilon then you assume that the origin of the ray is inside the plane. I was trying to point out that a test like this does only assures you that the origin of the ray is in that plane and it actually doesn't tell you anything about how it's realted to the triangle.

To sum it up, there is no fail-safe way to test if a ray is coplanar with a triangle and it's better as suggested before to not handle this case at all.

Consider this example:
A triangle is located in the xz-plane. The normal should be (0,1,0) but then due to numerical errors let say it's (0.0000001, 1, 0) (using floats).

The origin of a ray starts in (n, 0, 0) (where n is any real number) and it has a direction along the positive x-axis (1,0,0).

Since the plane goes throught the origin then d value is 0 so the plane equation is Pi = 0.0000001x + 1y + 0z. Now insert the origin of the ray. We have 0.0000001 * n + 1 * 0 + 0 * 0 = 0.0000001 * n.

Now in theory, we should have Pi = 0 for any value of n but reality we have an error that grows linearly with the distance between the ray and the triangle.

Finally insert for example the ray origin (-10000, 0, 0) in Pi. We get Pi = -0.001. The abs(-0.001) is surely greater than a predefined epsilon. However, we know that the ray is coplaner with the triangle we had from the start.

This example has only a small error along one axis. The situation is often much worse in reality.

What I was trying to say is that it's more intresting that your code can handle this situation without crashing then finding the intersection for a ray that is "coplanar" with the triangle.

/__fold


i think i'm with you now and have a good understanding about it. so this situation proves that it's best not to handle cases as intersection when ray is coplanar with the plane?

Share this post


Link to post
Share on other sites
Guest Anonymous Poster
Quote:
Original post by harmless
i think i'm with you now and have a good understanding about it. so this situation proves that it's best not to handle cases as intersection when ray is coplanar with the plane?


Well you can never be sure about weather the ray is coplanar to the triangle or not so it's pointless.

/__fold

Share this post


Link to post
Share on other sites

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