Ray Tracing: Rendering a Triangle

News (August, 31): We are working on Scratchapixel 3.0 at the moment (current version of 2). The idea is to make the project open source by storing the content of the website on GitHub as Markdown files. In practice, that means you and the rest of the community will be able to edit the content of the pages if you want to contribute (typos and bug fixes, rewording sentences). You will also be able to contribute by translating pages to different languages if you want to. Then when we publish the site we will translate the Markdown files to HTML. That means new design as well.

That's what we are busy with right now and why there won't be a lot of updates in the weeks to come. More news about SaP 3.0 soon.

We are looking for native Engxish (yes we know there's a typo here) speakers that will be willing to readproof a few lessons. If you are interested please get in touch on Discord, in the #scratchapixel3-0 channel. Also looking for at least one experienced full dev stack dev that would be willing to give us a hand with the next design.

Feel free to send us your requests, suggestions, etc. (on Discord) to help us improve the website.

And you can also donate). Donations go directly back into the development of the project. The more donation we get the more content you will get and the quicker we will be able to deliver it to you.

13 mns read.

Ray-triangle intersection: geometric solution

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.

In the previous paragraphs we learned how to calculate a plane's normal. Next, we need to find out the position of point P (for some illustrations we also used Phit), the point where the ray intersects the plane.

Step 1: Finding P

We know that P is somewhere on the ray defined by its origin \(O\) and its direction \(R\). We used \(D\) in the previous lesson but we will use the term \(R\) in this lesson to avoid confusion with the term \(D\) from the plane equation. The ray parametric equation is (equation 1):

$$P=O + tR.$$

Where \(t\) is distance from the ray origin \(O\) to P. To find P we need to find \(t\) (figure 1). 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 previous lesson for more information on this topic):

Corrected formula (Contribution: Boris the Brave).
$$ Ax + By + Cz + D = 0\\ D = -(Ax + By + Cz) $$

Where A, B, C can be seen as the components (or coordinates) of the normal to the plane (\(N_{plane}=(A, B, C)\)) 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's normal and we know that 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. Let's choose V0:

float D = dotProduct(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):

$$ \begin{array}{l} P = O + tR\\ Ax + By + Cz + D = 0\\ A * P_x + B * P_y + C * P_z + D = 0\\ A * (O_x + tR_x) + B * (O_y + tR_y) + C * (O_z + tR_z) + D = 0\\ A * O_x + B * O_y + C * O_z + A * tR_x + B * tR_y + C * tR_z + D = 0\\ t * (A * R_x + B * R_y + C * R_z) + A * O_x + B * O_y + C * O_z + D = 0\\ t = -{\dfrac{A * O_x + B * O_y + C * O_z + D}{A * R_x + B * R_y + C * R_z}}\\ t = -{\dfrac{ N(A,B,C) \cdot O + D}{N(A,B,C) \cdot R}} \end{array} $$
float t = - (dot(N, orig) + D) / dot(N, dir);

We now have computed \(t\) which we can use to compute the position of P:

Vec3f Phit = orig + t * dir;

There are two very important cases that 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 won't intersect (figure 2). For robustness, we need to handle that case if it happens. This is very simple. if the triangle and the ray are parallel, then 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 no possible intersection. This situation occurs when the normal of the triangle and the ray direction is 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\). Our code is not very robust because this term can potentially be 0 and we should always catch a possible division by 0. It 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 calculate \(t\), we will calculate the result of the term \(N \cdot R\) first and if the result is 0, the function will return 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. Equation 3 returns a valid result, even when the triangle is "behind" the ray. In that case, \(t\) is negative which causes the intersection point to be in the opposite direction to the ray direction. If we don't catch this "error" though, the triangle will show up in the final image which we don't want. Therefore we need to check for the sign of \(t\) before we can decide whether the intersection is valid. If \(t\) is lower than 0, the triangle is behind the ray's origin (with regards to the ray's direction) and is not visible. There is no intersection and we can return false again. If \(t\) is greater than 0, the triangle is "visible" to 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 intersect, 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: C and C' point in opposite directions.

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

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.

The solution to this problem is simple and is also called the "inside-outside" test (we have already used this term in the lesson on rasterization. In the context of rasterization, the test was applied for 2D triangles. We will use it here for 3D triangles. Imagine that you have a vector A which is aligned with the x-axis (figure 4). Let's imagine that this vector is aligned with one edge of our triangle (the edge defined by the two vertices v0-v01). B, the second edge, is defined by the vertices v0 and v2 of the triangles as shown in figure 4. Let's calculate the cross-product of these two vectors. As expected the result is a vector that is pointing in the same direction as the z-axis and as the normal of the triangle.

$$ \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} $$

Let's now 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).

$$ \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} $$

Because 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. You can see in figure 5, that 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 the 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.

Here is an example of the inside-outside test in pseudocode:

Vec3f edge0 = v1 - v0; Vec3f edge1 = v2 - v1; Vec3f edge2 = v0 - v2; Vec3f C0 = P - v0; Vec3f C1 = P - v1; Vec3f C2 = P - v2; if (dotProduct(N, crossProduct(edge0, C0)) > 0 && dotProduct(N, crossProduct(edge1, C1)) > 0 && dotProduct(N, crossProduct(edge2, C2)) > 0) return true; // P is inside the triangle
Note the similarities between this method and the method we used in the rasterization lesson to find if a pixel overlaps a (2D) triangle.

Let's write the complete ray-triangle intersection test routine 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 parallel, 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.

bool rayTriangleIntersect( const Vec3f &orig, const Vec3f &dir, const Vec3f &v0, const Vec3f &v1, const Vec3f &v2, float &t) { // compute plane's normal Vec3f v0v1 = v1 - v0; Vec3f v0v2 = v2 - v0; // no need to normalize Vec3f N = v0v1.crossProduct(v0v2); // N float area2 = N.length(); // Step 1: finding P // check if ray and plane are parallel. float NdotRayDirection = N.dotProduct(dir); if (fabs(NdotRayDirection) < kEpsilon) // almost 0 return false; // they are parallel so they don't intersect ! // compute d parameter using equation 2 float d = -N.dotProduct(v0); // compute t (equation 3) t = -(N.dotProduct(orig) + 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 Vec3f P = orig + t * dir; // Step 2: inside-outside test Vec3f C; // vector perpendicular to triangle's plane // edge 0 Vec3f edge0 = v1 - v0; Vec3f vp0 = P - v0; C = edge0.crossProduct(vp0); if (N.dotProduct(C) < 0) return false; // P is on the right side // edge 1 Vec3f edge1 = v2 - v1; Vec3f vp1 = P - v1; C = edge1.crossProduct(vp1); if (N.dotProduct(C) < 0) return false; // P is on the right side // edge 2 Vec3f edge2 = v0 - v2; Vec3f vp2 = P - v2; C = edge2.crossProduct(vp2); if (N.dotProduct(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 edges' 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 to test the other edges if one fails to pass 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 precomputed and stored in memory for each triangle of the scene.

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 the back. We can also compute what we call the barycentric coordinates. These coordinates are necessary for doing things such as applying a texture to the triangle.