Is my triangle intersection code correct?

Started by
12 comments, last by IsItSharp 9 years ago

I was wondering if my triangle intersection code is correct?


            Triangle.prototype.intersect = function(ray){
                var v1 = this.t1.sub(ray.direction);
                var v2 = this.t2.sub(ray.direction);
                var n1 = v1.cross(v2);
                n1 = n1.normalize();
                
                var d1 = -ray.origin.dot(n1);
                if((ray.direction.dot(n1) + d1) < 0)
                    return -1.0;
                return d1;
            }

My triangle is defined by 3 Vectors (t1,t2,t3).

Advertisement

I was wondering if my triangle intersection code is correct?

IMHO it isn't. It has at least 3 issues.

It seems that v1, v2 shall be some direction vectors that span the (infinite) plane in which the triangle lies. Then t1, t2, t3 are known to be on that plane, but ray.direction is not related to that plane. Hence computing v1, v2 must not use ray.direction. Now, with 3 positions, you can compute 6 difference vectors

t1 - t2, t1 - t3, t2 - t3 (and their negative variants)

from which you need 2, say

v1 := t1 - t3

v2 := t2 - t3

(This is your 1st issue: Your difference vectors are not correct.)

With the above you can describe the plane as an infinite set of points in space, reachable by any linear combination of v1 and v2 when starting at any known point on the plane. For example

P( a, b ) := t3 + v1 * a + v2 * b

Notice, however, that this works if and only if v1 and v2 are not co-linear itself, so you cannot find a scalar f so that v1 * f = v2. This mean, in other words, that t1, t2, and t3 must not lie on a straight line in space.

Another kind of definition for such a plane is to use not 2 direction vectors inside that plane, but instead the normal of the plane. Well, with 2 difference vectors you can compute the normal as their cross-product:

n1 := v1 x v2

The the plane can be given in the so-called Hesse Normal Form:

( x - t3 ) . n1 == 0

It means "any point x in space belongs to that plane if the distance to the plane is zero".

(This is your 2nd issue: You can use the Hesse Normal Form for computing whether / where the ray hits the plane but not whether / where the ray hits the triangle, because the information about the triangle is lost in that formula. The 3rd issue is that using the Hesse Normal Form in that way is not complete.)

Similar to the plane above, a ray can be understood as the set of points in space when starting at a known position (the ray.origin) and going for any distance along a direction (the ray.direction):

R( c ) := ray.origin + c * ray.direction

What you are now interested in is the set of points in space that are part of both the set of points of the plane and the set of points of the ray. Hence you want to look for any points with the condition

P( a, b ) == R( c )

Obviously you have to solve for the 3 unknowns a, b, and c. Luckily, the above condition is given in 3 dimensional space, hence giving you 3 equations, one for each dimension:

Px( a, b ) == Rx( c )

Py( a, b ) == Ry( c )

Pz( a, b ) == Rz( c )

For your it is sufficient to solve the above linear equation system for a and b. On the going you will see that the chance of a "division by zero" exists. If that occurs then the ray is parallel to the plane, hence there is no intersection (if the ray is distant from the plane) or infinite many intersections (if the ray is on the plane). You have to catch this, of course.

Considering here that the ray is not parallel to the plane, you get a unique solution for a and b, say a' and b'. However, so far you still just know that the ray hits the plane, but you want to know whether the ray hits the triangle. As said above, a and b in the plane formula allow you to reach any point on the plane. So you need to take a closer look at a' and b', especially whether they describe a point within the triangle. Because we have build v1 and v2 from the triangle's corners so that they describe 2 of its legs emanating from t3, any solution with

a' < 0 or a' > 1 or b' < 0 or b' > 1

cannot be part of the triangle. That gives the condition

0 <= a' <= 1 and 0 <= b' <= 1

Thinking further of the above limits, the figure described by them is a rhombus with the desired triangle is one half of (simply chose a' == 1 and try some b' between 0 and 1 to see that). That introduces an extended condition:

