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