BRL-CAD
tri_intersect.c
Go to the documentation of this file.
1 /* T R I _ I N T E R S E C T . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2011-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 
21 /** @file libgcv/tri_intersect.c
22  *
23  * Intersect 2 triangles using a modified Möller routine.
24  */
25 
26 #include "common.h"
27 
28 #include <math.h> /* ceil */
29 #include <string.h> /* memcpy */
30 
31 #include "vmath.h"
32 #include "bn.h"
33 #include "gcv.h"
34 #include "soup.h"
35 #include "tri_intersect.h"
36 
37 
38 /**********************************************************************/
39 /* stuff from the moller97 paper */
40 
41 
42 void
44  point_t VTX0, point_t VTX1, point_t VTX2,
45  fastf_t VV0, fastf_t VV1, fastf_t VV2,
46  fastf_t D0, fastf_t D1, fastf_t D2,
47  fastf_t *isect0, fastf_t *isect1,
48  point_t *isectpoint0, point_t *isectpoint1)
49 {
50  fastf_t tmp=D0/(D0-D1);
51  fastf_t diff[3];
52 
53  *isect0=VV0+(VV1-VV0)*tmp;
54  VSUB2(diff, VTX1, VTX0);
55  VSCALE(diff, diff, tmp);
56  VADD2(*isectpoint0, diff, VTX0);
57 
58  tmp=D0/(D0-D2);
59  *isect1=VV0+(VV2-VV0)*tmp;
60  VSUB2(diff, VTX2, VTX0);
61  VSCALE(diff, diff, tmp);
62  VADD2(*isectpoint1, VTX0, diff);
63 }
64 
65 
66 int
68  fastf_t VV0, fastf_t VV1, fastf_t VV2, fastf_t D0, fastf_t D1, fastf_t D2,
69  fastf_t D0D1, fastf_t D0D2, fastf_t *isect0, fastf_t *isect1,
70  point_t *isectpoint0, point_t *isectpoint1,
71  const struct bn_tol *tol)
72 {
73  if (D0D1>0.0f)
74  /* here we know that D0D2<=0.0 */
75  /* that is D0, D1 are on the same side, D2 on the other or on the plane */
76  gcv_fisect2(f->vert[2], f->vert[0], f->vert[1], VV2, VV0, VV1, D2, D0, D1, isect0, isect1, isectpoint0, isectpoint1);
77  else if (D0D2>0.0f)
78  /* here we know that d0d1<=0.0 */
79  gcv_fisect2(f->vert[1], f->vert[0], f->vert[2], VV1, VV0, VV2, D1, D0, D2, isect0, isect1, isectpoint0, isectpoint1);
80  else if (D1*D2>0.0f || !NEAR_ZERO(D0, tol->dist))
81  /* here we know that d0d1<=0.0 or that D0!=0.0 */
82  gcv_fisect2(f->vert[0], f->vert[1], f->vert[2], VV0, VV1, VV2, D0, D1, D2, isect0, isect1, isectpoint0, isectpoint1);
83  else if (!NEAR_ZERO(D1, tol->dist))
84  gcv_fisect2(f->vert[1], f->vert[0], f->vert[2], VV1, VV0, VV2, D1, D0, D2, isect0, isect1, isectpoint0, isectpoint1);
85  else if (!NEAR_ZERO(D2, tol->dist))
86  gcv_fisect2(f->vert[2], f->vert[0], f->vert[1], VV2, VV0, VV1, D2, D0, D1, isect0, isect1, isectpoint0, isectpoint1);
87  else
88  /* triangles are coplanar */
89  return 1;
90  return 0;
91 }
92 
93 
94 int
95 gcv_edge_edge_test(point_t V0, point_t U0, point_t U1, fastf_t Ax, fastf_t Ay, int i0, int i1)
96 {
97  fastf_t Bx, By, Cx, Cy, e, d, f;
98 
99  Bx = U0[i0] - U1[i0];
100  By = U0[i1] - U1[i1];
101  Cx = V0[i0] - U0[i0];
102  Cy = V0[i1] - U0[i1];
103  f = Ay * Bx - Ax * By;
104  d = By * Cx - Bx * Cy;
105  if ((f > 0 && d >= 0 && d <= f) || (f < 0 && d <= 0 && d >= f)) {
106  e = Ax * Cy - Ay * Cx;
107  if (f > 0) {
108  if (e >= 0 && e <= f)
109  return 1;
110  } else if (e <= 0 && e >= f)
111  return 1;
112  }
113  return 0;
114 }
115 
116 
117 int
118 gcv_edge_against_tri_edges(point_t V0, point_t V1, point_t U0, point_t U1, point_t U2, int i0, int i1)
119 {
120  fastf_t Ax, Ay;
121  Ax=V1[i0]-V0[i0];
122  Ay=V1[i1]-V0[i1];
123  /* test edge U0, U1 against V0, V1 */
124  if (gcv_edge_edge_test(V0, U0, U1, Ax, Ay, i0, i1)) return 1;
125  /* test edge U1, U2 against V0, V1 */
126  if (gcv_edge_edge_test(V0, U1, U2, Ax, Ay, i0, i1)) return 1;
127  /* test edge U2, U1 against V0, V1 */
128  if (gcv_edge_edge_test(V0, U2, U0, Ax, Ay, i0, i1)) return 1;
129  return 0;
130 }
131 
132 
133 int
134 gcv_point_in_tri(point_t V0, point_t U0, point_t U1, point_t U2, int i0, int i1)
135 {
136  fastf_t a, b, c, d0, d1, d2;
137  /* is T1 completely inside T2? */
138  /* check if V0 is inside tri(U0, U1, U2) */
139  a=U1[i1]-U0[i1];
140  b=-(U1[i0]-U0[i0]);
141  c=-a*U0[i0]-b*U0[i1];
142  d0=a*V0[i0]+b*V0[i1]+c;
143 
144  a=U2[i1]-U1[i1];
145  b=-(U2[i0]-U1[i0]);
146  c=-a*U1[i0]-b*U1[i1];
147  d1=a*V0[i0]+b*V0[i1]+c;
148 
149  a=U0[i1]-U2[i1];
150  b=-(U0[i0]-U2[i0]);
151  c=-a*U2[i0]-b*U2[i1];
152  d2=a*V0[i0]+b*V0[i1]+c;
153  if (d0*d1>0.0)
154  if (d0*d2>0.0) return 1;
155  return 0;
156 }
157 
158 
159 int
160 gcv_coplanar_tri_tri(vect_t N, vect_t V0, vect_t V1, vect_t V2, vect_t U0, vect_t U1, vect_t U2)
161 {
162  vect_t A;
163  short i0, i1;
164  /* first project onto an axis-aligned plane, that maximizes the area */
165  /* of the triangles, compute indices: i0, i1. */
166  A[0]=fabs(N[0]);
167  A[1]=fabs(N[1]);
168  A[2]=fabs(N[2]);
169  if (A[0]>A[1]) {
170  if (A[0]>A[2]) {
171  i0=1; /* A[0] is greatest */
172  i1=2;
173  } else {
174  i0=0; /* A[2] is greatest */
175  i1=1;
176  }
177  } else { /* A[0]<=A[1] */
178  if (A[2]>A[1]) {
179  i0=0; /* A[2] is greatest */
180  i1=1;
181  } else {
182  i0=0; /* A[1] is greatest */
183  i1=2;
184  }
185  }
186 
187  /* test all edges of triangle 1 against the edges of triangle 2 */
188  if (gcv_edge_against_tri_edges(V0, V1, U0, U1, U2, i0, i1))return 1;
189  if (gcv_edge_against_tri_edges(V1, V2, U0, U1, U2, i0, i1))return 1;
190  if (gcv_edge_against_tri_edges(V2, V0, U0, U1, U2, i0, i1))return 1;
191 
192  /* finally, test if tri1 is totally contained in tri2 or vice versa */
193  if (gcv_point_in_tri(V0, U0, U1, U2, i0, i1))return 1;
194  if (gcv_point_in_tri(U0, V0, V1, V2, i0, i1))return 1;
195 
196  return 0;
197 }
198 
199 
200 int
201 gcv_tri_tri_intersect_with_isectline(struct soup_s *UNUSED(left), struct soup_s *UNUSED(right), struct face_s *lf, struct face_s *rf, int *coplanar, point_t *isectpt, const struct bn_tol *tol)
202 {
203  vect_t D, isectpointA1={0}, isectpointA2={0}, isectpointB1={0}, isectpointB2={0};
204  fastf_t d1, d2, du0, du1, du2, dv0, dv1, dv2, du0du1, du0du2, dv0dv1, dv0dv2, vp0, vp1, vp2, up0, up1, up2, b, c, max, isect1[2]={0, 0}, isect2[2]={0, 0};
205  int i, smallest1=0, smallest2=0;
206 
207  /* compute plane equation of triangle(lf->vert[0], lf->vert[1], lf->vert[2]) */
208  d1=-VDOT(lf->plane, lf->vert[0]);
209  /* plane equation 1: lf->plane.X+d1=0 */
210 
211  /* put rf->vert[0], rf->vert[1], rf->vert[2] into plane equation 1 to compute signed distances to the plane*/
212  du0=VDOT(lf->plane, rf->vert[0])+d1;
213  du1=VDOT(lf->plane, rf->vert[1])+d1;
214  du2=VDOT(lf->plane, rf->vert[2])+d1;
215 
216  du0du1=du0*du1;
217  du0du2=du0*du2;
218 
219  if (du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */
220  return 0; /* no intersection occurs */
221 
222  /* compute plane of triangle (rf->vert[0], rf->vert[1], rf->vert[2]) */
223  d2=-VDOT(rf->plane, rf->vert[0]);
224  /* plane equation 2: rf->plane.X+d2=0 */
225 
226  /* put lf->vert[0], lf->vert[1], lf->vert[2] into plane equation 2 */
227  dv0=VDOT(rf->plane, lf->vert[0])+d2;
228  dv1=VDOT(rf->plane, lf->vert[1])+d2;
229  dv2=VDOT(rf->plane, lf->vert[2])+d2;
230 
231  dv0dv1=dv0*dv1;
232  dv0dv2=dv0*dv2;
233 
234  if (dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */
235  return 0; /* no intersection occurs */
236 
237  /* compute direction of intersection line */
238  VCROSS(D, lf->plane, rf->plane);
239 
240  /* compute and i to the largest component of D */
241  max=fabs(D[0]);
242  i=0;
243  b=fabs(D[1]);
244  c=fabs(D[2]);
245  if (b>max) max=b, i=1;
246  if (c>max) max=c, i=2;
247 
248  /* this is the simplified projection onto L*/
249  vp0=lf->vert[0][i];
250  vp1=lf->vert[1][i];
251  vp2=lf->vert[2][i];
252 
253  up0=rf->vert[0][i];
254  up1=rf->vert[1][i];
255  up2=rf->vert[2][i];
256 
257  /* compute interval for triangle 1 */
258  *coplanar=gcv_compute_intervals_isectline(lf, vp0, vp1, vp2, dv0, dv1, dv2, dv0dv1, dv0dv2, &isect1[0], &isect1[1], &isectpointA1, &isectpointA2, tol);
259  if (*coplanar)
260  return gcv_coplanar_tri_tri(lf->plane, lf->vert[0], lf->vert[1], lf->vert[2], rf->vert[0], rf->vert[1], rf->vert[2]);
261 
262  /* compute interval for triangle 2 */
263  gcv_compute_intervals_isectline(rf, up0, up1, up2, du0, du1, du2, du0du1, du0du2, &isect2[0], &isect2[1], &isectpointB1, &isectpointB2, tol);
264 
265  /* sort so that a<=b */
266  smallest1 = smallest2 = 0;
267 #define SORT2(a, b, smallest) if (a>b) { fastf_t _c; _c=a; a=b; b=_c; smallest=1; }
268  SORT2(isect1[0], isect1[1], smallest1);
269  SORT2(isect2[0], isect2[1], smallest2);
270 #undef SORT2
271 
272  if (isect1[1]<isect2[0] || isect2[1]<isect1[0])
273  return 0;
274 
275  /* at this point, we know that the triangles intersect */
276 
277  if (isect2[0] < isect1[0]) {
278  if (smallest1 == 0)
279  VMOVE(isectpt[0], isectpointA1);
280  else
281  VMOVE(isectpt[0], isectpointA2);
282  if (isect2[1] < isect1[1])
283  if (smallest2 == 0)
284  VMOVE(isectpt[1], isectpointB2);
285  else
286  VMOVE(isectpt[1], isectpointB1);
287  else if (smallest1 == 0)
288  VMOVE(isectpt[1], isectpointA2);
289  else
290  VMOVE(isectpt[1], isectpointA1);
291  } else {
292  if (smallest2 == 0)
293  VMOVE(isectpt[0], isectpointB1);
294  else
295  VMOVE(isectpt[0], isectpointB2);
296  if (isect2[1] > isect1[1])
297  if (smallest1 == 0)
298  VMOVE(isectpt[1], isectpointA2);
299  else
300  VMOVE(isectpt[1], isectpointA1);
301  else if (smallest2 == 0)
302  VMOVE(isectpt[1], isectpointB2);
303  else
304  VMOVE(isectpt[1], isectpointB1);
305  }
306  return 1;
307 }
308 
309 
310 /*
311  * Local Variables:
312  * tab-width: 8
313  * mode: C
314  * indent-tabs-mode: t
315  * c-file-style: "stroustrup"
316  * coding: utf-8
317  * End:
318  * ex: shiftwidth=4 tabstop=8
319  */
double dist
>= 0
Definition: tol.h:73
#define N
Definition: randmt.c:39
point_t vert[3]
Definition: soup.h:44
Header file for the BRL-CAD common definitions.
Definition: soup.h:43
#define SORT2(a, b, smallest)
plane_t plane
Definition: soup.h:45
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
void gcv_fisect2(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_intersect.c:43
int gcv_tri_tri_intersect_with_isectline(struct soup_s *left, struct soup_s *right, struct face_s *lf, struct face_s *rf, int *coplanar, point_t *isectpt, const struct bn_tol *tol)
int gcv_compute_intervals_isectline(struct face_s *f, 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, const struct bn_tol *tol)
Definition: tri_intersect.c:67
#define UNUSED(parameter)
Definition: common.h:239
Support for uniform tolerances.
Definition: tol.h:71
int gcv_coplanar_tri_tri(vect_t N, vect_t V0, vect_t V1, vect_t V2, vect_t U0, vect_t U1, vect_t U2)
int gcv_point_in_tri(point_t V0, point_t U0, point_t U1, point_t U2, int i0, int i1)
int gcv_edge_against_tri_edges(point_t V0, point_t V1, point_t U0, point_t U1, point_t U2, int i0, int i1)
int gcv_edge_edge_test(point_t V0, point_t U0, point_t U1, fastf_t Ax, fastf_t Ay, int i0, int i1)
Definition: tri_intersect.c:95
#define A
Definition: msr.c:51
double fastf_t
Definition: defines.h:300
Definition: soup.h:49