Blender  V3.3
GeomUtils.cpp
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
8 #include "GeomUtils.h"
9 
11 
12 // This internal procedure is defined below.
13 bool intersect2dSegPoly(Vec2r *seg, Vec2r *poly, unsigned n);
14 
15 bool intersect2dSeg2dArea(const Vec2r &min, const Vec2r &max, const Vec2r &A, const Vec2r &B)
16 {
17  Vec2r seg[2];
18  seg[0] = A;
19  seg[1] = B;
20 
21  Vec2r poly[5];
22  poly[0][0] = min[0];
23  poly[0][1] = min[1];
24  poly[1][0] = max[0];
25  poly[1][1] = min[1];
26  poly[2][0] = max[0];
27  poly[2][1] = max[1];
28  poly[3][0] = min[0];
29  poly[3][1] = max[1];
30  poly[4][0] = min[0];
31  poly[4][1] = min[1];
32 
33  return intersect2dSegPoly(seg, poly, 4);
34 }
35 
36 bool include2dSeg2dArea(const Vec2r &min, const Vec2r &max, const Vec2r &A, const Vec2r &B)
37 {
38  if ((((max[0] > A[0]) && (A[0] > min[0])) && ((max[0] > B[0]) && (B[0] > min[0]))) &&
39  (((max[1] > A[1]) && (A[1] > min[1])) && ((max[1] > B[1]) && (B[1] > min[1])))) {
40  return true;
41  }
42  return false;
43 }
44 
46  const Vec2r &p1, const Vec2r &p2, const Vec2r &p3, const Vec2r &p4, Vec2r &res)
47 {
48  real a1, a2, b1, b2, c1, c2; // Coefficients of line eqns
49  real r1, r2, r3, r4; // 'Sign' values
50  real denom, num; // Intermediate values
51 
52  // Compute a1, b1, c1, where line joining points p1 and p2 is "a1 x + b1 y + c1 = 0".
53  a1 = p2[1] - p1[1];
54  b1 = p1[0] - p2[0];
55  c1 = p2[0] * p1[1] - p1[0] * p2[1];
56 
57  // Compute r3 and r4.
58  r3 = a1 * p3[0] + b1 * p3[1] + c1;
59  r4 = a1 * p4[0] + b1 * p4[1] + c1;
60 
61  // Check signs of r3 and r4. If both point 3 and point 4 lie on same side of line 1,
62  // the line segments do not intersect.
63  if (r3 != 0 && r4 != 0 && r3 * r4 > 0.0) {
64  return DONT_INTERSECT;
65  }
66 
67  // Compute a2, b2, c2
68  a2 = p4[1] - p3[1];
69  b2 = p3[0] - p4[0];
70  c2 = p4[0] * p3[1] - p3[0] * p4[1];
71 
72  // Compute r1 and r2
73  r1 = a2 * p1[0] + b2 * p1[1] + c2;
74  r2 = a2 * p2[0] + b2 * p2[1] + c2;
75 
76  // Check signs of r1 and r2. If both point 1 and point 2 lie on same side of second line
77  // segment, the line segments do not intersect.
78  if (r1 != 0 && r2 != 0 && r1 * r2 > 0.0) {
79  return DONT_INTERSECT;
80  }
81 
82  // Line segments intersect: compute intersection point.
83  denom = a1 * b2 - a2 * b1;
84  if (fabs(denom) < M_EPSILON) {
85  return COLINEAR;
86  }
87 
88  num = b1 * c2 - b2 * c1;
89  res[0] = num / denom;
90 
91  num = a2 * c1 - a1 * c2;
92  res[1] = num / denom;
93 
94  return DO_INTERSECT;
95 }
96 
98  const Vec2r &p1, const Vec2r &p2, const Vec2r &p3, const Vec2r &p4, Vec2r &res)
99 {
100  real a1, a2, b1, b2, c1, c2; // Coefficients of line eqns
101  real denom, num; // Intermediate values
102 
103  // Compute a1, b1, c1, where line joining points p1 and p2 is "a1 x + b1 y + c1 = 0".
104  a1 = p2[1] - p1[1];
105  b1 = p1[0] - p2[0];
106  c1 = p2[0] * p1[1] - p1[0] * p2[1];
107 
108  // Compute a2, b2, c2
109  a2 = p4[1] - p3[1];
110  b2 = p3[0] - p4[0];
111  c2 = p4[0] * p3[1] - p3[0] * p4[1];
112 
113  // Line segments intersect: compute intersection point.
114  denom = a1 * b2 - a2 * b1;
115  if (fabs(denom) < M_EPSILON) {
116  return COLINEAR;
117  }
118 
119  num = b1 * c2 - b2 * c1;
120  res[0] = num / denom;
121 
122  num = a2 * c1 - a1 * c2;
123  res[1] = num / denom;
124 
125  return DO_INTERSECT;
126 }
127 
129  const Vec2r &p2,
130  const Vec2r &p3,
131  const Vec2r &p4,
132  real &t,
133  real &u,
134  real epsilon)
135 {
136  real a1, a2, b1, b2, c1, c2; // Coefficients of line eqns
137  real r1, r2, r3, r4; // 'Sign' values
138  real denom, num; // Intermediate values
139 
140  // Compute a1, b1, c1, where line joining points p1 and p2 is "a1 x + b1 y + c1 = 0".
141  a1 = p2[1] - p1[1];
142  b1 = p1[0] - p2[0];
143  c1 = p2[0] * p1[1] - p1[0] * p2[1];
144 
145  // Compute r3 and r4.
146  r3 = a1 * p3[0] + b1 * p3[1] + c1;
147  r4 = a1 * p4[0] + b1 * p4[1] + c1;
148 
149  // Check signs of r3 and r4. If both point 3 and point 4 lie on same side of line 1,
150  // the line segments do not intersect.
151  if (r3 != 0 && r4 != 0 && r3 * r4 > 0.0) {
152  return DONT_INTERSECT;
153  }
154 
155  // Compute a2, b2, c2
156  a2 = p4[1] - p3[1];
157  b2 = p3[0] - p4[0];
158  c2 = p4[0] * p3[1] - p3[0] * p4[1];
159 
160  // Compute r1 and r2
161  r1 = a2 * p1[0] + b2 * p1[1] + c2;
162  r2 = a2 * p2[0] + b2 * p2[1] + c2;
163 
164  // Check signs of r1 and r2. If both point 1 and point 2 lie on same side of second line
165  // segment, the line segments do not intersect.
166  if (r1 != 0 && r2 != 0 && r1 * r2 > 0.0) {
167  return DONT_INTERSECT;
168  }
169 
170  // Line segments intersect: compute intersection point.
171  denom = a1 * b2 - a2 * b1;
172  if (fabs(denom) < epsilon) {
173  return COLINEAR;
174  }
175 
176  real d1, e1;
177 
178  d1 = p1[1] - p3[1];
179  e1 = p1[0] - p3[0];
180 
181  num = -b2 * d1 - a2 * e1;
182  t = num / denom;
183 
184  num = -b1 * d1 - a1 * e1;
185  u = num / denom;
186 
187  return DO_INTERSECT;
188 }
189 
190 // AABB-triangle overlap test code by Tomas Akenine-Möller
191 // Function: int triBoxOverlap(real boxcenter[3], real boxhalfsize[3],real triverts[3][3]);
192 // History:
193 // 2001-03-05: released the code in its first version
194 // 2001-06-18: changed the order of the tests, faster
195 //
196 // Acknowledgement: Many thanks to Pierre Terdiman for suggestions and discussions on how to
197 // optimize code. Thanks to David Hunt for finding a ">="-bug!
198 
199 #define X 0
200 #define Y 1
201 #define Z 2
202 
203 #define FINDMINMAX(x0, x1, x2, min, max) \
204  { \
205  min = max = x0; \
206  if (x1 < min) { \
207  min = x1; \
208  } \
209  if (x1 > max) { \
210  max = x1; \
211  } \
212  if (x2 < min) { \
213  min = x2; \
214  } \
215  if (x2 > max) { \
216  max = x2; \
217  } \
218  } \
219  (void)0
220 
221 //======================== X-tests ========================//
222 #define AXISTEST_X01(a, b, fa, fb) \
223  { \
224  p0 = a * v0[Y] - b * v0[Z]; \
225  p2 = a * v2[Y] - b * v2[Z]; \
226  if (p0 < p2) { \
227  min = p0; \
228  max = p2; \
229  } \
230  else { \
231  min = p2; \
232  max = p0; \
233  } \
234  rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \
235  if (min > rad || max < -rad) { \
236  return 0; \
237  } \
238  } \
239  (void)0
240 
241 #define AXISTEST_X2(a, b, fa, fb) \
242  { \
243  p0 = a * v0[Y] - b * v0[Z]; \
244  p1 = a * v1[Y] - b * v1[Z]; \
245  if (p0 < p1) { \
246  min = p0; \
247  max = p1; \
248  } \
249  else { \
250  min = p1; \
251  max = p0; \
252  } \
253  rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z]; \
254  if (min > rad || max < -rad) { \
255  return 0; \
256  } \
257  } \
258  (void)0
259 
260 //======================== Y-tests ========================//
261 #define AXISTEST_Y02(a, b, fa, fb) \
262  { \
263  p0 = -a * v0[X] + b * v0[Z]; \
264  p2 = -a * v2[X] + b * v2[Z]; \
265  if (p0 < p2) { \
266  min = p0; \
267  max = p2; \
268  } \
269  else { \
270  min = p2; \
271  max = p0; \
272  } \
273  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \
274  if (min > rad || max < -rad) { \
275  return 0; \
276  } \
277  } \
278  (void)0
279 
280 #define AXISTEST_Y1(a, b, fa, fb) \
281  { \
282  p0 = -a * v0[X] + b * v0[Z]; \
283  p1 = -a * v1[X] + b * v1[Z]; \
284  if (p0 < p1) { \
285  min = p0; \
286  max = p1; \
287  } \
288  else { \
289  min = p1; \
290  max = p0; \
291  } \
292  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z]; \
293  if (min > rad || max < -rad) { \
294  return 0; \
295  } \
296  } \
297  (void)0
298 
299 //======================== Z-tests ========================//
300 #define AXISTEST_Z12(a, b, fa, fb) \
301  { \
302  p1 = a * v1[X] - b * v1[Y]; \
303  p2 = a * v2[X] - b * v2[Y]; \
304  if (p2 < p1) { \
305  min = p2; \
306  max = p1; \
307  } \
308  else { \
309  min = p1; \
310  max = p2; \
311  } \
312  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \
313  if (min > rad || max < -rad) { \
314  return 0; \
315  } \
316  } \
317  (void)0
318 
319 #define AXISTEST_Z0(a, b, fa, fb) \
320  { \
321  p0 = a * v0[X] - b * v0[Y]; \
322  p1 = a * v1[X] - b * v1[Y]; \
323  if (p0 < p1) { \
324  min = p0; \
325  max = p1; \
326  } \
327  else { \
328  min = p1; \
329  max = p0; \
330  } \
331  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y]; \
332  if (min > rad || max < -rad) { \
333  return 0; \
334  } \
335  } \
336  (void)0
337 
338 // This internal procedure is defined below.
339 bool overlapPlaneBox(Vec3r &normal, real d, Vec3r &maxbox);
340 
341 // Use separating axis theorem to test overlap between triangle and box need to test for overlap in
342 // these directions: 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle we
343 // do not even need to test these) 2) normal of the triangle 3) crossproduct(edge from tri,
344 // {x,y,z}-directin) this gives 3x3=9 more tests
345 bool overlapTriangleBox(Vec3r &boxcenter, Vec3r &boxhalfsize, Vec3r triverts[3])
346 {
347  Vec3r v0, v1, v2, normal, e0, e1, e2;
348  real min, max, d, p0, p1, p2, rad, fex, fey, fez;
349 
350  // This is the fastest branch on Sun
351  // move everything so that the boxcenter is in (0, 0, 0)
352  v0 = triverts[0] - boxcenter;
353  v1 = triverts[1] - boxcenter;
354  v2 = triverts[2] - boxcenter;
355 
356  // compute triangle edges
357  e0 = v1 - v0;
358  e1 = v2 - v1;
359  e2 = v0 - v2;
360 
361  // Bullet 3:
362  // Do the 9 tests first (this was faster)
363  fex = fabs(e0[X]);
364  fey = fabs(e0[Y]);
365  fez = fabs(e0[Z]);
366  AXISTEST_X01(e0[Z], e0[Y], fez, fey);
367  AXISTEST_Y02(e0[Z], e0[X], fez, fex);
368  AXISTEST_Z12(e0[Y], e0[X], fey, fex);
369 
370  fex = fabs(e1[X]);
371  fey = fabs(e1[Y]);
372  fez = fabs(e1[Z]);
373  AXISTEST_X01(e1[Z], e1[Y], fez, fey);
374  AXISTEST_Y02(e1[Z], e1[X], fez, fex);
375  AXISTEST_Z0(e1[Y], e1[X], fey, fex);
376 
377  fex = fabs(e2[X]);
378  fey = fabs(e2[Y]);
379  fez = fabs(e2[Z]);
380  AXISTEST_X2(e2[Z], e2[Y], fez, fey);
381  AXISTEST_Y1(e2[Z], e2[X], fez, fex);
382  AXISTEST_Z12(e2[Y], e2[X], fey, fex);
383 
384  // Bullet 1:
385  // first test overlap in the {x,y,z}-directions
386  // find min, max of the triangle each direction, and test for overlap in that direction -- this
387  // is equivalent to testing a minimal AABB around the triangle against the AABB
388 
389  // test in X-direction
390  FINDMINMAX(v0[X], v1[X], v2[X], min, max);
391  if (min > boxhalfsize[X] || max < -boxhalfsize[X]) {
392  return false;
393  }
394 
395  // test in Y-direction
396  FINDMINMAX(v0[Y], v1[Y], v2[Y], min, max);
397  if (min > boxhalfsize[Y] || max < -boxhalfsize[Y]) {
398  return false;
399  }
400 
401  // test in Z-direction
402  FINDMINMAX(v0[Z], v1[Z], v2[Z], min, max);
403  if (min > boxhalfsize[Z] || max < -boxhalfsize[Z]) {
404  return false;
405  }
406 
407  // Bullet 2:
408  // test if the box intersects the plane of the triangle
409  // compute plane equation of triangle: normal * x + d = 0
410  normal = e0 ^ e1;
411  d = -(normal * v0); // plane eq: normal.x + d = 0
412  if (!overlapPlaneBox(normal, d, boxhalfsize)) {
413  return false;
414  }
415 
416  return true; // box and triangle overlaps
417 }
418 
419 // Fast, Minimum Storage Ray-Triangle Intersection
420 //
421 // Tomas Möller
422 // Prosolvia Clarus AB
423 // Sweden
424 // <tompa@clarus.se>
425 //
426 // Ben Trumbore
427 // Cornell University
428 // Ithaca, New York
429 // <wbt@graphics.cornell.edu>
430 bool intersectRayTriangle(const Vec3r &orig,
431  const Vec3r &dir,
432  const Vec3r &v0,
433  const Vec3r &v1,
434  const Vec3r &v2,
435  real &t,
436  real &u,
437  real &v,
438  const real epsilon)
439 {
440  Vec3r edge1, edge2, tvec, pvec, qvec;
441  real det, inv_det;
442 
443  // find vectors for two edges sharing v0
444  edge1 = v1 - v0;
445  edge2 = v2 - v0;
446 
447  // begin calculating determinant - also used to calculate U parameter
448  pvec = dir ^ edge2;
449 
450  // if determinant is near zero, ray lies in plane of triangle
451  det = edge1 * pvec;
452 
453  // calculate distance from v0 to ray origin
454  tvec = orig - v0;
455  inv_det = 1.0 / det;
456 
457  qvec = tvec ^ edge1;
458 
459  if (det > epsilon) {
460  u = tvec * pvec;
461  if (u < 0.0 || u > det) {
462  return false;
463  }
464 
465  // calculate V parameter and test bounds
466  v = dir * qvec;
467  if (v < 0.0 || u + v > det) {
468  return false;
469  }
470  }
471  else if (det < -epsilon) {
472  // calculate U parameter and test bounds
473  u = tvec * pvec;
474  if (u > 0.0 || u < det) {
475  return false;
476  }
477 
478  // calculate V parameter and test bounds
479  v = dir * qvec;
480  if (v > 0.0 || u + v < det) {
481  return false;
482  }
483  }
484  else {
485  return false; // ray is parallel to the plane of the triangle
486  }
487 
488  u *= inv_det;
489  v *= inv_det;
490  t = (edge2 * qvec) * inv_det;
491 
492  return true;
493 }
494 
495 // Intersection between plane and ray, adapted from Graphics Gems, Didier Badouel
496 // The plane is represented by a set of points P implicitly defined as dot(norm, P) + d = 0.
497 // The ray is represented as r(t) = orig + dir * t.
499  const Vec3r &dir,
500  const Vec3r &norm,
501  const real d,
502  real &t,
503  const real epsilon)
504 {
505  real denom = norm * dir;
506 
507  if (fabs(denom) <= epsilon) { // plane and ray are parallel
508  if (fabs((norm * orig) + d) <= epsilon) {
509  return COINCIDENT; // plane and ray are coincident
510  }
511 
512  return COLINEAR;
513  }
514 
515  t = -(d + (norm * orig)) / denom;
516 
517  if (t < 0.0f) {
518  return DONT_INTERSECT;
519  }
520 
521  return DO_INTERSECT;
522 }
523 
524 bool intersectRayBBox(const Vec3r &orig,
525  const Vec3r &dir, // ray origin and direction
526  const Vec3r &boxMin,
527  const Vec3r &boxMax, // the bbox
528  real t0,
529  real t1,
530  real &tmin, // I0 = orig + tmin * dir is the first intersection
531  real &tmax, // I1 = orig + tmax * dir is the second intersection
532  real /*epsilon*/)
533 {
534  float tymin, tymax, tzmin, tzmax;
535  Vec3r inv_direction(1.0 / dir[0], 1.0 / dir[1], 1.0 / dir[2]);
536  int sign[3];
537  sign[0] = (inv_direction.x() < 0);
538  sign[1] = (inv_direction.y() < 0);
539  sign[2] = (inv_direction.z() < 0);
540 
541  Vec3r bounds[2];
542  bounds[0] = boxMin;
543  bounds[1] = boxMax;
544 
545  tmin = (bounds[sign[0]].x() - orig.x()) * inv_direction.x();
546  tmax = (bounds[1 - sign[0]].x() - orig.x()) * inv_direction.x();
547  tymin = (bounds[sign[1]].y() - orig.y()) * inv_direction.y();
548  tymax = (bounds[1 - sign[1]].y() - orig.y()) * inv_direction.y();
549  if ((tmin > tymax) || (tymin > tmax)) {
550  return false;
551  }
552  if (tymin > tmin) {
553  tmin = tymin;
554  }
555  if (tymax < tmax) {
556  tmax = tymax;
557  }
558  tzmin = (bounds[sign[2]].z() - orig.z()) * inv_direction.z();
559  tzmax = (bounds[1 - sign[2]].z() - orig.z()) * inv_direction.z();
560  if ((tmin > tzmax) || (tzmin > tmax)) {
561  return false;
562  }
563  if (tzmin > tmin) {
564  tmin = tzmin;
565  }
566  if (tzmax < tmax) {
567  tmax = tzmax;
568  }
569  return ((tmin < t1) && (tmax > t0));
570 }
571 
572 // Checks whether 3D points p lies inside or outside of the triangle ABC
573 bool includePointTriangle(const Vec3r &P, const Vec3r &A, const Vec3r &B, const Vec3r &C)
574 {
575  Vec3r AB(B - A);
576  Vec3r BC(C - B);
577  Vec3r CA(A - C);
578  Vec3r AP(P - A);
579  Vec3r BP(P - B);
580  Vec3r CP(P - C);
581 
582  Vec3r N(AB ^ BC); // triangle's normal
583 
584  N.normalize();
585 
586  Vec3r J(AB ^ AP), K(BC ^ BP), L(CA ^ CP);
587  J.normalize();
588  K.normalize();
589  L.normalize();
590 
591  if (J * N < 0) {
592  return false; // on the right of AB
593  }
594 
595  if (K * N < 0) {
596  return false; // on the right of BC
597  }
598 
599  if (L * N < 0) {
600  return false; // on the right of CA
601  }
602 
603  return true;
604 }
605 
606 void transformVertex(const Vec3r &vert, const Matrix44r &matrix, Vec3r &res)
607 {
608  HVec3r hvert(vert), res_tmp;
609  real scale;
610  for (unsigned int j = 0; j < 4; j++) {
611  scale = hvert[j];
612  for (unsigned int i = 0; i < 4; i++) {
613  res_tmp[i] += matrix(i, j) * scale;
614  }
615  }
616 
617  res[0] = res_tmp.x();
618  res[1] = res_tmp.y();
619  res[2] = res_tmp.z();
620 }
621 
622 void transformVertices(const vector<Vec3r> &vertices, const Matrix44r &trans, vector<Vec3r> &res)
623 {
624  size_t i;
625  res.resize(vertices.size());
626  for (i = 0; i < vertices.size(); i++) {
627  transformVertex(vertices[i], trans, res[i]);
628  }
629 }
630 
631 Vec3r rotateVector(const Matrix44r &mat, const Vec3r &v)
632 {
633  Vec3r res;
634  for (unsigned int i = 0; i < 3; i++) {
635  res[i] = 0;
636  for (unsigned int j = 0; j < 3; j++) {
637  res[i] += mat(i, j) * v[j];
638  }
639  }
640  res.normalize();
641  return res;
642 }
643 
644 // This internal procedure is defined below.
645 void fromCoordAToCoordB(const Vec3r &p, Vec3r &q, const real transform[4][4]);
646 
647 void fromWorldToCamera(const Vec3r &p, Vec3r &q, const real model_view_matrix[4][4])
648 {
649  fromCoordAToCoordB(p, q, model_view_matrix);
650 }
651 
652 void fromCameraToRetina(const Vec3r &p, Vec3r &q, const real projection_matrix[4][4])
653 {
654  fromCoordAToCoordB(p, q, projection_matrix);
655 }
656 
657 void fromRetinaToImage(const Vec3r &p, Vec3r &q, const int viewport[4])
658 {
659  // winX:
660  q[0] = viewport[0] + viewport[2] * (p[0] + 1.0) / 2.0;
661 
662  // winY:
663  q[1] = viewport[1] + viewport[3] * (p[1] + 1.0) / 2.0;
664 
665  // winZ:
666  q[2] = (p[2] + 1.0) / 2.0;
667 }
668 
669 void fromWorldToImage(const Vec3r &p,
670  Vec3r &q,
671  const real model_view_matrix[4][4],
672  const real projection_matrix[4][4],
673  const int viewport[4])
674 {
675  Vec3r p1, p2;
676  fromWorldToCamera(p, p1, model_view_matrix);
677  fromCameraToRetina(p1, p2, projection_matrix);
678  fromRetinaToImage(p2, q, viewport);
679  q[2] = p1[2];
680 }
681 
682 void fromWorldToImage(const Vec3r &p, Vec3r &q, const real transform[4][4], const int viewport[4])
683 {
685 
686  // winX:
687  q[0] = viewport[0] + viewport[2] * (q[0] + 1.0) / 2.0;
688 
689  // winY:
690  q[1] = viewport[1] + viewport[3] * (q[1] + 1.0) / 2.0;
691 }
692 
693 void fromImageToRetina(const Vec3r &p, Vec3r &q, const int viewport[4])
694 {
695  q = p;
696  q[0] = 2.0 * (q[0] - viewport[0]) / viewport[2] - 1.0;
697  q[1] = 2.0 * (q[1] - viewport[1]) / viewport[3] - 1.0;
698 }
699 
700 void fromRetinaToCamera(const Vec3r &p, Vec3r &q, real focal, const real projection_matrix[4][4])
701 {
702  if (projection_matrix[3][3] == 0.0) { // perspective
703  q[0] = (-p[0] * focal) / projection_matrix[0][0];
704  q[1] = (-p[1] * focal) / projection_matrix[1][1];
705  q[2] = focal;
706  }
707  else { // orthogonal
708  q[0] = p[0] / projection_matrix[0][0];
709  q[1] = p[1] / projection_matrix[1][1];
710  q[2] = focal;
711  }
712 }
713 
714 void fromCameraToWorld(const Vec3r &p, Vec3r &q, const real model_view_matrix[4][4])
715 {
716  real translation[3] = {
717  model_view_matrix[0][3],
718  model_view_matrix[1][3],
719  model_view_matrix[2][3],
720  };
721  for (unsigned short i = 0; i < 3; i++) {
722  q[i] = 0.0;
723  for (unsigned short j = 0; j < 3; j++) {
724  q[i] += model_view_matrix[j][i] * (p[j] - translation[j]);
725  }
726  }
727 }
728 
729 //
730 // Internal code
731 //
733 
734 // Copyright 2001, softSurfer (www.softsurfer.com)
735 // This code may be freely used and modified for any purpose providing that this copyright notice
736 // is included with it. SoftSurfer makes no warranty for this code, and cannot be held liable for
737 // any real or imagined damage resulting from its use. Users of this code must verify correctness
738 // for their application.
739 
740 #define PERP(u, v) ((u)[0] * (v)[1] - (u)[1] * (v)[0]) // 2D perp product
741 
742 inline bool intersect2dSegPoly(Vec2r *seg, Vec2r *poly, unsigned n)
743 {
744  if (seg[0] == seg[1]) {
745  return false;
746  }
747 
748  real tE = 0; // the maximum entering segment parameter
749  real tL = 1; // the minimum leaving segment parameter
750  real t, N, D; // intersect parameter t = N / D
751  Vec2r dseg = seg[1] - seg[0]; // the segment direction vector
752  Vec2r e; // edge vector
753 
754  for (unsigned int i = 0; i < n; i++) { // process polygon edge poly[i]poly[i+1]
755  e = poly[i + 1] - poly[i];
756  N = PERP(e, seg[0] - poly[i]);
757  D = -PERP(e, dseg);
758  if (fabs(D) < M_EPSILON) {
759  if (N < 0) {
760  return false;
761  }
762 
763  continue;
764  }
765 
766  t = N / D;
767  if (D < 0) { // segment seg is entering across this edge
768  if (t > tE) { // new max tE
769  tE = t;
770  if (tE > tL) { // seg enters after leaving polygon
771  return false;
772  }
773  }
774  }
775  else { // segment seg is leaving across this edge
776  if (t < tL) { // new min tL
777  tL = t;
778  if (tL < tE) { // seg leaves before entering polygon
779  return false;
780  }
781  }
782  }
783  }
784 
785  // tE <= tL implies that there is a valid intersection subsegment
786  return true;
787 }
788 
789 inline bool overlapPlaneBox(Vec3r &normal, real d, Vec3r &maxbox)
790 {
791  Vec3r vmin, vmax;
792 
793  for (unsigned int q = X; q <= Z; q++) {
794  if (normal[q] > 0.0f) {
795  vmin[q] = -maxbox[q];
796  vmax[q] = maxbox[q];
797  }
798  else {
799  vmin[q] = maxbox[q];
800  vmax[q] = -maxbox[q];
801  }
802  }
803  if ((normal * vmin) + d > 0.0f) {
804  return false;
805  }
806  if ((normal * vmax) + d >= 0.0f) {
807  return true;
808  }
809  return false;
810 }
811 
812 inline void fromCoordAToCoordB(const Vec3r &p, Vec3r &q, const real transform[4][4])
813 {
814  HVec3r hp(p);
815  HVec3r hq(0, 0, 0, 0);
816 
817  for (unsigned int i = 0; i < 4; i++) {
818  for (unsigned int j = 0; j < 4; j++) {
819  hq[i] += transform[i][j] * hp[j];
820  }
821  }
822 
823  if (hq[3] == 0) {
824  q = p;
825  return;
826  }
827 
828  for (unsigned int k = 0; k < 3; k++) {
829  q[k] = hq[k] / hq[3];
830  }
831 }
832 
833 } // namespace Freestyle::GeomUtils
#define D
#define K(key)
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble y2 _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat y2 _GL_VOID_RET _GL_VOID GLint GLint GLint y2 _GL_VOID_RET _GL_VOID GLshort GLshort GLshort y2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLuint *buffer _GL_VOID_RET _GL_VOID GLdouble t _GL_VOID_RET _GL_VOID GLfloat t _GL_VOID_RET _GL_VOID GLint t _GL_VOID_RET _GL_VOID GLshort t _GL_VOID_RET _GL_VOID GLdouble t
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble v1
#define X
Definition: GeomUtils.cpp:199
#define AXISTEST_Y02(a, b, fa, fb)
Definition: GeomUtils.cpp:261
#define FINDMINMAX(x0, x1, x2, min, max)
Definition: GeomUtils.cpp:203
#define AXISTEST_Z12(a, b, fa, fb)
Definition: GeomUtils.cpp:300
#define Z
Definition: GeomUtils.cpp:201
#define AXISTEST_Y1(a, b, fa, fb)
Definition: GeomUtils.cpp:280
#define Y
Definition: GeomUtils.cpp:200
#define AXISTEST_X2(a, b, fa, fb)
Definition: GeomUtils.cpp:241
#define AXISTEST_X01(a, b, fa, fb)
Definition: GeomUtils.cpp:222
#define PERP(u, v)
Definition: GeomUtils.cpp:740
#define AXISTEST_Z0(a, b, fa, fb)
Definition: GeomUtils.cpp:319
Various tools for geometry.
#define C
Definition: RandGen.cpp:25
ATTR_WARN_UNUSED_RESULT const BMVert * v2
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
ATTR_WARN_UNUSED_RESULT const BMVert * v
#define A
SIMD_FORCE_INLINE btVector3 transform(const btVector3 &point) const
static btDbvtVolume bounds(btDbvtNode **leaves, int count)
Definition: btDbvt.cpp:299
SIMD_FORCE_INLINE btScalar norm() const
Return the norm (length) of the vector.
Definition: btVector3.h:263
value_type z() const
Definition: VecMat.h:473
value_type x() const
Definition: VecMat.h:463
value_type y() const
Definition: VecMat.h:468
value_type x() const
Definition: VecMat.h:518
value_type z() const
Definition: VecMat.h:538
value_type y() const
Definition: VecMat.h:528
Vec< T, N > & normalize()
Definition: VecMat.h:105
IconTextureDrawCall normal
ccl_device_inline float2 fabs(const float2 &a)
Definition: math_float2.h:222
static float P(float k)
Definition: math_interp.c:25
#define N
#define B
#define L
void fromRetinaToImage(const Vec3r &p, Vec3r &q, const int viewport[4])
Definition: GeomUtils.cpp:657
void fromWorldToImage(const Vec3r &p, Vec3r &q, const real model_view_matrix[4][4], const real projection_matrix[4][4], const int viewport[4])
Definition: GeomUtils.cpp:669
void transformVertex(const Vec3r &vert, const Matrix44r &matrix, Vec3r &res)
Definition: GeomUtils.cpp:606
void fromWorldToCamera(const Vec3r &p, Vec3r &q, const real model_view_matrix[4][4])
Definition: GeomUtils.cpp:647
intersection_test intersect2dSeg2dSegParametric(const Vec2r &p1, const Vec2r &p2, const Vec2r &p3, const Vec2r &p4, real &t, real &u, real epsilon)
Definition: GeomUtils.cpp:128
bool intersect2dSeg2dArea(const Vec2r &min, const Vec2r &max, const Vec2r &A, const Vec2r &B)
Definition: GeomUtils.cpp:15
intersection_test intersect2dSeg2dSeg(const Vec2r &p1, const Vec2r &p2, const Vec2r &p3, const Vec2r &p4, Vec2r &res)
Definition: GeomUtils.cpp:45
bool intersect2dSegPoly(Vec2r *seg, Vec2r *poly, unsigned n)
Definition: GeomUtils.cpp:742
bool overlapTriangleBox(Vec3r &boxcenter, Vec3r &boxhalfsize, Vec3r triverts[3])
Definition: GeomUtils.cpp:345
bool overlapPlaneBox(Vec3r &normal, real d, Vec3r &maxbox)
Definition: GeomUtils.cpp:789
void fromCoordAToCoordB(const Vec3r &p, Vec3r &q, const real transform[4][4])
Definition: GeomUtils.cpp:812
void fromCameraToWorld(const Vec3r &p, Vec3r &q, const real model_view_matrix[4][4])
Definition: GeomUtils.cpp:714
bool intersectRayTriangle(const Vec3r &orig, const Vec3r &dir, const Vec3r &v0, const Vec3r &v1, const Vec3r &v2, real &t, real &u, real &v, const real epsilon)
Definition: GeomUtils.cpp:430
bool includePointTriangle(const Vec3r &P, const Vec3r &A, const Vec3r &B, const Vec3r &C)
Definition: GeomUtils.cpp:573
void fromCameraToRetina(const Vec3r &p, Vec3r &q, const real projection_matrix[4][4])
Definition: GeomUtils.cpp:652
void transformVertices(const vector< Vec3r > &vertices, const Matrix44r &trans, vector< Vec3r > &res)
Definition: GeomUtils.cpp:622
bool intersectRayBBox(const Vec3r &orig, const Vec3r &dir, const Vec3r &boxMin, const Vec3r &boxMax, real t0, real t1, real &tmin, real &tmax, real)
Definition: GeomUtils.cpp:524
bool include2dSeg2dArea(const Vec2r &min, const Vec2r &max, const Vec2r &A, const Vec2r &B)
Definition: GeomUtils.cpp:36
intersection_test intersectRayPlane(const Vec3r &orig, const Vec3r &dir, const Vec3r &norm, const real d, real &t, const real epsilon)
Definition: GeomUtils.cpp:498
void fromRetinaToCamera(const Vec3r &p, Vec3r &q, real focal, const real projection_matrix[4][4])
Definition: GeomUtils.cpp:700
void fromImageToRetina(const Vec3r &p, Vec3r &q, const int viewport[4])
Definition: GeomUtils.cpp:693
intersection_test intersect2dLine2dLine(const Vec2r &p1, const Vec2r &p2, const Vec2r &p3, const Vec2r &p4, Vec2r &res)
Definition: GeomUtils.cpp:97
Vec3r rotateVector(const Matrix44r &mat, const Vec3r &v)
Definition: GeomUtils.cpp:631
VecMat::Vec3< real > Vec3r
Definition: Geom.h:28
static const real M_EPSILON
Definition: Precision.h:15
double real
Definition: Precision.h:12
double sign(double arg)
Definition: utility.h:250
static double epsilon
#define min(a, b)
Definition: sort.c:35
float max