Carlos Ureña / Varios

Efficient ray-triangle intersection test

$$

1. Abstract

I introduce here a time-optimized ray-triangle intersection test, where time optimization is obtained at the cost of using a dual representation of the triangle. The algorithm is a variant of known algorithms, which precomputes as much values as possible. I show how this dual representation can be obtained from the usual representation, and vice-versa, the usual representation can be recovered from the dual.

2. Triangle and ray

The triangle is defined by three arbitrary 3D vertex points, we call them ${\mathbf{v}}_0,{\mathbf{v}}_1,{\mathbf{v}}_2$ (in no particular order), while the ray is defined by origin $\mathbf{o}$ and direction vector $\mathbf{d}$ (not necessarily normalized). In what follows, vertex ${\mathbf{v}}_0$ is called the triangle origin, while edges from the origin are called ${\mathbf{e}}_1={\mathbf{v}}_1-{\mathbf{v}}_0$ and ${\mathbf{e}}_2={\mathbf{v}}_2-{\mathbf{v}}_0$. The (normalized) normal of the triangle is $\mathbf{n}$, which can be obtained as: $$\mathbf{n}=\frac{{\mathbf{e}}_1\times {\mathbf{e}}_2}{\Vert {\mathbf{e}}_1\times {\mathbf{e}}_2\Vert }$$ all these vectors are represented by 3-tuples, holding their coordinates w.r.t world reference system , whose origin is $\mathbf{0}$.

3. Intersection point properties

Let us call $\mathbf{p}$ to the point of intersection between the line containing the ray and the plane $\mathrm{\Pi }$ containing the triangle (assuming they are not parallel), then obviously there exist three real values $t,c_1,c_2$ such that: $$\begin{array}{}\text{(1)}& \mathbf{p}=\mathbf{o}+t\mathbf{d}={\mathbf{v}}_0+c_1{\mathbf{e}}_1+c_2{\mathbf{e}}_2\end{array}$$ note that $t$ is the distance to the point of intersection from ray origin, and $(c_1,c_2)$ are the coordinates of $\mathbf{p}$ with respect to the coordinate frame defined by ${\mathbf{e}}_1,{\mathbf{e}}_2$ in plane $\mathrm{\Pi }$ (the pair of vectors $\{{\mathbf{e}}_1,{\mathbf{e}}_2\}$ form a basis for free vectors in plane $\mathrm{\Pi }$).

The problem of intersection computation is in fact the problem of computing $t,c_1,c_2$, because there is an intersection if and only if each of $c_1,c_2$ and $c_1+c_2$ are in the interval $[0,1]$, and $t>0$. In that case $\mathbf{p}$ is easily computed from $(\text{1})$.

Computation of $t$ can be done, as usual, by observing that vector $\mathbf{p}-{\mathbf{v}}_0$ is in $\mathrm{\Pi }$ thus perpendicular to $\mathbf{n}$, that is $(\mathbf{p}-{\mathbf{v}}_0)\cdot \mathbf{n}=0$. From this equality we can solve for $t$, we get: $$\begin{array}{}\text{(2)}& t=\frac{(\mathbf{o}-{\mathbf{v}}_0)\cdot \mathbf{n}}{-\mathbf{d}\cdot \mathbf{n}}=\frac{\mathbf{o}\cdot \mathbf{n}-k_3}{-\mathbf{d}\cdot \mathbf{n}}\text{where:}{k}_3={\mathbf{v}}_0\cdot \mathbf{n}\end{array}$$

Once $t$ is computed, we obtain $\mathbf{p}$ from ($\text{1}$). Then we can directly compute both $c_1$ and $c_2$ as follows: $$\begin{array}{}\text{(3)}& \begin{array}{l}{c}_1=(\mathbf{p}-{\mathbf{v}}_0)\cdot {\overline{\mathbf{e}}}_1=\mathbf{p}\cdot {\overline{\mathbf{e}}}_1-k_1\\ c_2=(\mathbf{p}-{\mathbf{v}}_0)\cdot {\overline{\mathbf{e}}}_2=\mathbf{p}\cdot {\overline{\mathbf{e}}}_2-k_2\end{array}\text{where:}\begin{array}{l}{k}_1={\mathbf{v}}_0\cdot {\overline{\mathbf{e}}}_1\\ k_2={\mathbf{v}}_0\cdot {\overline{\mathbf{e}}}_2\end{array}\end{array}$$ vectors ${\overline{\mathbf{e}}}_1$ and ${\overline{\mathbf{e}}}_2$ are called dual edges and are introduced in the next section.

4. Reciprocal basis

As stated, vectors $\{{\mathbf{e}}_1,{\mathbf{e}}_2\}$ form a basis for the space of all free vectors in plane $\mathrm{\Pi }$. Thus there exists another basis for that space which is called the reciprocal or dual basis, and which we call $\{{\overline{\mathbf{e}}}_1,{\overline{\mathbf{e}}}_2\}$. The following equalities hold: $$\begin{array}{}\text{(4)}& \begin{array}{cc}{\mathbf{e}}_1\cdot {\overline{\mathbf{e}}}_1=1& {\mathbf{e}}_1\cdot {\overline{\mathbf{e}}}_2=0\\ {\mathbf{e}}_2\cdot {\overline{\mathbf{e}}}_1=0& {\mathbf{e}}_2\cdot {\overline{\mathbf{e}}}_2=1\end{array}\end{array}$$

These equalities can be used to obtain explicit expressions for $\{{\overline{\mathbf{e}}}_1,{\overline{\mathbf{e}}}_2\}$ in terms of $\{{\mathbf{e}}_1,{\mathbf{e}}_2\}$. We know there exist four unknown real values $a,b,c$ and $d$ such that ${\overline{\mathbf{e}}}_1=a{\mathbf{e}}_1+b{\mathbf{e}}_2$ and ${\overline{\mathbf{e}}}_2=c{\mathbf{e}}_1+d{\mathbf{e}}_2$. By substituting these expression in ($\text{4}$) we arrive to a system of four equation and four unknowns, then by solving it we can write the dual basis in terms of the original one, as follows: $$\begin{array}{}\text{(5)}& \begin{array}{l}{\overline{\mathbf{e}}}_1=s[({\mathbf{e}}_2\cdot {\mathbf{e}}_2){\mathbf{e}}_1-({\mathbf{e}}_1\cdot {\mathbf{e}}_2){\mathbf{e}}_2]\\ {\overline{\mathbf{e}}}_2=s[({\mathbf{e}}_1\cdot {\mathbf{e}}_1){\mathbf{e}}_2-({\mathbf{e}}_1\cdot {\mathbf{e}}_2){\mathbf{e}}_1]\end{array}\text{where:}s=\frac{1}{({\mathbf{e}}_1\cdot {\mathbf{e}}_1)({\mathbf{e}}_2\cdot {\mathbf{e}}_2)-({\mathbf{e}}_1\cdot {\mathbf{e}}_2{)}^2}\end{array}$$ value $s$ must be finite when ${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ are not parallel neither null. By substituting ($\text{5}$) into ($\text{4}$) it can be verified that ($\text{4}$) holds.

5. Dual representation of triangles

The intersection algorithm must precompute and store, for each triangle, these values and vectors (which we call the dual representation of the triangle): $$\begin{array}{lcl}{\overline{\mathbf{e}}}_1=s[({\mathbf{e}}_2\cdot {\mathbf{e}}_2){\mathbf{e}}_1-({\mathbf{e}}_1\cdot {\mathbf{e}}_2){\mathbf{e}}_2]& & k_1={\mathbf{v}}_0\cdot {\overline{\mathbf{e}}}_1\\ {\overline{\mathbf{e}}}_2=s[({\mathbf{e}}_1\cdot {\mathbf{e}}_1){\mathbf{e}}_2-({\mathbf{e}}_1\cdot {\mathbf{e}}_2){\mathbf{e}}_1]& & k_2={\mathbf{v}}_0\cdot {\overline{\mathbf{e}}}_2\\ \mathbf{n}=({\mathbf{e}}_1\times {\mathbf{e}}_2)/\Vert {\mathbf{e}}_1\times {\mathbf{e}}_2\Vert & & k_3={\mathbf{v}}_0\cdot \mathbf{n}\end{array}$$ note that basis $\{{\overline{\mathbf{e}}}_1,{\overline{\mathbf{e}}}_2,\mathbf{n}\}$ is the dual basis of $\{{\mathbf{e}}_1,{\mathbf{e}}_2,\mathbf{n}\}$ (because, as $\mathbf{n}$ is normalized and perpendicular to $\mathrm{\Pi }$, it holds $\mathbf{n}\cdot \mathbf{n}=1$, and $\mathbf{n}\cdot {\overline{\mathbf{e}}}_i=\mathbf{n}\cdot {\mathbf{e}}_i=0$). Tuple $(k_1,k_2,k_3)$ holds the coordinates of ${\mathbf{v}}_0$ (as a free vector) w.r.t. basis $\{{\mathbf{e}}_1,{\mathbf{e}}_2,\mathbf{n}\}$.

For a nearly degenerated triangle (${\mathbf{e}}_1$ nearly parallel to ${\mathbf{e}}_2$, or either ${\mathbf{e}}_1$ or ${\mathbf{e}}_2$ nearly null), the absolute value of $1/s$ is near $0$. Thus, computing $1/s$ can be used for discarding nearly degenerated triangles.

6. The algorithm

The intersection algorithm can be written by using the dual representation, in such a way that it exits as soon as possible when there is no intersection:

1. Let $a=-\mathbf{d}\cdot \mathbf{n}$
2. If $\Vert a\Vert <ϵ$ then return (no intersection: line nearly parallel to plane)
3. Let $t=(\mathbf{o}\cdot \mathbf{n}-k_3)/a$
4. If $t\le 0$ then return (no intersection: line intersects plane in the negative half-line)
5. Let $\mathbf{p}=\mathbf{o}+t\mathbf{d}$
6. Let $c_1=\mathbf{p}\cdot {\overline{\mathbf{e}}}_1-k_1$
7. If $c_1<0$ or $1<c_1$ then return (no intersection: point not in triangle)
8. Let $c_2=\mathbf{p}\cdot {\overline{\mathbf{e}}}_2-k_2$
9. If $c_2<0$ or $1<c_2$ then return (no intersection: point not in triangle)
10. If $1<c_1+c_2$ then return (no intersection: point not in triangle)
11. Return $(\mathbf{p},t,c_1,c_2)$ (there is an intersection)

This algorithm requires at least two and at most four dot products. On the other hand, it requires storage for 12 real values per triangle.

7. Using homogeneous coordinates

In some case it can be efficient to use 4-tuple dot product routines implemented in hardware. In the algorithm above, this can be done for all the dots products involved, by extending the dual basis vectors with $k_1,k_2$ and $k_3$ as their homogeneous coordinates, then $t,c_1$ and $c_2$ can be computed as: $$\begin{array}{rcl}a& =& -{\mathbf{d}}^{\prime }\cdot {\mathbf{n}}^{\prime }\\ t& =& ({\mathbf{o}}^{\prime }\cdot {\mathbf{n}}^{\prime })/a\\ {\mathbf{p}}^{\prime }& =& {\mathbf{o}}^{\prime }+t{\mathbf{d}}^{\prime }\\ c_1& =& {\mathbf{p}}^{\prime }\cdot {\overline{\mathbf{e}}}_1^{\prime }\\ c_2& =& {\mathbf{p}}^{\prime }\cdot {\overline{\mathbf{e}}}_2^{\prime }\end{array}\text{where:}\{\begin{array}{rcl}{\mathbf{o}}^{\prime }& =& [\mathbf{o}:1]\\ {\mathbf{d}}^{\prime }& =& [\mathbf{d}:0]\\ {\mathbf{n}}^{\prime }& =& [\mathbf{n}:-k_3]\\ {\overline{\mathbf{e}}}_1^{\prime }& =& [{\overline{\mathbf{e}}}_1:-k_1]\\ {\overline{\mathbf{e}}}_2^{\prime }& =& [{\overline{\mathbf{e}}}_2:-k_2]\end{array}$$

8. Recovering original representation.

Vectors ${\mathbf{v}}_0$,${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ can be recovered from the dual representation. Edges $\{{\mathbf{e}}_1,{\mathbf{e}}_2\}$ can be easily obtained from $\{{\overline{\mathbf{e}}}_1,{\overline{\mathbf{e}}}_2\}$ by using the inverse version of ($\text{5}$), because equalities ($\text{4}$) are symmetric, that is: $$\begin{array}{}\text{(6)}& \begin{array}{l}{\mathbf{e}}_1=s^{\prime }[({\overline{\mathbf{e}}}_2\cdot {\overline{\mathbf{e}}}_2){\overline{\mathbf{e}}}_1-({\overline{\mathbf{e}}}_1\cdot {\overline{\mathbf{e}}}_2){\overline{\mathbf{e}}}_2]\\ {\mathbf{e}}_2=s^{\prime }[({\overline{\mathbf{e}}}_1\cdot {\overline{\mathbf{e}}}_1){\overline{\mathbf{e}}}_2-({\overline{\mathbf{e}}}_1\cdot {\overline{\mathbf{e}}}_2){\overline{\mathbf{e}}}_1]\end{array}\text{where:}{s}^{\prime }=\frac{1}{({\overline{\mathbf{e}}}_1\cdot {\overline{\mathbf{e}}}_1)({\overline{\mathbf{e}}}_2\cdot {\overline{\mathbf{e}}}_2)-({\overline{\mathbf{e}}}_1\cdot {\overline{\mathbf{e}}}_2{)}^2}\end{array}$$ Vertex ${\mathbf{v}}_0$ can be computed from its coordinates $(k_1,k_2,k_3)$ (w.r.t basis $\{{\mathbf{e}}_1,{\mathbf{e}}_2,\mathbf{n}\}$), namely: $${\mathbf{v}}_0=k_1{\mathbf{e}}_1+k_2{\mathbf{e}}_2+k_3\mathbf{n}$$

9. See also:

• Ray-triangle intersection by Dan Sunday (GeomAlgorithms) -- in fact this uses dual vectors, but (implicitly) computes then for every intersection.

(this page uses MathJax)