Intersect triangle with segment

Hi!
Basic task - to find the intersection point between a line segment with corner points S1 and S2 and a triangle with vertices A,B,C.

To do this I wrote a function that calculates partial coefficients or “weights” of each triangle’ vertices at the intersection point (Pint), so it can be found as Pint=WaA+WbB+Wc*C.

I intuitively did it through volumes, though I haven’t seen such algorithm anywhere, so I would like a confirmation if this is correct (it seems so…). Here is a working code of this function:

//Macros for calculating determinant of 3x3 matrix composed from 3 vectors
#define Determinant3xV(a,b,c) (a[0]*b[1]*c[2]+a[1]*b[2]*c[0]+a[2]*b[0]*c[1]-a[2]*b[1]*c[0]-a[1]*b[0]*c[2]-a[0]*b[2]*c[1])

//------------------------------------------------------------------------------

//Function checks if specified line' segment S1-S2 crosses triangle with given
//vertexes (A,B,C), and writes each triangle' vertex' attribute-part-value (at
//the intersection point) to PartsABC[3] variable
bool PolarizeTriangleWithSegment(const float S1[3], const float S2[3],
                                 const float A[3], const float B[3], const float C[3],
                                 float *PartsABC){
 //Calculate vectors
 float R[3] = {S2[0]-S1[0], S2[1]-S1[1], S2[2]-S1[2]},
       S1A[3] = {A[0]-S1[0], A[1]-S1[1], A[2]-S1[2]},
       S1B[3] = {B[0]-S1[0], B[1]-S1[1], B[2]-S1[2]},
       S1C[3] = {C[0]-S1[0], C[1]-S1[1], C[2]-S1[2]},
       S2A[3] = {A[0]-S2[0], A[1]-S2[1], A[2]-S2[2]},
       S2B[3] = {B[0]-S2[0], B[1]-S2[1], B[2]-S2[2]},
       S2C[3] = {C[0]-S2[0], C[1]-S2[1], C[2]-S2[2]};
 //Calculate 6X volumes
 float S1ABC = Determinant3xV(S1A,S1B,S1C);
 float S2BAC = Determinant3xV(S2B,S2A,S2C);
  if(S1ABC*S2BAC<0.f)return false; //No intersection with segment
 float S1ABCS2 = S1ABC + S2BAC;
  if(S1ABCS2==0.f)return false;    //Triangle is degenerate or S1==S2
 //Calculate volumes' parts
 PartsABC[0] = Determinant3xV(S1B,S1C,R)/S1ABCS2;
  if(PartsABC[0]<0.f)return false; //No intersection with segment
 PartsABC[1] = Determinant3xV(S1C,S1A,R)/S1ABCS2;
  if(PartsABC[1]<0.f)return false; //No intersection with segment
 PartsABC[2] = Determinant3xV(S1A,S1B,R)/S1ABCS2;
  if(PartsABC[2]<0.f)return false; //No intersection with segment
 //Weights are calculated; Pint = PartsABC[0]*A + PartsABC[1]*B + PartsABC[2]*C
 return true;
}

Is this correct?

I’m not sure whether that specific formulation is correct.

I start from the fact that the intersection of a triangle and a line segment, both expressed in barycentric coordinates, is given by


[ xa xb xc  0  0 -1  0  0 ]   [ a ]   [  0 ]
[ ya yb yc  0  0  0 -1  0 ] * [ b ] = [  0 ]
[ za zb zc  0  0  0  0 -1 ]   [ c ]   [  0 ]
[  1  1  1  0  0  0  0  0 ]   [ s ]   [  1 ]
[  0  0  0 xs xt -1  0  0 ]   [ t ]   [  0 ]
[  0  0  0 ys yt  0 -1  0 ]   [ x ]   [  0 ]
[  0  0  0 zs zt  0  0 -1 ]   [ y ]   [  0 ]
[  0  0  0  1  1  0  0  0 ]   [ z ]   [  1 ]
.

Multiplying both sides by the inverse gives:


[ a ]   [ xa xb xc  0  0 -1  0  0 ]-1   [  0 ]
[ b ] = [ ya yb yc  0  0  0 -1  0 ]   * [  0 ]
[ c ]   [ za zb zc  0  0  0  0 -1 ]     [  0 ]
[ s ]   [  1  1  1  0  0  0  0  0 ]     [  1 ]
[ t ]   [  0  0  0 xs xt -1  0  0 ]     [  0 ]
[ x ]   [  0  0  0 ys yt  0 -1  0 ]     [  0 ]
[ y ]   [  0  0  0 zs zt  0  0 -1 ]     [  0 ]
[ z ]   [  0  0  0  1  1  0  0  0 ]     [  1 ]
.

Here, [xyz][abc] are the vertices of the triangle, [xyz][st] are the vertices of the line segment, [abc] are the barycentric coordinates of the triangle, [st] are the barycentric coordinates of the line segment, and [xyz] are the spatial coordinates.

The first three rows of the matrix express the mapping between the triangle’s barycentric coordinates and spatial coordinates, while the fourth expresses the constraint that a+b+c=1. The last four rows are similar but for the line segment.

If all five barycentric coordinates [abcst] are positive, then the segment intersects the triangle. If any of [abc] are negative, the intersection between the segment’s line and the triangle’s plane is outside the triangle. If any of [st] are negative, the intersection is outside the segment. If the matrix has a zero determinant (and thus no inverse), the line through the segment and the plane of the triangle are parallel, and thus have no intersection.

Given the structure of the matrix, calculating its inverse via Cramer’s rule will result in may terms disappearing due to multiplication by zero, many terms being simplified due to multiplication by one, with the remaining non-trivial terms consisting of determinants of 3x3 matrices. And you’re only interested in the fourth and eighth columns of the inverse, as the others end up being multiplied by zero.

So determining a closed-form expression for the solution is essentially a matter of symbolically calculating the inverse of the above matrix. The determinant of an NxN matrix is the dot product between any one row or column of the matrix and the determinants of the (N-1)x(N-1) submatrices formed from the remaining rows or columns. Naturally, the first choice would be a row or column consisting entirely of zeros, followed by a row or column consisting entirely of zeros or units (i.e. 0, 1 or -1). In some cases you’ll have a choice between several such rows or columns, leading to different formulations of equal complexity.

Wow… That is too much for me, I guess… I solved the problem geometrically, through relativity of the volumes (which is calculated as a determinant of 3x3 matrix composed of three vectors forming the frustum). Anyway, if there is a code for calculating an intersection point, which is proven to be right, then it would be enough for me to see it gives the same result as mine. Could anyone share the code of finding the intersection point so I can run an output comparison test?

You can find the intersection of a line segment with the triangle’s plane using:


// triangle's vertices are A,B,C
// segment's vertices are S,T
vec3 N = cross(B - A, C - A);
float k0 = dot(N, A);
float ks = dot(N, S) - k0;
float kt = dot(N, T) - k0;
vec3 result = (kt * S - ks * T) / (kt - ks);
//

ks and kt must have opposite sign, otherwise both endpoints of the segment are on the same side of the plane. If kt-ks is close to zero, the line is almost parallel to the plane, so the result is ill-conditioned.

This doesn’t tell you whether the point of intersection lies inside the triangle, though.

GLM has glm::intersectLineTriangle().

Thank you, GClements, indeed the algorithm gives the same result for intersection point. I manually checked conditions - no_intersection_conditions reported correctly (“false” returned). I rewrote the code to GLSL-style, so here it is, just for a sake of topic completion:


//Function checks if specified line' segment S-T crosses triangle with given
//vertexes (A,B,C), and writes each triangle' vertex' attribute-part-value (at
//the intersection point) to WeightsABC[3] variable
bool IntersectTriangleWithSegment(const vec3& A, const vec3& B, const vec3& C,
                                 const vec3& S, const vec3& T,
                                 float *WeightsABC){
 //Calculate vectors
 vec3 R=T-S,
      SA=A-S, SB=B-S, SC=C-S,
      TA=A-T, TB=B-T, TC=C-T;
 //Calculate 6X volumes
 float SABC = determinant(mat3(SA,SB,SC));
 float TBAC = determinant(mat3(TB,TA,TC));
  if(SABC*TBAC<0.f)return false; //No intersection with segment
 float SABCT = SABC + TBAC;
  if(SABCT==0.f)return false;    //Triangle is degenerate or S==T
 //Calculate volumes' parts
 WeightsABC[0] = determinant(mat3(SB,TC,R))/SABCT;
 WeightsABC[1] = determinant(mat3(SC,SA,R))/SABCT;
 WeightsABC[2] = determinant(mat3(SA,SB,R))/SABCT;
 //Check if the intersection point lies outside the triangle
 if(WeightsABC[0]<0.f)return false;
 if(WeightsABC[1]<0.f)return false;
 if(WeightsABC[2]<0.f)return false;
 //Weights are calculated; Pint=WeightsABC[0]*A+WeightsABC[1]*B+WeightsABC[2]*C
 return true;
}

Though I see it can be optimized a lot, this way it just looks clearer.
So, the advantage it gives is that it allows to calculate not only the coordinates of intersection point, but also a value of interpolated attribute, like a color, or texture coordinate at the intersection point.

Another possible approach is:


[ x ]   [ xa xb xc ]   [ a ]
[ y ] = [ ya yb yc ] * [ b ]
[ z ]   [ za zb zc ]   [ c ]

[ x ]   [ xs xt ]   [ s ]
[ y ] = [ ys yt ] * [ t ]
[ z ]   [ zs zt ]

[ xa xb xc ]   [ a ]   [ xs xt ]   [ s ]
[ ya yb yc ] * [ b ] = [ ys yt ] * [ t ]
[ za zb zc ]   [ c ]   [ zs zt ]        

[ xa xb xc -xs -xt ]   [ a ]   [ 0 ]
[ ya yb yc -ys -yt ] * [ b ] = [ 0 ]
[ za zb zc -zs -zt ]   [ c ]   [ 0 ]
                       [ s ]   
                       [ t ]   


Adding the constraints a+b+c=1 and s+t=1 gives 5 equations in 5 unknowns:


[ xa xb xc -xs -xt ]   [ a ]   [ 0 ]
[ ya yb yc -ys -yt ] * [ b ] = [ 0 ]
[ za zb zc -zs -zt ]   [ c ]   [ 0 ]
[  1  1  1   0   0 ]   [ s ]   [ 1 ]
[  0  0  0   1   1 ]   [ t ]   [ 1 ]

[ a ]   [ xa xb xc -xs -xt ]-1   [ 0 ]
[ b ] = [ ya yb yc -ys -yt ]   * [ 0 ]
[ c ]   [ za zb zc -zs -zt ]     [ 0 ]
[ s ]   [  1  1  1   0   0 ]     [ 1 ]
[ t ]   [  0  0  0   1   1 ]     [ 1 ]


Expanding out the fourth and fifth columns of the inverse results in each of the 5 variables being the sum of 4 determinants of 3x3 matrices whose columns are taken from A,B,C,-S,-T, divided by the determinant of the above matrix (itself a sum of 6 such determinants).

Once a,b,c,s,t are calculated, x,y,z can be calculated as either aA+bB+cC or sS+t*T.

The determinant of a 3x3 matrix with columns A,B,C is the triple product A.BxC (= B.CxA = C.AxB). This gives rise to the following formulation:


N = (A-C)×(B-C) // = (B-A)×(C-A) = (C-B)×(A-B) = (B×C + C×A + A×B) - normal vector
det = (S-T).N // determinant of 5x5 matrix
a = ((S-T).B×C + (B-C).S×T)/det
b = ((S-T).C×A + (C-A).S×T)/det
c = ((S-T).A×B + (A-B).S×T)/det
det_ABC = (A×B).C // determinant of 3x3 matrix with columns A,B,C
s = ( det_ABC - T.N)/det
t = (-det_ABC + S.N)/det
XYZ = s*S + t*T // = a*A + b*B + c*C
.