Ray Triangle Collision

From GPWiki
Jump to: navigation, search

This is a simple method using a transformation matrix.


Creating the transformation matrix

Now, the triangle is made of vertices like so:


Ray Triangle Collision triangle.png 
A = \begin{bmatrix}
A_x \\
A_y \\
A_z
\end{bmatrix}
,
B = \begin{bmatrix}
B_x \\
B_y \\
B_z
\end{bmatrix}
,
C = \begin{bmatrix}
C_x \\
C_y \\
C_z
\end{bmatrix}

For our purpose, we want a world -> local transformation matrix that uses AB as the X axis, AC as the Y axis, and the triangle's normal as the Z axis.

Ray Triangle Collision axes.png 
\begin{matrix}
a & = & B - A \\
b & = & C - A \\
c & = & a \times b
\end{matrix}

Note that the cross product is calculated as follows:


c = 
\begin{bmatrix}
a_y b_z - a_z b_y \\
a_z b_x - a_x b_z \\
a_x b_y - a_y b_x
\end{bmatrix}

We'll need that later

These 4 vectors are placed into a matrix:


T = 
\begin{bmatrix}
a_x & b_x & c_x & A_x \\
a_y & b_y & c_y & A_y \\
a_z & b_z & c_z & A_z \\
0 & 0 & 0 & 1
\end{bmatrix}

This is the matrix to convert from {a,b,c} coordinates to {x,y,z} coordinates. We want to do the reverse, so we invert the matrix.

Knowing a few things about this matrix, we can invert the upper-left 3x3 matrix, and use it to transform the translation.


M = 
\begin{bmatrix}
a_x & b_x & c_x \\
a_y & b_y & c_y \\
a_z & b_z & c_z
\end{bmatrix}


The method I used is that shown in Inversion of 3 x 3 matrices


M^{-1} = {1 \over det(M)} adj(M)

where det(M) is the determinant of M, and adj(M) is the adjugate of M


First, we create the adjugate matrix:


adj(M) = 
\begin{bmatrix}
b_y c_z - b_z c_y & b_z c_x - b_x c_z & b_x c_y - b_y c_x \\
c_y a_z - c_z a_y & c_z a_x - c_x a_z & c_x a_y - c_y a_x \\
a_y b_z - a_z b_y & a_z b_x - a_x b_z & a_x b_y - a_y b_x
\end{bmatrix}


Notice that row 3 of the matrix was already calculated when we calculated c. So, we can remove those redundant operations:


adj(M) =
\begin{bmatrix}
b_y c_z - b_z c_y & b_z c_x - b_x c_z & b_x c_y - b_y c_x \\
c_y a_z - c_z a_y & c_z a_x - c_x a_z & c_x a_y - c_y a_x \\
c_x & c_y & c_z
\end{bmatrix}


Now, we calculate the determinant.


\begin{matrix}
det(M) & = & a_y b_z c_x + a_z b_x c_y + a_x b_y c_z - a_z b_y c_x - a_x b_z c_y - a_y b_x c_z \\
       & = & c_x (a_y b_z - a_z b_y) + c_y (a_z b_x - a_x b_z) + c_z (a_x b_y - a_y b_x) \\
       & = & c_x^2 + c_y^2 + c_z^2
\end{matrix}


We now divide the matrix by the determinant:


M^{-1} = {1 \over {c_x^2 + c_y^2 + c_z^2}} \begin{bmatrix}
b_y c_z - b_z c_y & b_z c_x - b_x c_z & b_x c_y - b_y c_x \\
c_y a_z - c_z a_y & c_z a_x - c_x a_z & c_x a_y - a_x c_y \\
c_x & c_y & c_z
\end{bmatrix}


Now, we transform the translation:


\begin{matrix}
O & = & -M^{-1} A \\
  & = & -
