Wayverb
tri_cube_c.h
1 #pragma once
2 
3 // from
4 // http://www.realtimerendering.com/resources/GraphicsGems/gemsiii/triangleCube.c
5 
6 #include <math.h>
7 
8 /* this version of SIGN3 shows some numerical instability, and is improved
9  * by using the uncommented macro that follows, and a different test with it */
10 #ifdef OLD_TEST
11 #define SIGN3(A) \
12  (((A).x < 0) ? 4 : 0 | ((A).y < 0) ? 2 : 0 | ((A).z < 0) ? 1 : 0)
13 #else
14 #define EPS 10e-5
15 #define SIGN3(A) \
16  ((((A).x < EPS) ? 4 : 0) | (((A).x > -EPS) ? 32 : 0) | \
17  (((A).y < EPS) ? 2 : 0) | (((A).y > -EPS) ? 16 : 0) | \
18  (((A).z < EPS) ? 1 : 0) | (((A).z > -EPS) ? 8 : 0))
19 #endif
20 
21 #define CROSS(A, B, C) \
22  { \
23  (C).x = (A).y * (B).z - (A).z * (B).y; \
24  (C).y = -(A).x * (B).z + (A).z * (B).x; \
25  (C).z = (A).x * (B).y - (A).y * (B).x; \
26  }
27 #define SUB(A, B, C) \
28  { \
29  (C).x = (A).x - (B).x; \
30  (C).y = (A).y - (B).y; \
31  (C).z = (A).z - (B).z; \
32  }
33 #define LERP(A, B, C) ((B) + (A) * ((C) - (B)))
34 #define MIN3(a, b, c) \
35  ((((a) < (b)) && ((a) < (c))) ? (a) : (((b) < (c)) ? (b) : (c)))
36 #define MAX3(a, b, c) \
37  ((((a) > (b)) && ((a) > (c))) ? (a) : (((b) > (c)) ? (b) : (c)))
38 #define INSIDE 0
39 #define OUTSIDE 1
40 
41 typedef struct {
42  float x;
43  float y;
44  float z;
45 } Point3;
46 
47 typedef struct {
48  Point3 v1; /* Vertex1 */
49  Point3 v2; /* Vertex2 */
50  Point3 v3; /* Vertex3 */
51 } Triangle3;
52 
53 /*___________________________________________________________________________*/
54 
55 /* Which of the six face-plane(s) is point P outside of? */
56 
57 long face_plane(Point3 p) {
58  long outcode;
59 
60  outcode = 0;
61  if (p.x > .5)
62  outcode |= 0x01;
63  if (p.x < -.5)
64  outcode |= 0x02;
65  if (p.y > .5)
66  outcode |= 0x04;
67  if (p.y < -.5)
68  outcode |= 0x08;
69  if (p.z > .5)
70  outcode |= 0x10;
71  if (p.z < -.5)
72  outcode |= 0x20;
73  return (outcode);
74 }
75 
76 /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
77 
78 /* Which of the twelve edge plane(s) is point P outside of? */
79 
80 long bevel_2d(Point3 p) {
81  long outcode;
82 
83  outcode = 0;
84  if (p.x + p.y > 1.0)
85  outcode |= 0x001;
86  if (p.x - p.y > 1.0)
87  outcode |= 0x002;
88  if (-p.x + p.y > 1.0)
89  outcode |= 0x004;
90  if (-p.x - p.y > 1.0)
91  outcode |= 0x008;
92  if (p.x + p.z > 1.0)
93  outcode |= 0x010;
94  if (p.x - p.z > 1.0)
95  outcode |= 0x020;
96  if (-p.x + p.z > 1.0)
97  outcode |= 0x040;
98  if (-p.x - p.z > 1.0)
99  outcode |= 0x080;
100  if (p.y + p.z > 1.0)
101  outcode |= 0x100;
102  if (p.y - p.z > 1.0)
103  outcode |= 0x200;
104  if (-p.y + p.z > 1.0)
105  outcode |= 0x400;
106  if (-p.y - p.z > 1.0)
107  outcode |= 0x800;
108  return (outcode);
109 }
110 
111 /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
112 
113 /* Which of the eight corner plane(s) is point P outside of? */
114 
115 long bevel_3d(Point3 p) {
116  long outcode;
117 
118  outcode = 0;
119  if ((p.x + p.y + p.z) > 1.5)
120  outcode |= 0x01;
121  if ((p.x + p.y - p.z) > 1.5)
122  outcode |= 0x02;
123  if ((p.x - p.y + p.z) > 1.5)
124  outcode |= 0x04;
125  if ((p.x - p.y - p.z) > 1.5)
126  outcode |= 0x08;
127  if ((-p.x + p.y + p.z) > 1.5)
128  outcode |= 0x10;
129  if ((-p.x + p.y - p.z) > 1.5)
130  outcode |= 0x20;
131  if ((-p.x - p.y + p.z) > 1.5)
132  outcode |= 0x40;
133  if ((-p.x - p.y - p.z) > 1.5)
134  outcode |= 0x80;
135  return (outcode);
136 }
137 
138 /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
139 
140 /* Test the point "alpha" of the way from P1 to P2 */
141 /* See if it is on a face of the cube */
142 /* Consider only faces in "mask" */
143 
144 long check_point(Point3 p1, Point3 p2, float alpha, long mask) {
145  Point3 plane_point;
146 
147  plane_point.x = LERP(alpha, p1.x, p2.x);
148  plane_point.y = LERP(alpha, p1.y, p2.y);
149  plane_point.z = LERP(alpha, p1.z, p2.z);
150  return (face_plane(plane_point) & mask);
151 }
152 
153 /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
154 
155 /* Compute intersection of P1 --> P2 line segment with face planes */
156 /* Then test intersection point to see if it is on cube face */
157 /* Consider only face planes in "outcode_diff" */
158 /* Note: Zero bits in "outcode_diff" means face line is outside of */
159 
160 long check_line(Point3 p1, Point3 p2, long outcode_diff) {
161  if ((0x01 & outcode_diff) != 0)
162  if (check_point(p1, p2, (.5 - p1.x) / (p2.x - p1.x), 0x3e) == INSIDE)
163  return (INSIDE);
164  if ((0x02 & outcode_diff) != 0)
165  if (check_point(p1, p2, (-.5 - p1.x) / (p2.x - p1.x), 0x3d) == INSIDE)
166  return (INSIDE);
167  if ((0x04 & outcode_diff) != 0)
168  if (check_point(p1, p2, (.5 - p1.y) / (p2.y - p1.y), 0x3b) == INSIDE)
169  return (INSIDE);
170  if ((0x08 & outcode_diff) != 0)
171  if (check_point(p1, p2, (-.5 - p1.y) / (p2.y - p1.y), 0x37) == INSIDE)
172  return (INSIDE);
173  if ((0x10 & outcode_diff) != 0)
174  if (check_point(p1, p2, (.5 - p1.z) / (p2.z - p1.z), 0x2f) == INSIDE)
175  return (INSIDE);
176  if ((0x20 & outcode_diff) != 0)
177  if (check_point(p1, p2, (-.5 - p1.z) / (p2.z - p1.z), 0x1f) == INSIDE)
178  return (INSIDE);
179  return (OUTSIDE);
180 }
181 
182 /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
183 
184 /* Test if 3D point is inside 3D triangle */
185 
186 long point_triangle_intersection(Point3 p, Triangle3 t) {
187  long sign12, sign23, sign31;
188  Point3 vect12, vect23, vect31, vect1h, vect2h, vect3h;
189  Point3 cross12_1p, cross23_2p, cross31_3p;
190 
191  /* First, a quick bounding-box test: */
192  /* If P is outside triangle bbox, there cannot be an intersection. */
193 
194  if (p.x > MAX3(t.v1.x, t.v2.x, t.v3.x))
195  return (OUTSIDE);
196  if (p.y > MAX3(t.v1.y, t.v2.y, t.v3.y))
197  return (OUTSIDE);
198  if (p.z > MAX3(t.v1.z, t.v2.z, t.v3.z))
199  return (OUTSIDE);
200  if (p.x < MIN3(t.v1.x, t.v2.x, t.v3.x))
201  return (OUTSIDE);
202  if (p.y < MIN3(t.v1.y, t.v2.y, t.v3.y))
203  return (OUTSIDE);
204  if (p.z < MIN3(t.v1.z, t.v2.z, t.v3.z))
205  return (OUTSIDE);
206 
207  /* For each triangle side, make a vector out of it by subtracting vertexes;
208  */
209  /* make another vector from one vertex to point P. */
210  /* The crossproduct of these two vectors is orthogonal to both and the */
211  /* signs of its X,Y,Z components indicate whether P was to the inside or */
212  /* to the outside of this triangle side. */
213 
214  SUB(t.v1, t.v2, vect12)
215  SUB(t.v1, p, vect1h);
216  CROSS(vect12, vect1h, cross12_1p)
217  sign12 = SIGN3(
218  cross12_1p); /* Extract X,Y,Z signs as 0..7 or 0...63 integer */
219 
220  SUB(t.v2, t.v3, vect23)
221  SUB(t.v2, p, vect2h);
222  CROSS(vect23, vect2h, cross23_2p)
223  sign23 = SIGN3(cross23_2p);
224 
225  SUB(t.v3, t.v1, vect31)
226  SUB(t.v3, p, vect3h);
227  CROSS(vect31, vect3h, cross31_3p)
228  sign31 = SIGN3(cross31_3p);
229 
230 /* If all three crossproduct vectors agree in their component signs, */
231 /* then the point must be inside all three. */
232 /* P cannot be OUTSIDE all three sides simultaneously. */
233 
234 /* this is the old test; with the revised SIGN3() macro, the test
235  * needs to be revised. */
236 #ifdef OLD_TEST
237  if ((sign12 == sign23) && (sign23 == sign31))
238  return (INSIDE);
239  else
240  return (OUTSIDE);
241 #else
242  return ((sign12 & sign23 & sign31) == 0) ? OUTSIDE : INSIDE;
243 #endif
244 }
245 
246 /*. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . */
247 
248 /**********************************************/
249 /* This is the main algorithm procedure. */
250 /* Triangle t is compared with a unit cube, */
251 /* centered on the origin. */
252 /* It returns INSIDE (0) or OUTSIDE(1) if t */
253 /* intersects or does not intersect the cube. */
254 /**********************************************/
255 
256 long t_c_intersection(Triangle3 t) {
257  long v1_test, v2_test, v3_test;
258  float d;
259  float denom;
260  Point3 vect12, vect13, norm;
261  Point3 hitpp, hitpn, hitnp, hitnn;
262 
263  /* First compare all three vertexes with all six face-planes */
264  /* If any vertex is inside the cube, return immediately! */
265 
266  if ((v1_test = face_plane(t.v1)) == INSIDE)
267  return (INSIDE);
268  if ((v2_test = face_plane(t.v2)) == INSIDE)
269  return (INSIDE);
270  if ((v3_test = face_plane(t.v3)) == INSIDE)
271  return (INSIDE);
272 
273  /* If all three vertexes were outside of one or more face-planes, */
274  /* return immediately with a trivial rejection! */
275 
276  if ((v1_test & v2_test & v3_test) != 0)
277  return (OUTSIDE);
278 
279  /* Now do the same trivial rejection test for the 12 edge planes */
280 
281  v1_test |= bevel_2d(t.v1) << 8;
282  v2_test |= bevel_2d(t.v2) << 8;
283  v3_test |= bevel_2d(t.v3) << 8;
284  if ((v1_test & v2_test & v3_test) != 0)
285  return (OUTSIDE);
286 
287  /* Now do the same trivial rejection test for the 8 corner planes */
288 
289  v1_test |= bevel_3d(t.v1) << 24;
290  v2_test |= bevel_3d(t.v2) << 24;
291  v3_test |= bevel_3d(t.v3) << 24;
292  if ((v1_test & v2_test & v3_test) != 0)
293  return (OUTSIDE);
294 
295  /* If vertex 1 and 2, as a pair, cannot be trivially rejected */
296  /* by the above tests, then see if the v1-->v2 triangle edge */
297  /* intersects the cube. Do the same for v1-->v3 and v2-->v3. */
298  /* Pass to the intersection algorithm the "OR" of the outcode */
299  /* bits, so that only those cube faces which are spanned by */
300  /* each triangle edge need be tested. */
301 
302  if ((v1_test & v2_test) == 0)
303  if (check_line(t.v1, t.v2, v1_test | v2_test) == INSIDE)
304  return (INSIDE);
305  if ((v1_test & v3_test) == 0)
306  if (check_line(t.v1, t.v3, v1_test | v3_test) == INSIDE)
307  return (INSIDE);
308  if ((v2_test & v3_test) == 0)
309  if (check_line(t.v2, t.v3, v2_test | v3_test) == INSIDE)
310  return (INSIDE);
311 
312  /* By now, we know that the triangle is not off to any side, */
313  /* and that its sides do not penetrate the cube. We must now */
314  /* test for the cube intersecting the interior of the triangle. */
315  /* We do this by looking for intersections between the cube */
316  /* diagonals and the triangle...first finding the intersection */
317  /* of the four diagonals with the plane of the triangle, and */
318  /* then if that intersection is inside the cube, pursuing */
319  /* whether the intersection point is inside the triangle itself. */
320 
321  /* To find plane of the triangle, first perform crossproduct on */
322  /* two triangle side vectors to compute the normal vector. */
323 
324  SUB(t.v1, t.v2, vect12);
325  SUB(t.v1, t.v3, vect13);
326  CROSS(vect12, vect13, norm)
327 
328  /* The normal vector "norm" X,Y,Z components are the coefficients */
329  /* of the triangles AX + BY + CZ + D = 0 plane equation. If we */
330  /* solve the plane equation for X=Y=Z (a diagonal), we get */
331  /* -D/(A+B+C) as a metric of the distance from cube center to the */
332  /* diagonal/plane intersection. If this is between -0.5 and 0.5, */
333  /* the intersection is inside the cube. If so, we continue by */
334  /* doing a point/triangle intersection. */
335  /* Do this for all four diagonals. */
336 
337  d = norm.x * t.v1.x + norm.y * t.v1.y + norm.z * t.v1.z;
338 
339  /* if one of the diagonals is parallel to the plane, the other will
340  * intersect the plane */
341  if (fabs(denom = (norm.x + norm.y + norm.z)) > EPS)
342  /* skip parallel diagonals to the plane; division by 0 can occur */
343  {
344  hitpp.x = hitpp.y = hitpp.z = d / denom;
345  if (fabs(hitpp.x) <= 0.5)
346  if (point_triangle_intersection(hitpp, t) == INSIDE)
347  return (INSIDE);
348  }
349  if (fabs(denom = (norm.x + norm.y - norm.z)) > EPS) {
350  hitpn.z = -(hitpn.x = hitpn.y = d / denom);
351  if (fabs(hitpn.x) <= 0.5)
352  if (point_triangle_intersection(hitpn, t) == INSIDE)
353  return (INSIDE);
354  }
355  if (fabs(denom = (norm.x - norm.y + norm.z)) > EPS) {
356  hitnp.y = -(hitnp.x = hitnp.z = d / denom);
357  if (fabs(hitnp.x) <= 0.5)
358  if (point_triangle_intersection(hitnp, t) == INSIDE)
359  return (INSIDE);
360  }
361  if (fabs(denom = (norm.x - norm.y - norm.z)) > EPS) {
362  hitnn.y = hitnn.z = -(hitnn.x = d / denom);
363  if (fabs(hitnn.x) <= 0.5)
364  if (point_triangle_intersection(hitnn, t) == INSIDE)
365  return (INSIDE);
366  }
367 
368  /* No edge touched the cube; no cube diagonal touched the triangle. */
369  /* We're done...there was no intersection. */
370 
371  return (OUTSIDE);
372 }
Point3
Definition: tri_cube_c.h:41
Triangle3
Definition: tri_cube_c.h:47