$
\newcommand{\vp}{\mathbf{p}}
\newcommand{\vo}{\mathbf{o}}
\newcommand{\vd}{\mathbf{d}}
\newcommand{\ve}{\mathbf{e}}
\newcommand{\ved}{\bar{\mathbf{e}}}
\newcommand{\vv}{\mathbf{v}}
\newcommand{\vn}{\mathbf{n}}
$
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 $\vv_0,\vv_1,\vv_2$ (in no particular order),
while the ray is defined by origin $\vo$ and direction vector $\vd$ (not necessarily normalized).
In what follows, vertex $\vv_0$ is called the
triangle origin, while edges from the origin are called
$\ve_1=\vv_1-\vv_0$ and $\ve_2=\vv_2-\vv_0$. The (normalized) normal of the triangle is $\vn$, which
can be obtained as:
$$
\vn ~~=~~\frac{\ve_1\times\ve_2}{\|\ve_1\times\ve_2\|}
$$
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 $\vp$ to the point of intersection between the line containing the ray
and the plane $\Pi$ containing the triangle (assuming
they are not parallel), then obviously there exist three real values $t,c_1,c_2$ such that:
$$
\begin{equation}\label{eq:p}
\vp~~=~~ \vo+t\vd ~~=~~ \vv_0\,+\,c_1\ve_1 \,+\, c_2\ve_2
\end{equation}
$$
note that $t$ is the distance to the point of intersection from ray origin, and $(c_1,c_2)$ are
the coordinates of $\vp$ with respect to the coordinate frame defined by $\ve_1,\ve_2$ in plane $\Pi$
(the pair of vectors $\{\ve_1,\ve_2\}$ form a basis for free vectors in plane $\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 $\vp$ is easily computed from $(\ref{eq:p})$.
Computation of $t$ can be done, as usual, by observing that vector $\vp-\vv_0$ is in $\Pi$ thus perpendicular
to $\vn$, that is $(\vp-\vv_0)\cdot\vn=0$. From this equality we can solve for
$t$, we get:
$$
\begin{equation}
t ~~=~~ \frac{(\vo-\vv_0)\cdot\vn}{-\vd\cdot\vn} ~~=~~ \frac{\vo\cdot\vn-k_3}{-\vd\cdot\vn}
~~~~~~~~~~~~\mbox{where:} ~~~k_3\,=\,\vv_0\cdot\vn
\end{equation}
$$
Once $t$ is computed, we obtain $\vp$ from ($\ref{eq:p}$).
Then we can directly compute both $c_1$ and $c_2$ as follows:
$$
\begin{equation}
\begin{array}{l}
c_1~~=~~ (\vp-\vv_0)\cdot\ved_1 ~~=~~ \vp\cdot\ved_1\,-\,k_1 \\
c_2~~=~~ (\vp-\vv_0)\cdot\ved_2 ~~=~~ \vp\cdot\ved_2\,-\,k_2
\end{array}
~~~~~~~~~~~\mbox{where:}~~~~~~
\begin{array}{l}
k_1~~=~~ \vv_0\cdot\ved_1\\
k_2~~=~~ \vv_0\cdot\ved_2
\end{array}
\end{equation}
$$
vectors $\ved_1$ and $\ved_2$ are called
dual edges and are introduced in the next section.
4. Reciprocal basis
As stated, vectors $\{\ve_1,\ve_2\}$ form a
basis for the space of all free vectors in plane $\Pi$. Thus there exists
another basis for that space which is called the
reciprocal or
dual basis, and which we call
$\{\ved_1,\ved_2\}$. The following equalities hold:
$$
\begin{equation}\label{eq:dualbasisprop}
\begin{array}{ccc}
\ve_1\cdot\ved_1 ~=~ 1 & ~~~~~~ \ve_1\cdot\ved_2 ~=~ 0 \\
\ve_2\cdot\ved_1 ~=~ 0 & ~~~~~~ \ve_2\cdot\ved_2 ~=~ 1 \\
\end{array}
\end{equation}
$$
These equalities can be used to obtain explicit expressions for $\{\ved_1,\ved_2\}$ in terms
of $\{\ve_1,\ve_2\}$. We know there exist four unknown real values $a,b,c$ and $d$ such that
$\ved_1=a\ve_1+b\ve_2$ and $\ved_2=c\ve_1+d\ve_2$. By substituting these expression in ($\ref{eq:dualbasisprop}$)
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{equation}\label{eq:dualbasis}
\begin{array}{l}
\ved_1 ~~=~~ s\,\left[ (\ve_2\cdot\ve_2)\ve_1 - (\ve_1\cdot\ve_2)\ve_2 \right] \\
\ved_2 ~~=~~ s\,\left[ (\ve_1\cdot\ve_1)\ve_2 - (\ve_1\cdot\ve_2)\ve_1 \right]
\end{array}
~~~~~~\mbox{where:} ~~~
s ~~=~~ \frac{1}{(\ve_1\cdot\ve_1)(\ve_2\cdot\ve_2) \,-\, (\ve_1\cdot\ve_2)^2}
\end{equation}
$$
value $s$ must be finite when $\ve_1$ and $\ve_2$ are not parallel neither null. By substituting ($\ref{eq:dualbasis}$) into
($\ref{eq:dualbasisprop}$) it can be verified that ($\ref{eq:dualbasisprop}$) 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}
\ved_1 ~=~ s\,\left[ (\ve_2\cdot\ve_2)\ve_1 - (\ve_1\cdot\ve_2)\ve_2 \right] & ~~~~~~~~~~~~~~~~~~~~~~~~~~& k_1~=~ \vv_0\cdot\ved_1 \\
\ved_2 ~=~ s\,\left[ (\ve_1\cdot\ve_1)\ve_2 - (\ve_1\cdot\ve_2)\ve_1 \right] & ~~~~~~~~~~~~~~~~~~~~~~~~~~& k_2~=~ \vv_0\cdot\ved_2 \\
\vn ~~=~~(\ve_1\times\ve_2)/\|\ve_1\times\ve_2\| & ~~~~~~~~~~~~~~~~~~~~~~~~~ & k_3 ~=~ \vv_0\cdot\vn
\end{array}
$$
note that basis $\{\ved_1,\ved_2,\vn\}$ is the dual basis of $\{\ve_1,\ve_2,\vn\}$ (because, as $\vn$ is normalized and perpendicular to $\Pi$, it holds
$\vn\cdot\vn=1$, and $\vn\cdot\ved_i=\vn\cdot\ve_i=0$). Tuple $(k_1,k_2,k_3)$ holds the coordinates of $\vv_0$ (as a free vector) w.r.t. basis $\{\ve_1,\ve_2,\vn\}$.
For a nearly degenerated triangle ($\ve_1$ nearly parallel to $\ve_2$, or either $\ve_1$ or $\ve_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:
- Let $a\,=\,-\vd\cdot\vn$
- If $\|a\|<\epsilon$ then return (no intersection: line nearly parallel to plane)
- Let $t\,=~(\vo\cdot\vn-k_3)/a$
- If $t\leq 0$ then return (no intersection: line intersects plane in the negative half-line)
- Let $\vp\,=~\vo+t\vd$
- Let $c_1\,=~\vp\cdot\ved_1\,-\,k_1$
- If $c_1<0$ or $1< c_1$ then return (no intersection: point not in triangle)
- Let $c_2\,=~\vp\cdot\ved_2\,-\,k_2$
- If $c_2<0$ or $1< c_2$ then return (no intersection: point not in triangle)
- If $1< c_1+c_2$ then return (no intersection: point not in triangle)
- Return $(\vp,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 & = & -\vd'\cdot\vn' \\
t & = & (\vo'\cdot\vn')/a \\
\vp' & = & \vo'+t\vd' \\
c_1 & = & \vp'\cdot\ved'_1 \\
c_2 & = & \vp'\cdot\ved'_2
\end{array}
~~~~~~~~~\mbox{where:}~~~~
\left\{\begin{array}{rcl}
\vo' & = & [\vo:1] \\
\vd' & = & [\vd:0] \\
\vn' & = & [\vn:-k_3] \\
\ved'_1 & = & [\ved_1:-k_1] \\
\ved'_2 & = & [\ved_2:-k_2]
\end{array}\right.
$$
8. Recovering original representation.
Vectors $\vv_0$,$\ve_1$ and $\ve_2$ can be recovered from the dual representation.
Edges $\{\ve_1,\ve_2\}$ can be easily obtained from $\{\ved_1,\ved_2\}$
by using the inverse version of (\ref{eq:dualbasis}), because equalities (\ref{eq:dualbasisprop}) are
symmetric, that is:
$$
\begin{equation}\label{eq:dualbasisdual}
\begin{array}{l}
\ve_1 ~~=~~ s'\,\left[ (\ved_2\cdot\ved_2)\ved_1 - (\ved_1\cdot\ved_2)\ved_2 \right] \\
\ve_2 ~~=~~ s'\,\left[ (\ved_1\cdot\ved_1)\ved_2 - (\ved_1\cdot\ved_2)\ved_1 \right]
\end{array}
~~~~~~~~~~\mbox{where:}~~
s' ~=~ \frac{1}{(\ved_1\cdot\ved_1)(\ved_2\cdot\ved_2) \,-\, (\ved_1\cdot\ved_2)^2}
\end{equation}
$$
Vertex $\vv_0$ can be computed from its coordinates $(k_1,k_2,k_3)$ (w.r.t basis $\{\ve_1,\ve_2,\vn\}$), namely:
$$
\vv_0 ~=~ k_1\ve_1 \,+\, k_2\ve_2 \,+\, k_3\vn
$$
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)