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));
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
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
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);
The following shows the major compoenents of the intersection computation along with the resulting point of intersection
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);
PtBeta
Pt