\begin{bmatrix}
(M^{-1})_{x,x} & (M^{-1})_{x,y} & (M^{-1})_{x,z} \\
(M^{-1})_{y,x} & (M^{-1})_{y,y} & (M^{-1})_{y,z} \\
(M^{-1})_{z,x} & (M^{-1})_{z,y} & (M^{-1})_{z,z}
\end{bmatrix}
\begin{bmatrix}
A_x \\ A_y \\ A_z
\end{bmatrix}
\\
  & = & -
\begin{bmatrix}
(M^{-1})_{x,x} A_x + (M^{-1})_{x,y} A_y + (M^{-1})_{x,z} A_z \\
(M^{-1})_{y,x} A_x + (M^{-1})_{y,y} A_y + (M^{-1})_{y,z} A_z \\
(M^{-1})_{z,x} A_x + (M^{-1})_{z,y} A_y + (M^{-1})_{z,z} A_z
\end{bmatrix}
\end{matrix}


We can now create the {x,y,z} -> {a,b,c} transformation matrix R:


R = 
\begin{bmatrix}
(M^{-1})_{x,x} & (M^{-1})_{x,y} & (M^{-1})_{x,z} & O_x \\
(M^{-1})_{y,x} & (M^{-1})_{y,y} & (M^{-1})_{y,z} & O_y \\
(M^{-1})_{z,x} & (M^{-1})_{z,y} & (M^{-1})_{z,z} & O_z \\
0 & 0 & 0 & 1
\end{bmatrix}


Phew.

Here's some C code to do that:

