Intersections.h
1 
2 #pragma once
3 
4 #define VF_EPSILON 0.000001
5 #define VF_CROSS(dest,v1,v2) dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
6 #define VF_DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
7 #define VF_SUB(dest,v1,v2) dest[0]=v1[0]-v2[0]; dest[1]=v1[1]-v2[1]; dest[2]=v1[2]-v2[2];
8 
9 namespace VoxelFarm
10 {
11  bool intersectPointPlane(
12  const double orig[3], const double dir[3],
13  const double vert0[3], const double vert1[3], const double vert2[3],
14  double *x, double *y, double *z);
15 
16  int planeBoxOverlap(const double normal[3], double d, const double maxbox[3]);
17 
18  int intersectTriangleAABB(const double boxcenter[3], const double boxhalfsize[3], const double triverts[3][3]);
19 
20  int intersectTriangle(
21  const double orig[3], const double dir[3],
22  const double vert0[3], const double vert1[3], const double vert2[3],
23  double *t, double *u, double *v);
24 
25  int intersectPointTriangle(
26  const double orig[3], const double dir[3],
27  const double vert0[3], const double vert1[3], const double vert2[3],
28  double *x, double *y, double *z,
29  double *t, double *u, double *v);
30 }
31 
33 // inline function implementations
35 
36 #define VF_X 0
37 #define VF_Y 1
38 #define VF_Z 2
39 
40 #define VF_FINDMINMAX(x0,x1,x2,min,max) \
41  min = max = x0; \
42  if (x1<min) min = x1; \
43  if (x1>max) max = x1; \
44  if (x2<min) min = x2; \
45  if (x2>max) max = x2;
46 
47 
48 /*======================== X-tests ========================*/
49 #define VF_AXISTEST_X01(a, b, fa, fb) \
50  p0 = a*v0[VF_Y] - b*v0[VF_Z]; \
51  p2 = a*v2[VF_Y] - b*v2[VF_Z]; \
52  if (p0<p2) { min = p0; max = p2; } \
53  else { min = p2; max = p0; } \
54  rad = fa * boxhalfsize[VF_Y] + fb * boxhalfsize[VF_Z]; \
55  if (min>rad || max < -rad) return 0;
56 
57 #define VF_AXISTEST_X2(a, b, fa, fb) \
58  p0 = a*v0[VF_Y] - b*v0[VF_Z]; \
59  p1 = a*v1[VF_Y] - b*v1[VF_Z]; \
60  if (p0<p1) { min = p0; max = p1; } \
61  else { min = p1; max = p0; } \
62  rad = fa * boxhalfsize[VF_Y] + fb * boxhalfsize[VF_Z]; \
63  if (min>rad || max < -rad) return 0;
64 
65 /*======================== Y-tests ========================*/
66 #define VF_AXISTEST_Y02(a, b, fa, fb) \
67  p0 = -a*v0[VF_X] + b*v0[VF_Z]; \
68  p2 = -a*v2[VF_X] + b*v2[VF_Z]; \
69  if (p0<p2) { min = p0; max = p2; } \
70  else { min = p2; max = p0; } \
71  rad = fa * boxhalfsize[VF_X] + fb * boxhalfsize[VF_Z]; \
72  if (min>rad || max < -rad) return 0;
73 
74 #define VF_AXISTEST_Y1(a, b, fa, fb) \
75  p0 = -a*v0[VF_X] + b*v0[VF_Z]; \
76  p1 = -a*v1[VF_X] + b*v1[VF_Z]; \
77  if (p0<p1) { min = p0; max = p1; } \
78  else { min = p1; max = p0; } \
79  rad = fa * boxhalfsize[VF_X] + fb * boxhalfsize[VF_Z]; \
80  if (min>rad || max < -rad) return 0;
81 
82 /*======================== Z-tests ========================*/
83 
84 #define VF_AXISTEST_Z12(a, b, fa, fb) \
85  p1 = a*v1[VF_X] - b*v1[VF_Y]; \
86  p2 = a*v2[VF_X] - b*v2[VF_Y]; \
87  if (p2<p1) { min = p2; max = p1; } \
88  else { min = p1; max = p2; } \
89  rad = fa * boxhalfsize[VF_X] + fb * boxhalfsize[VF_Y]; \
90  if (min>rad || max < -rad) return 0;
91 
92 #define VF_AXISTEST_Z0(a, b, fa, fb) \
93  p0 = a*v0[VF_X] - b*v0[VF_Y]; \
94  p1 = a*v1[VF_X] - b*v1[VF_Y]; \
95  if (p0<p1) { min = p0; max = p1; } \
96  else { min = p1; max = p0; } \
97  rad = fa * boxhalfsize[VF_X] + fb * boxhalfsize[VF_Y]; \
98  if (min>rad || max < -rad) return 0;
99 
100 
101 inline bool VoxelFarm::intersectPointPlane(
102  const double orig[3], const double dir[3],
103  const double vert0[3], const double vert1[3], const double vert2[3],
104  double *x, double *y, double *z)
105 {
106  double u[3], v[3], n[3], w0[3];
107  double r, a, b;
108 
109  // get triangle edge vectors and plane normal
110  VF_SUB(u, vert1, vert0);
111  VF_SUB(v, vert2, vert0);
112  VF_CROSS(n, u, v);
113 
114  VF_SUB(w0, orig, vert0);
115  b = VF_DOT(n, dir);
116 
117  if (b > -VF_EPSILON && b < VF_EPSILON)
118  {
119  return false;
120  }
121 
122  a = -VF_DOT(n, w0);
123  r = a / b;
124  if (r < 0.0)
125  {
126  return false;
127  }
128 
129  *x = orig[0] + r*dir[0];
130  *y = orig[1] + r*dir[1];
131  *z = orig[2] + r*dir[2];
132 
133  return true;
134 }
135 
136 inline int VoxelFarm::planeBoxOverlap(const double normal[3], double d, const double maxbox[3])
137 {
138  int q;
139  double vmin[3], vmax[3];
140  for (q = VF_X; q <= VF_Z; q++)
141  {
142  if (normal[q] > 0.0f)
143  {
144  vmin[q] = -maxbox[q];
145  vmax[q] = maxbox[q];
146  }
147  else
148  {
149  vmin[q] = maxbox[q];
150  vmax[q] = -maxbox[q];
151  }
152  }
153  if (VF_DOT(normal, vmin) + d > 0.0f)
154  {
155  return 0;
156  }
157  if (VF_DOT(normal, vmax) + d >= 0.0f)
158  {
159  return 1;
160  }
161 
162  return 0;
163 }
164 
165 inline int VoxelFarm::intersectTriangleAABB(const double boxcenter[3], const double boxhalfsize[3], const double triverts[3][3])
166 {
167  double v0[3], v1[3], v2[3];
168  double min, max, d, p0, p1, p2, rad, fex, fey, fez;
169  double normal[3], e0[3], e1[3], e2[3];
170 
171  VF_SUB(v0, triverts[0], boxcenter);
172  VF_SUB(v1, triverts[1], boxcenter);
173  VF_SUB(v2, triverts[2], boxcenter);
174 
175  VF_SUB(e0, v1, v0);
176  VF_SUB(e1, v2, v1);
177  VF_SUB(e2, v0, v2);
178 
179  fex = fabs(e0[VF_X]);
180  fey = fabs(e0[VF_Y]);
181  fez = fabs(e0[VF_Z]);
182  VF_AXISTEST_X01(e0[VF_Z], e0[VF_Y], fez, fey);
183  VF_AXISTEST_Y02(e0[VF_Z], e0[VF_X], fez, fex);
184  VF_AXISTEST_Z12(e0[VF_Y], e0[VF_X], fey, fex);
185 
186  fex = fabs(e1[VF_X]);
187  fey = fabs(e1[VF_Y]);
188  fez = fabs(e1[VF_Z]);
189  VF_AXISTEST_X01(e1[VF_Z], e1[VF_Y], fez, fey);
190  VF_AXISTEST_Y02(e1[VF_Z], e1[VF_X], fez, fex);
191  VF_AXISTEST_Z0(e1[VF_Y], e1[VF_X], fey, fex);
192 
193  fex = fabs(e2[VF_X]);
194  fey = fabs(e2[VF_Y]);
195  fez = fabs(e2[VF_Z]);
196  VF_AXISTEST_X2(e2[VF_Z], e2[VF_Y], fez, fey);
197  VF_AXISTEST_Y1(e2[VF_Z], e2[VF_X], fez, fex);
198  VF_AXISTEST_Z12(e2[VF_Y], e2[VF_X], fey, fex);
199 
200  VF_FINDMINMAX(v0[VF_X], v1[VF_X], v2[VF_X], min, max);
201  if (min > boxhalfsize[VF_X] || max<-boxhalfsize[VF_X])
202  {
203  return 0;
204  }
205 
206  VF_FINDMINMAX(v0[VF_Y], v1[VF_Y], v2[VF_Y], min, max);
207  if (min>boxhalfsize[VF_Y] || max<-boxhalfsize[VF_Y])
208  {
209  return 0;
210  }
211 
212  VF_FINDMINMAX(v0[VF_Z], v1[VF_Z], v2[VF_Z], min, max);
213  if (min>boxhalfsize[VF_Z] || max < -boxhalfsize[VF_Z])
214  {
215  return 0;
216  }
217 
218  VF_CROSS(normal, e0, e1);
219  d = -VF_DOT(normal, v0);
220  if (!planeBoxOverlap(normal, d, boxhalfsize))
221  {
222  return 0;
223  }
224 
225  return 1;
226 }
227 
228 inline int VoxelFarm::intersectTriangle(
229  const double orig[3], const double dir[3],
230  const double vert0[3], const double vert1[3], const double vert2[3],
231  double *t, double *u, double *v)
232 {
233  double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
234  double det, inv_det;
235  VF_SUB(edge1, vert1, vert0);
236  VF_SUB(edge2, vert2, vert0);
237  VF_CROSS(pvec, dir, edge2);
238  det = VF_DOT(edge1, pvec);
239  if (det > -VF_EPSILON && det < VF_EPSILON)
240  {
241  return 0;
242  }
243  inv_det = 1.0 / det;
244  VF_SUB(tvec, orig, vert0);
245  *u = VF_DOT(tvec, pvec) * inv_det;
246  if (*u < 0.0 || *u > 1.0)
247  {
248  return 0;
249  }
250  VF_CROSS(qvec, tvec, edge1);
251  *v = VF_DOT(dir, qvec) * inv_det;
252  if (*v < 0.0 || *u + *v > 1.0)
253  {
254  return 0;
255  }
256  *t = VF_DOT(edge2, qvec) * inv_det;
257  return 1;
258 }
259 
260 inline int VoxelFarm::intersectPointTriangle(
261  const double orig[3], const double dir[3],
262  const double vert0[3], const double vert1[3], const double vert2[3],
263  double *x, double *y, double *z,
264  double *t, double *u, double *v)
265 {
266  double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
267  double det, inv_det;
268  VF_SUB(edge1, vert1, vert0);
269  VF_SUB(edge2, vert2, vert0);
270  VF_CROSS(pvec, dir, edge2);
271  det = VF_DOT(edge1, pvec);
272  if (det > -VF_EPSILON && det < VF_EPSILON)
273  {
274  return 0;
275  }
276  inv_det = 1.0 / det;
277  VF_SUB(tvec, orig, vert0);
278  *u = VF_DOT(tvec, pvec) * inv_det;
279  if (*u < 0.0 || *u > 1.0)
280  {
281  return 0;
282  }
283  VF_CROSS(qvec, tvec, edge1);
284  *v = VF_DOT(dir, qvec) * inv_det;
285  if (*v < 0.0 || *u + *v > 1.0)
286  {
287  return 0;
288  }
289  *t = VF_DOT(edge2, qvec) * inv_det;
290 
291  //intersection point
292  *x = orig[0] + dir[0] * (*t);
293  *y = orig[1] + dir[1] * (*t);
294  *z = orig[2] + dir[2] * (*t);
295 
296  return 1;
297 }
298 
299 #undef VF_FINDMINMAX
300 #undef VF_X
301 #undef VF_Y
302 #undef VF_Z
303 
304 #undef VF_AXISTEST_X01
305 #undef VF_AXISTEST_X2
306 #undef VF_AXISTEST_Y02
307 #undef VF_AXISTEST_Y1
308 #undef VF_AXISTEST_Z12
309 #undef VF_AXISTEST_Z0
Contains all classes and functions for the VoxelFarm engine.