Ray-Triangle Intersection: Geometric Solution

Conventional Approach To The Ray-Triangle Intersection Test

In the previous paragraphs we have learned how to compute the plane's normal (which is the same as the triangle's normal). Next what we need to find out is the position of point P (for some illustrations we also used Phit), the point of intersection between the plane and the ray.

Figure 1: intersection of a ray and a triangle. The triangle lies in a plane. The value \(t\) is the distance from the ray origin to the intersection point.

Step 1: Finding P

We know that P is somewhere on the ray defined by its origin O and its direction D. Consequently we can write the following equation (1):

$$\scriptsize P=O + tD$$

Where t is distance from the ray origin O to P. To find P we need to find t (figure xx). What else do we know? We know the plane's normal which we have already computed and the plane equation (2) which is (check the chapter on the ray-plane intersection in lesson 7 for more information on this topic):

$$\scriptsize Ax + By + Cz + D = 0$$

Where A, B, C can be seen as the components (or coordinates) of the normal to the plane (\(\small N_{plane}=(A, B, C)\normalsize\)) and D is the distance from the origin (0, 0, 0) to the plane (if we trace a line from the origin to the plane, parallel to the plane's normal). The variables x, y and z are the coordinates of any point that lies on this plane.

We know the plane normal and we know than any of the three triangle's vertices (v0, v1, v2) lie in the plane. It is therefore possible to compute D. Any of the three vertices can be chosen but we use the first one v0 by convention:

float D = vecDot(N, v0); // or if you want to compute the dot product directly float D = N.x * v0.x + N.y * v0.y + N.z * v0.z;

We also know that the point P which is the intersection point of the ray and the plane lies in the plane. Consequently we can substitute P (from equation 1) to (x, y, z) in equation 2 and solve for t (equation 3):

$$\scriptsize\begin{array}{l} P = O + tD\\ Ax + By + Cz + D = 0\\ A * P_x + B * P_y + C * P_z + D = 0\\ A * (O_x + tD_x) + B * (O_y + tD_y) + C * (O_z + tD_z) + D = 0\\ A * O_x + B * O_y + C * O_z + A * tD_x + B * tD_y + C * tD_z + D = 0\\ t * (A * D_x + B * D_y + C * D_z) + A * O_x + B * O_y + C * O_z + D = 0\\ t = -{{A * O_x + B * O_y + C * O_z + D} \over {A * D_x + B * D_y + C * D_z}}\\ t = -{{ N(A,B,C) \cdot O + D} \over {N(A,B,C) \cdot D}} \end{array}\normalsize$$

float t = - (dot(triNormal, rayOrigin) + D) / dot(triNormal, rayDirection);

We have computed t which we can use to compute the position of P:

Point pointOnTriPlane; pointOnTriPlane = rayOrigin + t * rayDirection;

There are two very important cases though we need to look at before we check if the point is inside the triangle.

The Ray And The Triangle Are Parallel

If the ray and the plane are parallel they wont intersect (figure 2). For robustness we need to handle that case if it happens. This is very simple. if the triangle is parallel to the ray direction, it means that the triangle's normal and the ray's direction should be perpendicular.

Figure 2: several situations can occur. The ray can intersect the triangle or miss it. If the ray is parallel to the triangle there is not possible intersection. This situation occurs when the normal of the triangle and the ray direction are perpendicular (and the dot product of these two vectors is 0).

We have learned that the dot product of two perpendicular vectors is 0. If you look at the denominator of equation 3 (the term below the line) we do already compute the dot product between the triangle's normal N and the ray direction D. In fact our code is not very robust because this term can potentially be 0 and we should always catch a possible division by 0. In turns out that when this term is equal to 0 the ray is parallel to the triangle. Since they don't intersect in that case there's therefore no need to compute t anyway. Conclusion: before we compute t in our code, we will compute the result of the term \(\small N \cdot D \normalsize \) first and if the result is 0, the function will return the value false (no intersection).

The Triangle Is "Behind" The Ray

Figure 3: If a triangle is "behind" the ray, it shouldn't be visible. Whenever the value of t computed with equation 3 is lower than 0, the intersection point is located behind the origin of the ray and should be discarded. There is no intersection in that case.

So far we have assumed that the triangle was always in front of the ray. But what happens if the triangle is behind the ray with the ray still looking in the same direction? Normally the triangle shouldn't be visible. Even when the triangle is "behind" the ray, we do get a valid value with equation 3. In that case, t is negative which causes the intersection point to be in the opposition direction than the ray's direction. If we don't catch this "error" though, the triangle will show up in the final frame which is wrong (geometry behind the camera shouldn't be visible). The conclusion of all this is that we need to check for the sign of t before we can decide if the intersection point P is valid or not. If t is lower than 0, the triangle is behind the ray's origin (in regards to the ray's direction) and is not visible. There's no intersection and we can return from the function. If t is greater than 0, the triangle is visible for that ray and we can proceed to the next step.