… and a' + b' <= 1

EDIT: Other ways of solution may exist as well. The above way is the straight forward one.

P( a, b ) := t3 + v1 * a + v2 * b

How can solve this equation with two unknown values?

IsitSharp: that's not an equation, it's a function definition. What you are to solve is this set of equations:
Px( a, b ) - Rx( c ) = 0
Py( a, b ) - Ry( c ) = 0
Pz( a, b ) - Rz( c ) = 0
Which if I write ray.origin and ray.direction as r0, r1 respectively,

t3x + v1x * a + v2x * b - r0x - r1x * c = 0

(repeat for y, z in place of x)

Which I think gives you the matrix...


[[ v1x, v2x, -r1x ]    [[a],    [[r0x - t3x],
 [ v1y, v2y, -r1y ]  *  [b],  =  [r0y - t3y],
 [ v1z, v2z, -r1z ]]    [c]]     [r0z - t3z]]

(Edit: Haegarr pointed out the negations I forgot to include)

... but it's been seven years since I last did this by hand. I remember how frustrating it was when I first tried to figure this out. I used http://www.mathsisfun.com/algebra/systems-linear-equations-matrices.html as my reference just now.

RIP GameDev.net: launched 2 unusably-broken forum engines in as many years, and now has ceased operating as a forum at all, happy to remain naught but an advertising platform with an attached social media presense, headed by a staff who by their own admission have no idea what their userbase wants or expects.Here's to the good times; shame they exist in the past.

Which I think gives you the matrix...

… almost, except that all of v1 and v2 need to be negated, or else all the other coefficients r0, r1, and t3 need to be negated.

t3x + v1x * a + v2x * b - r0x - r1x * c = 0

So basically if i solve this for a i have something like this:

a = [TERM]

then it would be:

t3x + v1x * [TERM] + v2x * b - r0x - r1x * c = 0

If so:

is there some kind of mathematical function solver out there? I feel no pleasure if i have to do this on my own :(


is there some kind of mathematical function solver out there?

I almost hate to write it, but … Google knows! Try "linear system solver" and you'll get plenty of hits, for example http://wims.unice.fr/wims/wims.cgi.

Is this code maybe correct?


            Triangle.prototype.intersect = function(ray){
                var EPSILON = 0.000001;
                var edge1 = this.t2.sub(this.t1);
                var edge2 = this.t3.sub(this.t1);
                
                var pvec = ray.direction.cross(edge2);
                
                var det = edge1.dot(pvec);
                
                if(det > -EPSILON && det < EPSILON)
                    return -1.0;
                var inv_det = 1.0 / det;
                
                var tvec = ray.origin.sub(this.t1);
                
                var u = tvec.dot(pvec) * inv_det;
                if( u < 0.0 || u > 1.0)
                    return -1.0;
                
                var qvec = tvec.cross(edge1);
                var v = ray.direction.dot(qvec) * inv_det;
                if(v < 0.0 || v > 1.0)
                    return -1.0;
                
                var t = edge2.dot(qvec) * inv_det;
            }


Which I think gives you the matrix...

… almost, except that all of v1 and v2 need to be negated, or else all the other coefficients r0, r1, and t3 need to be negated.

Derp, forgot t3 and r0 were moving back to the other side of the equality, and r1 was negative. Fix'd, thanks.

RIP GameDev.net: launched 2 unusably-broken forum engines in as many years, and now has ceased operating as a forum at all, happy to remain naught but an advertising platform with an attached social media presense, headed by a staff who by their own admission have no idea what their userbase wants or expects.Here's to the good times; shame they exist in the past.

There is a second method to calculate intersection. It's more intuitive and it has his advantages in some cases.
You can calculate intersection with triangle plane if you already have abcd coefficients, then project triangle and intersection point on one of the coordinate planes, calculate (u,v) and check if point is within triangle.

This topic is closed to new replies.

Advertisement