Developer Documentation
Algorithms.hh
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#ifndef GEO_ALGORITHMS_HH
45#define GEO_ALGORITHMS_HH
46
47
48//== INCLUDES =================================================================
49
50#include <cfloat>
51#include <ACG/Math/VectorT.hh>
52#include <vector>
53#include <iostream>
54
55#include "../Math/Matrix3x3T.hh"
56
57
58namespace ACG {
59namespace Geometry {
60
61
62//== 3D STUFF =================================================================
63
64
65
67template<typename Scalar>
68bool
69circumCenter( const VectorT<Scalar,3>& _v0,
70 const VectorT<Scalar,3>& _v1,
71 const VectorT<Scalar,3>& _v2,
72 const VectorT<Scalar,3>& _v3,
73 VectorT<Scalar,3>& _result );
74
75
77template<typename Scalar>
78Scalar
80 const VectorT<Scalar,3>& _v1,
81 const VectorT<Scalar,3>& _v2,
82 const VectorT<Scalar,3>& _v3 )
83{
85 return circumCenter(_v0, _v1, _v2, _v3, cc) ? (cc-_v0).sqrnorm() : FLT_MAX;
86}
87
88
90template<typename Scalar>
91Scalar
93 const VectorT<Scalar,3>& _v1,
94 const VectorT<Scalar,3>& _v2,
95 const VectorT<Scalar,3>& _v3 )
96{
97 return sqrt(circumRadiusSquared(_v0, _v1, _v2, _v3));
98}
99
112template<typename Scalar>
113bool edgeConvexPolygonIntersection(std::vector<VectorT<Scalar,3> > _polygon_points,
116 VectorT<Scalar,3> &_result);
117
118
133template<typename Scalar>
134bool
136 const VectorT<Scalar,3>& _v1,
137 VectorT<Scalar,3>& _axis,
138 Scalar& _angle,
139 bool _degree = true);
140
141
150template <typename Scalar>
153
154
171template<typename Vec>
172bool
173triangleIntersection( const Vec& _o,
174 const Vec& _dir,
175 const Vec& _v0,
176 const Vec& _v1,
177 const Vec& _v2,
178 typename Vec::value_type& _t,
179 typename Vec::value_type& _u,
180 typename Vec::value_type& _v );
181
182
195template<typename Vec>
196bool
197axisAlignedBBIntersection( const Vec& _o,
198 const Vec& _dir,
199 const Vec& _bbmin,
200 const Vec& _bbmax,
201 typename Vec::value_type& _t0,
202 typename Vec::value_type& _t1 );
203
204
205//== 2D STUFF =================================================================
206
208template<typename Scalar>
209Scalar
211 const VectorT<Scalar,2>& _v0,
212 const VectorT<Scalar,2>& _v1 )
213{
214 VectorT<Scalar,2> d1(_p-_v0), d2(_v1-_v0);
215 return (d1[0]*d2[1]-d1[1]*d2[0]);
216}
217
218
220template<typename Scalar>
221bool
223 const VectorT<Scalar,2>& _v1,
224 const VectorT<Scalar,2>& _v2 )
225{
226 return ( pointLineOrientation(_v0, _v1, _v2) < Scalar(0.0) );
227}
228
229
231template<typename Scalar>
232bool
234 const VectorT<Scalar,2>& _v1,
235 const VectorT<Scalar,2>& _v2 )
236{
237 return ( pointLineOrientation(_v0, _v1, _v2) > Scalar(0.0) );
238}
239
240
242template<typename Scalar>
243bool
245 const VectorT<Scalar,2>& _v1,
246 const VectorT<Scalar,2>& _v2,
247 const VectorT<Scalar,2>& _v3,
248 Scalar& _t1,
249 Scalar& _t2 );
250
251
252//===========================================================================
255//===========================================================================
256
257
266template<class Vec>
267typename Vec::value_type
268distPointLineSquared( const Vec& _p,
269 const Vec& _v0,
270 const Vec& _v1,
271 Vec* _min_v=0);
272
273
274
286template<class Vec>
287typename Vec::value_type
288distPointLine( const Vec& _p,
289 const Vec& _v0,
290 const Vec& _v1,
291 Vec* _min_v=0 )
292{ return sqrt(distPointLineSquared(_p, _v0, _v1, _min_v)); }
293
294
296template <class Vec>
297typename Vec::value_type
298distPointTriangleSquared( const Vec& _p,
299 const Vec& _v0,
300 const Vec& _v1,
301 const Vec& _v2,
302 Vec& _nearestPoint );
303
315template <class Vec>
316typename Vec::value_type
317distPointTriangleSquaredStable( const Vec& _p,
318 const Vec& _v0,
319 const Vec& _v1,
320 const Vec& _v2,
321 Vec& _nearestPoint );
322
324template <class Vec>
325typename Vec::value_type
326distPointTriangle( const Vec& _p,
327 const Vec& _v0,
328 const Vec& _v1,
329 const Vec& _v2,
330 Vec& _nearestPoint )
331{ return sqrt(distPointTriangleSquared(_p, _v0, _v1, _v2, _nearestPoint)); }
332
343template <class Vec>
344typename Vec::value_type
346 const Vec& _v0,
347 const Vec& _v1,
348 const Vec& _v2,
349 Vec& _nearestPoint )
350{ return sqrt(distPointTriangleSquaredStable(_p, _v0, _v1, _v2, _nearestPoint)); }
351
360template < typename VectorT , typename ValueT >
361inline
362ValueT
363distPointPlane(const VectorT& _porigin,
364 const VectorT& _pnormal,
365 const VectorT& _point);
366
367
370//===========================================================================
373//===========================================================================
374
376template<typename Scalar>
377Scalar
379 const VectorT<Scalar,3>& _v01,
380 const VectorT<Scalar,3>& _v10,
381 const VectorT<Scalar,3>& _v11,
382 VectorT<Scalar,3>* _min_v0=0,
383 VectorT<Scalar,3>* _min_v1=0,
384 bool _fastApprox=false );
385
386
388template<typename Scalar>
389Scalar
391 const VectorT<Scalar,3>& _v01,
392 const VectorT<Scalar,3>& _v10,
393 const VectorT<Scalar,3>& _v11,
394 VectorT<Scalar,3>* _min_v0=0,
395 VectorT<Scalar,3>* _min_v1=0 )
396{
397 return sqrt(distLineLineSquared(_v00, _v01, _v10, _v11,
398 _min_v0, _min_v1));
399}
400
404//===========================================================================
407//===========================================================================
408
409
416template < typename VectorT >
417VectorT projectToEdge(const VectorT& _start ,
418 const VectorT& _end ,
419 const VectorT& _point );
420
421
428template < typename VectorT >
429inline
431projectToPlane(const VectorT& _porigin,
432 const VectorT& _pnormal,
433 const VectorT& _point);
434
437//===========================================================================
440//===========================================================================
441
447template<typename Scalar>
448bool
450 const VectorT<Scalar,2>& _u,
451 const VectorT<Scalar,2>& _v,
452 const VectorT<Scalar,2>& _w,
453 VectorT<Scalar,3>& _result );
454
455
457template<typename Scalar>
458bool
460 const VectorT<Scalar,2>& _u,
461 const VectorT<Scalar,2>& _v,
462 const VectorT<Scalar,2>& _w )
463{
465 if (baryCoord(_p, _u, _v, _w, bary))
466 return ((bary[0]>0.0) && (bary[1]>0.0) && (bary[2]>0.0));
467 else {
468 std::cerr << "point in triangle error\n";
469 return false;
470 }
471}
472
473template<typename Scalar>
474bool
476 const VectorT<Scalar,2>& _v1,
477 const VectorT<Scalar,2>& _v2,
478 VectorT<Scalar,2>& _result );
479
482//===========================================================================
485//===========================================================================
486
489template<typename Scalar>
490bool
492 const VectorT<Scalar,3>& _u,
493 const VectorT<Scalar,3>& _v,
494 const VectorT<Scalar,3>& _w,
495 VectorT<Scalar,3>& _result );
496
497
499template<typename Scalar>
500bool
502 const VectorT<Scalar,3>& _v1,
503 const VectorT<Scalar,3>& _v2,
504 VectorT<Scalar,3>& _center,
505 Scalar& _radius);
506
507
509template<typename Scalar>
510Scalar
512 const VectorT<Scalar,3>& _v1,
513 const VectorT<Scalar,3>& _v2 );
514
515
517template<typename Scalar>
518Scalar
520 const VectorT<Scalar,3>& _v1,
521 const VectorT<Scalar,3>& _v2 )
522{
523 return sqrt(minRadiusSquared(_v0, _v1, _v2));
524}
525
526
528template<typename Scalar>
529bool
531 const VectorT<Scalar,3>& _v1,
532 const VectorT<Scalar,3>& _v2,
533 VectorT<Scalar,3>& _result );
534
535
537template<typename Scalar>
538Scalar
540 const VectorT<Scalar,3>& _v1,
541 const VectorT<Scalar,3>& _v2 );
542
543
545template<typename Scalar>
546Scalar
548 const VectorT<Scalar,3>& _v1,
549 const VectorT<Scalar,3>& _v2 )
550{
551 return sqrt(circumRadiusSquared(_v0, _v1, _v2));
552}
553
561template<class VectorT>
562int isObtuse(const VectorT& _p0,
563 const VectorT& _p1,
564 const VectorT& _p2 );
565
568//===========================================================================
571//===========================================================================
572
573
580template <class Vec>
581typename Vec::value_type
582triangleAreaSquared( const Vec& _v0,
583 const Vec& _v1,
584 const Vec& _v2 );
585
586
593template <class Vec>
594typename Vec::value_type
595triangleArea( const Vec& _v0,
596 const Vec& _v1,
597 const Vec& _v2 )
598{
599 return sqrt(triangleAreaSquared(_v0,_v1,_v2));
600}
601
602
609template <typename Scalar, int N>
610Scalar
612 const VectorT<Scalar, N>& _v1,
613 const VectorT<Scalar, N>& _v2 );
614
621template <typename Scalar, int N>
622Scalar
623roundness( const VectorT<Scalar, N>& _v0,
624 const VectorT<Scalar, N>& _v1,
625 const VectorT<Scalar, N>& _v2 );
626
629template<typename Vector>
630Vector closestPointLineSegment(Vector x, Vector p1, Vector p2) {
631 const auto delta = ((p2-p1)|(x-p1)) / (p2-p1).sqrnorm();
632 //std::cout << "\x1b[32mdelta = " << delta << "\x1b[0m" << std::endl;
633 if (delta <= 0) {
634 //std::cout << "\x1b[43mdelta <= 0\x1b[0m" << std::endl;
635 return p1;
636 } else if (delta >= 1) {
637 //std::cout << "\x1b[43mdelta >= 1\x1b[0m" << std::endl;
638 return p2;
639 } else if (delta != delta) { // p1 = p2 or x = p1
640 //std::cout << "\x1b[43mdelta != delta\x1b[0m" << std::endl;
641 return (x-p1).sqrnorm() < (x-p2).sqrnorm() ? p1 : p2;
642 } else {
643 //std::cout << "\x1b[43mdelta \\in [0, 1]\x1b[0m" << std::endl;
644 return (1 - delta) * p1 + delta * p2;
645 }
646};
647
648template<typename Vector>
649Vector closestPointTri(Vector p, Vector a, Vector b, Vector c) {
650 constexpr double thresh = 1e-8;
651
652 const auto n = ((b - a) % (c - a)); // normalization unnecessary
653
654 if ((b-a).sqrnorm() < thresh || (c-a).sqrnorm() < thresh || n.sqrnorm() < thresh) {
655 //std::cout << "\x1b[42mDegenerate case.\x1b[0m" << std::endl;
656 // Degenerate triangle. Find distance to longest segment.
657 std::array<ACG::Vec3d, 2> max_segment = {a, b};
658 double max_len = (b-a).sqrnorm();
659 if ((c-a).sqrnorm() > max_len)
660 max_segment = {a, c};
661 if ((c-b).sqrnorm() > max_len)
662 max_segment = {b, c};
663
664 // closestPointLineSegment is stable, even if the segment is super short
665 return closestPointLineSegment(p, max_segment[0], max_segment[1]);
666 }
667
668 const auto abd = Matrix3x3d::fromColumns(a-c, b-c, n).inverse() * (p - c);
669 const bool alpha = (abd[0] >= 0.0),
670 beta = (abd[1] >= 0.0),
671 gamma = (1.0-abd[0]-abd[1] >= 0.0);
672
673 if (alpha && beta && gamma) {
674 //std::cout << "\x1b[42mInside case.\x1b[0m" << std::endl;
675 // Inside triangle.
676 return abd[0] * a + abd[1] * b + (1.0 - abd[0] - abd[1]) * c;
677 } else if (!alpha) {
678 //std::cout << "\x1b[42m!alpha case.\x1b[0m" << std::endl;
679 // Closest to line segment (b, c).
680 return closestPointLineSegment(p, b, c);
681 } else if (!beta) {
682 //std::cout << "\x1b[42m!beta case.\x1b[0m" << std::endl;
683 // Closest to line segment (a, c).
684 return closestPointLineSegment(p, a, c);
685 } else if (!gamma) {
686 //std::cout << "\x1b[42m!gamma case.\x1b[0m" << std::endl;
687 // Closest to line segment (a, b).
688 return closestPointLineSegment(p, a, b);
689 } else {
690 throw std::logic_error("This cannot happen.");
691 }
692}
693
694//=============================================================================
695} // namespace Geometry
696} // namespace ACG
697//=============================================================================
698#if defined(INCLUDE_TEMPLATES) && !defined(GEO_ALGORITHMS_C)
699#define GEO_ALGORITHMS_TEMPLATES
700#include "Algorithms.cc"
701#endif
702//=============================================================================
703#endif // GEO_ALGORITHMS_HH defined
704//=============================================================================
705
decltype(std::declval< S >() *std::declval< S >()) sqrnorm() const
compute squared euclidean norm
Definition: Vector11T.hh:422
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
Scalar minRadius(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, const VectorT< Scalar, 3 > &_v2)
return radius of min. enclosing circle of triangle (_v0,_v1,_v2)
Definition: Algorithms.hh:519
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
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
Scalar distLineLine(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=0, VectorT< Scalar, 3 > *_min_v1=0)
distance of lines (_v00, _v01) and (_v10, _v11)
Definition: Algorithms.hh:390
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
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
bool isCW(const VectorT< Scalar, 2 > &_v0, const VectorT< Scalar, 2 > &_v1, const VectorT< Scalar, 2 > &_v2)
are 3 vertices in clockwise order? in 2D
Definition: Algorithms.hh:233
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 isInTriangle(const VectorT< Scalar, 2 > &_p, const VectorT< Scalar, 2 > &_u, const VectorT< Scalar, 2 > &_v, const VectorT< Scalar, 2 > &_w)
is point _p in triangle (_v0,_v1,_v2) in 2D
Definition: Algorithms.hh:459
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
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
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
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
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
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
VectorT projectToPlane(const VectorT &_porigin, const VectorT &_pnormal, const VectorT &_point)
projects a point to a plane
Definition: Algorithms.cc:894
Scalar pointLineOrientation(const VectorT< Scalar, 2 > &_p, const VectorT< Scalar, 2 > &_v0, const VectorT< Scalar, 2 > &_v1)
orientation of point _p w.r.t. line through _v0,_v1 in 2D
Definition: Algorithms.hh:210
VectorT< Scalar, 3 > perpendicular(const VectorT< Scalar, 3 > &v)
find a vector that's perpendicular to _v
Definition: Algorithms.cc:1152
Vec::value_type triangleArea(const Vec &_v0, const Vec &_v1, const Vec &_v2)
return area of triangle (_v0, _v1, _v2)
Definition: Algorithms.hh:595
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 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
bool isCCW(const VectorT< Scalar, 2 > &_v0, const VectorT< Scalar, 2 > &_v1, const VectorT< Scalar, 2 > &_v2)
are 3 vertices in counterclockwise order? in 2D
Definition: Algorithms.hh:222
Vec::value_type distPointTriangle(const Vec &_p, const Vec &_v0, const Vec &_v1, const Vec &_v2, Vec &_nearestPoint)
distance from point _p to triangle (_v0, _v1, _v2)
Definition: Algorithms.hh:326
Vec::value_type distPointLine(const Vec &_p, const Vec &_v0, const Vec &_v1, Vec *_min_v=0)
Compute distance from point to line segment.
Definition: Algorithms.hh:288
ValueT distPointPlane(const VectorT &_porigin, const VectorT &_pnormal, const VectorT &_point)
Checks the distance from a point to a plane.
Definition: Algorithms.cc:862
VectorT projectToEdge(const VectorT &_start, const VectorT &_end, const VectorT &_point)
Definition: Algorithms.cc:874
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
Scalar circumRadius(const VectorT< Scalar, 3 > &_v0, const VectorT< Scalar, 3 > &_v1, const VectorT< Scalar, 3 > &_v2, const VectorT< Scalar, 3 > &_v3)
return radius of circumcircle of tetrahedron (_v0,_v1,_v2,_v3)
Definition: Algorithms.hh:92
int isObtuse(const VectorT &_p0, const VectorT &_p1, const VectorT &_p2)
Definition: Algorithms.cc:1003
Vec::value_type distPointTriangleStable(const Vec &_p, const Vec &_v0, const Vec &_v1, const Vec &_v2, Vec &_nearestPoint)
distance from point _p to triangle (_v0, _v1, _v2)
Definition: Algorithms.hh:345
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
Namespace providing different geometric functions concerning angles.