Step 2: Is P Inside Or Outside The Triangle?

Now that we have found the point P which is the point where the ray and the plane intersection we still have to find out if P is inside the triangle (in which case the ray intersects the triangle) or if P is outside (in which case the rays misses the triangle). Figure 2 illustrates these two cases.

Figure 4:

The solution to this problem is simple and is called the "inside-outside" test. A simple demonstration will help to understand how this test works. Imagine that you have a vector A which is aligned with the x-axis (figure 4). Lets imagine that this vector is actually aligned with one edge of our triangle (the edge defined by the two vertices v0v1). Now, the second edge B, is defined by the vertices v0 and v2 of the triangles as showed in figure 4. Lets compute the cross product of these two vectors. As expected the result is a vector which is pointing in the same direction as the z-axis and the normal of the triangle.

$$\scriptsize\begin{array}{l} A=(1, 0, 0)\\ B=(1, 1, 0)\\ C_x = A_y B_z - A_z B_y = 0\\ C_y = A_z B_x - A_x B_z = 0\\ C_z = A_x B_y - A_y B_x = 1 * 1 - 0 * 1 = 1\\ C = (0, 0, 1) \end{array}$$

Now lets imagine that rather than having the coordinates (1, 1, 0) the vertex v2 has the coordinates (1, -1, 0). In other words we have mirrored its position about the x-axis. If we now compute the cross product AxB' we get the result C'=(0, 0, -1).

$$\scriptsize\begin{array}{l} A=(1, 0, 0)\\ B=(1, -1, 0)\\ C_x = A_y B_z - A_z B_y = 0\\ C_y = A_z B_x - A_x B_z = 0\\ C_z = A_x B_y - A_y B_x = 1 * -1 - 0 * 1 = -1\\ C = (0, 0, -1) \end{array}$$

Figure 5: if P is on the left side of A, the dot product N.C is positive. If P is on the right side (P'), N.C' is negative. The vector C is computed from v0 and P (C=P-v0).

Because the C and N are pointing in the same direction their dot product returns a value greater than 0 (positive). However because C' and N are pointing in opposite directions, their dot product returns a value lower than 0 (negative). What does that test tell us? We know that the point where the ray intersects the triangle and the triangle are in the same plane. We also know from the test we have just made that if a point P which is in the triangle's plane (such as the vertex v2 or the intersection point) is on the left side of vector A, then the dot product between the triangle's normal and vector C is positive (C is the result of the cross product between A and B. In this scenario, A = (v1 - v0) and B  = (P - v0)). However if P is on the right side of A (as with v2') then this dot product is negative. And as you can see in figure 5, the point P is inside the triangle when it is located on the left side of A. To apply the technique we have just described to the ray-triangle intersection problem, we need to repeat the left/right test for each edge of the triangle. If for each one of triangle's edges we find that point P is on the left side of vector C (where C is defined as v1-v0, v2-v1 and v0-v2 respectively for each edge of the triangle), then we know for sure that P is inside the triangle. If the test fails for any of the triangle edges, then P lies outside the triangle's boundaries. This process is illustrated in figure 6.

Figure 6: to find out if P is inside the triangle we can test if the dot product of the vector along the edge and the vector defined by the first vertex of the tested edge and P is positive (meaning P is on the left side of the edge). If P is on the left of all three edges, then P is inside the triangle.

In pseudocode we can write:

vector edge0 = v1 - v0; vector edge1 = v2 - v1; vector edge2 = v0 - v2; vector C0 = P - v0; vector C1 = P - v1; vector C2 = P - v2; if (dot(N, cross(edge0, C0)) > 0 && dot(N, cross(edge1, C1)) > 0 && dot(N, cross(edge2, C2)) > 0) return true; // P is inside the triangle

Lets write a complete ray-triangle intersection test in pseudo code. First we will compute the triangle's normal, then test if the ray and the triangle are parallel. If they are the intersection test fails. If they are not, we compute \(t\) from which we can compute the intersection point P. If the inside-out test succeeds (we test if P is on the left side of each one of the triangle's edges) then the ray intersects the triangle and P is inside the triangle's boundaries. The test is successful.

