API tools faq
paste
mrDIMAS

Triangle-triangle Intersection

mrDIMAS
Jan 8th, 2017
355
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 9.69 KB | None | 0 0
raw download clone embed print report
  1. /* Triangle/triangle intersection test routine,
  2.  * by Tomas Moller, 1997.
  3.  * See article "A Fast Triangle-Triangle Intersection Test",
  4.  * Journal of Graphics Tools, 2(2), 1997
  5.  *
  6.  * int tri_tri_intersect(float V0[3],float V1[3],float V2[3],
  7.  *                         float U0[3],float U1[3],float U2[3])
  8.  *
  9.  * parameters: vertices of triangle 1: V0,V1,V2
  10.  *             vertices of triangle 2: U0,U1,U2
  11.  * result    : returns 1 if the triangles intersect, otherwise 0
  12.  *
  13.  */
  14.  
  15. #include <math.h>
  16.  
  17.  
  18. /* if USE_EPSILON_TEST is true then we do a check:
  19.          if |dv|<EPSILON then dv=0.0;
  20.    else no check is done (which is less robust)
  21. */
  22. #define USE_EPSILON_TEST TRUE  
  23. #define EPSILON 0.000001
  24.  
  25.  
  26. /* some macros */
  27. #define CROSS(dest,v1,v2)                      \
  28.               dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
  29.               dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
  30.               dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
  31.  
  32. #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
  33.  
  34. #define SUB(dest,v1,v2)          \
  35.             dest[0]=v1[0]-v2[0]; \
  36.             dest[1]=v1[1]-v2[1]; \
  37.             dest[2]=v1[2]-v2[2];
  38.  
  39. /* sort so that a<=b */
  40. #define SORT(a,b)       \
  41.              if(a>b)    \
  42.              {          \
  43.                float c; \
  44.                c=a;     \
  45.                a=b;     \
  46.                b=c;     \
  47.              }
  48.  
  49. #define ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1) \
  50.               isect0=VV0+(VV1-VV0)*D0/(D0-D1);    \
  51.               isect1=VV0+(VV2-VV0)*D0/(D0-D2);
  52.  
  53.  
  54. #define COMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,isect0,isect1) \
  55.   if(D0D1>0.0f)                                         \
  56.   {                                                     \
  57.     /* here we know that D0D2<=0.0 */                   \
  58.     /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \
  59.     ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1);          \
  60.   }                                                     \
  61.   else if(D0D2>0.0f)                                    \
  62.   {                                                     \
  63.     /* here we know that d0d1<=0.0 */                   \
  64.     ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1);          \
  65.   }                                                     \
  66.   else if(D1*D2>0.0f || D0!=0.0f)                       \
  67.   {                                                     \
  68.     /* here we know that d0d1<=0.0 or that D0!=0.0 */   \
  69.     ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1);          \
  70.   }                                                     \
  71.   else if(D1!=0.0f)                                     \
  72.   {                                                     \
  73.     ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1);          \
  74.   }                                                     \
  75.   else if(D2!=0.0f)                                     \
  76.   {                                                     \
  77.     ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1);          \
  78.   }                                                     \
  79.   else                                                  \
  80.   {                                                     \
  81.     /* triangles are coplanar */                        \
  82.     return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2);      \
  83.   }
  84.  
  85.  
  86.  
  87. /* this edge to edge test is based on Franlin Antonio's gem:
  88.    "Faster Line Segment Intersection", in Graphics Gems III,
  89.    pp. 199-202 */
  90. #define EDGE_EDGE_TEST(V0,U0,U1)                      \
  91.   Bx=U0[i0]-U1[i0];                                   \
  92.   By=U0[i1]-U1[i1];                                   \
  93.   Cx=V0[i0]-U0[i0];                                   \
  94.   Cy=V0[i1]-U0[i1];                                   \
  95.   f=Ay*Bx-Ax*By;                                      \
  96.   d=By*Cx-Bx*Cy;                                      \
  97.   if((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f))  \
  98.   {                                                   \
  99.     e=Ax*Cy-Ay*Cx;                                    \
  100.     if(f>0)                                           \
  101.     {                                                 \
  102.       if(e>=0 && e<=f) return 1;                      \
  103.     }                                                 \
  104.     else                                              \
  105.     {                                                 \
  106.       if(e<=0 && e>=f) return 1;                      \
  107.     }                                                 \
  108.   }                                
  109.  
  110. #define EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2) \
  111. {                                              \
  112.   float Ax,Ay,Bx,By,Cx,Cy,e,d,f;               \
  113.   Ax=V1[i0]-V0[i0];                            \
  114.   Ay=V1[i1]-V0[i1];                            \
  115.   /* test edge U0,U1 against V0,V1 */          \
  116.   EDGE_EDGE_TEST(V0,U0,U1);                    \
  117.   /* test edge U1,U2 against V0,V1 */          \
  118.   EDGE_EDGE_TEST(V0,U1,U2);                    \
  119.   /* test edge U2,U1 against V0,V1 */          \
  120.   EDGE_EDGE_TEST(V0,U2,U0);                    \
  121. }
  122.  
  123. #define POINT_IN_TRI(V0,U0,U1,U2)           \
  124. {                                           \
  125.   float a,b,c,d0,d1,d2;                     \
  126.   /* is T1 completly inside T2? */          \
  127.   /* check if V0 is inside tri(U0,U1,U2) */ \
  128.   a=U1[i1]-U0[i1];                          \
  129.   b=-(U1[i0]-U0[i0]);                       \
  130.   c=-a*U0[i0]-b*U0[i1];                     \
  131.   d0=a*V0[i0]+b*V0[i1]+c;                   \
  132.                                             \
  133.   a=U2[i1]-U1[i1];                          \
  134.   b=-(U2[i0]-U1[i0]);                       \
  135.   c=-a*U1[i0]-b*U1[i1];                     \
  136.   d1=a*V0[i0]+b*V0[i1]+c;                   \
  137.                                             \
  138.   a=U0[i1]-U2[i1];                          \
  139.   b=-(U0[i0]-U2[i0]);                       \
  140.   c=-a*U2[i0]-b*U2[i1];                     \
  141.   d2=a*V0[i0]+b*V0[i1]+c;                   \
  142.   if(d0*d1>0.0)                             \
  143.   {                                         \
  144.     if(d0*d2>0.0) return 1;                 \
  145.   }                                         \
  146. }
  147.  
  148. int coplanar_tri_tri(float N[3],float V0[3],float V1[3],float V2[3],
  149.                      float U0[3],float U1[3],float U2[3])
  150. {
  151.    float A[3];
  152.    short i0,i1;
  153.    /* first project onto an axis-aligned plane, that maximizes the area */
  154.    /* of the triangles, compute indices: i0,i1. */
  155.    A[0]=fabs(N[0]);
  156.    A[1]=fabs(N[1]);
  157.    A[2]=fabs(N[2]);
  158.    if(A[0]>A[1])
  159.    {
  160.       if(A[0]>A[2])  
  161.       {
  162.           i0=1;      /* A[0] is greatest */
  163.           i1=2;
  164.       }
  165.       else
  166.       {
  167.           i0=0;      /* A[2] is greatest */
  168.           i1=1;
  169.       }
  170.    }
  171.    else   /* A[0]<=A[1] */
  172.    {
  173.       if(A[2]>A[1])
  174.       {
  175.           i0=0;      /* A[2] is greatest */
  176.           i1=1;                                          
  177.       }
  178.       else
  179.       {
  180.           i0=0;      /* A[1] is greatest */
  181.           i1=2;
  182.       }
  183.     }              
  184.                
  185.     /* test all edges of triangle 1 against the edges of triangle 2 */
  186.     EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2);
  187.     EDGE_AGAINST_TRI_EDGES(V1,V2,U0,U1,U2);
  188.     EDGE_AGAINST_TRI_EDGES(V2,V0,U0,U1,U2);
  189.                
  190.     /* finally, test if tri1 is totally contained in tri2 or vice versa */
  191.     POINT_IN_TRI(V0,U0,U1,U2);
  192.     POINT_IN_TRI(U0,V0,V1,V2);
  193.  
  194.     return 0;
  195. }
  196.  
  197.  
  198. int tri_tri_intersect(float V0[3],float V1[3],float V2[3],
  199.                       float U0[3],float U1[3],float U2[3])
  200. {
  201.   float E1[3],E2[3];
  202.   float N1[3],N2[3],d1,d2;
  203.   float du0,du1,du2,dv0,dv1,dv2;
  204.   float D[3];
  205.   float isect1[2], isect2[2];
  206.   float du0du1,du0du2,dv0dv1,dv0dv2;
  207.   short index;
  208.   float vp0,vp1,vp2;
  209.   float up0,up1,up2;
  210.   float b,c,max;
  211.  
  212.   /* compute plane equation of triangle(V0,V1,V2) */
  213.   SUB(E1,V1,V0);
  214.   SUB(E2,V2,V0);
  215.   CROSS(N1,E1,E2);
  216.   d1=-DOT(N1,V0);
  217.   /* plane equation 1: N1.X+d1=0 */
  218.  
  219.   /* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/
  220.   du0=DOT(N1,U0)+d1;
  221.   du1=DOT(N1,U1)+d1;
  222.   du2=DOT(N1,U2)+d1;
  223.  
  224.   /* coplanarity robustness check */
  225. #if USE_EPSILON_TEST==TRUE
  226.   if(fabs(du0)<EPSILON) du0=0.0;
  227.   if(fabs(du1)<EPSILON) du1=0.0;
  228.   if(fabs(du2)<EPSILON) du2=0.0;
  229. #endif
  230.   du0du1=du0*du1;
  231.   du0du2=du0*du2;
  232.  
  233.   if(du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */
  234.     return 0;                    /* no intersection occurs */
  235.  
  236.   /* compute plane of triangle (U0,U1,U2) */
  237.   SUB(E1,U1,U0);
  238.   SUB(E2,U2,U0);
  239.   CROSS(N2,E1,E2);
  240.   d2=-DOT(N2,U0);
  241.   /* plane equation 2: N2.X+d2=0 */
  242.  
  243.   /* put V0,V1,V2 into plane equation 2 */
  244.   dv0=DOT(N2,V0)+d2;
  245.   dv1=DOT(N2,V1)+d2;
  246.   dv2=DOT(N2,V2)+d2;
  247.  
  248. #if USE_EPSILON_TEST==TRUE
  249.   if(fabs(dv0)<EPSILON) dv0=0.0;
  250.   if(fabs(dv1)<EPSILON) dv1=0.0;
  251.   if(fabs(dv2)<EPSILON) dv2=0.0;
  252. #endif
  253.  
  254.   dv0dv1=dv0*dv1;
  255.   dv0dv2=dv0*dv2;
  256.        
  257.   if(dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */
  258.     return 0;                    /* no intersection occurs */
  259.  
  260.   /* compute direction of intersection line */
  261.   CROSS(D,N1,N2);
  262.  
  263.   /* compute and index to the largest component of D */
  264.   max=fabs(D[0]);
  265.   index=0;
  266.   b=fabs(D[1]);
  267.   c=fabs(D[2]);
  268.   if(b>max) max=b,index=1;
  269.   if(c>max) max=c,index=2;
  270.  
  271.         /* this is the simplified projection onto L*/
  272.         vp0=V0[index];
  273.         vp1=V1[index];
  274.         vp2=V2[index];
  275.  
  276.         up0=U0[index];
  277.         up1=U1[index];
  278.         up2=U2[index];
  279.  
  280.   /* compute interval for triangle 1 */
  281.   COMPUTE_INTERVALS(vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2,isect1[0],isect1[1]);
  282.  
  283.   /* compute interval for triangle 2 */
  284.   COMPUTE_INTERVALS(up0,up1,up2,du0,du1,du2,du0du1,du0du2,isect2[0],isect2[1]);
  285.  
  286.   SORT(isect1[0],isect1[1]);
  287.   SORT(isect2[0],isect2[1]);
  288.  
  289.   if(isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0;
  290.   return 1;
  291. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!