Ray Triangle Collision

From GPWiki

This is a simple method using a transformation matrix.

Creating the transformation matrix

Now, the triangle is made of vertices like so:

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.

Note that the cross product is calculated as follows:

We'll need that later

These 4 vectors are placed into a matrix:

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.

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

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

First, we create the adjugate matrix:

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

Now, we calculate the determinant.

We now divide the matrix by the determinant:

Now, we transform the translation:

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

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:

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:

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.

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

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:

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.

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.