Here is a solution to the ray-triangle intersection problem that does not require explicitly intersecting the ray with a full 3D plane. This approach instead takes advantage of a simple property of triangles, namely that any point in a triangle may be expressed as some offset from a key vertex in directions determined by vectors derived from this key vertex and the remaining two vertices.
Ross Beveridge, September 20, 2018
In this example the actual scalar variables are being made explicit, hence the next setup block is rather large in so much as it declares all the necessary symbolic scalar variables. We will use the symbols more below and offer a fuller explanation, but for the moment, note the conventions match the PowerPoint presentation with a triangle defined by 3D points A, B and C.
latex.matrix_delimiters("[", "]")
var('ax', 'ay', 'az', 'bx', 'by', 'bz', 'cx', 'cy', 'cz');
var('lx', 'ly', 'lz', 'dx', 'dy', 'dz');
var('t', 'beta', 'gamma');
Av = vector(SR, 3, ('ax', 'ay', 'az'));
Bv = vector(SR, 3, ('bx', 'by', 'bz'));
Cv = vector(SR, 3, ('cx', 'cy', 'cz'));
Lv = vector(SR, 3, ('lx', 'ly', 'lz'));
Dv = vector(SR, 3, ('dx', 'dy', 'dz'));
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);
The ray is parameterized in the standard manner, by a single value t. The ray originates at a point L and proceeds in a direction defined by a nunit length vector D.
Rt = L + t * D;
pretty_print(LatexExpr("R(t) = L + t D ="), L, "+", t, "*", D, "=", Rt);
The triangle is defined by three points: A, B and C. Any position on the infinite plane containing A, B and C may be reached using two parameters, one point (A) and two vectors (B-A) and (C-A). Do notice that early in the semester the importance of understanding what is meant by a Point versus a Vector in geometric terms is key to understanding this development.
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) \;\;"))
pretty_print(tstr , "=", Pbg);
In the previous lecture we saw how in 2D to solve for the t and s values where tow parameterized lines crossed. One key step is to create a single equation where the point on one object equals a point on the other object. The same approach is taken here, we need to crreate a single linear algebraic equation as follows:
pretty_print(Rt, " = ", Pbg);
Unless the geometry is very badly conditioned, namely the ray is parallel to the plane of the triangle, there will exist a combination of t value for the ray along with alpha and gamma for the plane such that the equation above is true. The next step is to do some re-writing of the eqaution to place it into a form of a standard linear algebraic equation of a kind we are used to solving.
pretty_print(Rt, " + ", -1 * Pbg, "=", matrix(SR, 3,1, 0));
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);
The next question, how most efficiently to solve for the parameters beta, gamma and t is complicated. One option is to find a simple linear algebra library that includes matrix inverse.
Another option is to work out Cramer's Rule. This approach involves solving for the determinants of four matrices, the original M matrix and then three derived by substituting Y into successive columns.
Below is the brute force application of Cramer's Rule in order to solve for the parameters. While this is a good start, keep in mind when actually implementing this approach that there is a lot of repitition and common terms should be merged. The most obvious example is the determinant of the matrix M which should only be computed once.
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];
pretty_print("M = ", MM);
pretty_print(LatexExpr("M_{1} = \;"), MMs1, LatexExpr(", \;\; M_{2} = \;"), MMs2, LatexExpr(", \;\; M_{3} = \;"), MMs3);
detM = MM.determinant();
detM1 = MMs1.determinant();
detM2 = MMs2.determinant();
detM3 = MMs3.determinant();
sbeta = detM1 / detM; sgamma = detM2 / detM; stval = detM3 / detM;
pretty_print(beta, LatexExpr("\;=\; \\frac{|M_{1}|}{|M|}"));
pretty_print(gamma, LatexExpr("\;=\; \\frac{|M_{2}|}{|M|}"));
pretty_print(t, LatexExpr("\;=\; \\frac{|M_{3}|}{|M|}"));
pretty_print(beta, LatexExpr("\;=\;"), sbeta);
pretty_print(gamma, LatexExpr("\;=\;"), sgamma);
pretty_print(t, LatexExpr("\;=\;"), stval);