---

guar_intersect.cpp

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#include "stdafx.h"

#include <iostream>

#include "gul_types.h"
#include "gul_vector.h"
#include "guar_exact.h"
#include "guar_intersect.h"
#include "gul_io.h"

namespace guar {

using gul::Ptr;
using gul::point;
using gul::point2;
using gul::triangle;
using gul::triangle2;
using gul::cross_product;
using gul::test;
using gul::compare;
using std::cout;

/* ----------------------------------------------------------------------
  calculate the intersection point of two lines (in 3D). It must exist,
  or the result will be wrong.
----------------------------------------------------------------------- */
bool RegularIntersectLines(
  const point<rational>& A,     // point on line 1
  const point<rational>& B,     // direction vector of line 1
  const point<rational>& a,     // point on line 2
  const point<rational>& b,     // direction vector of line 2
  rational *lambda,  // parameter value of intersect. point (for line 1)
  rational *mu       // parameter value of intersect. point (for line 2)
)
{  
  rational A0,A1,A2,B0,B1,B2;
  rational a0,a1,a2,b0,b1,b2;
  rational N1,N2,N3;
  
  // more convenient notation
  A0 = A.x;  A1 = A.y;  A2 = A.z;
  B0 = B.x;  B1 = B.y;  B2 = B.z;
  a0 = a.x;  a1 = a.y;  a2 = a.z;
  b0 = b.x;  b1 = b.y;  b2 = b.z;
  
  N1 = B1*b0 - b1*B0;

  if( test(N1) ) // calculation without z (for exact arithmetic the
  {              // goodness of N1 doesn't matter)
    *lambda = (b0*(a1-A1) - b1*(a0-A0)) / N1; 
    *mu = (B0*(a1-A1) - B1*(a0-A0)) / N1; 
  }
  else
  {
    N2 = B2*b0 - b2*B0; 

    if( test(N2) )  // calculation without y (for exact arithmetic the
    {               // goodness of N2 doesn't matter)
      *lambda = (b0*(a2-A2) - b2*(a0-A0)) / N2;
      *mu = (B0*(a2-A2) - B2*(a0-A0)) / N2;
    }
    else
    {
      N3 = B2*b1 - b2*B1;

      if( test(N3) ) // calculation without y (for exact arithmetic the 
      {              // goodness of N3 doesn't matter)
        *lambda = (b1*(a2-A2) - b2*(a1-A1)) / N3;
        *mu = (B1*(a2-A2) - B2*(a1-A1)) / N3;
      }
      else
      {
        return false;  // no intersection point or lines colinear
      }      
    }
  }  
  return true;
}

/*-----------------------------------------------------------------------
  Estimate the (u,v) parameter values of a point P for a NURBS surface F,
  assuming that A = F(0,0), B = F(1,0), C = F(0,1)
------------------------------------------------------------------------*/
template< class T >
GULAPI void BarycentricCoordinates(
      const point<T>& P1, const point<T>& P2,
      const point<T>& P3, const point<T>& Pin,
      T *w, T *u, T *v )
{
  point2<T> W,U,V,P;
  point2<T> WU,WV,UV,WP,UP;
  point<T> UVW,PVW,PWU,PUV,P12,P13,N;
  T aUVW,aPVW,aPWU,aPUV;
  rational one_neg(ULong(1),-1);

  P12 = P2 - P1;
  P13 = P3 - P1;
  N = cross_product( P12, P13 );

  // calculate barycentric coordinates in xy-plane
  if( test(N.z) )
  {
    W.x = P1.x;
    W.y = P1.y;

    U.x = P2.x;
    U.y = P2.y;

    V.x = P3.x;
    V.y = P3.y;

    WU.x = P12.x;
    WU.y = P12.y;

    WV.x = P13.x;
    WV.y = P13.y;

    P.x = Pin.x;
    P.y = Pin.y;
  }
  else if( test(N.y) ) // calculate barycentric coordinates in xz-plane
  {
    W.x = P1.x;
    W.y = P1.z;

    U.x = P2.x;
    U.y = P2.z;

    V.x = P3.x;
    V.y = P3.z;

    WU.x = P12.x;
    WU.y = P12.z;

    WV.x = P13.x;
    WV.y = P13.z;

    P.x = Pin.x;
    P.y = Pin.z;
  }
  else // calculate barycentric coordinates in yz-plane
  {
    gul::Assert<InternalError>( ndebug || test(N.x) );

    W.x = P1.y;
    W.y = P1.z;

    U.x = P2.y;
    U.y = P2.z;

    V.x = P3.y;
    V.y = P3.z;

    WU.x = P12.y;
    WU.y = P12.z;

    WV.x = P13.y;
    WV.y = P13.z;

    P.x = Pin.y;
    P.y = Pin.z;
  }

  UVW = cross_product(WU,WV);
  aUVW = UVW.z;
  if( test(aUVW) < 0 ) aUVW = aUVW * one_neg;

  WP = P-W;
  PVW = cross_product(WP,WV);
  aPVW = PVW.z;
  if( test(aPVW) < 0 ) aPVW = aPVW * one_neg;

  *u = aPVW/aUVW;

  PWU = cross_product(WP,WU);
  aPWU = PWU.z;
  if( test(aPWU) < 0 ) aPWU = aPWU *one_neg;
   
  *v = aPWU/aUVW;

  UV = V-U;
  UP = P-U;    

  PUV = cross_product(UV,UP);
  aPUV = PUV.z;
  if( test(aPUV) < 0 ) aPUV = aPUV * one_neg;

  *w = aPUV/aUVW;
}
// template instantiation
template GULAPI void BarycentricCoordinates(
      const point<rational>& P1, const point<rational>& P2, 
      const point<rational>& P3, const point<rational>& Pin, 
      rational *w, rational *u, rational *v );

/*-----------------------------------------------------------------------
  Estimate the (u,v) parameter values of a point P for a NURBS surface F,
  assuming that A = F(0,0), B = F(1,0), C = F(0,1)
------------------------------------------------------------------------*/
template< class T >
GULAPI void BarycentricCoordinates(
      const point<T>& P1, const point<T>& P2,
      const point<T>& P3, const point<T>& P12,
      const point<T>& P13, const point<T>& P23,
      const point<T>& N, const point<T>& Pin, 
       T *w, T *u, T *v )
{
  point2<T> W,U,V,P;
  point2<T> WU,WV,UV,WP,UP;
  point<T> UVW,PVW,PWU,PUV;
  T aUVW,aPVW,aPWU,aPUV;
  rational one_neg(ULong(1),-1);

  // calculate barycentric coordinates in xy-plane
  if( test(N.z) )
  {
    W.x = P1.x;
    W.y = P1.y;

    U.x = P2.x;
    U.y = P2.y;

    V.x = P3.x;
    V.y = P3.y;

    WU.x = P12.x;
    WU.y = P12.y;

    WV.x = P13.x;
    WV.y = P13.y;

    UV.x = P23.x;
    UV.y = P23.y;

    P.x = Pin.x;
    P.y = Pin.y;
  }
  else if( test(N.y) ) // calculate barycentric coordinates in xz-plane
  {
    W.x = P1.x;
    W.y = P1.z;

    U.x = P2.x;
    U.y = P2.z;

    V.x = P3.x;
    V.y = P3.z;

    WU.x = P12.x;
    WU.y = P12.z;

    WV.x = P13.x;
    WV.y = P13.z;

    UV.x = P23.x;
    UV.y = P23.z;

    P.x = Pin.x;
    P.y = Pin.z;
  }
  else // calculate barycentric coordinates in yz-plane
  {
    gul::Assert<InternalError>( ndebug || test(N.x) );

    W.x = P1.y;
    W.y = P1.z;

    U.x = P2.y;
    U.y = P2.z;

    V.x = P3.y;
    V.y = P3.z;

    WU.x = P12.y;
    WU.y = P12.z;

    WV.x = P13.y;
    WV.y = P13.z;

    UV.x = P23.y;
    UV.y = P23.z;

    P.x = Pin.y;
    P.y = Pin.z;
  }

  /*
  cout << "N = (" << N << ")\n";
  cout << "P12 = (" << P12 << "), ";
  cout << "WU = (" << WU << ")\n";
  cout << "P13 = (" << P13 << "), ";
  cout << "WV = (" << WV << ")\n";
  */

  UVW = cross_product(WU,WV);
  aUVW = UVW.z;
  if( test(aUVW) < 0 ) aUVW = aUVW * one_neg;

  WP = P-W;
  PVW = cross_product(WP,WV);
  aPVW = PVW.z;
  if( test(aPVW) < 0 ) aPVW = aPVW * one_neg;

  *u = aPVW/aUVW;

  PWU = cross_product(WP,WU);
  aPWU = PWU.z;
  if( test(aPWU) < 0 ) aPWU = aPWU *one_neg;
   
  *v = aPWU/aUVW;

  UP = P-U;    
  PUV = cross_product(UV,UP);
  aPUV = PUV.z;
  if( test(aPUV) < 0 ) aPUV = aPUV * one_neg;

  *w = aPUV/aUVW;
}
// template instantiation
template GULAPI void BarycentricCoordinates(
      const point<rational>& P1, const point<rational>& P2, 
      const point<rational>& P3, const point<rational>& P12, 
      const point<rational>& P13, const point<rational>& P23, 
      const point<rational>& N, const point<rational>& Pin, 
      rational *w, rational *u, rational *v );


/*----------------------------------------------------------------
  calculates the endpoints of the intersection line segment of two
  triangles (the cases with both triangles being coplanar is not 
  handled).
  it also checks if the two halves of the triangles left and right
  of the intersection line are lying above or below the other
  triangle
----------------------------------------------------------------*/
GULAPI int IntersectTriangles( 
           const triangle<rational>& tri0, 
           const triangle<rational>& tri1,
           Ptr<point<rational> >& retP )
{
  point<rational>  v11,v12,v13,n,v21,v22,v23,g,t;
  point<rational>  P1,P2,P3;
  rational A,B,C,R,a,b,c,r,N;
  rational one(ULong(1)),zero;
  rational lambda,mu,l1,r1,l2,r2,l;
  bool result,init;
  int dim;

  // calculate plane equation for 1.triangle
  // (Ax + By + Cz = R)

  v11 = tri0.P2 - tri0.P1;       // 1. direction vector
  v12 = tri0.P3 - tri0.P1;       // 2. direction vector
  v13 = tri0.P3 - tri0.P2;       // (needed later)
  n = cross_product( v11, v12 ); // normal vector of the plane
  A = n.x;  B = n.y;  C = n.z;
  R = n * tri0.P1;               // right side of plane equation

  // calculate plane equation for 2.triangle
  // (ax + by + cz = r)
  v21 = tri1.P2 - tri1.P1;       // 1. direction vector
  v22 = tri1.P3 - tri1.P1;       // 2. direction vector
  v23 = tri1.P3 - tri1.P2;       // (needed later)
  n = cross_product( v21, v22 ); // normal vector of the plane
  a = n.x;  b = n.y;  c = n.z;
  r = n * tri1.P1;               // right side of plane equation

  // calculate intersection line
  // (G : x = g + lambda * t)

  // if intersection line is not parallel to the X-axis,
  // set g.x = 0 and t.x = 1
  // ==> g.y, g.z and t.x, t.z
   
  N = B*c - b*C;    // determinant
  dim = 1;          // assume intersection has dimension 1

  if( test(N) )     // if not parallel to the X-axis 
  {
    g.x = zero;
    g.y = (R*c - r*C)/N;
    g.z = (B*r - b*R)/N;
    t.x = one;
    t.y = (C*a - c*A)/N;
    t.z = (A*b - a*B)/N;    
  }
  else
  {
    N = A*c - a*C;        
    
    if( test(N) )   // if not parallel to the Y-axis
    {
      g.x = (R*c - r*C) / N;
      g.y = zero;
      g.z = (A*r - a*R) / N;
      t.x = (C*b - c*B) / N;
      t.y = one;
      t.z = (B*a - b*A) / N;
    }
    else
    {
      N = A*b - a*B;
          
      if( test(N) )   // if not parallel to the Z-axis
      {
        g.x = (R*b - r*B) / N;
        g.y = (A*r - a*R) / N;
        g.z = zero;
        t.x = (B*c - b*C) / N;
        t.y = (C*a - c*A) / N;
        t.z = one;      
      }
      else        // normal vectors of the two planes are linear dependent
      {
        // check if one point of the second triangle lies within the plane
        // of the first triangle
        if( test(A*tri1.P1.x + B*tri1.P1.y + C*tri1.P1.z - R) )
          return 0; // no intersection
        else
          dim = 2;
      }  
    }
  }

  if( dim == 1 )
  {
    // intersect intersection line of the planes with the border of the first
    // triangle. the result is an interval [l1,r1] in the parametric domain
    // of the intersection line.

    // convenient notation
    P1 = tri0.P1;
    P2 = tri0.P2;
    P3 = tri0.P3;

    // intersection line: X = g + mu * t
    // line through triangle edge P1P2:  X = P1 + lambda * v11

    init = false;
    result = RegularIntersectLines( P1, v11, g, t, &lambda, &mu );
    if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
    {
      init = true;
      l1 = mu;
      r1 = mu;
    }      

    // line through triangle edge P1P3:  X = P1 + lambda * v12

    result = RegularIntersectLines( P1, v12, g, t, &lambda, &mu );
    if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
    {
      if( init )
      {
        if( compare(mu,l1) < 0 ) l1 = mu;
        else if( compare(mu,r1) > 0 ) r1 = mu;
      }
      else
      {
        init = true;
        l1 = mu;
        r1 = mu;
      }      
    } 

    // line through triangle edge P2P3:  X = P2 + lambda * v13

    result = RegularIntersectLines( P2, v13, g, t, &lambda, &mu );
    if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
    {
      if( init )
      {
        if( compare(mu,l1) < 0 ) l1 = mu;
        else if( compare(mu,r1) > 0 ) r1 = mu;
      }
      else
      {
        init = true;
        l1 = mu;
        r1 = mu;
      }      
    } 

    if( !init )
      return 0;
    cout << "l1 = " << l1.dump() << ", r1 = " << r1.dump() << "\n";

    // intersect intersection line of the planes with the border of the second
    // triangle. the result is an interval [l2,r2] in the parametric domain
    // of the intersection line.

    // convenient notation
    P1 = tri1.P1;
    P2 = tri1.P2;
    P3 = tri1.P3;

    // intersection line: X = g + mu * t
    // line through triangle edge P1P2:  X = P1 + lambda * v21

    init = false;
    result = RegularIntersectLines( P1, v21, g, t, &lambda, &mu );
    if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
    {
      init = true;
      l2 = mu;
      r2 = mu;
    }      
  
    // line through triangle edge P1P3:  X = P1 + lambda * v22

    result = RegularIntersectLines( P1, v22, g, t, &lambda, &mu );
    if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
    {
      if( init )
      {
        if( compare(mu,l2) < 0 ) l2 = mu;
        else if( compare(mu,r2) > 0 ) r2 = mu;
      }
      else
      {
        init = true;
        l2 = mu;
        r2 = mu;
      }      
    } 

    // line through triangle edge P2P3:  X = P2 + lambda * v23

    result = RegularIntersectLines( P2, v23, g, t, &lambda, &mu );
    if( result && (test(lambda) >= 0) && (compare(lambda,one)<=0) )
    {
      if( init )
      {
        if( compare(mu,l2) < 0 ) l2 = mu;
        else if( compare(mu,r2) > 0 ) r2 = mu;
      }
      else
      {
        init = true;
        l2 = mu;
        r2 = mu;
      }      
    } 

    if( !init )
      return 0;
    cout << "l2 = " << l2.dump() << ", r2 = " << r2.dump() << "\n";

    // intersection of the intervals [l1,r1] and [l2,r2]  delivers the 
    // resulting segment of the intersection line 

    if(compare(l1,l2)>0) l = l1; else l = l2;
    if(compare(r1,r2)<0) r = r1; else r = r2;
    if( compare(r,l) < 0 )
      return 0;

    // calculate the endpoints of the intersection segment
    retP[0] = g + l * t;
    retP[1] = g + r * t;
    return 2;
  }

  return 0;
}

/*----------------------------------------------------------------
  calculates the endpoints of the intersection line segment of two
  triangles (the cases with both triangles being coplanar or just
  touching in one point is not handled).
  it also checks if the two halves of the triangles left and right
  of the intersection line are lying above or below the other
  triangle
----------------------------------------------------------------*/
GULAPI bool IntersectTrianglesUV( 
    const triangle<rational>& tri1, const triangle2<rational>& tri1uv, 
    const triangle<rational>& tri2, const triangle2<rational>& tri2uv,
    point2<rational> *S1uv, int *S1flag,
    point2<rational> *S2uv, int *S2flag,
    point<rational> *S )
{
  point<rational>  v11,v12,v13,n,N,v21,v22,v23,g,t;
  point<rational>  P1,P2,P3;
  rational A,B,C,R,a,b,c,r,T;
  rational one(ULong(1)),zero;
  rational lambda,mu,l1,r1,l2,r2,l;
  bool result,init;
  int dim;
  bool inter1_12, inter1_13, inter1_23;
  bool inter2_12, inter2_13, inter2_23;
  int  mid1,mid2,t0,t1,sign,side;
  point<rational> Pc1,Pc2;
  rational u,v,w;
  point2<rational> Pc1uv, Pc2uv;

  // calculate plane equation for 1.triangle
  // (Ax + By + Cz = R)

  v11 = tri1.P2 - tri1.P1;       // 1. direction vector
  v12 = tri1.P3 - tri1.P1;       // 2. direction vector
  v13 = tri1.P3 - tri1.P2;       // (needed later)
  N = cross_product( v11, v12 ); // normal vector of the plane
  A = N.x;  B = N.y;  C = N.z;
  R = N * tri1.P1;               // right side of plane equation
  // cout << "N = (" << N << ")\n";

  // calculate plane equation for 2.triangle
  // (ax + by + cz = r)
  v21 = tri2.P2 - tri2.P1;       // 1. direction vector
  v22 = tri2.P3 - tri2.P1;       // 2. direction vector
  v23 = tri2.P3 - tri2.P2;       // (needed later)
  n = cross_product( v21, v22 ); // normal vector of the plane
  a = n.x;  b = n.y;  c = n.z;
  r = n * tri2.P1;               // right side of plane equation

  // calculate intersection line
  // (G : x = g + lambda * t)

  // if intersection line is not parallel to the X-axis,
  // set g.x = 0 and t.x = 1
  // ==> g.y, g.z and t.x, t.z
   
  T = B*c - b*C;    // determinant
  dim = 1;          // assume intersection has dimension 1

  if( test(T) )     // if not parallel to the X-axis 
  {
    g.x = zero;
    g.y = (R*c - r*C)/T;
    g.z = (B*r - b*R)/T;
    t.x = one;
    t.y = (C*a - c*A)/T;
    t.z = (A*b - a*B)/T;    
  }
  else
  {
    T = A*c - a*C;        
    
    if( test(T) )   // if not parallel to the Y-axis
    {
      g.x = (R*c - r*C) / T;
      g.y = zero;
      g.z = (A*r - a*R) / T;
      t.x = (C*b - c*B) / T;
      t.y = one;
      t.z = (B*a - b*A) / T;
    }
    else
    {
      T = A*b - a*B;
          
      if( test(T) )   // if not parallel to the Z-axis
      {
        g.x = (R*b - r*B) / T;
        g.y = (A*r - a*R) / T;
        g.z = zero;
        t.x = (B*c - b*C) / T;
        t.y = (C*a - c*A) / T;
        t.z = one;      
      }
      else        // normal vectors of the two planes are linear dependent
      {
        // check if one point of the second triangle lies within the plane
        // of the first triangle
        if( test(A*tri2.P1.x + B*tri2.P1.y + C*tri2.P1.z - R) )
          return false; // no intersection
        else
          dim = 2;
      }  
    }
  }

  if( dim == 1 )
  {
    // intersect intersection line of the planes with the border of the first
    // triangle. the result is an interval [l1,r1] in the parametric domain
    // of the intersection line.

    mid1 = false; 
    inter1_12 = false;
    inter1_13 = false;
    inter1_23 = false;

    // convenient notation
    P1 = tri1.P1;
    P2 = tri1.P2;
    P3 = tri1.P3;

    // intersection line: X = g + mu * t
    // line through triangle edge P1P2:  X = P1 + lambda * v11

    init = false;
    result = RegularIntersectLines( P1, v11, g, t, &lambda, &mu );
    if( result )
    {
      t0 = test(lambda);
      if( t0 >= 0 )
      {
        t1 = compare(lambda,one);
        if( t1 <= 0 )
        {
          inter1_12 = true;

          init = true;
          l1 = mu;
          r1 = mu;

          if( (t0 > 0) && (t1 < 0) )
          {
            mid1 = true;
            Pc1 = tri1.P1;  // will be used to classify the two halves of
                            // tri 0 left and right of the intersection line
                            // as lying above or below tri2
          }
        }
      }
    }

    // line through triangle edge P1P3:  X = P1 + lambda * v12

    result = RegularIntersectLines( P1, v12, g, t, &lambda, &mu );
    if( result )
    {
      t0 = test(lambda);
      if( t0 >= 0 )
      {
        t1 = compare(lambda,one);
        if( t1 <= 0 )
        {
          inter1_13 = true;

          if( init )
          {
            if( compare(mu,l1) < 0 ) l1 = mu;
            else if( compare(mu,r1) > 0 ) r1 = mu;
          }
          else
          {
            init = true;
            l1 = mu;
            r1 = mu;
          }

          if( (t0 > 0) && (t1 < 0) )
          { 
            mid1 = true;
            Pc1 = tri1.P1;  // will be used to classify the two halves of
                            // tri 0 left and right of the intersection line
                            // as lying above or below tri2
          }
        }
      }
    }

    // line through triangle edge P2P3:  X = P2 + lambda * v13

    result = RegularIntersectLines( P2, v13, g, t, &lambda, &mu );
    if( result )
    {
      t0 = test(lambda);
      if( t0 >= 0 )
      {
        t1 = compare(lambda,one);
        if( t1 <= 0 )
        {
          inter1_23 = true;

          if( init )
          {
            if( compare(mu,l1) < 0 ) l1 = mu;
            else if( compare(mu,r1) > 0 ) r1 = mu;
          }
          else
          {
            init = true;
            l1 = mu;
            r1 = mu;
          }

          if( (t0 > 0) && (t1 < 0) )
          { 
            mid1 = true;
            Pc1 = tri1.P2;  // will be used to classify the two halves of
                            // tri 0 left and right of the intersection line
                            // as lying above or below tri2
          }
        }
      }
    }

    if( !init )
      return 0;
    // cout << "l1 = " << l1.dump() << ", r1 = " << r1.dump() << "\n";

    // intersect intersection line of the planes with the border of the second
    // triangle. the result is an interval [l2,r2] in the parametric domain
    // of the intersection line.

    mid2 = false;
    inter2_12 = false;
    inter2_13 = false;
    inter2_23 = false;

    // convenient notation
    P1 = tri2.P1;
    P2 = tri2.P2;
    P3 = tri2.P3;

    // intersection line: X = g + mu * t
    // line through triangle edge P1P2:  X = P1 + lambda * v21

    init = false;
    result = RegularIntersectLines( P1, v21, g, t, &lambda, &mu );
    if( result )
    {
      t0 = test(lambda);
      if( t0 >= 0 )
      {
        t1 = compare(lambda,one);
        if( t1 <= 0 )
        {
          inter2_12 = true;

          init = true;
          l2 = mu;
          r2 = mu;

          if( (t0 > 0) && (t1 < 0) )
          {
            mid2 = true;
            Pc2 = tri2.P1;   // will be used to classify the two halves of
                            // tri 0 left and right of the intersection line
                            // as lying above or below tri2
          }
        }
      }
    }
  
    // line through triangle edge P1P3:  X = P1 + lambda * v22

    result = RegularIntersectLines( P1, v22, g, t, &lambda, &mu );
    if( result )
    {
      t0 = test(lambda);
      if( t0 >= 0 )
      {
        t1 = compare(lambda,one);
        if( t1 <= 0 )
        {
          inter2_13 = true;

          if( init )
          {
            if( compare(mu,l2) < 0 ) l2 = mu;
            else if( compare(mu,r2) > 0 ) r2 = mu;
          }
          else
          {
            init = true;
            l2 = mu;
            r2 = mu;
          }

          if( (t0 > 0) && (t1 < 0) )
          {
            mid2 = true;
            Pc2 = tri2.P1;  // will be used to classify the two halves of
                            // tri 0 left and right of the intersection line
                            // as lying above or below tri2
          }
        }
      }
    }

    // line through triangle edge P2P3:  X = P2 + lambda * v23

    result = RegularIntersectLines( P2, v23, g, t, &lambda, &mu );
    if( result )
    {
      t0 = test(lambda);
      if( t0 >= 0 )
      {
        t1 = compare(lambda,one);
        if( t1 <= 0 )
        {
          inter2_23 = true;

          if( init )
          {
            if( compare(mu,l2) < 0 ) l2 = mu;
            else if( compare(mu,r2) > 0 ) r2 = mu;
          }
          else
          {
            init = true;
            l2 = mu;
            r2 = mu;
          }

          if( (t0 > 0) && (t1 < 0) )
          {
            mid2 = true;
            Pc2 = tri2.P2;  // will be used to classify the two halves of
                            // tri 0 left and right of the intersection line
                            // as lying above or below tri2
          }
        }
      }
    }

    if( !init )
      return 0;
    // cout << "l2 = " << l2.dump() << ", r2 = " << r2.dump() << "\n";

    // intersection of the intervals [l1,r1] and [l2,r2]  delivers the 
    // resulting segment of the intersection line 

    if(compare(l1,l2)>0) l = l1; else l = l2;
    if(compare(r1,r2)<0) r = r1; else r = r2;
    t1 = compare(r,l);
    if( t1 < 0 ) return false;
    else if( t1 == 0 ) // single point intersection
    {
      // retP[0] = g + l * t;
      return false;
    }

    // calculate the endpoints of the intersection segment
    S[0] = g + l * t;
    S[1] = g + r * t;

    // in triangle coordinates of triangle 1
    BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
                            S[0], &u, &v, &w );
    S1uv[0] = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;

    BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
                            S[1], &u, &v, &w );
    S1uv[1] = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;

