Ray-Casting algorithm in Ray-triangle intersection

Question

Currently I started studying about ray-casting when I came across this following problem based on ray-triangle intersection. The problem was:

You are provided with a triangle with vertices ,(x1,y1,z1), (x2,y2,z2) and (x3,y3,z3). A ray with origin (a1,b1,c1) and direction (a2,b2,c2) is also given. Your task is to find:

1)Whether or not the ray intersects the triangle

2) If the ray intersects the triangle, what's the point of intersection and also find the distance of theat point from the origin of the ray

Examples:

1) For (x1,y1,z1)=(-2,2,6), (x2,y2,z2)=(2,2,6), (x3,y3,z3)=(0,-4,6) , (a1,b1,c1)=(1,0,0),(a2,b2,c2)=(-0.2,0,1) the answer is: coordinates of intersection:(-0.2,0,6) and the distance of the origin of the ray from the point of intersection is 6.12

2) For (x1,y1,z1)=(-2,2,1), (x2,y2,z2)=(2,2,1), (x3,y3,z3)=(0,-4,1) , (a1,b1,c1)=(0,0,0),(a2,b2,c2)=(0,0,1) the answer is: coordinates of intersection:(0,0,1) and the distance of the origin of the ray from the point of intersection is 1

3) For (x1,y1,z1)=(-10,-2.3,0), (x2,y2,z2)=(4.4,20.3,9.5), (x3,y3,z3)=(9.8,-10,0) , (a1,b1,c1)=(0,0,0),(a2,b2,c2)=(0.68,-1.14,1.82) the answer is: coordinates of intersection:(0.67, -1.12, 1.79) and the distace of the origin of the ray from the point of intersection is 2.22

All I need is to understand the general equations to solve this question. I need to know how the equations are formed and how the problem is solved. This seems a pretty hard question for me.

Answer 1

The way this is done in most raytracers and raycasters is actually quite simple. The intersection distance (from ray origin) is calculated first, and it is typically calculated for several objects for the same ray, so that the objects can be handled in order of increasing distance (unless fully reflected). Then, the intersection point on the relevant surface(s) is parametrized if needed, for example, for shading purposes.

The direction vector for the ray is always an unit vector (length 1), so that it can be parametrized as $$\vec{p}(t) = \vec{p}_0 + t \hat{d}$$ with $\lVert\hat{d}\rVert = 1$. In your case, $$\vec{p}_0 = (a_1, \; b_1, \; c_1)$$ and $$\hat{d} = \left(\frac{a_2}{\sqrt{a_2^2 + b_2^2 + c_2^2}}, \; \frac{b_2}{\sqrt{a_2^2 + b_2^2 + c_2^2}}, \; \frac{c_2}{\sqrt{a_2^2 + b_2^2 + c_2^2}}\right)$$ If $\hat{d}$ were not an unit vector, then $t$ would be the distance in units of the length of the direction vector.

For each planar surface, the four components of the plane are stored separately; the equation for the plane being $$x n_x + y n_y + n_z z = n_d$$ where $\hat{n} = (n_x, n_y, n_z)$ is the plane normal (typically pointing outwards for one-sided surfaces), preferably an unit vector to make the math simpler, and so that $n_d$ is the minimum distance from origin. In your case, $$\vec{n} = (x_2 - x_1, y_2 - y_1, z_2 - z_1) \times (x_3 - x_1, y_3 - y_1, z_3 - z_1)$$ where $\times$ represents a vector cross product; the plane unit normal is then $$\hat{n} = \frac{\vec{n}}{\lVert\vec{n}\rVert} = \frac{\vec{n}}{\sqrt{\vec{n} \cdot \vec{n}}}$$ and $n_d$ is the dot product between $\hat{n}$ and any point on the plane, usually calculated using one of $$\begin{array}{rl} n_d \; =& \hat{n} \cdot (x_1 - a_1, y_1 - b_1, z_1 - c_1) \\ \; =& \hat{n} \cdot (x_2 - a_1, y_2 - b_1, z_2 - c_1) \\ \; =& \hat{n} \cdot (x_3 - a_1, y_3 - b_1, z_3 - c_1) \end{array}$$ Note that with numerical calculations using floating-point values the three may not be exactly the same value due to rounding.