// // Compute the intersection of ray with a triangle using geometric solution // Input: // points v0, v1, v2 are the triangle's vertices // rayOrig and rayDir are the ray's origin (point) and the ray's direction // Return: // return true is the ray intersects the triangle, false otherwise // bool intersectTriangle(point v0, point v1, point v2, point rayOrig, vector rayDir) { // compute plane's normal vector A, B; A = v1 - v0; B = v2 - v0; // no need to normalize vector N = cross(A, B); // N // // Step 1: finding P // // check if ray and plane are parallel ? float NdotRayDirection = Dot(N, rayDir); if (NdotRayDirection == 0) return false; // they are parallel so they don't intersect ! // compute d parameter using equation 2 float d = dot(N, v0); // compute t (equation 3) float t = -(Dot(N, rayOrig) + d) / NdotRayDirection; // check if the triangle is in behind the ray if (t < 0) return false; // the triangle is behind // compute the intersection point using equation 1 point P = rayOrig + t * rayDir; // // Step 2: inside-outside test // vector C; // vector perpendicular to triangle's plane // edge 0 vector edge0 = v1 - v0; vector VP0 = P - v0; C = cross(edge0, VP0); if (dot(N, C) < 0) return false; // P is on the right side // edge 1 vector edge1 = v2 - v1; vector VP1 = P - v1; C = cross(edge0, VP1); if (dot(N, C) < 0) return false; // P is on the right side // edge 2 vector edge2 = v0 - v2; vector VP2 = P - v2; C = cross(CA, CP); if (dot(N, C) < 0) return false; // P is on the right side; return true; // this ray hits the triangle }

The "inside-outside" technique we have just described works for any convex polygon. Repeat the technique we have used for triangles for each edge of the polygon. Compute the cross product of the vector defined by the two edge's vertices and the vector defined by the first edge's vertex and the point. Compute the dot product of the resulting vector and the polygon's normal. The sign of the resulting dot product determines if the point is on the right or left side of that edge. Iterate through each edge of the polygon. There's no need testing the other edges if one fails passing the test.

Note that this technique can be made faster if the triangle's normal as well as the value D from the plane equation are pre-computed and stored in memory for each triangle of the scene.

Source Code

In the following code we implemented the method described in this chapter to compute the intersection of a ray with a triangle. The constructor of the Triangle class set the position of a default triangle which we can use for quick testings. Change the values of these vertices if you wich to experiment with different triangle shapes. You could also play with the object-to-world matrix to move, rotate or scale this default triangle (transform the triangle's input vertices). The intersect method (line 29) now takes an IsectData as one of the method's parameters (instead a t). This structure (lines 1-5) can be used to store useful information about the intersection (only t for now but we will extend it as we progress in this lesson).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
struct IsectData { float t; }; class Object { public: Object(const Matrix44f &o2w) : objectToWorld(o2w) { worldToObject = objectToWorld.inverse(); std::cerr << objectToWorld << "\n" << worldToObject << std::endl; color = Vec3f(drand48(), drand48(), drand48()); } virtual bool intersect(const Ray<float> &r, IsectData &isectData) const = 0; Matrix44f objectToWorld, worldToObject; Vec3f color; }; class Triangle : public Object { public: Triangle(const Matrix44f &o2w) : Object(o2w) { v0 = Vec3f(-1, -1, 0); v1 = Vec3f( 1, -1, 0); v2 = Vec3f( 0, 1, 0); } bool intersect(const Ray<float> &r, IsectData &isectData) const { Vec3f v0v1 = v1 - v0; Vec3f v0v2 = v2 - v0; Vec3f N = cross(v0v1, v0v2); float nDotRay = dot(N, r.dir); if (dot(N, r.dir) == 0) return false; // ray parallel to triangle float d = dot(N, v0); float t = -(dot(N, r.orig) + d) / nDotRay; // inside-out test Vec3f Phit = r(t); // inside-out test edge0 Vec3f v0p = Phit - v0; float v = dot(N, cross(v0v1, v0p)); if (v < 0) return false; // P outside triangle // inside-out test edge1 Vec3f v1p = Phit - v1; Vec3f v1v2 = v2 - v1; float w = dot(N, cross(v1v2, v1p)); if (w < 0) return false; // P outside triangle // inside-out test edge2 Vec3f v2p = Phit - v2; Vec3f v2v0 = v0 - v2; float u = dot(N, cross(v2v0, v2p)); if (u < 0) return false; // P outside triangle isectData.t = t; return true; } float d; Vec3f v0, v1, v2; };

The download section of this lesson contain the complete souce code which you need to compile with the files from lesson 8 (use instructions for lesson 8 to compile the code).

What's Next?

In this chapter we have presented a technique to compute the ray-triangle intersection test using simple geometry. However there is much more to the ray-triangle intersection test which we haven't considered yet such as whether the ray hits the triangle from the front or from the back. We can also compute what we call the barycentric coordinates of the hit point, which are the two-dimensional coordinates of P in regards to the first vertex of the triangle. These coordinates are necessary for doing things such as applying a texture to the triangle. The next chapter will look into these two topics.

Chapter 3 of 6