# Is my triangle intersection code correct?

## Recommended Posts

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).

##### Share on other sites

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

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.

Edited by haegarr

##### Share on other sites

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

How can solve this equation with two unknown values?

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

Edited by Wyrframe

##### Share on other sites

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.

Edited by haegarr

##### Share on other sites

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 :(

##### Share on other sites

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.

##### Share on other sites

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;
}


##### Share on other sites

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.

##### Share on other sites

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.

##### Share on other sites

In fact there is a lot of triangle intersection code algorithms, out of all, most interesting are:

which is just a standard either single-side or double-side barycentric test. It doesn't need any additional storage and in general is very fast on GPUs and CPUs.

I find it very fast on the CPU. It needs more storage compared to MT test (plus for shading you will need original vertex coordinates), so I find it to behave a bit worse (on some of the GPUs) compared to MT.

Out of others - Plucker test is one of the most robust, Badouel (projection) test behave very fast on SIMD SSE on CPU, and so on. I'm also linking to (a bit older) good summary of algorithms and their comparsions in terms of speed - http://gggj.ujaen.es/docs/grapp14_jimenez.pdf

I hope this can help. I could link you to my implementation of gpu-raytracer on GitHub, but that repo is heavily WIP and I haven't pushed for few weeks (I have modifications locally), plus the setup is a bit harder as of now (but I'm working on it!). https://github.com/Zgragselus/OpenTracer ... there is obsolete CPU code with variation of barycentric test and woop test in OpenCL (which is used in runtime).

##### Share on other sites

Hi Vilem,

i tried to convert your version of triangle intersection from here: https://github.com/Zgragselus/OpenTracer/blob/master/OpenTracer/Math/Intersection/Intersection.h to my code:

            Triangle.prototype.intersect = function(ray){
var e1 = this.t2.sub(this.t1);
var e2 = this.t3.sub(this.t1);

var pvec = ray.direction.cross(e2);
var det = e1.cross(pvec);
if(det > -0.001 && det < 0.001)
return -1.0;

var inv_det = 1.0 / det;
var tvec = ray.origin.sub(this.t1);

var bx = tvec.dot(pvec) * inv_det;
if(bx < 0.0 || bx > 1.0)
return -1.0;

var qvec = tvec.cross(e1);
var by = ray.direction.dot(qvec) * inv_det;
if(by < 0.0 || (bx+by) > 1.0)
return -1.0;

var d = e2.dot(qvec) * inv_det; //distance?

var bz = 1.0 - bx - by;
var bw = 0.0;

return d;
}


For my algorithm i need the distance from ray.origin to the intersection point. Is d the right value to return?

##### Share on other sites

Yes, 'd' is the distance from camera to hitpoint.

Also, I'll do a little shameless self-promotion; some time ago I've written an article about intersection tests and their derivations here on GD - http://www.gamedev.net/page/resources/_/technical/math-and-physics/intersection-math-algorithms-learn-to-derive-r3033, where I focused mainly to show the way how to derive some different intersection and give out an idea how to derive your own (of course it is not always as trivial and in some cases not even solvable (surfaces defined by differential functions or such) ~ just using numeric approximation) .. for any of the primitives, the derivation is straight-forward (which is good!).

Also note, please don't try to run the code from that repo, it won't work (sadly I don't have time to keep it up to date right now, I've finally managed to pull some parts of the code that were license-protected by the company I've been working with ~ so I'm free to release it, yet I can't find time to re-format it into some useful library).

##### Share on other sites

Is this code maybe correct?

            Triangle.prototype.intersect = function(ray){
var edge1 = this.t2.sub(this.t1);
var edge2 = this.t3.sub(this.t1);

var s1 = ray.direction.cross(edge2);
var divisor = s1.dot(edge1);
if(divisor == 0.0)
return -1.0;

var invDivisor = 1 / divisor;
var distance = ray.origin.sub(this.t1);
var barycCoord1 = distance.dot(s1) * invDivisor;
if(barycCoord1 < 0.0 || barycCoord1 > 1.0)
return -1.0;

var s2 = distance.cross(edge1);
var barycCoord2 = ray.direction.dot(s2) * invDivisor;
if(barycCoord2 < 0.0 || (barycCoord1 + barycCoord2) > 1.0)
return -1.0;
var intersectionDistance = edge2.dot(s2) * invDivisor;
return intersectionDistance;

}

I think there is an error somehwere:

If i use the triangles as a light i get this:

(i tried to build 4 walls around the sphere i don't know if my coordinates are correct).

Edited by IsItSharp

## Create an account

Register a new account

• ### Forum Statistics

• Total Topics
628383
• Total Posts
2982368

• 10
• 9
• 15
• 24
• 11