Part of the Khronos Group
OpenGL.org

The Industry's Foundation for High Performance Graphics

from games to virtual reality, mobile phones to supercomputers

Results 1 to 6 of 6

Thread: Intersect triangle with segment

  1. #1
    Junior Member Regular Contributor
    Join Date
    Dec 2010
    Location
    Oakville, ON, CA
    Posts
    165

    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=Wa*A+Wb*B+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:

    Code :
    //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?

  2. #2
    Senior Member OpenGL Guru
    Join Date
    Jun 2013
    Posts
    2,402
    Quote Originally Posted by Yandersen View Post
    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
    Code :
    [ 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:
    Code :
    [ 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.

  3. #3
    Junior Member Regular Contributor
    Join Date
    Dec 2010
    Location
    Oakville, ON, CA
    Posts
    165
    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?

  4. #4
    Senior Member OpenGL Guru
    Join Date
    Jun 2013
    Posts
    2,402
    You can find the intersection of a line segment with the triangle's plane using:
    Code :
    // 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().

  5. #5
    Junior Member Regular Contributor
    Join Date
    Dec 2010
    Location
    Oakville, ON, CA
    Posts
    165
    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:
    Code :
    //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.
    Last edited by Yandersen; 08-26-2015 at 11:27 AM.

  6. #6
    Senior Member OpenGL Guru
    Join Date
    Jun 2013
    Posts
    2,402
    Another possible approach is:
    Code :
    [ 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:
    Code :
    [ 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 a*A+b*B+c*C or s*S+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:

    Code :
    N = (A-C)(B-C) // = (B-A)(C-A) = (C-B)(A-B) = (BC + CA + AB) - normal vector
    det = (S-T).N // determinant of 5x5 matrix
    a = ((S-T).BC + (B-C).ST)/det
    b = ((S-T).CA + (C-A).ST)/det
    c = ((S-T).AB + (A-B).ST)/det
    det_ABC = (AB).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
    .

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •