public static int intersectTriangles(Vec4[] v,Vec4[] u,Vec4[] intersectionVertices){
  double e1x=v[1].x - v[0].x;
  double e1y=v[1].y - v[0].y;
  double e1z=v[1].z - v[0].z;
  double e2x=v[2].x - v[0].x;
  double e2y=v[2].y - v[0].y;
  double e2z=v[2].z - v[0].z;
  Vec4 n1=new Vec4(e1y * e2z - e1z * e2y,e1z * e2x - e1x * e2z,e1x * e2y - e1y * e2x);
  double d1=-n1.dot3(v[0]);
  double du0=n1.dot3(u[0]) + d1;
  double du1=n1.dot3(u[1]) + d1;
  double du2=n1.dot3(u[2]) + d1;
  if (Math.abs(du0) < EPSILON)   du0=0;
  if (Math.abs(du1) < EPSILON)   du1=0;
  if (Math.abs(du2) < EPSILON)   du2=0;
  double du0du1=du0 * du1;
  double du0du2=du0 * du2;
  if (du0du1 > 0 && du0du2 > 0)   return -1;
  e1x=u[1].x - u[0].x;
  e1y=u[1].y - u[0].y;
  e1z=u[1].z - u[0].z;
  e2x=u[2].x - u[0].x;
  e2y=u[2].y - u[0].y;
  e2z=u[2].z - u[0].z;
  Vec4 n2=new Vec4(e1y * e2z - e1z * e2y,e1z * e2x - e1x * e2z,e1x * e2y - e1y * e2x);
  double d2=-n2.dot3(u[0]);
  double dv0=n2.dot3(v[0]) + d2;
  double dv1=n2.dot3(v[1]) + d2;
  double dv2=n2.dot3(v[2]) + d2;
  if (Math.abs(dv0) < EPSILON)   dv0=0;
  if (Math.abs(dv1) < EPSILON)   dv1=0;
  if (Math.abs(dv2) < EPSILON)   dv2=0;
  double dv0dv1=dv0 * dv1;
  double dv0dv2=dv0 * dv2;
  if (dv0dv1 > 0 && dv0dv2 > 0)   return -1;
  Vec4 ld=n1.cross3(n2);
  double max=Math.abs(ld.x);
  int index=0;
  double b=Math.abs(ld.y);
  double c=Math.abs(ld.z);
  if (b > max) {
    max=b;
    index=1;
  }
  if (c > max) {
    index=2;
  }
  double vp0=v[0].x;
  double vp1=v[1].x;
  double vp2=v[2].x;
  double up0=u[0].x;
  double up1=u[1].x;
  double up2=u[2].x;
  if (index == 1) {
    vp0=v[0].y;
    vp1=v[1].y;
    vp2=v[2].y;
    up0=u[0].y;
    up1=u[1].y;
    up2=u[2].y;
  }
 else   if (index == 2) {
    vp0=v[0].z;
    vp1=v[1].z;
    vp2=v[2].z;
    up0=u[0].z;
    up1=u[1].z;
    up2=u[2].z;
  }
  TriangleIntersection isectA=compute_intervals_isectline(v,vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2);
  if (isectA == null)   return coplanarTriangles(n1,v,u) ? 0 : -1;
  int smallest1=0;
  if (isectA.s0 > isectA.s1) {
    double cc=isectA.s0;
    isectA.s0=isectA.s1;
    isectA.s1=cc;
    smallest1=1;
  }
  TriangleIntersection isectB=compute_intervals_isectline(u,up0,up1,up2,du0,du1,du2,du0du1,du0du2);
  int smallest2=0;
  if (isectB.s0 > isectB.s1) {
    double cc=isectB.s0;
    isectB.s0=isectB.s1;
    isectB.s1=cc;
    smallest2=1;
  }
  if (isectA.s1 < isectB.s0 || isectB.s1 < isectA.s0)   return -1;
  if (isectB.s0 < isectA.s0) {
    if (smallest1 == 0)     intersectionVertices[0]=isectA.p0;
 else     intersectionVertices[0]=isectA.p1;
    if (isectB.s1 < isectA.s1) {
      if (smallest2 == 0)       intersectionVertices[1]=isectB.p1;
 else       intersectionVertices[1]=isectB.p0;
    }
 else {
      if (smallest1 == 0)       intersectionVertices[1]=isectA.p1;
 else       intersectionVertices[1]=isectA.p0;
    }
  }
 else {
    if (smallest2 == 0)     intersectionVertices[0]=isectB.p0;
 else     intersectionVertices[0]=isectB.p1;
    if (isectB.s1 > isectA.s1) {
      if (smallest1 == 0)       intersectionVertices[1]=isectA.p1;
 else       intersectionVertices[1]=isectA.p0;
    }
 else {
      if (smallest2 == 0)       intersectionVertices[1]=isectB.p1;
 else       intersectionVertices[1]=isectB.p0;
    }
  }
  return 1;
}