    // calculate the edge/side pairs of the intersection segment in tri1
    if( mid1 )
    {
      sign = test(n * Pc1 - r);

      // this is a bit of a overkill since the coordinates of 'Pc1' in the 
      // parametric domain are known
      BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
                              Pc1, &u, &v, &w );
      Pc1uv = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;

      side = test( cross_product(S1uv[1]-S1uv[0],Pc1uv-S1uv[0]).z );

      gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
      if( side < 0 ) side = 0; 
      else if( side > 0 ) side = 1;
      S1flag[side] = sign;
      S1flag[!side] = -sign;
    }
    else
    {
      if( inter1_12 && inter1_13 )
        Pc1 = tri1.P1;
      else if( inter1_12 && inter1_23 )
        Pc1 = tri1.P2;
      else
      {
        gul::Assert<InternalError>( ndebug || (inter1_13 && inter1_23));
        Pc1 = tri1.P3;
      }

      sign = test(n * Pc1 - r);

      // this is a bit of a overkill since the coordinates of 'Pc1' in the 
      // parametric domain are known
      BarycentricCoordinates( tri1.P1, tri1.P2, tri1.P3, v11, v12, v13, N, 
                              Pc1, &u, &v, &w );
      Pc1uv = u*tri1uv.P1 + v*tri1uv.P2 + w*tri1uv.P3;

      side = test( cross_product(S1uv[1]-S1uv[0],Pc1uv-S1uv[0]).z );

      gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
      S1flag[side] = sign;
      S1flag[!side] = 0;  // undetermined
    }

    // in triangle coordinates of triangle 2
    BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
                            S[0], &u, &v, &w );
    S2uv[0] = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;

    BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
                            S[1], &u, &v, &w );
    S2uv[1] = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;

    // calculate the edge/side pairs of the intersection segment in tri2
    if( mid2 )
    {
      sign = test(N * Pc2 - R);

      // this is a bit of a overkill since the coordinates of 'Pc1' in the 
      // parametric domain are known
      BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
                              Pc2, &u, &v, &w );
      Pc2uv = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;

      side = test( cross_product(S2uv[1]-S2uv[0],Pc2uv-S2uv[0]).z );

      gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
      if( side < 0 ) side = 0; 
      else if( side > 0 ) side = 1;
      S2flag[side] = sign;
      S2flag[!side] = -sign;
    }
    else
    {
      if( inter2_12 && inter2_13 )
        Pc2 = tri2.P1;
      else if( inter2_12 && inter2_23 )
        Pc2 = tri2.P2;
      else
      {
        gul::Assert<InternalError>( ndebug || (inter2_13 && inter2_23));
        Pc2 = tri1.P3;
      }

      sign = test(N * Pc2 - R);

      // this is a bit of a overkill since the coordinates of 'Pc1' in the 
      // parametric domain are known
      BarycentricCoordinates( tri2.P1, tri2.P2, tri2.P3, v21, v22, v23, n, 
                              Pc2, &u, &v, &w );
      Pc2uv = u*tri2uv.P1 + v*tri2uv.P2 + w*tri2uv.P3;

      side = test( cross_product(S2uv[1]-S2uv[0],Pc2uv-S2uv[0]).z );

      gul::Assert<InternalError>( ndebug || (side != 0) && (sign != 0));
      S2flag[side] = sign;
      S2flag[!side] = 0;  // undetermined
    }

    // two intersection points
    return true;
  }

  return false;
}

}

---