When $\vec{p}_0$, $\hat{d}$, $\hat{n}$, and $n_d$ are known, we can easily solve the (signed) distance $t$ at which the ray intersects the plane: $$t = \frac{n_d - \vec{p}_0 \cdot \hat{n}}{\hat{d} \cdot \hat{n}}$$ Note that if $\hat{d}$ and $\hat{n}$ are perpendicular, the plane is parallel to the ray, and the denominator is zero. If $t \lt 0$, the intersection occurs "before" the ray starts, and if $t = 0$, at ray origin.

For $t \ge 0$, the point where the ray intersects the plane is $$\vec{p}(t) = \vec{p}_0 + t \hat{n}$$

The next step is to determine whether the point is within the triangle or not. There are several different methods to do this, depending on whether you need the 2D coordinates for the intersection point on the plane (for shading purposes or similar), or if you only are interested whether the point is inside the triangle or not.

The typical method is to find the planar cartesian coordinates $(u,v)$ for the intersection point. Typically, we choose $(0,0)$ to be the vertex at the clockwise end of the longest edge, $(w,0)$ at the counterclockwise end of the longest edge, and $(g,h)$ to be the third vertex for a triangle, with $h \gt 0$. Also, $0 \le w \le g$. Below, I shall assume $\vec{v}_1$ is $(0,0)$, $\vec{v}_2$ is $(w,0)$, and $\vec{v}_3$ is $(g,h)$ -- but note that this assumes the edge between the first two is the longest one; relabel/renumber the vertices if necessary.