double **create_collision_matrix (double *A, double *B, double *C){

 double a[3], b[3], c[3], in_det;
 static double R[4][4], Rp[4];
 /* Fill in Rp pointers */
 Rp[0] = R[0];
 Rp[1] = R[1];
 Rp[2] = R[2];
 Rp[3] = R[3];
 /* a = B - A */
 a[0] = B[0] - A[0]; 
 a[1] = B[1] - A[1]; 
 a[2] = B[2] - A[2];
 /* b = C - B */
 b[0] = C[0] - A[0];
 b[1] = C[1] - A[1];
 b[2] = C[2] - A[2];
 /* c = a × b */
 c[0] = a[1] * b[2] - a[2] * b[1];
 c[1] = a[2] * b[0] - a[0] * b[2];
 c[2] = a[0] * b[1] - a[1] * b[0];
 /* M^(-1) = (1/det(M)) * adj(M) */
 in_det = 1 / (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
 R[0][0] = (b[1] * c[2] - b[2] * c[1]) * in_det;
 R[0][1] = (b[2] * c[0] - b[0] * c[2]) * in_det;
 R[0][2] = (b[0] * c[1] - b[1] * c[0]) * in_det;
 R[1][0] = (c[1] * a[2] - c[2] * a[1]) * in_det;
 R[1][1] = (c[2] * a[0] - c[0] * a[2]) * in_det;
 R[1][2] = (c[0] * a[1] - c[1] * a[0]) * in_det;
 R[2][0] = (c[0]) * in_det;
 R[2][1] = (c[1]) * in_det;
 R[2][2] = (c[2]) * in_det;
 
 /* O = M^(-1) * A */
 R[0][3] = -(R[0][0] * A[0] + R[0][1] * A[1] + R[0][2] * A[2]);
 R[1][3] = -(R[1][0] * A[0] + R[1][1] * A[1] + R[1][2] * A[2]);
 R[2][3] = -(R[2][0] * A[0] + R[2][1] * A[1] + R[2][2] * A[2]);
 /* fill in last row of 4x4 matrix */
 R[3][0] = R[3][1] = R[3][2] = 0;
 R[3][3] = 1;
 
 return Rp;

}

What if the triangle needs to be rotated?

If we store the collision matrices in the model, then we need a quick way of transforming them.

We know that the transformation matrix T is used to transform from {a,b,c} coordinates to {x,y,z} local coordinates. If the local->world transformation matrix is W, then transforming from {a,b,c} coordinates to {x,y,z} world coordinates is like so:


\begin{matrix}
P & = & W (T p) \\
  & = & (W T) p
\end{matrix}


We know that

(AB)-1 = B-1A-1

so we can calculate W-1, and calculate the new matrix by T-1W-1.

So, we calculate the inverse of W for the entire object, and apply the matrix R (which is the inverse of matrix T) to it.

If we can assume that W has no scaling or skewing, then the linear transformation part of the W-1 is the transpose of W, and it's very simple:


\begin{bmatrix}
x_x & x_y & x_z & x \\
y_x & y_y & y_z & y \\
z_x & z_y & z_z & z \\
0 & 0 & 0 & 1
\end{bmatrix}
^{-1} =
\begin{bmatrix}
x_x & y_x & z_x & -(x_x x + y_x y + z_x z) \\
x_y & y_y & z_y & -(x_y x + y_y y + z_y z) \\
x_z & y_z & z_z & -(x_z x + y_z y + z_z z) \\
0 & 0 & 0 & 1
\end{bmatrix}


If we cannot assume this, then we use the standard matrix inversion on the 3x3 linear transformation matrix, and apply it to the translation part of the matrix.

This can all be done when the object transformation matrix is constructed.

The matrix multiplication can then be done on a triangle-by-triangle basis.


Using the transformation matrix

We use a simple step-by-step transformation on both ends of the ray.

We'll use the ray j->k.


Ray Triangle Collision line.png


First, we do the z axis of the {x,y,z} -> {a,b,c} transformation.


\begin{matrix}
j'_z & = & R_{z,x} j_x + R_{z,y} j_y + R_{z,z} j_z + R_{z,w} \\
k'_z & = & R_{z,x} k_x + R_{z,y} k_y + R_{z,z} k_z + R_{z,w}
\end{matrix}

If j'z and k'z are the same sign, then the ray does not intersect with the triangle, and we can go home.

If both are zero, which is statistically very unlikely, we can probably go home.

If they are of opposite sign, or one of them is zero, then the ray might intersect with the triangle, and we must continue.


We then do the x axis.

We calculate the x value of the point of intersection:


i'_x = j'_x + j'_z {k'_x - j'_x \over j'_z - k'_z}


If i'x is over 1 or under 0, then the ray does not intersect with the triangle. Otherwise, we do the same with the y axis.

If i'y is over 1 or under 0, or i'x + i'y is over 1, then the ray does not intersect with the triangle, and we can go home.

If we have reached this point, then we have found an intersection between the ray and the triangle.


Ray Triangle Collision collision.png


Here's some code to do this:

double *collide (double **R, double *j, double *k){

 double J[3], K[3];
 static double i[2];
 J[2] = R[2][0] * j[0] + R[2][1] * j[1] + R[2][2] * j[2] + R[2][3];
 K[2] = R[2][0] * k[0] + R[2][1] * k[1] + R[2][2] * k[2] + R[2][3];
 if (J[2] * K[2] >= 0){
   return NULL;
 }
 J[0] = R[0][0] * j[0] + R[0][1] * j[1] + R[0][2] * j[2] + R[0][3];
 K[0] = R[0][0] * k[0] + R[0][1] * k[1] + R[0][2] * k[2] + R[0][3];
 i[0] = J[0] + J[2] * ((K[0] - J[0]) / (J[2] - K[2]));
 if (i[0] < 0 || i[0] > 1){
   return NULL;
 }
 
 J[1] = R[1][0] * j[0] + R[1][1] * j[1] + R[1][2] * j[2] + R[1][3];
 K[1] = R[1][0] * k[0] + R[1][1] * k[1] + R[1][2] * k[2] + R[1][3];
 i[1] = J[1] + J[2] * ((K[1] - J[1]) / (J[2] - K[2]));
 if (i[1] < 0 || i[1] > 1 || i[0] + i[1] > 1){
   return NULL;
 }
 return i;

}


If we wanted to, we could use i'x and i'y to do further detection - e.g. if we have a texture that has holes in it, we could determine whether the ray passes through one of those holes.


See also

Fast, minimum storage ray-triangle intersection.