Developer Documentation
Algorithms.cc
1 /*===========================================================================*\
2  * *
3  * OpenFlipper *
4  * Copyright (c) 2001-2015, RWTH-Aachen University *
5  * Department of Computer Graphics and Multimedia *
6  * All rights reserved. *
7  * www.openflipper.org *
8  * *
9  *---------------------------------------------------------------------------*
10  * This file is part of OpenFlipper. *
11  *---------------------------------------------------------------------------*
12  * *
13  * Redistribution and use in source and binary forms, with or without *
14  * modification, are permitted provided that the following conditions *
15  * are met: *
16  * *
17  * 1. Redistributions of source code must retain the above copyright notice, *
18  * this list of conditions and the following disclaimer. *
19  * *
20  * 2. Redistributions in binary form must reproduce the above copyright *
21  * notice, this list of conditions and the following disclaimer in the *
22  * documentation and/or other materials provided with the distribution. *
23  * *
24  * 3. Neither the name of the copyright holder nor the names of its *
25  * contributors may be used to endorse or promote products derived from *
26  * this software without specific prior written permission. *
27  * *
28  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31  * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32  * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39  * *
40 \*===========================================================================*/
41 
42 
43 //=============================================================================
44 //
45 // CLASS Geometry - IMPLEMENTATION
46 //
47 //=============================================================================
48 
49 #define GEO_ALGORITHMS_C
50 
51 //== INCLUDES =================================================================
52 
53 #include "Algorithms.hh"
54 #include <ACG/Utils/NumLimitsT.hh>
55 #include <ACG/Utils/VSToolsT.hh>
56 #include <ACG/Math/GLMatrixT.hh>
57 
58 
59 #ifdef max
60 # undef max
61 #endif
62 
63 #ifdef min
64 # undef min
65 #endif
66 
67 
68 //----------------------------------------------------------------------------
69 
70 
71 namespace ACG {
72 namespace Geometry {
73 
74 
75 //== IMPLEMENTATION ==========================================================
76 
77 
78 //== 3D STUFF ================================================================
79 
80 
81 template<typename Scalar>
82 bool
84  const VectorT<Scalar,3> & _u,
85  const VectorT<Scalar,3> & _v,
86  const VectorT<Scalar,3> & _w,
87  VectorT<Scalar,3> & _result )
88 {
89  VectorT<Scalar,3> vu = _v - _u,
90  wu = _w - _u,
91  pu = _p - _u;
92 
93 
94  // find largest absolute coodinate of normal
95  Scalar nx = vu[1]*wu[2] - vu[2]*wu[1],
96  ny = vu[2]*wu[0] - vu[0]*wu[2],
97  nz = vu[0]*wu[1] - vu[1]*wu[0],
98  ax = fabs(nx),
99  ay = fabs(ny),
100  az = fabs(nz);
101 
102 
103  unsigned char max_coord;
104 
105  if ( ax > ay ) {
106  if ( ax > az ) {
107  max_coord = 0;
108  }
109  else {
110  max_coord = 2;
111  }
112  }
113  else {
114  if ( ay > az ) {
115  max_coord = 1;
116  }
117  else {
118  max_coord = 2;
119  }
120  }
121 
122 
123  // solve 2D problem
124  switch (max_coord) {
125 
126  case 0:
127  {
128  if (1.0+ax == 1.0) return false;
129  _result[1] = 1.0 + (pu[1]*wu[2]-pu[2]*wu[1])/nx - 1.0;
130  _result[2] = 1.0 + (vu[1]*pu[2]-vu[2]*pu[1])/nx - 1.0;
131  _result[0] = 1.0 - _result[1] - _result[2];
132  }
133  break;
134 
135  case 1:
136  {
137  if (1.0+ay == 1.0) return false;
138  _result[1] = 1.0 + (pu[2]*wu[0]-pu[0]*wu[2])/ny - 1.0;
139  _result[2] = 1.0 + (vu[2]*pu[0]-vu[0]*pu[2])/ny - 1.0;
140  _result[0] = 1.0 - _result[1] - _result[2];
141  }
142  break;
143 
144  case 2:
145  {
146  if (1.0+az == 1.0) return false;
147  _result[1] = 1.0 + (pu[0]*wu[1]-pu[1]*wu[0])/nz - 1.0;
148  _result[2] = 1.0 + (vu[0]*pu[1]-vu[1]*pu[0])/nz - 1.0;
149  _result[0] = 1.0 - _result[1] - _result[2];
150  }
151  break;
152  }
153 
154  return true;
155 }
156 
157 
158 //-----------------------------------------------------------------------------
159 
160 
161 template <class Vec>
162 typename Vec::value_type
163 triangleAreaSquared( const Vec& _v0,
164  const Vec& _v1,
165  const Vec& _v2 )
166 {
167  return 0.25 * ((_v1-_v0)%(_v2-_v0)).sqrnorm();
168 }
169 
170 //-----------------------------------------------------------------------------
171 template<typename Scalar>
172 bool edgeConvexPolygonIntersection(std::vector<VectorT<Scalar,3> > _polygon_points,
173  VectorT<Scalar,3> _v0,
174  VectorT<Scalar,3> _v1,
175  VectorT<Scalar,3> &_result)
176 {
177  if(_polygon_points.size() < 3)
178  {
179  return false;
180  }
181 
182  // compute center of gravity
183  VectorT<Scalar,3> cog(0.0);
184  for(size_t i = 0; i<_polygon_points.size(); i++)
185  {
186  cog += _polygon_points[i];
187  }
188  cog /= ((Scalar)_polygon_points.size());
189 
190  //get face normal
191  VectorT<Scalar,3> n = ( _polygon_points[0] - cog )%(_polygon_points[1] - cog );
192 
193  size_t c = 1;
194  while((std::fabs(n[0])<1e-30) & (std::fabs(n[1])<1e-30) & (std::fabs(n[2])<1e-30))
195  {
196  n = ( _polygon_points[c] - cog )%(_polygon_points[(c+1)%_polygon_points.size()] - cog );
197  c++;
198  if(c == _polygon_points.size()+1)
199  {
200  std::cerr << "warning: no valid normal found \n";
201  return false;
202  }
203  }
204  n = n.normalize();
205 
206  if(n.norm() <= 0.01)
207  {
208  std::cerr << "warning: no valid normal found normal"<< n<<" norm "<<n.norm()<< " \n";
209  }
210 
211  //compute rotation to xy plane
212  VectorT<Scalar,3> z(0.0, 0.0, 1.0);
213  VectorT<Scalar,3> axis(0.0, 0.0, 0.0);
214  Scalar angle = 0.0;
215  bool rotation = rotationOfTwoVectors(n, z, axis, angle, true);
216 
218  R.identity();
219 
220  //set to 0.0 if isnan in rotation function or if not set at all
221  if((angle != 0.0) && rotation)
222  {
223  R.rotate(angle, axis);
224  }
225 
226  //rotate all points to xy plane
227  std::vector<VectorT<Scalar,2> > rotated_points;
228  for(size_t i = 0; i<_polygon_points.size(); i++)
229  {
230  VectorT<Scalar,3> new_point_3 = _polygon_points[i] - cog;
231  VectorT<Scalar,4> new_point_4(new_point_3[0], new_point_3[1], new_point_3[2], 0);
232  new_point_4 = R*new_point_4;
233  rotated_points.push_back(VectorT<Scalar,2>(new_point_4[0], new_point_4[1]));
234  }
235 
236  //compute ray plane intersection
237  Scalar d = n|cog;
238  if((n|(_v1 - _v0)) == 0.0)
239  {
240  // if((n|_v0)-d <= 0.00000001)
241  // {
242  // std::cerr << __FUNCTION__ << " parallel to face "<< d<<"\n";
243  // _result = _v0;
244  // return true;
245  // }
246  return false;
247  }
248 
249  Scalar alpha = (d - (n|_v0))/(n|(_v1 - _v0));
250  _result = _v0 + alpha*(_v1 - _v0);
251 
252  //returns false if not on edge _v0, _v1
253  if(alpha > 1.0 || alpha < 0.0)
254  {
255  return false;
256  }
257 
258  VectorT<Scalar,4> rotated_result(_result[0] - cog[0], _result[1] - cog[1], _result[2] - cog[2], 0);
259  rotated_result = R*rotated_result;
260  //std::cerr <<" angle "<< angle <<" normal "<< n <<"result "<<_result<<" alpha "<<alpha<< " rot res "<< rotated_result<<"in plane: "<<((_result|n) - d)<<"\n";
261  VectorT<Scalar,2> intersect(rotated_result[0], rotated_result[1]);
262 
263  //compare normal direction to intersection
264  size_t points_count = rotated_points.size();
265  for(size_t i = 0; i<points_count; i++)
266  {
267  VectorT<Scalar,2> e = rotated_points[((i+1)%points_count)] - rotated_points[i];
268  VectorT<Scalar,2> n_e(e[1], -e[0]);
269  VectorT<Scalar,2> mid_e = 0.5 * (rotated_points[((i+1)%points_count)] + rotated_points[i]);
270  VectorT<Scalar,2> cmp0 = - mid_e;
271  VectorT<Scalar,2> cmp1 = intersect - mid_e;
272  int sgn0 = ((n_e|cmp0) < 0 )? -1 : ((n_e|cmp0) > 0);
273  int sgn1 = ((n_e|cmp1) < 0 )? -1 : ((n_e|cmp1) > 0);
274 
275  if(sgn1 != sgn0 && sgn1 != 0)
276  {
277  return false;
278  }
279  }
280  return true;
281 
282 }
283 
284 
285 //-----------------------------------------------------------------------------
286 
287 
288 template<class Vec>
289 typename Vec::value_type
290 distPointLineSquared( const Vec& _p,
291  const Vec& _v0,
292  const Vec& _v1,
293  Vec* _min_v )
294 {
295  Vec d1(_p-_v0), d2(_v1-_v0), min_v(_v0);
296  typename Vec::value_type t = (d1|d2)/ d2.sqrnorm();
297 
298  if (t > 1.0) d1 = _p - (min_v = _v1);
299  else if (t > 0.0) d1 = _p - (min_v = _v0 + d2*t);
300 
301  if (_min_v) *_min_v = min_v;
302  return d1.sqrnorm();
303 }
304 
305  //-----------------------------------------------------------------------------
306 
307 template <class Vec>
308 typename Vec::value_type
309 distPointTriangleSquared( const Vec& _p,
310  const Vec& _v0,
311  const Vec& _v1,
312  const Vec& _v2,
313  Vec& _nearestPoint,
314  bool _stable)
315 {
316  Vec v0v1 = _v1 - _v0;
317  Vec v0v2 = _v2 - _v0;
318  Vec n = v0v1 % v0v2; // not normalized !
319  typename Vec::value_type d = n.sqrnorm();
320 
321 
322  // Check if the triangle is degenerated
323  if (d < FLT_MIN && d > -FLT_MIN) {
324  if (_stable) {
325  const double l0 = v0v1.sqrnorm();
326  const double l1 = v0v2.sqrnorm();
327  const double l2 = (_v2 - _v1).sqrnorm();
328  if (l0 > l1 && l0 > l2) {
329  return distPointLineSquared(_p, _v0, _v1, &_nearestPoint);
330  }
331  else if (l1 > l0 && l1 > l2) {
332  return distPointLineSquared(_p, _v0, _v2, &_nearestPoint);
333  }
334  else {
335  return distPointLineSquared(_p, _v1, _v2, &_nearestPoint);
336  }
337  } else {
338  std::cerr << "distPointTriangleSquared: Degenerated triangle !\n";
339  return -1.0;
340  }
341  }
342  typename Vec::value_type invD = typename Vec::value_type(1.0) / d;
343 
344 
345  // these are not needed for every point, should still perform
346  // better with many points against one triangle
347  Vec v1v2 = _v2 - _v1;
348  typename Vec::value_type inv_v0v2_2 = typename Vec::value_type(1.0) / v0v2.sqrnorm();
349  typename Vec::value_type inv_v0v1_2 = typename Vec::value_type(1.0) / v0v1.sqrnorm();
350  typename Vec::value_type inv_v1v2_2 = typename Vec::value_type(1.0) / v1v2.sqrnorm();
351 
352 
353  Vec v0p = _p - _v0;
354  Vec t = v0p % n;
355  typename Vec::value_type s01, s02, s12;
356  typename Vec::value_type a = (t | v0v2) * -invD;
357  typename Vec::value_type b = (t | v0v1) * invD;
358 
359 
360  if (a < 0)
361  {
362  // Calculate the distance to an edge or a corner vertex
363  s02 = ( v0v2 | v0p ) * inv_v0v2_2;
364  if (s02 < 0.0)
365  {
366  s01 = ( v0v1 | v0p ) * inv_v0v1_2;
367  if (s01 <= 0.0) {
368  v0p = _v0;
369  } else if (s01 >= 1.0) {
370  v0p = _v1;
371  } else {
372  v0p = _v0 + v0v1 * s01;
373  }
374  } else if (s02 > 1.0) {
375  s12 = ( v1v2 | ( _p - _v1 )) * inv_v1v2_2;
376  if (s12 >= 1.0) {
377  v0p = _v2;
378  } else if (s12 <= 0.0) {
379  v0p = _v1;
380  } else {
381  v0p = _v1 + v1v2 * s12;
382  }
383  } else {
384  v0p = _v0 + v0v2 * s02;
385  }
386  } else if (b < 0.0) {
387  // Calculate the distance to an edge or a corner vertex
388  s01 = ( v0v1 | v0p ) * inv_v0v1_2;
389  if (s01 < 0.0)
390  {
391  s02 = ( v0v2 | v0p ) * inv_v0v2_2;
392  if (s02 <= 0.0) {
393  v0p = _v0;
394  } else if (s02 >= 1.0) {
395  v0p = _v2;
396  } else {
397  v0p = _v0 + v0v2 * s02;
398  }
399  } else if (s01 > 1.0) {
400  s12 = ( v1v2 | ( _p - _v1 )) * inv_v1v2_2;
401  if (s12 >= 1.0) {
402  v0p = _v2;
403  } else if (s12 <= 0.0) {
404  v0p = _v1;
405  } else {
406  v0p = _v1 + v1v2 * s12;
407  }
408  } else {
409  v0p = _v0 + v0v1 * s01;
410  }
411  } else if (a+b > 1.0) {
412  // Calculate the distance to an edge or a corner vertex
413  s12 = ( v1v2 | ( _p - _v1 )) * inv_v1v2_2;
414  if (s12 >= 1.0) {
415  s02 = ( v0v2 | v0p ) * inv_v0v2_2;
416  if (s02 <= 0.0) {
417  v0p = _v0;
418  } else if (s02 >= 1.0) {
419  v0p = _v2;
420  } else {
421  v0p = _v0 + v0v2*s02;
422  }
423  } else if (s12 <= 0.0) {
424  s01 = ( v0v1 | v0p ) * inv_v0v1_2;
425  if (s01 <= 0.0) {
426  v0p = _v0;
427  } else if (s01 >= 1.0) {
428  v0p = _v1;
429  } else {
430  v0p = _v0 + v0v1 * s01;
431  }
432  } else {
433  v0p = _v1 + v1v2 * s12;
434  }
435  } else {
436  // Calculate the distance to an interior point of the triangle
437  _nearestPoint = _p - n*((n|v0p) * invD);
438  return (_nearestPoint - _p).sqrnorm();
439  }
440 
441  _nearestPoint = v0p;
442 
443  return (_nearestPoint - _p).sqrnorm();
444 }
445 
446 //-----------------------------------------------------------------------------
447 
448 
449 template <class Vec>
450 typename Vec::value_type
451 distPointTriangleSquared( const Vec& _p,
452  const Vec& _v0,
453  const Vec& _v1,
454  const Vec& _v2,
455  Vec& _nearestPoint )
456 {
457  return distPointTriangleSquared(_p, _v0, _v1, _v2, _nearestPoint, false);
458 }
459 
460 
461 //-----------------------------------------------------------------------------
462 
463 
464 template <class Vec>
465 typename Vec::value_type
467  const Vec& _v0,
468  const Vec& _v1,
469  const Vec& _v2,
470  Vec& _nearestPoint )
471 {
472  return distPointTriangleSquared(_p, _v0, _v1, _v2, _nearestPoint, true);
473 }
474 
475 
476 //-----------------------------------------------------------------------------
477 
478 
479 
480 //
481 // Modified code of Dave Eberly (www.magic-software.com)
482 //
483 
484 template<typename Scalar>
485 Scalar
487  const VectorT<Scalar,3>& _v01,
488  const VectorT<Scalar,3>& _v10,
489  const VectorT<Scalar,3>& _v11,
490  VectorT<Scalar,3>* _min_v0,
491  VectorT<Scalar,3>* _min_v1,
492  bool _fastApprox )
493 {
494  VectorT<Scalar,3> kDiff = _v00 - _v10;
495  VectorT<Scalar,3> d0 = _v01-_v00;
496  VectorT<Scalar,3> d1 = _v11-_v10;
497 
498  Scalar fA00 = d0.sqrnorm();
499  Scalar fA01 = -(d0|d1);
500  Scalar fA11 = d1.sqrnorm();
501  Scalar fB0 = (kDiff|d0);
502  Scalar fC = kDiff.sqrnorm();
503  Scalar fDet = fabs(fA00*fA11-fA01*fA01);
504  Scalar fB1, fS, fT, fSqrDist, fTmp;
505 
506 
507  if ( fDet >= FLT_MIN )
508  {
509  // line segments are not parallel
510  fB1 = -(kDiff|d1);
511  fS = fA01*fB1-fA11*fB0;
512  fT = fA01*fB0-fA00*fB1;
513 
514 
515  // conservative approximation only?
516  if (_fastApprox)
517  {
518  if ( fS < 0.0 ) fS = 0.0;
519  else if ( fS > fDet ) fS = fDet;
520  if ( fT < 0.0 ) fT = 0.0;
521  else if ( fT > fDet ) fT = fDet;
522  }
523 
524 
525  if ( fS >= 0.0 )
526  {
527  if ( fS <= fDet )
528  {
529  if ( fT >= 0.0 )
530  {
531  if ( fT <= fDet ) // region 0 (interior)
532  {
533  // minimum at two interior points of 3D lines
534  Scalar fInvDet = 1.0/fDet;
535  fS *= fInvDet;
536  fT *= fInvDet;
537  fSqrDist = fS*(fA00*fS+fA01*fT+2.0*fB0) +
538  fT*(fA01*fS+fA11*fT+2.0*fB1)+fC;
539  }
540  else // region 3 (side)
541  {
542  fT = 1.0;
543  fTmp = fA01+fB0;
544  if ( fTmp >= 0.0 )
545  {
546  fS = 0.0;
547  fSqrDist = fA11+2.0*fB1+fC;
548  }
549  else if ( -fTmp >= fA00 )
550  {
551  fS = 1.0;
552  fSqrDist = fA00+fA11+fC+2.0*(fB1+fTmp);
553  }
554  else
555  {
556  fS = -fTmp/fA00;
557  fSqrDist = fTmp*fS+fA11+2.0*fB1+fC;
558  }
559  }
560  }
561  else // region 7 (side)
562  {
563  fT = 0.0;
564  if ( fB0 >= 0.0 )
565  {
566  fS = 0.0;
567  fSqrDist = fC;
568  }
569  else if ( -fB0 >= fA00 )
570  {
571  fS = 1.0;
572  fSqrDist = fA00+2.0*fB0+fC;
573  }
574  else
575  {
576  fS = -fB0/fA00;
577  fSqrDist = fB0*fS+fC;
578  }
579  }
580  }
581  else
582  {
583  if ( fT >= 0.0 )
584  {
585  if ( fT <= fDet ) // region 1 (side)
586  {
587  fS = 1.0;
588  fTmp = fA01+fB1;
589  if ( fTmp >= 0.0 )
590  {
591  fT = 0.0;
592  fSqrDist = fA00+2.0*fB0+fC;
593  }
594  else if ( -fTmp >= fA11 )
595  {
596  fT = 1.0;
597  fSqrDist = fA00+fA11+fC+2.0*(fB0+fTmp);
598  }
599  else
600  {
601  fT = -fTmp/fA11;
602  fSqrDist = fTmp*fT+fA00+2.0*fB0+fC;
603  }
604  }
605  else // region 2 (corner)
606  {
607  fTmp = fA01+fB0;
608  if ( -fTmp <= fA00 )
609  {
610  fT = 1.0;
611  if ( fTmp >= 0.0 )
612  {
613  fS = 0.0;
614  fSqrDist = fA11+2.0*fB1+fC;
615  }
616  else
617  {
618  fS = -fTmp/fA00;
619  fSqrDist = fTmp*fS+fA11+2.0*fB1+fC;
620  }
621  }
622  else
623  {
624  fS = 1.0;
625  fTmp = fA01+fB1;
626  if ( fTmp >= 0.0 )
627  {
628  fT = 0.0;
629  fSqrDist = fA00+2.0*fB0+fC;
630  }
631  else if ( -fTmp >= fA11 )
632  {
633  fT = 1.0;
634  fSqrDist = fA00+fA11+fC+2.0*(fB0+fTmp);
635  }
636  else
637  {
638  fT = -fTmp/fA11;
639  fSqrDist = fTmp*fT+fA00+2.0*fB0+fC;
640  }
641  }
642  }
643  }
644  else // region 8 (corner)
645  {
646  if ( -fB0 < fA00 )
647  {
648  fT = 0.0;
649  if ( fB0 >= 0.0 )
650  {
651  fS = 0.0;
652  fSqrDist = fC;
653  }
654  else
655  {
656  fS = -fB0/fA00;
657  fSqrDist = fB0*fS+fC;
658  }
659  }
660  else
661  {
662  fS = 1.0;
663  fTmp = fA01+fB1;
664  if ( fTmp >= 0.0 )
665  {
666  fT = 0.0;
667  fSqrDist = fA00+2.0*fB0+fC;
668  }
669  else if ( -fTmp >= fA11 )
670  {
671  fT = 1.0;
672  fSqrDist = fA00+fA11+fC+2.0*(fB0+fTmp);
673  }
674  else
675  {
676  fT = -fTmp/fA11;
677  fSqrDist = fTmp*fT+fA00+2.0*fB0+fC;
678  }
679  }
680  }
681  }
682  }
683  else
684  {
685  if ( fT >= 0.0 )
686  {
687  if ( fT <= fDet ) // region 5 (side)
688  {
689  fS = 0.0;
690  if ( fB1 >= 0.0 )
691  {
692  fT = 0.0;
693  fSqrDist = fC;
694  }
695  else if ( -fB1 >= fA11 )
696  {
697  fT = 1.0;
698  fSqrDist = fA11+2.0*fB1+fC;
699  }
700  else
701  {
702  fT = -fB1/fA11;
703  fSqrDist = fB1*fT+fC;
704  }
705  }
706  else // region 4 (corner)
707  {
708  fTmp = fA01+fB0;
709  if ( fTmp < 0.0 )
710  {
711  fT = 1.0;
712  if ( -fTmp >= fA00 )
713  {
714  fS = 1.0;
715  fSqrDist = fA00+fA11+fC+2.0*(fB1+fTmp);
716  }
717  else
718  {
719  fS = -fTmp/fA00;
720  fSqrDist = fTmp*fS+fA11+2.0*fB1+fC;
721  }
722  }
723  else
724  {
725  fS = 0.0;
726  if ( fB1 >= 0.0 )
727  {
728  fT = 0.0;
729  fSqrDist = fC;
730  }
731  else if ( -fB1 >= fA11 )
732  {
733  fT = 1.0;
734  fSqrDist = fA11+2.0*fB1+fC;
735  }
736  else
737  {
738  fT = -fB1/fA11;
739  fSqrDist = fB1*fT+fC;
740  }
741  }
742  }
743  }
744  else // region 6 (corner)
745  {
746  if ( fB0 < 0.0 )
747  {
748  fT = 0.0;
749  if ( -fB0 >= fA00 )
750  {
751  fS = 1.0;
752  fSqrDist = fA00+2.0*fB0+fC;
753  }
754  else
755  {
756  fS = -fB0/fA00;
757  fSqrDist = fB0*fS+fC;
758  }
759  }
760  else
761  {
762  fS = 0.0;
763  if ( fB1 >= 0.0 )
764  {
765  fT = 0.0;
766  fSqrDist = fC;
767  }
768  else if ( -fB1 >= fA11 )
769  {
770  fT = 1.0;
771  fSqrDist = fA11+2.0*fB1+fC;
772  }
773  else
774  {
775  fT = -fB1/fA11;
776  fSqrDist = fB1*fT+fC;
777  }
778  }
779  }
780  }
781  }
782  else
783  {
784  // line segments are parallel
785  if ( fA01 > 0.0 )
786  {
787  // direction vectors form an obtuse angle
788  if ( fB0 >= 0.0 )
789  {
790  fS = 0.0;
791  fT = 0.0;
792  fSqrDist = fC;
793  }
794  else if ( -fB0 <= fA00 )
795  {
796  fS = -fB0/fA00;
797  fT = 0.0;
798  fSqrDist = fB0*fS+fC;
799  }
800  else
801  {
802  fB1 = -(kDiff|d1);
803  fS = 1.0;
804  fTmp = fA00+fB0;
805  if ( -fTmp >= fA01 )
806  {
807  fT = 1.0;
808  fSqrDist = fA00+fA11+fC+2.0*(fA01+fB0+fB1);
809  }
810  else
811  {
812  fT = -fTmp/fA01;
813  fSqrDist = fA00+2.0*fB0+fC+fT*(fA11*fT+2.0*(fA01+fB1));
814  }
815  }
816  }
817  else
818  {
819  // direction vectors form an acute angle
820  if ( -fB0 >= fA00 )
821  {
822  fS = 1.0;
823  fT = 0.0;
824  fSqrDist = fA00+2.0*fB0+fC;
825  }
826  else if ( fB0 <= 0.0 )
827  {
828  fS = -fB0/fA00;
829  fT = 0.0;
830  fSqrDist = fB0*fS+fC;
831  }
832  else
833  {
834  fB1 = -(kDiff|d1);
835  fS = 0.0;
836  if ( fB0 >= -fA01 )
837  {
838  fT = 1.0;
839  fSqrDist = fA11+2.0*fB1+fC;
840  }
841  else
842  {
843  fT = -fB0/fA01;
844  fSqrDist = fC+fT*(2.0*fB1+fA11*fT);
845  }
846  }
847  }
848  }
849 
850 
851  if (_min_v0) *_min_v0 = _v00 + fS*d0;
852  if (_min_v1) *_min_v1 = _v10 + fT*d1;
853 
854  return fabs(fSqrDist);
855 }
856 
857 //-----------------------------------------------------------------------------
858 
859 template < typename VectorT , typename ValueT >
860 inline
861 ValueT
862 distPointPlane(const VectorT& _porigin,
863  const VectorT& _pnormal,
864  const VectorT& _point)
865 {
866  assert( fabs(_pnormal.norm()) > 0.0000000001) ;
867  return( ( (_point - _porigin) | _pnormal ) );
868 }
869 
870 
871 //-----------------------------------------------------------------------------
872 
873 template < typename VectorT >
874 VectorT projectToEdge(const VectorT& _start , const VectorT& _end , const VectorT& _point )
875 {
876  VectorT direction = ( _end - _start ).normalize();
877  assert( fabs(direction.norm()) > 0.0000000001) ;
878  const VectorT projected_point = ( ( _point - _start ) | direction ) * direction + _start;
879 
880  if ( ( ( projected_point - _start ) | direction ) > ( ( _end - _start ) | direction ) )
881  return _end;
882 
883  if ( ( ( projected_point - _start ) | direction ) < 0 )
884  return _start;
885 
886  return projected_point;
887 }
888 
889 //-----------------------------------------------------------------------------
890 
891 template < typename VectorT >
892 inline
893 VectorT
894 projectToPlane(const VectorT& _porigin,
895  const VectorT& _pnormal,
896  const VectorT& _point)
897 {
898  return (_point - _pnormal * distPointPlane< VectorT , double >( _porigin , _pnormal , _point ) );
899 }
900 
901 //-----------------------------------------------------------------------------
902 
903 
904 template<typename Scalar>
905 bool
907  const VectorT<Scalar,3>& _v1,
908  const VectorT<Scalar,3>& _v2,
909  VectorT<Scalar,3>& _result )
910 {
911  VectorT<Scalar,3> a(_v1-_v2),
912  b(_v2-_v0),
913  c(_v0-_v1),
914  G((_v0+_v1+_v2)/3.0);
915 
916  Scalar d0 = -(b|c),
917  d1 = -(c|a),
918  d2 = -(a|b),
919  e0 = d1*d2,
920  e1 = d2*d0,
921  e2 = d0*d1,
922  e = e0+e1+e2;
923 
924  if (e <= NumLimitsT<Scalar>::min()) return false;
925 
926  VectorT<Scalar,3> H((e0*_v0 + e1*_v1 + e2*_v2)/e);
927 
928  _result = (Scalar)0.5 * ((Scalar)3.0*G - H);
929 
930  return true;
931 }
932 
933 
934 
935 //-----------------------------------------------------------------------------
936 
937 
938 template<typename Scalar>
939 Scalar
941  const VectorT<Scalar,3>& _v1,
942  const VectorT<Scalar,3>& _v2 )
943 {
944  VectorT<Scalar,3> v0v1(_v1-_v0),
945  v0v2(_v2-_v0),
946  v1v2(_v2-_v1);
947 
948  Scalar denom = 4.0*((v0v1%v0v2).sqrnorm());
949  if (denom < NumLimitsT<Scalar>::min() &&
950  denom > -NumLimitsT<Scalar>::min())
951  return NumLimitsT<Scalar>::max();
952 
953  return ( v0v1.sqrnorm() *
954  v0v2.sqrnorm() *
955  v1v2.sqrnorm() /
956  denom );
957 }
958 
959 //-----------------------------------------------------------------------------
960 
961 template<typename Scalar>
962 bool
964  const VectorT<Scalar,3>& _v1,
965  VectorT<Scalar,3>& _axis,
966  Scalar& _angle,
967  bool _degree ) {
968 
969  // Copy axes
970  VectorT<Scalar,3> v0 = _v0;
971  VectorT<Scalar,3> v1 = _v1;
972 
973  // Normalize axes
974  v0.normalize();
975  v1.normalize();
976 
977  // Get rotation axis
978  _axis = (v0 % v1).normalize();
979 
980  // Is nan?
981  if ( std::isnan(_axis[0]) || std::isnan(_axis[1]) || std::isnan(_axis[2]) ) {
982  return false;
983  }
984 
985  // Get rotation angle (in radiant)
986  _angle = acos(v0 | v1);
987 
988  // Is nan?
989  if ( std::isnan(_angle) )
990  _angle = 0.0;
991 
992  // Convert to degree
993  if(_degree) {
994  _angle *= 180.0 / M_PI;
995  }
996 
997  return true;
998 }
999 
1000 //-----------------------------------------------------------------------------
1001 
1002 template<class VectorT>
1003 int isObtuse(const VectorT& _p0,
1004  const VectorT& _p1,
1005  const VectorT& _p2 )
1006 {
1007  const double a0 = ((_p1-_p0)|(_p2-_p0));
1008 
1009  if ( a0<0.0) return 1;
1010  else
1011  {
1012  const double a1 = ((_p0-_p1)|(_p2-_p1));
1013 
1014  if (a1 < 0.0 ) return 2;
1015  else
1016  {
1017  const double a2 = ((_p0-_p2)|(_p1-_p2));
1018 
1019  if (a2 < 0.0 ) return 3;
1020  else return 0;
1021  }
1022  }
1023 }
1024 
1025 //-----------------------------------------------------------------------------
1026 
1027 
1028 template<typename Scalar>
1029 bool
1031  const VectorT<Scalar,3>& _v1,
1032  const VectorT<Scalar,3>& _v2,
1033  VectorT<Scalar,3>& _center,
1034  Scalar& _radius)
1035 {
1036  VectorT<Scalar,3> a(_v1-_v2),
1037  b(_v2-_v0),
1038  c(_v0-_v1);
1039 
1040  Scalar d0 = -(b|c),
1041  d1 = -(c|a),
1042  d2 = -(a|b);
1043 
1044 
1045  // obtuse angle ?
1046  if (d2 < NumLimitsT<Scalar>::min())
1047  {
1048  _center = (_v0+_v1)*0.5;
1049  _radius = 0.5 * c.norm();
1050  return true;
1051  }
1052  if (d0 < NumLimitsT<Scalar>::min())
1053  {
1054  _center = (_v1+_v2)*0.5;
1055  _radius = 0.5 * a.norm();
1056  return true;
1057  }
1058  if (d1 < NumLimitsT<Scalar>::min())
1059  {
1060  _center = (_v2+_v0)*0.5;
1061  _radius = 0.5 * b.norm();
1062  return true;
1063  }
1064 
1065 
1066  // acute angle
1067  VectorT<Scalar,3> G((_v0+_v1+_v2)/3.0);
1068 
1069  Scalar e0 = d1*d2,
1070  e1 = d2*d0,
1071  e2 = d0*d1,
1072  e = e0+e1+e2;
1073 
1074  if ( e <= NumLimitsT<Scalar>::min() ) return false;
1075 
1076  VectorT<Scalar,3> H((e0*_v0 + e1*_v1 + e2*_v2)/e);
1077 
1078  _center = (Scalar)0.5 * ((Scalar)3.0*G - H);
1079  _radius = (_center-_v0).norm();
1080 
1081  return true;
1082 }
1083 
1084 
1085 //-----------------------------------------------------------------------------
1086 
1087 
1088 template<typename Scalar>
1089 Scalar
1091  const VectorT<Scalar,3>& _v1,
1092  const VectorT<Scalar,3>& _v2 )
1093 {
1094  VectorT<Scalar,3> v0v1(_v1-_v0),
1095  v0v2(_v2-_v0),
1096  v1v2(_v2-_v1);
1097 
1098  Scalar denom = 4.0*((v0v1%v0v2).sqrnorm());
1099  if (denom < NumLimitsT<Scalar>::min() &&
1100  denom > -NumLimitsT<Scalar>::min())
1101  return NumLimitsT<Scalar>::max();
1102 
1103  Scalar l0 = v0v1.sqrnorm(),
1104  l1 = v0v2.sqrnorm(),
1105  l2 = v1v2.sqrnorm(),
1106  l3 = l0*l1*l2/denom;
1107 
1108  return std::max(std::max(l0, l1), std::max(l2, l3));
1109 }
1110 
1111 
1112 //-----------------------------------------------------------------------------
1113 
1114 
1115 template<typename Scalar>
1116 bool
1118  const VectorT<Scalar,3>& _v1,
1119  const VectorT<Scalar,3>& _v2,
1120  const VectorT<Scalar,3>& _v3,
1121  VectorT<Scalar,3>& _result )
1122 {
1124  v0v1(_v1-_v0),
1125  v0v2(_v2-_v0),
1126  v0v3(_v3-_v0);
1127 
1128  Scalar denom = ((v0v1[0]*v0v2[1]*v0v3[2] +
1129  v0v1[1]*v0v2[2]*v0v3[0] +
1130  v0v1[2]*v0v2[0]*v0v3[1]) -
1131  (v0v1[0]*v0v2[2]*v0v3[1] +
1132  v0v1[1]*v0v2[0]*v0v3[2] +
1133  v0v1[2]*v0v2[1]*v0v3[0])) * 2.0;
1134 
1135  if (denom < NumLimitsT<Scalar>::min() &&
1136  denom > -NumLimitsT<Scalar>::min()) return false;
1137 
1138 
1139  _result = _v0 + (( v0v3.sqrnorm()*(v0v1%v0v2) +
1140  v0v2.sqrnorm()*(v0v3%v0v1) +
1141  v0v1.sqrnorm()*(v0v2%v0v3) ) / denom);
1142 
1143  return true;
1144 }
1145 
1146 
1147 //-----------------------------------------------------------------------------
1148 
1149 
1150 template <typename Scalar>
1153 {
1154  if (fabs(v[0]) < fabs(v[1])) {
1155  if (fabs(v[0]) < fabs(v[2]))
1156  return VectorT<Scalar, 3>(Scalar(1.0) - v[0] * v[0], -v[0] * v[1], -v[0] * v[2]).normalize();
1157  } else {
1158  if (fabs(v[1]) < fabs(v[2]))
1159  return VectorT<Scalar, 3>(-v[1] * v[0], Scalar(1.0) - v[1] * v[1], -v[1] * v[2]).normalize();
1160  }
1161 
1162  return VectorT<Scalar, 3>(-v[2] * v[0], -v[2] * v[1], Scalar(1.0) - v[2] * v[2]).normalize();
1163 }
1164 
1165 
1166 
1167 //== 2D STUFF ================================================================
1168 
1169 
1170 
1171 template<typename Scalar>
1172 bool
1174  const VectorT<Scalar,2> & _u,
1175  const VectorT<Scalar,2> & _v,
1176  const VectorT<Scalar,2> & _w,
1177  VectorT<Scalar,3> & _result )
1178 {
1179  Scalar v0(_v[0]-_u[0]), v1(_v[1]-_u[1]),
1180  w0(_w[0]-_u[0]), w1(_w[1]-_u[1]),
1181  p0(_p[0]-_u[0]), p1(_p[1]-_u[1]),
1182  denom = v0*w1-v1*w0;
1183 
1184  if (1.0+fabs(denom) == 1.0) {
1185  std::cerr << "degen tri: ("
1186  << _u << "), ("
1187  << _v << "), ("
1188  << _w << ")\n";
1189  return false;
1190  }
1191 
1192  _result[1] = 1.0 + ((p0*w1-p1*w0)/denom) - 1.0;
1193  _result[2] = 1.0 + ((v0*p1-v1*p0)/denom) - 1.0;
1194  _result[0] = 1.0 - _result[1] - _result[2];
1195 
1196  return true;
1197 }
1198 
1199 
1200 //-----------------------------------------------------------------------------
1201 
1202 
1203 template<typename Scalar>
1204 bool
1206  const VectorT<Scalar,2>& _v1,
1207  const VectorT<Scalar,2>& _v2,
1208  const VectorT<Scalar,2>& _v3,
1209  Scalar& _t1,
1210  Scalar& _t2 )
1211 {
1212  _t1 = ((_v0[1]-_v2[1])*(_v3[0]-_v2[0])-(_v0[0]-_v2[0])*(_v3[1]-_v2[1])) /
1213  ((_v1[0]-_v0[0])*(_v3[1]-_v2[1])-(_v1[1]-_v0[1])*(_v3[0]-_v2[0]));
1214 
1215  _t2 = ((_v0[1]-_v2[1])*(_v1[0]-_v0[0])-(_v0[0]-_v2[0])*(_v1[1]-_v0[1])) /
1216  ((_v1[0]-_v0[0])*(_v3[1]-_v2[1])-(_v1[1]-_v0[1])*(_v3[0]-_v2[0]));
1217 
1218  return ((_t1>0.0) && (_t1<1.0) && (_t2>0.0) && (_t2<1.0));
1219 }
1220 
1221 
1222 //-----------------------------------------------------------------------------
1223 
1224 template<typename Scalar>
1225 bool
1227  const VectorT<Scalar,2>& _v1,
1228  const VectorT<Scalar,2>& _v2,
1229  VectorT<Scalar,2>& _result )
1230 {
1231  Scalar x0(_v0[0]), y0(_v0[1]), xy0(x0*x0+y0*y0),
1232  x1(_v1[0]), y1(_v1[1]), xy1(x1*x1+y1*y1),
1233  x2(_v2[0]), y2(_v2[1]), xy2(x2*x2+y2*y2),
1234  a(x0*y1 - x0*y2 - x1*y0 + x1*y2 + x2*y0 - x2*y1),
1235  b(xy0*y1 - xy0*y2 - xy1*y0 + xy1*y2 + xy2*y0 - xy2*y1),
1236  c(xy0*x1 - xy0*x2 - xy1*x0 + xy1*x2 + xy2*x0 - xy2*x1);
1237 
1238  if (Scalar(1.0)+a == Scalar(1.0)) {
1239  std::cerr << "circumcircle: colinear points\n";
1240  return false;
1241  }
1242 
1243  _result[0] = 0.5*b/a;
1244  _result[1] = -0.5*c/a;
1245 
1246  return true;
1247 }
1248 
1249 
1250 
1251 //== N-DIM STUFF ==============================================================
1252 
1253 
1254 template <typename Scalar, int N>
1255 Scalar
1257  const VectorT<Scalar, N>& _v1,
1258  const VectorT<Scalar, N>& _v2 )
1259 {
1260  VectorT<Scalar,3> d0 = _v0 - _v1;
1261  VectorT<Scalar,3> d1 = _v1 - _v2;
1262 
1263  // finds the max squared edge length
1264  Scalar l2, maxl2 = d0.sqrnorm();
1265  if ((l2=d1.sqrnorm()) > maxl2)
1266  maxl2 = l2;
1267  // keep searching for the max squared edge length
1268  d1 = _v2 - _v0;
1269  if ((l2=d1.sqrnorm()) > maxl2)
1270  maxl2 = l2;
1271 
1272  // squared area of the parallelogram spanned by d0 and d1
1273  Scalar a2 = (d0 % d1).sqrnorm();
1274 
1275  // the area of the triangle would be
1276  // sqrt(a2)/2 or length * height / 2
1277  // aspect ratio = length / height
1278  // = length * length / (2*area)
1279  // = length * length / sqrt(a2)
1280 
1281  // returns the length of the longest edge
1282  // divided by its corresponding height
1283  return sqrt( (maxl2 * maxl2) / a2 );
1284 }
1285 
1286 
1287 //-----------------------------------------------------------------------------
1288 
1289 
1290 template <typename Scalar, int N>
1291 Scalar
1293  const VectorT<Scalar, N>& _v1,
1294  const VectorT<Scalar, N>& _v2 )
1295 {
1296  const double FOUR_ROOT3 = 6.928203230275509;
1297 
1298  double area = 0.5*((_v1-_v0)%(_v2-_v0)).norm();
1299 
1300  return (Scalar) (FOUR_ROOT3 * area / ((_v0-_v1).sqrnorm() +
1301  (_v1-_v2).sqrnorm() +
1302  (_v2-_v0).sqrnorm() ));
1303 }
1304 
1305 template<typename Vec>
1306 bool
1307 triangleIntersection( const Vec& _o,
1308  const Vec& _dir,
1309  const Vec& _v0,
1310  const Vec& _v1,
1311  const Vec& _v2,
1312  typename Vec::value_type& _t,
1313  typename Vec::value_type& _u,
1314  typename Vec::value_type& _v )
1315 {
1316  //This code effectively replicates the method described by Moeller et al. in "Fast, Minimum Storage Ray-Triangle Intersection".
1317  Vec edge1, edge2, tvec, pvec, qvec;
1318  typename Vec::value_type det, inv_det;
1319 
1320  //find vectors for two edges sharing v0
1321  edge1 = _v1-_v0;
1322  edge2 = _v2-_v0;
1323 
1324  //begin calculating determinant - also used to calculate u parameter
1325  pvec = _dir % edge2;
1326 
1327  //if determinant is near zero, the ray lies in plane of triangle
1328  det = edge1 | pvec;
1329 
1330  static const typename Vec::value_type EPSILON = std::numeric_limits<typename Vec::value_type>::epsilon() * 1e2;
1331  if (det > -EPSILON && det < EPSILON) {
1332  return false;
1333  }
1334  inv_det = typename Vec::value_type(1.0) / det;
1335 
1336  //calculate distance from vert0 to ray origin
1337  tvec = _o - _v0;
1338 
1339  //calculate U parameter and test bounds
1340  _u = (tvec | pvec) * inv_det;
1341  if (_u < 0.0 || _u > 1.0)
1342  return false;
1343 
1344  //prepare to test V parameter
1345  qvec = tvec % edge1;
1346 
1347  //calculate V parameter and test bounds
1348  _v = (_dir | qvec) * inv_det;
1349  if (_v < 0.0 || _u + _v > 1.0)
1350  return false;
1351 
1352  //Intersection found! Calculate t and exit...
1353  _t = (edge2 | qvec) * inv_det;
1354  return true;
1355 }
1356 
1357 template<typename Vec>
1358 bool
1360  const Vec& _dir,
1361  const Vec& _bbmin,
1362  const Vec& _bbmax,
1363  typename Vec::value_type& tmin,
1364  typename Vec::value_type& tmax )
1365 {
1366  /*
1367  * Ray-box intersection using IEEE numerical properties to ensure that the
1368  * test is both robust and efficient, as described in:
1369  *
1370  * Amy Williams, Steve Barrus, R. Keith Morley, and Peter Shirley
1371  * "An Efficient and Robust Ray-Box Intersection Algorithm"
1372  * Journal of graphics tools, 10(1):49-54, 2005
1373  *
1374  */
1375  typename Vec::value_type tymin, tymax, tzmin, tzmax;
1376  Vec inv_dir;
1377 
1378  inv_dir[0] = 1/_dir[0];
1379  inv_dir[1] = 1/_dir[1];
1380  inv_dir[2] = 1/_dir[2];
1381 
1382  if (inv_dir[0] >= 0) {
1383  tmin = (_bbmin[0] - _o[0]) * inv_dir[0];
1384  tmax = (_bbmax[0] - _o[0]) * inv_dir[0];
1385  }
1386  else {
1387  tmin = (_bbmax[0] - _o[0]) * inv_dir[0];
1388  tmax = (_bbmin[0] - _o[0]) * inv_dir[0];
1389  }
1390 
1391  if (inv_dir[1] >= 0) {
1392  tymin = (_bbmin[1] - _o[1]) * inv_dir[1];
1393  tymax = (_bbmax[1] - _o[1]) * inv_dir[1];
1394  }
1395  else {
1396  tymin = (_bbmax[1] - _o[1]) * inv_dir[1];
1397  tymax = (_bbmin[1] - _o[1]) * inv_dir[1];
1398  }
1399 
1400  if ( (tmin > tymax) || (tymin > tmax) )
1401  return false;
1402  if (tymin > tmin)
1403  tmin = tymin;
1404  if (tymax < tmax)
1405  tmax = tymax;
1406 
1407  if (inv_dir[2] >= 0) {
1408  tzmin = (_bbmin[2] - _o[2]) * inv_dir[2];
1409  tzmax = (_bbmax[2] - _o[2]) * inv_dir[2];
1410  }
1411  else {
1412  tzmin = (_bbmax[2] - _o[2]) * inv_dir[2];
1413  tzmax = (_bbmin[2] - _o[2]) * inv_dir[2];
1414  }
1415 
1416  if ( (tmin > tzmax) || (tzmin > tmax) )
1417  return false;
1418  if (tzmin > tmin)
1419  tmin = tzmin;
1420  if (tzmax < tmax)
1421  tmax = tzmax;
1422 
1423  return true;
1424 }
1425 
1426 //=============================================================================
1427 } // namespace Geometry
1428 } // namespace ACG
1429 //=============================================================================
Namespace providing different geometric functions concerning angles.
VectorT projectToEdge(const VectorT &_start, const VectorT &_end, const VectorT &_point)
Definition: Algorithms.cc:874
auto norm() const -> decltype(std::sqrt(std::declval< VectorT< S, DIM >>().sqrnorm()))
compute euclidean norm
Definition: Vector11T.hh:409
decltype(std::declval< S >() *std::declval< S >()) sqrnorm() const
compute squared euclidean norm
Definition: Vector11T.hh:397
ValueT distPointPlane(const VectorT &_porigin, const VectorT &_pnormal, const VectorT &_point)
Checks the distance from a point to a plane.
Definition: Algorithms.cc:862
Scalar distLineLineSquared(const VectorT< Scalar, 3 > &_v00, const VectorT< Scalar, 3 > &_v01, const VectorT< Scalar, 3 > &_v10, const VectorT< Scalar, 3 > &_v11, VectorT< Scalar, 3 > *_min_v0, VectorT< Scalar, 3 > *_min_v1, bool _fastApprox)
squared distance of lines (_v00, _v01) and (_v10, _v11)
Definition: Algorithms.cc:486
Scalar roundness(const VectorT< Scalar, N > &_v0, const VectorT< Scalar, N > &_v1, const VectorT< Scalar, N > &_v2)
return roundness of triangle: 1=equilateral, 0=colinear
Definition: Algorithms.cc:1292
int isObtuse(const VectorT &_p0, const VectorT &_p1, const VectorT &_p2)
Definition: Algorithms.cc:1003
Vec::value_type triangleAreaSquared(const Vec &_v0, const Vec &_v1, const Vec &_v2)
return squared area of triangle (_v0, _v1, _v2)
Definition: Algorithms.cc:163
bool circumCenter(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, const VectorT< Scalar, 3 > &_v2, VectorT< Scalar, 3 > &_result)
return circumcenter of triangle (_v0,_v1,_v2)
Definition: Algorithms.cc:906
auto normalize() -> decltype(*this/=std::declval< VectorT< S, DIM >>().norm())
Definition: Vector11T.hh:429
bool edgeConvexPolygonIntersection(std::vector< VectorT< Scalar, 3 > > _polygon_points, VectorT< Scalar, 3 > _v0, VectorT< Scalar, 3 > _v1, VectorT< Scalar, 3 > &_result)
Get intersection point of a ray and a convex polygon.
Definition: Algorithms.cc:172
VectorT projectToPlane(const VectorT &_porigin, const VectorT &_pnormal, const VectorT &_point)
projects a point to a plane
Definition: Algorithms.cc:894
Vec::value_type distPointTriangleSquaredStable(const Vec &_p, const Vec &_v0, const Vec &_v1, const Vec &_v2, Vec &_nearestPoint)
squared distance from point _p to triangle (_v0, _v1, _v2)
Definition: Algorithms.cc:466
Vec::value_type distPointLineSquared(const Vec &_p, const Vec &_v0, const Vec &_v1, Vec *_min_v)
squared distance from point _p to line segment (_v0,_v1)
Definition: Algorithms.cc:290
Scalar aspectRatio(const VectorT< Scalar, N > &_v0, const VectorT< Scalar, N > &_v1, const VectorT< Scalar, N > &_v2)
return aspect ratio (length/height) of triangle
Definition: Algorithms.cc:1256
VectorT< Scalar, 3 > perpendicular(const VectorT< Scalar, 3 > &v)
find a vector that&#39;s perpendicular to _v
Definition: Algorithms.cc:1152
4x4 matrix implementing OpenGL commands.
Definition: GLMatrixT.hh:79
bool rotationOfTwoVectors(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, VectorT< Scalar, 3 > &_axis, Scalar &_angle, bool _degree)
Get rotation axis and signed angle of rotation between two vectors.
Definition: Algorithms.cc:963
bool axisAlignedBBIntersection(const Vec &_o, const Vec &_dir, const Vec &_bbmin, const Vec &_bbmax, typename Vec::value_type &tmin, typename Vec::value_type &tmax)
Intersect a ray and an axis aligned bounding box.
Definition: Algorithms.cc:1359
bool minSphere(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, const VectorT< Scalar, 3 > &_v2, VectorT< Scalar, 3 > &_center, Scalar &_radius)
construct min. enclosing sphere
Definition: Algorithms.cc:1030
Scalar circumRadiusSquared(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, const VectorT< Scalar, 3 > &_v2)
return squared radius of circumcircle of triangle (_v0,_v1,_v2)
Definition: Algorithms.cc:940
static Scalar min()
Return the smallest absolte value a scalar type can store.
Definition: NumLimitsT.hh:96
bool baryCoord(const VectorT< Scalar, 3 > &_p, const VectorT< Scalar, 3 > &_u, const VectorT< Scalar, 3 > &_v, const VectorT< Scalar, 3 > &_w, VectorT< Scalar, 3 > &_result)
Definition: Algorithms.cc:83
void identity()
setup an identity matrix
Scalar minRadiusSquared(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, const VectorT< Scalar, 3 > &_v2)
return squared radius of min. enclosing circle of triangle (_v0,_v1,_v2)
Definition: Algorithms.cc:1090
void rotate(Scalar angle, Scalar x, Scalar y, Scalar z, MultiplyFrom _mult_from=MULT_FROM_RIGHT)
bool triangleIntersection(const Vec &_o, const Vec &_dir, const Vec &_v0, const Vec &_v1, const Vec &_v2, typename Vec::value_type &_t, typename Vec::value_type &_u, typename Vec::value_type &_v)
Intersect a ray and a triangle.
Definition: Algorithms.cc:1307
bool lineIntersection(const VectorT< Scalar, 2 > &_v0, const VectorT< Scalar, 2 > &_v1, const VectorT< Scalar, 2 > &_v2, const VectorT< Scalar, 2 > &_v3, Scalar &_t1, Scalar &_t2)
intersect two line segments (_v0,_v1) and (_v2,_v3)
Definition: Algorithms.cc:1205