Ray Triangle Intersection Example

Here is worked numerical example of Ray Triangle intersection with a 3D visualization of the process and result.

Ross Beveridge, September 20, 2018

The actual 3D position of the 3 triangle corners as well as the ray starting point and direction are given here at the head of the file. Initially the example is setup using these values:

Av = vector(SR, 3, (3,0,0)); Bv = vector(SR, 3, (0,3,0)); Cv = vector(SR, 3, (0,0,3)); Lv = vector(SR, 3, (0,0,0)); Dv = vector(SR, 3, (1,1,1));

To further explore consider these alternatives:

Av = vector(SR, 3, (6,0,0)); Bv = vector(SR, 3, (0,6,0)); Cv = vector(SR, 3, (0,0,6)); Lv = vector(SR, 3, (1,1,1)); Dv = vector(SR, 3, (1,1,1));

Av = vector(SR, 3, (6,0,0)); Bv = vector(SR, 3, (0,6,0)); Cv = vector(SR, 3, (0,0,6)); Lv = vector(SR, 3, (1,1,1)); Dv = vector(SR, 3, (1,0,0));

Av = vector(SR, 3, (6,0,0)); Bv = vector(SR, 3, (0,6,0)); Cv = vector(SR, 3, (0,0,6)); Lv = vector(SR, 3, (0,0,0)); Dv = vector(SR, 3, (0,1,1));

In [52]:
latex.matrix_delimiters("[", "]")
var('t', 'beta', 'gamma');
# What follows are specific values for triangle and ray specification
Av = vector(SR, 3, (3,0,0));
Bv = vector(SR, 3, (0,3,0));
Cv = vector(SR, 3, (0,0,3)); 
Lv = vector(SR, 3, (0,0,0));
Dv = vector(SR, 3, (1,1,1));

Dv = Dv / Dv.norm();
A = matrix(SR, 3,1, Av); B = matrix(SR, 3,1, Bv); C = matrix(SR, 3,1, Cv);
L = matrix(SR, 3,1, Lv); D = matrix(SR, 3,1, Dv);

Here is the parameterized ray as well as the parameterized plane containing the Triangle

In [53]:
Rt = L + t * D;
pretty_print(LatexExpr("R(t) = L + t D ="), L, "+", t, "*", D, "=", Rt);
tstr = LatexExpr("P(") + latex(beta) + "," + latex(gamma) + LatexExpr(") =");
Pbg =  A + beta * (B - A) + gamma * (C - A);
pretty_print(tstr , LatexExpr("A + " + latex(beta) + " (B-A) + \gamma (C-A) \;\;"), "=", Pbg);

Here is the solution to the intersection between the Ray and Triangle

In [54]:
XV = matrix(SR, 3,1, (beta, gamma, t));
YV = matrix(SR, 3,1, (Av - Lv));
MM = matrix(SR, 3,3, [(Av-Bv), (Av-Cv), Dv]); 
MM = MM.transpose();
pretty_print(LatexExpr("M X = Y"), " ... " , MM, "*", XV, "=", YV);
MMs1 = copy(MM); MMs2 = copy(MM); MMs3 = copy(MM);
for i in range(3):
    MMs1[i,0] = YV[i,0];
    MMs2[i,1] = YV[i,0];
    MMs3[i,2] = YV[i,0];
detM  = MM.determinant();
detM1 = MMs1.determinant();
detM2 = MMs2.determinant();
detM3 = MMs3.determinant();
sbeta = detM1 / detM; sgamma = detM2 / detM; stval = detM3 / detM;
pretty_print(beta, "=" , sbeta, ",  ", gamma, "=" , sgamma, ",   ", t, "=" , stval);

Intersection Shown in 3D

The following shows the major compoenents of the intersection computation along with the resulting point of intersection

In [55]:
bs     = 7.0;
Pt     = Lv + stval * Dv;
Ptm    = Lv + bs * Dv;
PtBeta = Av + sbeta * (Bv - Av); 
contents = [
  polygon3d([Av, Bv, Cv], color=Color('#006633'), alpha=0.35), 
  arrow3d(Lv, (Lv + Dv), width=3),
  point(Pt, color='red', size=16),
  line3d([Lv, Ptm], thickness=5, color = Color('#f3ba29')),
  line3d([Av, PtBeta], thickness=5, color = Color('#ff4000')),
  line3d([PtBeta, Pt], thickness=5, color = Color('#40FF00')),
  line3d([(0,0,0),(bs, 0,  0)],thickness=5,color='red'),
  line3d([(0,0,0),(0, bs,  0)],thickness=5,color='green'),
  line3d([(0,0,0),(0,  0, bs)],thickness=5,color='blue')
  ];
show(sum(contents), aspect_ratio=(1,1,1), figsize=8);
In [56]:
PtBeta
Out[56]:
(2, 1, 0)
In [57]:
Pt
Out[57]:
(1, 1, 1)