Developer Documentation
Loading...
Searching...
No Matches
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
71namespace ACG {
72namespace Geometry {
73
74
75//== IMPLEMENTATION ==========================================================
76
77
78//== 3D STUFF ================================================================
79
80
81template<typename Scalar>
82bool
83baryCoord( const VectorT<Scalar,3> & _p,
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
161template <class Vec>
162typename Vec::value_type
163triangleAreaSquared( const Vec& _v0,
164 const Vec& _v1,
165 const Vec& _v2 )
166{
167 return 0.25 * ((_v1-_v0)%(_v2-_v0)).sqrnorm();
168}
169
170//-----------------------------------------------------------------------------
171template<typename Scalar>
172bool 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
217 GLMatrixT<Scalar> R;
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
288template<class Vec>
289typename Vec::value_type
290distPointLineSquared( 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
307template <class Vec>
308typename Vec::value_type
309distPointTriangleSquared( 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
449template <class Vec>
450typename Vec::value_type
451distPointTriangleSquared( 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
464template <class Vec>
465typename Vec::value_type
466distPointTriangleSquaredStable( const Vec& _p,
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
484template<typename Scalar>
485Scalar
486distLineLineSquared( const VectorT<Scalar,3>& _v00,
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
859template < typename VectorT , typename ValueT >
860inline
861ValueT
862distPointPlane(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
873template < typename VectorT >
874VectorT 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
891template < typename VectorT >
892inline
893VectorT
894projectToPlane(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
904template<typename Scalar>
905bool
906circumCenter( const VectorT<Scalar,3>& _v0,
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
938template<typename Scalar>
939Scalar
940circumRadiusSquared( const VectorT<Scalar,3>& _v0,
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
961template<typename Scalar>
962bool
963rotationOfTwoVectors( const VectorT<Scalar,3>& _v0,
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
1002template<class VectorT>
1003int 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
1028template<typename Scalar>
1029bool
1030minSphere(const VectorT<Scalar,3>& _v0,
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
1088template<typename Scalar>
1089Scalar
1090minRadiusSquared( const VectorT<Scalar,3>& _v0,
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
1115template<typename Scalar>
1116bool
1117circumCenter( const VectorT<Scalar,3>& _v0,
1118 const VectorT<Scalar,3>& _v1,
1119 const VectorT<Scalar,3>& _v2,
1120 const VectorT<Scalar,3>& _v3,
1121 VectorT<Scalar,3>& _result )
1122{
1123 VectorT<Scalar,3>
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
1150template <typename Scalar>
1151VectorT<Scalar,3>
1152perpendicular( const VectorT<Scalar,3>& v )
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
1171template<typename Scalar>
1172bool
1173baryCoord( const VectorT<Scalar,2> & _p,
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
1203template<typename Scalar>
1204bool
1205lineIntersection( const VectorT<Scalar,2>& _v0,
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
1224template<typename Scalar>
1225bool
1226circumCenter( const VectorT<Scalar,2>& _v0,
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
1254template <typename Scalar, int N>
1255Scalar
1256aspectRatio( const VectorT<Scalar, N>& _v0,
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
1290template <typename Scalar, int N>
1291Scalar
1292roundness( const VectorT<Scalar, N>& _v0,
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
1305template<typename Vec>
1306bool
1307triangleIntersection( 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
1357template<typename Vec>
1358bool
1359axisAlignedBBIntersection( const Vec& _o,
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//=============================================================================
ACG::GLMatrixT
4x4 matrix implementing OpenGL commands.
Definition GLMatrixT.hh:80
ACG::GLMatrixT::rotate
void rotate(Scalar angle, Scalar x, Scalar y, Scalar z, MultiplyFrom _mult_from=MULT_FROM_RIGHT)
Definition GLMatrixT_impl.hh:155
ACG::Matrix4x4T::identity
void identity()
setup an identity matrix
Definition Matrix4x4T_impl.hh:250
ACG::NumLimitsT::min
static Scalar min()
Return the smallest absolte value a scalar type can store.
Definition NumLimitsT.hh:95
OpenMesh::VectorT::sqrnorm
decltype(std::declval< S >() *std::declval< S >()) sqrnorm() const
compute squared euclidean norm
Definition Vector11T.hh:422
OpenMesh::VectorT::normalize
auto normalize() -> decltype(*this/=std::declval< VectorT< S, DIM > >().norm())
Definition Vector11T.hh:454
OpenMesh::VectorT::norm
auto norm() const -> decltype(std::sqrt(std::declval< VectorT< S, DIM > >().sqrnorm()))
compute euclidean norm
Definition Vector11T.hh:434
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:1090
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:906
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:1256
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:83
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:172
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:163
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:290
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:466
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:1292
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:963
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:486
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:1307
ACG::Geometry::projectToPlane
VectorT projectToPlane(const VectorT &_porigin, const VectorT &_pnormal, const VectorT &_point)
projects a point to a plane
Definition Algorithms.cc:894
ACG::Geometry::perpendicular
VectorT< Scalar, 3 > perpendicular(const VectorT< Scalar, 3 > &v)
find a vector that's perpendicular to _v
Definition Algorithms.cc:1152
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:1359
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:1205
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:862
ACG::Geometry::projectToEdge
VectorT projectToEdge(const VectorT &_start, const VectorT &_end, const VectorT &_point)
Definition Algorithms.cc:874
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:940
ACG::Geometry::isObtuse
int isObtuse(const VectorT &_p0, const VectorT &_p1, const VectorT &_p2)
Definition Algorithms.cc:1003
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:1030
ACG
Namespace providing different geometric functions concerning angles.
Definition BaseObjectData.hh:68