public static int intersectTriangles ( Vec4 [ ] v , Vec4 [ ] u , Vec4 [ ] intersectionVertices ) { double e1x = v [ 1 ] . x - v [ 0 ] . x ; double e1y = v [ 1 ] . y - v [ 0 ] . y ; double e1z = v [ 1 ] . z - v [ 0 ] . z ; double e2x = v [ 2 ] . x - v [ 0 ] . x ; double e2y = v [ 2 ] . y - v [ 0 ] . y ; double e2z = v [ 2 ] . z - v [ 0 ] . z ; Vec4 n1 = new Vec4 ( e1y * e2z - e1z * e2y , e1z * e2x - e1x * e2z , e1x * e2y - e1y * e2x ) ; double d1 = - n1 . dot3 ( v [ 0 ] ) ; double du0 = n1 . dot3 ( u [ 0 ] ) + d1 ; double du1 = n1 . dot3 ( u [ 1 ] ) + d1 ; double du2 = n1 . dot3 ( u [ 2 ] ) + d1 ; if ( Math . abs ( du0 ) < EPSILON ) du0 = 0 ; if ( Math . abs ( du1 ) < EPSILON ) du1 = 0 ; if ( Math . abs ( du2 ) < EPSILON ) du2 = 0 ; double du0du1 = du0 * du1 ; double du0du2 = du0 * du2 ; if ( du0du1 > 0 && du0du2 > 0 ) return - 1 ; e1x = u [ 1 ] . x - u [ 0 ] . x ; e1y = u [ 1 ] . y - u [ 0 ] . y ; e1z = u [ 1 ] . z - u [ 0 ] . z ; e2x = u [ 2 ] . x - u [ 0 ] . x ; e2y = u [ 2 ] . y - u [ 0 ] . y ; e2z = u [ 2 ] . z - u [ 0 ] . z ; Vec4 n2 = new Vec4 ( e1y * e2z - e1z * e2y , e1z * e2x - e1x * e2z , e1x * e2y - e1y * e2x ) ; double d2 = - n2 . dot3 ( u [ 0 ] ) ; double dv0 = n2 . dot3 ( v [ 0 ] ) + d2 ; double dv1 = n2 . dot3 ( v [ 1 ] ) + d2 ; double dv2 = n2 . dot3 ( v [ 2 ] ) + d2 ; if ( Math . abs ( dv0 ) < EPSILON ) dv0 = 0 ; if ( Math . abs ( dv1 ) < EPSILON ) dv1 = 0 ; if ( Math . abs ( dv2 ) < EPSILON ) dv2 = 0 ; double dv0dv1 = dv0 * dv1 ; double dv0dv2 = dv0 * dv2 ; if ( dv0dv1 > 0 && dv0dv2 > 0 ) return - 1 ; Vec4 ld = n1 . cross3 ( n2 ) ; double max = Math . abs ( ld . x ) ; int index = 0 ; double b = Math . abs ( ld . y ) ; double c = Math . abs ( ld . z ) ; if ( b > max ) { max = b ; index = 1 ; } if ( c > max ) { index = 2 ; } double vp0 = v [ 0 ] . x ; double vp1 = v [ 1 ] . x ; double vp2 = v [ 2 ] . x ; double up0 = u [ 0 ] . x ; double up1 = u [ 1 ] . x ; double up2 = u [ 2 ] . x ; if ( index == 1 ) { vp0 = v [ 0 ] . y ; vp1 = v [ 1 ] . y ; vp2 = v [ 2 ] . y ; up0 = u [ 0 ] . y ; up1 = u [ 1 ] . y ; up2 = u [ 2 ] . y ; } else if ( index == 2 ) { vp0 = v [ 0 ] . z ; vp1 = v [ 1 ] . z ; vp2 = v [ 2 ] . z ; up0 = u [ 0 ] . z ; up1 = u [ 1 ] . z ; up2 = u [ 2 ] . z ; } TriangleIntersection isectA = compute_intervals_isectline ( v , vp0 , vp1 , vp2 , dv0 , dv1 , dv2 , dv0dv1 , dv0dv2 ) ; if ( isectA == null ) return coplanarTriangles ( n1 , v , u ) ? 0 : - 1 ; int smallest1 = 0 ; if ( isectA . s0 > isectA . s1 ) { double cc = isectA . s0 ; isectA . s0 = isectA . s1 ; isectA . s1 = cc ; smallest1 = 1 ; } TriangleIntersection isectB = compute_intervals_isectline ( u , up0 , up1 , up2 , du0 , du1 , du2 , du0du1 , du0du2 ) ; int smallest2 = 0 ; if ( isectB . s0 > isectB . s1 ) { double cc = isectB . s0 ; isectB . s0 = isectB . s1 ; isectB . s1 = cc ; smallest2 = 1 ; } if ( isectA . s1 < isectB . s0 || isectB . s1 < isectA . s0 ) return - 1 ; if ( isectB . s0 < isectA . s0 ) { if ( smallest1 == 0 ) intersectionVertices [ 0 ] = isectA . p0 ; else intersectionVertices [ 0 ] = isectA . p1 ; if ( isectB . s1 < isectA . s1 ) { if ( smallest2 == 0 ) intersectionVertices [ 1 ] = isectB . p1 ; else intersectionVertices [ 1 ] = isectB . p0 ; } else { if ( smallest1 == 0 ) intersectionVertices [ 1 ] = isectA . p1 ; else intersectionVertices [ 1 ] = isectA . p0 ; } } else { if ( smallest2 == 0 ) intersectionVertices [ 0 ] = isectB . p0 ; else intersectionVertices [ 0 ] = isectB . p1 ; if ( isectB . s1 > isectA . s1 ) { if ( smallest1 == 0 ) intersectionVertices [ 1 ] = isectA . p1 ; else intersectionVertices [ 1 ] = isectA . p0 ; } else { if ( smallest2 == 0 ) intersectionVertices [ 1 ] = isectB . p1 ; else intersectionVertices [ 1 ] = isectB . p0 ; } } return 1 ; }
intersects two triangles and returns their intersection vertices .
