Scratchapixel: Intersection

In Ray-Tracing Rendering Technique, a very common and important method is called . It will be used to judge whether the ray will intersect with a shape (like Triangle, Sphere, Plane, Disk, etc.). Even though, we just need to care about Triangle in actual use (we choose to convert complex object into triangle mesh and compute the intersection of a ray with all these triangle). In this post, I’ll try to summary the algorithm to realise for these different shapes. And these algorithm are learned from .Scratchapixel.

1. Represent a Ray

At first, we need to define the ray, so that we can use it to compute the intersection.

And it’s a really simple definition, a point and a direction. Just as:

// the minimum requirement for defining a ray is a position and a direction
Vec3f orig;  //ray origin 
Vec3f dir;  //ray direction (normalized)

And along this ray, we can define the point as: .

class Ray 
{ 
public: 
    Ray(), orig(0), dir(0,0,-1) {} 
    Ray(const Vec3f &o, const Vec3f &d) : orig(o), dir(d) {} 
    // etc.
    ... 
    Vec3f orig; 
    Vec3f dir; 
}; 

2. Finding the hit point in a plane

A plane can be represent as . And can be seen as the components of the normal to the plane () and is the distance from the origin to the plane.

Now we have a ray , and a plane . We can use them to compute the interscted point of the ray and the plane:

$$\begin{align*} & Ax + By + Cz + D = 0 \\ & A(P_x) + B(P_y) + C(P_z) + D = 0 \\ & A(O_x + t*R_x) + B(O_y + t*R_y) + C(O_z + t*R_z) + D = 0 \\ & A*O_x + B*O_y + C*O_z + A*t*R_x + B*t*R_y + C*t*R_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 = - \frac{A*O_x + B*O_y + C*O_z + D}{A*R_x + B*R_y + C*R_z} \\ & t = - \frac{N(A,B,C)\cdot O + D}{N(A,B,C)\cdot R} \end{align*}$$

We get the hit point .

3. Excluding special cases

Before we try to compute the intersection, we need to exclue some special cases, so that the algorithm could be robust enough. If the ray is parallel to the plane, then there will be no intersection actually. And the same as the situation that the plane is behind origin of the ray.

3.1 Parallel

We can use the dot product to judge whether the ray and plane are parallel. If the normal of the plane and the ray are perpendicular. Then their dot product should be zero, and it means the situation of parallel happens.

3.2 Behind

We have computed the value of in to represent the hit point. And we can use the value of to judege. If , it means the plane is behind origin of the ray.

4. Triangle

4.1 Geometric Solution

We can use the position of the hit point to judge whether the ray will intersect inside a triangle or not with just a simple process. As the following picture shows, there are two points and .

We know that the cross product follows the right-hand theorem. It means the direction after cross product will be determined by the order in the two sides of cross product. Such as and will have opposite direction. lies in the left side of , and lies in the right side of . Based on this, if , and has the same direction as the normal of the triangle, it means the point should be inside the triangle.

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

4.2 Möller-Trumbore Solution

Barycentric Coordinates

Barycentric coordinates can be used to express the position of any point located on the triangle plane. It can be inferred as:

$$\begin{align*} P = &\ A + u*AB + v*AC \\ = &\ A + u*(B - A) + v*(C - A) \\ = &\ A + u*B - uA + v*C - v*A \\ = &\ (1-u-v)*A + u*B + v*C \\ = &\ w*A + u*B + v*C \end{align*}$$

And we can know that . If and , it means the point is inside the triangle . We can use the following formulas to compute the barycentric coordinates:

$$\begin{align*} & u = \frac{TriangleCAP_{Area}}{TriangleABC_{Area}} = \frac{||CA \times CP||/2}{||AB \times AC||/2} = \frac{||CA \times CP||}{||AB \times AC||}\\ & v = \frac{TriangleABP_{Area}}{TriangleABC_{Area}} = \frac{||AB \times AP||/2}{||AB \times AC||/2} = \frac{||AB \times AP||}{||AB \times AC||}\\ & w = \frac{TriangleBCP_{Area}}{TriangleABC_{Area}} = \frac{||BC \times BP||/2}{||AB \times AC||/2} = \frac{||BC \times BP||}{||AB \times AC||} \end{align*}$$

Möller-Trumbore Algorithm

Now we use Barycentric coordinates to represent the hit point , so we use the following equations to compute the intersection:

$$\begin{aligned} & \begin{cases} \begin{aligned} P = &\ A + u*AB + v*AC \\ = &\ A + u*(B - A) + v*(C - A) \\ P = &\ O + t*R \end{aligned} \end{cases} \\ & so, \\ & O + t*R = A + u*(B - A) + v*(C - A) \\ & O - A = -t*R + u*(B - A) + v*(C - A) \\ & [-R \quad (B-A) \quad (C-A)] \begin{bmatrix} t \\ u \\ v \end{bmatrix} = O - A \end{aligned}$$

We now know the vector , and we need to compute the value of . There is a technique named Cramer’s rule in mathematics, and it will help us to compute faster in this situation. With the help of Cramer’s rule, we can simplify the above equations, such as:

$$\begin{aligned} t = \frac{[|O-A \quad (B-A) \quad (C-A)|]}{[|-R \quad (B-A) \quad (C-A)|]} \end{aligned}$$

So,

$$\begin{aligned} \begin{bmatrix} t \\ u \\ v \end{bmatrix} & = \frac{1}{[|-R \quad (B-A) \quad (C-A)|]} \begin{bmatrix} |(O-A) \quad (B-A) \quad (C-A)| \\ |-R \quad (O-A) \quad (C-A)| \\ |-R \quad (B-A) \quad (O-A)| \end{bmatrix} \\ & = \frac{1}{(R \times (C-A)) \cdot (B-A)} \begin{bmatrix} ((O-A) \times (B-A)) \cdot (C-A) \\ (R \times (C-A)) \cdot (O-A) \\ ((O-A)) \times (B-A)) \cdot R \end{bmatrix} \\ & = \frac{1}{(R \times E_2) \cdot E_1} \begin{bmatrix} (T \times E_1) \cdot E_2 \\ (R \times E_2) \cdot T \\ (T \times E_1) \cdot R \end{bmatrix} \\ & = \frac{1}{P \cdot E_1} \begin{bmatrix} Q \cdot E_2 \\ P \cdot T \\ Q \cdot R \end{bmatrix} \end{aligned}$$

where , , , and .

Now we can use the above equations to compute . Then we can judge the hit point inside the triangle or not.

5. Sphere

Using the parametric form of ray and sphere, it’s simple to get:

$$\begin{aligned} & \begin{cases} \begin{aligned} P = &\ O + t*D \\ R^2 = &\ |P - C|^2 \end{aligned} \end{cases}, so \\ & |O + t*D - C|^2 - R^2 = 0 \\ & D^2t^2 + 2D(O-C)t + (O-C)^2 = 0 \\ \end{aligned}$$

And this is just a quadratic function, with . We can use the following equations to get result of :

$$\begin{align} & t = \frac{-b\ \pm\ \sqrt{b^2-4ac}}{2a} \\ & \triangle = b^2-4ac \end{align}$$

• When , there is two intersections, like situation ;
• When , there is only one intersection, like situation ;
• When , there is no intersection, like situation .

Notice that we should exlude the situation , because some intersections are actually behind the origin of the ray. And if , it means this intersection is behind the origin of the ray.

Besides, when we really implement the above algorithm, we will use the following equations instead to be more stable:

$$\begin{align} & q = -\frac{1}{2}(b+sign(b)\sqrt{\triangle}) \\ & t_1 = \frac{q}{a} \\ & t_2 = \frac{c}{q} \end{align}$$

6. Plane and Disk

It’s easy to find intersection for Plane and Disk. Just like we find hit point before. We exlude the situation of parallel, then we can compute the intersection for Plane.

And we have one more step for Disk. We just need to compute the distance between the hit point and the center point of the disk, then we compare this distance and the radius of Disk.

7. Box