The planar unit vectors are also often precalculated for each plane, and stored separately: $$\hat{e}_u = \frac{\vec{v}_2 - \vec{v}_1}{\left\lVert \vec{v}_2 - \vec{v}_1 \right\rVert}$$ The second unit vector must be perpendicular to the first, which can easily be done assuming the three vertex vectors are not collinear: $$\vec{e}_v = (\vec{v}_3 - \vec{v}_1) - \hat{e}_u \left( \hat{e}_u \cdot (\vec{v}_3 - \vec{v}_1 \right )$$ and normalizing the result: $$\hat{e}_v = \frac{\vec{e}_v}{\left\lVert\vec{e}_v\right\rVert}$$

The planar coordinates for the three vertices are then $(0,0)$ for $\vec{v}_1$, $(w,0)$ for $\vec{v}_2$ where $$w = \left\lVert \vec{v}_2 - \vec{v}_1 \right\rVert = \hat{e}_u \cdot \left( \vec{v}_2 - \vec{v}_1 \right)$$ and $(g,h)$ for $\vec{v}_3$ where $$\begin{cases} g = \hat{e}_u \cdot \left( \vec{v}_3 - \vec{v}_1 \right ) \\ h = \hat{e}_v \cdot \left( \vec{v}_3 - \vec{v}_1 \right ) \end{cases}$$ Note that the slopes $k_1 = h/g$ and $k_2 = h/(w-g)$ may be useful later, and will not change unless the shape or size of the triangle changes.

Similarly, the coordinates $(u,v)$ for the intersection point $\vec{p}$ are $$\begin{cases} u = \hat{e}_u \cdot \left( \vec{p} - \vec{v}_1 \right ) \\ v = \hat{e}_v \cdot \left( \vec{p} - \vec{v}_1 \right ) \end{cases}$$

If the edge 12 is the longest one in the triangle, then $0 \le g \le w$, $0 \le h$, and it is trivial to check if $(u,v)$ is within the triangle: $$\begin{array}{} \text{if } u \lt 0 \text{ then } \vec{p} \text{ is outside } \\ \text{if } u \gt w \text{ then } \vec{p} \text{ is outside } \\ \text{if } v \lt 0 \text{ then } \vec{p} \text{ is outside } \\ \text{if } v \gt h \text{ then } \vec{p} \text{ is outside } \\ \text{if } u \le g \text{ and } v \gt u \frac{h}{g} \text{ then } \vec{p} \text{ is outside } \\ \text{if } u \ge g \text{ and } v \gt (w - u)\frac{h}{w - g} \text{ then } \vec{p} \text{ is outside } \\ \text{else } \vec{p} \text{ is inside } \end{array}$$

For other planar polygons, and for generic texture support for triangles, you can instead store the 2D $(u,v)$ coordinates for each vertex, and rely on a 2D point-in-polygon test instead.

Although the above are pretty annoying to calculate by hand, if you write the functions needed to compute these in the general case, you'll find that the resulting code needs relatively few basic arithmetic operations per ray per triangle, especially if you precalculate the plane normal, plane distance from origin, the $u$ and $v$ planar unit normals, and the two slopes.

If you program in C, I recommend you start with types

typedef struct {
    double x;
    double y;
    double z;
} Vec3D;

static Vec3D vec3d(const double x, const double y, const double z)
{
    Vec3D result = { x, y, z };
    return result;
}

typedef struct {
    Vec3D  vertex[3];
    Vec3D  unit_u;
    Vec3D  unit_v;
    double u1, u2, v2; /* w, g, h */
} Triangle;

Triangle triangle(Vec3D v1, Vec3D v2, Vec3D v3)
{
    const double dd12 = (v2.x - v1.x)*(v2.x - v1.x)
                      + (v2.y - v1.y)*(v2.y - v1.y)
                      + (v2.z - v1.z)*(v2.z - v1.z);
    const double dd23 = (v3.x - v2.x)*(v3.x - v2.x)
                      + (v3.y - v2.y)*(v3.y - v2.y)
                      + (v3.z - v2.z)*(v3.z - v2.z);
    const double dd31 = (v1.x - v3.x)*(v1.x - v3.x)
                      + (v1.y - v3.y)*(v1.y - v3.y)
                      + (v1.z - v3.z)*(v1.z - v3.z);
    Triangle result;

    if (dd12 >= dd23 && dd12 >= dd31) {
        result.vertex[0] = v1;
        result.vertex[1] = v2;
        result.vertex[3] = v3;
    } else
    if (dd23 >= dd12 && dd23 >= dd31) {
        result.vertex[0] = v2;
        result.vertex[1] = v3;
        result.vertex[2] = v1;
    } else {
        result.vertex[0] = v3;
        result.vertex[1] = v1;
        result.vertex[2] = v2;
    }

    /* TODO: precalculate rest of Triangle result.
             result.vertex[0] to result.vertex[1]
             is the longest edge in the triangle;
             use result.vertex[0..2] instead of v1/v2/v3.
    */

    return result;
}

Answer 2

The ray can be parameterized by the equation $\mathbf{y}(t) = (a_1,b_1,c_1) + t(a_2,b_2,c_2)$ for all $t\geq 0$. The boundaries of the triangle may be parameterized by

\begin{align} \mathbf{s}_1(\tau) &= (x_1,y_1,z_1) + \tau(x_2-x_1,y_2-y_1,z_2-z_1)\\ \mathbf{s}_2(\tau) &= (x_2,y_2,z_2) + \tau(x_3 - x_2,y_3-y_2,z_3-z_2)\\ \mathbf{s}_3(\tau) &= (x_3,y_3,z_3) + \tau(x_1 - x_3,y_1 - y_3,z_1-z_3) \end{align} for all $\tau\in[0,1]$. The triangle lies in a plane with normal vector

$$ \mathbf{n} = (x_2-x_1,y_2-y_1,z_2-z_1)\times(x_2 - x_3,y_2 - y_3,z_2-z_3) $$

If $(a_2,b_2,c_2) \cdot \mathbf{n} \neq 0$, then the line intersects the plane for some $t\in\mathbb{R}$. The intersection point with the ray, if it exists, is given by solving $\mathbf{n}\cdot(\mathbf{y}(t) - (x_1,y_1,z_1)) = 0$ for $t\geq 0$. Lastly, you have to check if the intersection point lies on the triangle. To do this let $\mathbf{y}$ be the intersection point and suppose it does not lie on a boundary of the triangle (if it does then you are done). Define the line $\mathbf{x}(t) = \mathbf{y} + t(\mathbf{s}_1(0.5) - \mathbf{y})$ for all $t\geq 0$. If $\mathbf{x}(t) = \mathbf{s}_2(\tau)$ or $\mathbf{x}(t) = \mathbf{s}_3(\tau)$ have a solution for $t \geq 0$ and $\tau\in[0,1]$, then the point is exterior to the triangle, otherwise it is interior. The distance from the intersection point to the origin of the ray is simply $\|\mathbf{y} - (a_1,b_1,c_1)\|$.