BRL-CAD
tri_tri.c
Go to the documentation of this file.
1 /* T R I _ T R I . C
2  * BRL-CAD
3  *
4  * Published in 2012 by the United States Government.
5  * This work is in the public domain.
6  *
7  */
8 /* Triangle/triangle intersection test routine,
9  * by Tomas Moller, 1997.
10  * See article "A Fast Triangle-Triangle Intersection Test",
11  * Journal of Graphics Tools, 2(2), 1997
12  * updated: 2001-06-20 (added line of intersection)
13  *
14  * Updated June 1999: removed the divisions -- a little faster now!
15  * Updated October 1999: added {} to CROSS and SUB macros
16  *
17  * Calculate whether two coplanar triangles intersect:
18  *
19  * int bn_tri_tri_isect_coplanar(point_t V0, point_t V1, point_t V2,
20  * point_t U0, point_t U1, point_t U2,
21  * int area_flag)
22  * parameters: vertices of triangle 1: V0, V1, V2
23  * vertices of triangle 2: U0, U1, U2
24  * flag to tell function to require non-zero area: area_flag
25  * result : returns 1 if the triangles intersect, otherwise 0
26  *
27  * Calculate whether two triangles intersect:
28  *
29  * int bn_tri_tri_isect(point_t V0, point_t V1, point_t V2,
30  * point_t U0, point_t U1, point_t U2)
31  * parameters: vertices of triangle 1: V0, V1, V2
32  * vertices of triangle 2: U0, U1, U2
33  * result : returns 1 if the triangles intersect, otherwise 0
34  *
35  * This version computes the line of intersection as well (if they
36  * are not coplanar):
37  *
38  * int bn_tri_tri_isect_with_line(point_t V0, point_t V1, point_t V2,
39  * point_t U0, point_t U1, point_t U2,
40  * int *coplanar, point_t *isectpt1,
41  * point_t *isectpt2);
42  *
43  * coplanar returns whether the tris are coplanar
44  * isectpt1, isectpt2 are the endpoints of the line of intersection
45  *
46  * License: public domain (from Moller's collection of public domain code at
47  * http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/code/)
48  *
49  * The changes made for BRL-CAD integration were to use the point_t
50  * data type instead of fastf_t arrays and use vmath's vector macros
51  * instead of the locally defined versions. The function names were
52  * changed to bn_tri_tri_isect and bn_tri_tri_isect_with_line.
53  * A number of minor changes were made for C89 compatibility. BRL-CAD's
54  * NEAR_ZERO macro was used in place of exact floating point comparisons.
55  *
56  */
57 
58 #include "common.h"
59 #include <math.h>
60 #include "vmath.h"
61 #include "bn/plane_struct.h"
62 #include "bn/plane_calc.h"
63 #include "bn/tol.h"
64 #include "bn/tri_tri.h"
65 
66 /* if USE_EPSILON_TEST is true then we do a check:
67  if |dv|<EPSILON then dv=0.0;
68  else no check is done (which is less robust)
69 */
70 #define USE_EPSILON_TEST 1
71 #define EPSILON 0.000001
72 
73 /* sort so that a<=b */
74 #define SORT(a, b) \
75  if (a>b) { \
76  fastf_t c_tmp; \
77  c_tmp=a; \
78  a=b; \
79  b=c_tmp; \
80  }
81 
82 #define SORT2(a, b, smallest) \
83  if (a>b) { \
84  fastf_t c_tmp; \
85  c_tmp=a; \
86  a=b; \
87  b=c_tmp; \
88  smallest=1; \
89  } \
90  else smallest=0;
91 
92 
93 /* this edge to edge test is based on Franklin Antonio's gem:
94  "Faster Line Segment Intersection", in Graphics Gems III,
95  pp. 199-202 */
96 #define EDGE_EDGE_TEST(V0, U0, U1) \
97  Bx=U0[i0]-U1[i0]; \
98  By=U0[i1]-U1[i1]; \
99  Cx=V0[i0]-U0[i0]; \
100  Cy=V0[i1]-U0[i1]; \
101  f=Ay*Bx-Ax*By; \
102  d=By*Cx-Bx*Cy; \
103  if ((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f)) { \
104  e=Ax*Cy-Ay*Cx; \
105  if (f>0) { \
106  if (e>=0 && e<=f) return 1; \
107  } else { \
108  if (e<=0 && e>=f) return 1; \
109  } \
110  }
111 
112 #define EDGE_EDGE_TEST_AREA(V0, U0, U1) \
113  Bx=U0[i0]-U1[i0]; \
114  By=U0[i1]-U1[i1]; \
115  Cx=V0[i0]-U0[i0]; \
116  Cy=V0[i1]-U0[i1]; \
117  f=Ay*Bx-Ax*By; \
118  d=By*Cx-Bx*Cy; \
119  if ((f>0 && d>0 && d<f && !NEAR_EQUAL(d, f, EPSILON)) || (f<0 && d<0 && d>f && !NEAR_EQUAL(d, f, EPSILON))) { \
120  e=Ax*Cy-Ay*Cx; \
121  if (f>EPSILON) { \
122  if (e>EPSILON && e<f && !NEAR_EQUAL(e, f, EPSILON)) return 1; \
123  } else { \
124  if (e<EPSILON && e>f && !NEAR_EQUAL(e, f, EPSILON)) return 1; \
125  } \
126  }
127 
128 
129 #define EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2) \
130  { \
131  fastf_t Ax, Ay, Bx, By, Cx, Cy, e, d, f; \
132  Ax=V1[i0]-V0[i0]; \
133  Ay=V1[i1]-V0[i1]; \
134  /* test edge U0, U1 against V0, V1 */ \
135  EDGE_EDGE_TEST(V0, U0, U1); \
136  /* test edge U1, U2 against V0, V1 */ \
137  EDGE_EDGE_TEST(V0, U1, U2); \
138  /* test edge U2, U1 against V0, V1 */ \
139  EDGE_EDGE_TEST(V0, U2, U0); \
140  }
141 
142 
143 #define EDGE_AGAINST_TRI_EDGES_AREA(V0, V1, U0, U1, U2) \
144  { \
145  fastf_t Ax, Ay, Bx, By, Cx, Cy, e, d, f; \
146  Ax=V1[i0]-V0[i0]; \
147  Ay=V1[i1]-V0[i1]; \
148  /* test edge U0, U1 against V0, V1 */ \
149  EDGE_EDGE_TEST_AREA(V0, U0, U1); \
150  /* test edge U1, U2 against V0, V1 */ \
151  EDGE_EDGE_TEST_AREA(V0, U1, U2); \
152  /* test edge U2, U1 against V0, V1 */ \
153  EDGE_EDGE_TEST_AREA(V0, U2, U0); \
154  }
155 
156 
157 #define POINT_IN_TRI(V0, U0, U1, U2) \
158  { \
159  fastf_t a, b, c, d0, d1, d2; \
160  /* is T1 completely inside T2? */ \
161  /* check if V0 is inside tri(U0, U1, U2) */ \
162  a=U1[i1]-U0[i1]; \
163  b=-(U1[i0]-U0[i0]); \
164  c=-a*U0[i0]-b*U0[i1]; \
165  d0=a*V0[i0]+b*V0[i1]+c; \
166  \
167  a=U2[i1]-U1[i1]; \
168  b=-(U2[i0]-U1[i0]); \
169  c=-a*U1[i0]-b*U1[i1]; \
170  d1=a*V0[i0]+b*V0[i1]+c; \
171  \
172  a=U0[i1]-U2[i1]; \
173  b=-(U0[i0]-U2[i0]); \
174  c=-a*U2[i0]-b*U2[i1]; \
175  d2=a*V0[i0]+b*V0[i1]+c; \
176  if (d0*d1>0.0) { \
177  if (d0*d2>0.0) return 1; \
178  } \
179  }
180 
181 
182 int bn_tri_tri_isect_coplanar(point_t V0, point_t V1, point_t V2,
183  point_t U0, point_t U1, point_t U2, int area_flag)
184 {
185  int ret;
186  fastf_t A[3];
187  short i0, i1;
188  point_t E1, E2, N;
189  plane_t P1, P2;
190  static const struct bn_tol tol = {
191  BN_TOL_MAGIC, EPSILON, EPSILON*EPSILON, 1e-6, 1-1e-6
192  };
193 
194  /* compute plane of triangle (V0, V1, V2) */
195  ret = bn_mk_plane_3pts(P1, V0, V1, V2, &tol);
196  if (ret) return -1;
197  /* compute plane of triangle (U0, U1, U2) */
198  ret = bn_mk_plane_3pts(P2, U0, U1, U2, &tol);
199  if (ret) return -1;
200  /* verify that triangles are coplanar */
201  if (bn_coplanar(P1, P2, &tol) <= 0) return -1;
202 
203  /* first project onto an axis-aligned plane, that maximizes the area */
204  /* of the triangles, compute indices: i0, i1. */
205  VSUB2(E1, V1, V0);
206  VSUB2(E2, V2, V0);
207  VCROSS(N, E1, E2);
208 
209  A[0]=fabs(N[0]);
210  A[1]=fabs(N[1]);
211  A[2]=fabs(N[2]);
212  if (A[0]>A[1]) {
213  if (A[0]>A[2]) {
214  i0=1; /* A[0] is greatest */
215  i1=2;
216  } else {
217  i0=0; /* A[2] is greatest */
218  i1=1;
219  }
220  } else {
221  /* A[0]<=A[1] */
222  if (A[2]>A[1]) {
223  i0=0; /* A[2] is greatest */
224  i1=1;
225  } else {
226  i0=0; /* A[1] is greatest */
227  i1=2;
228  }
229  }
230 
231  /* test all edges of triangle 1 against the edges of triangle 2 */
232  if (!area_flag) {
233  EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2);
234  EDGE_AGAINST_TRI_EDGES(V1, V2, U0, U1, U2);
235  EDGE_AGAINST_TRI_EDGES(V2, V0, U0, U1, U2);
236  } else {
237  EDGE_AGAINST_TRI_EDGES_AREA(V0, V1, U0, U1, U2);
238  EDGE_AGAINST_TRI_EDGES_AREA(V1, V2, U0, U1, U2);
239  EDGE_AGAINST_TRI_EDGES_AREA(V2, V0, U0, U1, U2);
240  }
241 
242  /* finally, test if tri1 is totally contained in tri2 or vice versa */
243  POINT_IN_TRI(V0, U0, U1, U2);
244  POINT_IN_TRI(U0, V0, V1, V2);
245 
246  return 0;
247 }
248 
249 
250 /* Internal coplanar test used by more general routines. Separate from
251  * the external function because this version of the test can assume
252  * coplanar and does not need to test for area-only intersections. */
253 int coplanar_tri_tri(point_t N, point_t V0, point_t V1, point_t V2,
254  point_t U0, point_t U1, point_t U2)
255 {
256  fastf_t A[3];
257  short i0, i1;
258 
259  /* first project onto an axis-aligned plane, that maximizes the area */
260  /* of the triangles, compute indices: i0, i1. */
261  A[0]=fabs(N[0]);
262  A[1]=fabs(N[1]);
263  A[2]=fabs(N[2]);
264  if (A[0]>A[1]) {
265  if (A[0]>A[2]) {
266  i0=1; /* A[0] is greatest */
267  i1=2;
268  } else {
269  i0=0; /* A[2] is greatest */
270  i1=1;
271  }
272  } else {
273  /* A[0]<=A[1] */
274  if (A[2]>A[1]) {
275  i0=0; /* A[2] is greatest */
276  i1=1;
277  } else {
278  i0=0; /* A[1] is greatest */
279  i1=2;
280  }
281  }
282 
283  /* test all edges of triangle 1 against the edges of triangle 2 */
284  EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2);
285  EDGE_AGAINST_TRI_EDGES(V1, V2, U0, U1, U2);
286  EDGE_AGAINST_TRI_EDGES(V2, V0, U0, U1, U2);
287 
288  /* finally, test if tri1 is totally contained in tri2 or vice versa */
289  POINT_IN_TRI(V0, U0, U1, U2);
290  POINT_IN_TRI(U0, V0, V1, V2);
291 
292  return 0;
293 }
294 
295 
296 #define NEWCOMPUTE_INTERVALS(VV0, VV1, VV2, D0, D1, D2, D0D1, D0D2, A, B, C, X0, X1) \
297  { \
298  if (D0D1>0.0f) { \
299  /* here we know that D0D2<=0.0 */ \
300  /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \
301  A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \
302  } else if (D0D2>0.0f) { \
303  /* here we know that d0d1<=0.0 */ \
304  A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \
305  } else if (D1*D2>0.0f || !NEAR_ZERO(D0, SMALL_FASTF)) { \
306  /* here we know that d0d1<=0.0 or that D0!=0.0 */ \
307  A=VV0; B=(VV1-VV0)*D0; C=(VV2-VV0)*D0; X0=D0-D1; X1=D0-D2; \
308  } else if (!NEAR_ZERO(D1, SMALL_FASTF)) { \
309  A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \
310  } else if (!NEAR_ZERO(D2, SMALL_FASTF)) { \
311  A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \
312  } else { \
313  /* triangles are coplanar */ \
314  return coplanar_tri_tri(N1, V0, V1, V2, U0, U1, U2); \
315  } \
316  }
317 
318 
319 int bn_tri_tri_isect(point_t V0, point_t V1, point_t V2,
320  point_t U0, point_t U1, point_t U2)
321 {
322  point_t E1, E2;
323  point_t N1, N2;
324  fastf_t d1, d2;
325  fastf_t du0, du1, du2, dv0, dv1, dv2;
326  fastf_t D[3];
327  fastf_t isect1[2], isect2[2];
328  fastf_t du0du1, du0du2, dv0dv1, dv0dv2;
329  short index;
330  fastf_t vp0, vp1, vp2;
331  fastf_t up0, up1, up2;
332  fastf_t bb, cc, max;
333  /* parameters for triangle 1 interval computation */
334  fastf_t a, b, c, tri_x0, tri_x1;
335  /* parameters for triangle 2 interval computation */
336  fastf_t d, e, f, tri_y0, tri_y1;
337 
338  fastf_t xx, yy, xxyy, tmp;
339 
340  /* compute plane equation of triangle(V0, V1, V2) */
341  VSUB2(E1, V1, V0);
342  VSUB2(E2, V2, V0);
343  VCROSS(N1, E1, E2);
344  d1=-VDOT(N1, V0);
345  /* plane equation 1: N1.X+d1=0 */
346 
347  /* put U0, U1, U2 into plane equation 1 to compute signed distances to the plane*/
348  du0=VDOT(N1, U0)+d1;
349  du1=VDOT(N1, U1)+d1;
350  du2=VDOT(N1, U2)+d1;
351 
352  /* coplanarity robustness check */
353 #if USE_EPSILON_TEST
354  if (fabs(du0)<EPSILON) du0=0.0;
355  if (fabs(du1)<EPSILON) du1=0.0;
356  if (fabs(du2)<EPSILON) du2=0.0;
357 #endif
358  du0du1=du0*du1;
359  du0du2=du0*du2;
360 
361  if (du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */
362  return 0; /* no intersection occurs */
363 
364  /* compute plane of triangle (U0, U1, U2) */
365  VSUB2(E1, U1, U0);
366  VSUB2(E2, U2, U0);
367  VCROSS(N2, E1, E2);
368  d2=-VDOT(N2, U0);
369  /* plane equation 2: N2.X+d2=0 */
370 
371  /* put V0, V1, V2 into plane equation 2 */
372  dv0=VDOT(N2, V0)+d2;
373  dv1=VDOT(N2, V1)+d2;
374  dv2=VDOT(N2, V2)+d2;
375 
376 #if USE_EPSILON_TEST
377  if (fabs(dv0)<EPSILON) dv0=0.0;
378  if (fabs(dv1)<EPSILON) dv1=0.0;
379  if (fabs(dv2)<EPSILON) dv2=0.0;
380 #endif
381 
382  dv0dv1=dv0*dv1;
383  dv0dv2=dv0*dv2;
384 
385  if (dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */
386  return 0; /* no intersection occurs */
387 
388  /* compute direction of intersection line */
389  VCROSS(D, N1, N2);
390 
391  /* compute and index to the largest component of D */
392  max=(float)fabs(D[0]);
393  index=0;
394  bb=(float)fabs(D[1]);
395  cc=(float)fabs(D[2]);
396  if (bb>max) max=bb, index=1;
397  if (cc>max) index=2;
398 
399  /* this is the simplified projection onto L*/
400  vp0=V0[index];
401  vp1=V1[index];
402  vp2=V2[index];
403 
404  up0=U0[index];
405  up1=U1[index];
406  up2=U2[index];
407 
408  /* compute interval for triangle 1 */
409  NEWCOMPUTE_INTERVALS(vp0, vp1, vp2, dv0, dv1, dv2, dv0dv1, dv0dv2, a, b, c, tri_x0, tri_x1);
410 
411  /* compute interval for triangle 2 */
412  NEWCOMPUTE_INTERVALS(up0, up1, up2, du0, du1, du2, du0du1, du0du2, d, e, f, tri_y0, tri_y1);
413 
414  xx=tri_x0*tri_x1;
415  yy=tri_y0*tri_y1;
416  xxyy=xx*yy;
417 
418  tmp=a*xxyy;
419  isect1[0]=tmp+b*tri_x1*yy;
420  isect1[1]=tmp+c*tri_x0*yy;
421 
422  tmp=d*xxyy;
423  isect2[0]=tmp+e*xx*tri_y1;
424  isect2[1]=tmp+f*xx*tri_y0;
425 
426  SORT(isect1[0], isect1[1]);
427  SORT(isect2[0], isect2[1]);
428 
429  if (isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0;
430  return 1;
431 }
432 
433 
434 void calc_isect2(point_t VTX0, point_t VTX1, point_t VTX2, fastf_t VV0, fastf_t VV1,
435  fastf_t VV2, fastf_t D0, fastf_t D1, fastf_t D2, fastf_t *isect0,
436  fastf_t *isect1, point_t isectpoint0, point_t isectpoint1)
437 {
438  fastf_t tmp=D0/(D0-D1);
439  point_t diff;
440  *isect0=VV0+(VV1-VV0)*tmp;
441  VSUB2(diff, VTX1, VTX0);
442  VSCALE(diff, diff, tmp);
443  VADD2(isectpoint0, diff, VTX0);
444  tmp=D0/(D0-D2);
445  *isect1=VV0+(VV2-VV0)*tmp;
446  VSUB2(diff, VTX2, VTX0);
447  VSCALE(diff, diff, tmp);
448  VADD2(isectpoint1, VTX0, diff);
449 }
450 
451 
452 int compute_intervals_isectline(point_t VERT0, point_t VERT1, point_t VERT2,
453  fastf_t VV0, fastf_t VV1, fastf_t VV2, fastf_t D0, fastf_t D1, fastf_t D2,
454  fastf_t D0D1, fastf_t D0D2, fastf_t *isect0, fastf_t *isect1,
455  point_t isectpoint0, point_t isectpoint1)
456 {
457  if (D0D1>0.0f) {
458  /* here we know that D0D2<=0.0 */
459  /* that is D0, D1 are on the same side, D2 on the other or on the plane */
460  calc_isect2(VERT2, VERT0, VERT1, VV2, VV0, VV1, D2, D0, D1, isect0, isect1, isectpoint0, isectpoint1);
461  } else if (D0D2>0.0f) {
462  /* here we know that d0d1<=0.0 */
463  calc_isect2(VERT1, VERT0, VERT2, VV1, VV0, VV2, D1, D0, D2, isect0, isect1, isectpoint0, isectpoint1);
464  } else if (D1*D2>0.0f || !NEAR_ZERO(D0, SMALL_FASTF)) {
465  /* here we know that d0d1<=0.0 or that D0!=0.0 */
466  calc_isect2(VERT0, VERT1, VERT2, VV0, VV1, VV2, D0, D1, D2, isect0, isect1, isectpoint0, isectpoint1);
467  } else if (!NEAR_ZERO(D1, SMALL_FASTF)) {
468  calc_isect2(VERT1, VERT0, VERT2, VV1, VV0, VV2, D1, D0, D2, isect0, isect1, isectpoint0, isectpoint1);
469  } else if (!NEAR_ZERO(D2, SMALL_FASTF)) {
470  calc_isect2(VERT2, VERT0, VERT1, VV2, VV0, VV1, D2, D0, D1, isect0, isect1, isectpoint0, isectpoint1);
471  } else {
472  /* triangles are coplanar */
473  return 1;
474  }
475  return 0;
476 }
477 
478 
479 int bn_tri_tri_isect_with_line(point_t V0, point_t V1, point_t V2,
480  point_t U0, point_t U1, point_t U2,
481  int *coplanar, point_t *isectpt1, point_t *isectpt2)
482 {
483  point_t E1, E2;
484  point_t N1, N2;
485  fastf_t d1, d2;
486  fastf_t du0, du1, du2, dv0, dv1, dv2;
487  fastf_t D[3];
488  fastf_t isect1[2] = {0, 0};
489  fastf_t isect2[2] = {0, 0};
490  point_t isectpointA1 = VINIT_ZERO;
491  point_t isectpointA2 = VINIT_ZERO;
492  point_t isectpointB1 = VINIT_ZERO;
493  point_t isectpointB2 = VINIT_ZERO;
494  fastf_t du0du1, du0du2, dv0dv1, dv0dv2;
495  short index;
496  fastf_t vp0, vp1, vp2;
497  fastf_t up0, up1, up2;
498  fastf_t b, c, max;
499  int smallest1, smallest2;
500 
501  /* compute plane equation of triangle(V0, V1, V2) */
502  VSUB2(E1, V1, V0);
503  VSUB2(E2, V2, V0);
504  VCROSS(N1, E1, E2);
505  d1=-VDOT(N1, V0);
506  /* plane equation 1: N1.X+d1=0 */
507 
508  /* put U0, U1, U2 into plane equation 1 to compute signed distances to the plane*/
509  du0=VDOT(N1, U0)+d1;
510  du1=VDOT(N1, U1)+d1;
511  du2=VDOT(N1, U2)+d1;
512 
513  /* coplanarity robustness check */
514 #if USE_EPSILON_TEST
515  if (fabs(du0)<EPSILON) du0=0.0;
516  if (fabs(du1)<EPSILON) du1=0.0;
517  if (fabs(du2)<EPSILON) du2=0.0;
518 #endif
519  du0du1=du0*du1;
520  du0du2=du0*du2;
521 
522  if (du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */
523  return 0; /* no intersection occurs */
524 
525  /* compute plane of triangle (U0, U1, U2) */
526  VSUB2(E1, U1, U0);
527  VSUB2(E2, U2, U0);
528  VCROSS(N2, E1, E2);
529  d2=-VDOT(N2, U0);
530  /* plane equation 2: N2.X+d2=0 */
531 
532  /* put V0, V1, V2 into plane equation 2 */
533  dv0=VDOT(N2, V0)+d2;
534  dv1=VDOT(N2, V1)+d2;
535  dv2=VDOT(N2, V2)+d2;
536 
537 #if USE_EPSILON_TEST
538  if (fabs(dv0)<EPSILON) dv0=0.0;
539  if (fabs(dv1)<EPSILON) dv1=0.0;
540  if (fabs(dv2)<EPSILON) dv2=0.0;
541 #endif
542 
543  dv0dv1=dv0*dv1;
544  dv0dv2=dv0*dv2;
545 
546  if (dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */
547  return 0; /* no intersection occurs */
548 
549  /* compute direction of intersection line */
550  VCROSS(D, N1, N2);
551 
552  /* compute and index to the largest component of D */
553  max=fabs(D[0]);
554  index=0;
555  b=fabs(D[1]);
556  c=fabs(D[2]);
557  if (b>max) max=b, index=1;
558  if (c>max) index=2;
559 
560  /* this is the simplified projection onto L*/
561  vp0=V0[index];
562  vp1=V1[index];
563  vp2=V2[index];
564 
565  up0=U0[index];
566  up1=U1[index];
567  up2=U2[index];
568 
569  /* compute interval for triangle 1 */
570  *coplanar=compute_intervals_isectline(V0, V1, V2, vp0, vp1, vp2, dv0, dv1, dv2,
571  dv0dv1, dv0dv2, &isect1[0], &isect1[1], isectpointA1, isectpointA2);
572  if (*coplanar) return coplanar_tri_tri(N1, V0, V1, V2, U0, U1, U2);
573 
574 
575  /* compute interval for triangle 2 */
576  compute_intervals_isectline(U0, U1, U2, up0, up1, up2, du0, du1, du2,
577  du0du1, du0du2, &isect2[0], &isect2[1], isectpointB1, isectpointB2);
578 
579  SORT2(isect1[0], isect1[1], smallest1);
580  SORT2(isect2[0], isect2[1], smallest2);
581 
582  if (isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0;
583 
584  /* at this point, we know that the triangles intersect */
585 
586  if (isect2[0]<isect1[0]) {
587  if (smallest1==0) { VMOVE(*isectpt1, isectpointA1); } else { VMOVE(*isectpt1, isectpointA2); }
588 
589  if (isect2[1]<isect1[1]) {
590  if (smallest2==0) { VMOVE(*isectpt2, isectpointB2); } else { VMOVE(*isectpt2, isectpointB1); }
591  } else {
592  if (smallest1==0) { VMOVE(*isectpt2, isectpointA2); } else { VMOVE(*isectpt2, isectpointA1); }
593  }
594  } else {
595  if (smallest2==0) { VMOVE(*isectpt1, isectpointB1); } else { VMOVE(*isectpt1, isectpointB2); }
596 
597  if (isect2[1]>isect1[1]) {
598  if (smallest1==0) { VMOVE(*isectpt2, isectpointA2); } else { VMOVE(*isectpt2, isectpointA1); }
599  } else {
600  if (smallest2==0) { VMOVE(*isectpt2, isectpointB2); } else { VMOVE(*isectpt2, isectpointB1); }
601  }
602  }
603  return 1;
604 }
605 
606 
607 /*
608  * Local Variables:
609  * tab-width: 8
610  * mode: C
611  * indent-tabs-mode: t
612  * c-file-style: "stroustrup"
613  * End:
614  * ex: shiftwidth=4 tabstop=8
615  */
int compute_intervals_isectline(point_t VERT0, point_t VERT1, point_t VERT2, fastf_t VV0, fastf_t VV1, fastf_t VV2, fastf_t D0, fastf_t D1, fastf_t D2, fastf_t D0D1, fastf_t D0D2, fastf_t *isect0, fastf_t *isect1, point_t isectpoint0, point_t isectpoint1)
Definition: tri_tri.c:452
int bn_coplanar(const plane_t a, const plane_t b, const struct bn_tol *tol)
Test if two planes are identical. If so, their dot products will be either +1 or -1, with the distance from the origin equal in magnitude.
#define EDGE_AGAINST_TRI_EDGES(V0, V1, U0, U1, U2)
Definition: tri_tri.c:129
void calc_isect2(point_t VTX0, point_t VTX1, point_t VTX2, fastf_t VV0, fastf_t VV1, fastf_t VV2, fastf_t D0, fastf_t D1, fastf_t D2, fastf_t *isect0, fastf_t *isect1, point_t isectpoint0, point_t isectpoint1)
Definition: tri_tri.c:434
#define N
Definition: randmt.c:39
#define BN_TOL_MAGIC
Definition: magic.h:74
#define SMALL_FASTF
Definition: defines.h:342
Header file for the BRL-CAD common definitions.
#define EPSILON
Definition: tri_tri.c:71
int bn_tri_tri_isect(point_t V0, point_t V1, point_t V2, point_t U0, point_t U1, point_t U2)
Definition: tri_tri.c:319
#define EDGE_AGAINST_TRI_EDGES_AREA(V0, V1, U0, U1, U2)
Definition: tri_tri.c:143
#define SORT(a, b)
Definition: tri_tri.c:74
#define POINT_IN_TRI(V0, U0, U1, U2)
Definition: tri_tri.c:157
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define SORT2(a, b, smallest)
Definition: tri_tri.c:82
Support for uniform tolerances.
Definition: tol.h:71
#define NEWCOMPUTE_INTERVALS(VV0, VV1, VV2, D0, D1, D2, D0D1, D0D2, A, B, C, X0, X1)
Definition: tri_tri.c:296
int bn_tri_tri_isect_coplanar(point_t V0, point_t V1, point_t V2, point_t U0, point_t U1, point_t U2, int area_flag)
Tomas Möller's triangle/triangle intersection routines from the article.
Definition: tri_tri.c:182
#define A
Definition: msr.c:51
int coplanar_tri_tri(point_t N, point_t V0, point_t V1, point_t V2, point_t U0, point_t U1, point_t U2)
Definition: tri_tri.c:253
int bn_tri_tri_isect_with_line(point_t V0, point_t V1, point_t V2, point_t U0, point_t U1, point_t U2, int *coplanar, point_t *isectpt1, point_t *isectpt2)
Definition: tri_tri.c:479
double fastf_t
Definition: defines.h:300
int bn_mk_plane_3pts(plane_t plane, const point_t a, const point_t b, const point_t c, const struct bn_tol *tol)