OGS
AnalyticalGeometry.cpp
Go to the documentation of this file.
1
15#include "AnalyticalGeometry.h"
16
17#include <Eigen/Geometry>
18#include <algorithm>
19#include <cmath>
20#include <limits>
21
22#include "BaseLib/StringTools.h"
24#include "PointVec.h"
25#include "Polyline.h"
26#include "predicates.h"
27
29{
31 MathLib::Point3d const& c)
32{
33 return orient2d(const_cast<double*>(a.data()),
34 const_cast<double*>(b.data()),
35 const_cast<double*>(c.data()));
36}
37
39 MathLib::Point3d const& b,
40 MathLib::Point3d const& c)
41{
42 return orient2dfast(const_cast<double*>(a.data()),
43 const_cast<double*>(b.data()),
44 const_cast<double*>(c.data()));
45}
46} // namespace ExactPredicates
47
48namespace GeoLib
49{
51 MathLib::Point3d const& p1,
52 MathLib::Point3d const& p2)
53{
54 double const orientation = ExactPredicates::getOrientation2d(p0, p1, p2);
55 if (orientation > 0)
56 {
57 return CCW;
58 }
59 if (orientation < 0)
60 {
61 return CW;
62 }
63 return COLLINEAR;
64}
65
67 MathLib::Point3d const& p1,
68 MathLib::Point3d const& p2)
69{
70 double const orientation =
72 if (orientation > 0)
73 {
74 return CCW;
75 }
76 if (orientation < 0)
77 {
78 return CW;
79 }
80 return COLLINEAR;
81}
82
83bool parallel(Eigen::Vector3d v, Eigen::Vector3d w)
84{
85 const double eps(std::numeric_limits<double>::epsilon());
86 double const eps_squared = eps * eps;
87
88 // check degenerated cases
89 if (v.squaredNorm() < eps_squared)
90 {
91 return false;
92 }
93
94 if (w.squaredNorm() < eps_squared)
95 {
96 return false;
97 }
98
99 v.normalize();
100 w.normalize();
101
102 bool parallel(true);
103 if (std::abs(v[0] - w[0]) > eps)
104 {
105 parallel = false;
106 }
107 if (std::abs(v[1] - w[1]) > eps)
108 {
109 parallel = false;
110 }
111 if (std::abs(v[2] - w[2]) > eps)
112 {
113 parallel = false;
114 }
115
116 if (!parallel)
117 {
118 parallel = true;
119 // change sense of direction of v_normalised
120 v *= -1.0;
121 // check again
122 if (std::abs(v[0] - w[0]) > eps)
123 {
124 parallel = false;
125 }
126 if (std::abs(v[1] - w[1]) > eps)
127 {
128 parallel = false;
129 }
130 if (std::abs(v[2] - w[2]) > eps)
131 {
132 parallel = false;
133 }
134 }
135
136 return parallel;
137}
138
140 GeoLib::LineSegment const& s1,
141 GeoLib::Point& s)
142{
143 GeoLib::Point const& pa{s0.getBeginPoint()};
144 GeoLib::Point const& pb{s0.getEndPoint()};
145 GeoLib::Point const& pc{s1.getBeginPoint()};
146 GeoLib::Point const& pd{s1.getEndPoint()};
147
148 if (!isCoplanar(pa, pb, pc, pd))
149 {
150 return false;
151 }
152
153 auto const& a = s0.getBeginPoint().asEigenVector3d();
154 auto const& b = s0.getEndPoint().asEigenVector3d();
155 auto const& c = s1.getBeginPoint().asEigenVector3d();
156 auto const& d = s1.getEndPoint().asEigenVector3d();
157
158 Eigen::Vector3d const v = b - a;
159 Eigen::Vector3d const w = d - c;
160 Eigen::Vector3d const qp = c - a;
161 Eigen::Vector3d const pq = a - c;
162
163 double const eps = std::numeric_limits<double>::epsilon();
164 double const squared_eps = eps * eps;
165 // handle special cases here to avoid computing intersection numerical
166 if (qp.squaredNorm() < squared_eps || (d - a).squaredNorm() < squared_eps)
167 {
168 s = pa;
169 return true;
170 }
171 if ((c - b).squaredNorm() < squared_eps ||
172 (d - b).squaredNorm() < squared_eps)
173 {
174 s = pb;
175 return true;
176 }
177
178 auto isLineSegmentIntersectingAB =
179 [&v](Eigen::Vector3d const& ap, std::size_t i)
180 {
181 // check if p is located at v=(a,b): (ap = t*v, t in [0,1])
182 return 0.0 <= ap[i] / v[i] && ap[i] / v[i] <= 1.0;
183 };
184
185 if (parallel(v, w))
186 { // original line segments (a,b) and (c,d) are parallel
187 if (parallel(pq, v))
188 { // line segment (a,b) and (a,c) are also parallel
189 // Here it is already checked that the line segments (a,b) and (c,d)
190 // are parallel. At this point it is also known that the line
191 // segment (a,c) is also parallel to (a,b). In that case it is
192 // possible to express c as c(t) = a + t * (b-a) (analog for the
193 // point d). Since the evaluation of all three coordinate equations
194 // (x,y,z) have to lead to the same solution for the parameter t it
195 // is sufficient to evaluate t only once.
196
197 // Search id of coordinate with largest absolute value which is will
198 // be used in the subsequent computations. This prevents division by
199 // zero in case the line segments are parallel to one of the
200 // coordinate axis.
201 std::size_t i_max(std::abs(v[0]) <= std::abs(v[1]) ? 1 : 0);
202 i_max = std::abs(v[i_max]) <= std::abs(v[2]) ? 2 : i_max;
203 if (isLineSegmentIntersectingAB(qp, i_max))
204 {
205 s = pc;
206 return true;
207 }
208 Eigen::Vector3d const ad = d - a;
209 if (isLineSegmentIntersectingAB(ad, i_max))
210 {
211 s = pd;
212 return true;
213 }
214 return false;
215 }
216 return false;
217 }
218
219 // general case
220 const double sqr_len_v(v.squaredNorm());
221 const double sqr_len_w(w.squaredNorm());
222
223 Eigen::Matrix2d mat;
224 mat(0, 0) = sqr_len_v;
225 mat(0, 1) = -v.dot(w);
226 mat(1, 1) = sqr_len_w;
227 mat(1, 0) = mat(0, 1);
228
229 Eigen::Vector2d rhs{v.dot(qp), w.dot(pq)};
230
231 rhs = mat.partialPivLu().solve(rhs);
232
233 // no theory for the following tolerances, determined by testing
234 // lower tolerance: little bit smaller than zero
235 const double l(-1.0 * std::numeric_limits<float>::epsilon());
236 // upper tolerance a little bit greater than one
237 const double u(1.0 + std::numeric_limits<float>::epsilon());
238 if (rhs[0] < l || u < rhs[0] || rhs[1] < l || u < rhs[1])
239 {
240 return false;
241 }
242
243 // compute points along line segments with minimal distance
244 GeoLib::Point const p0(a[0] + rhs[0] * v[0], a[1] + rhs[0] * v[1],
245 a[2] + rhs[0] * v[2]);
246 GeoLib::Point const p1(c[0] + rhs[1] * w[0], c[1] + rhs[1] * w[1],
247 c[2] + rhs[1] * w[2]);
248
249 double const min_dist(std::sqrt(MathLib::sqrDist(p0, p1)));
250 double const min_seg_len(
251 std::min(std::sqrt(sqr_len_v), std::sqrt(sqr_len_w)));
252 if (min_dist < min_seg_len * 1e-6)
253 {
254 s[0] = 0.5 * (p0[0] + p1[0]);
255 s[1] = 0.5 * (p0[1] + p1[1]);
256 s[2] = 0.5 * (p0[2] + p1[2]);
257 return true;
258 }
259
260 return false;
261}
262
266 GeoLib::Point& intersection_pnt)
267{
268 std::size_t const n_segs(ply->getNumberOfSegments());
269 // Neighbouring segments always intersects at a common vertex. The algorithm
270 // checks for intersections of non-neighbouring segments.
271 for (seg_it0 = ply->begin(); seg_it0 != ply->end() - 2; ++seg_it0)
272 {
273 seg_it1 = seg_it0 + 2;
274 std::size_t const seg_num_0 = seg_it0.getSegmentNumber();
275 for (; seg_it1 != ply->end(); ++seg_it1)
276 {
277 // Do not check first and last segment, because they are
278 // neighboured.
279 if (!(seg_num_0 == 0 && seg_it1.getSegmentNumber() == n_segs - 1))
280 {
281 if (lineSegmentIntersect(*seg_it0, *seg_it1, intersection_pnt))
282 {
283 return true;
284 }
285 }
286 }
287 }
288 return false;
289}
290
291void rotatePoints(Eigen::Matrix3d const& rot_mat,
292 std::vector<GeoLib::Point*>& pnts)
293{
294 rotatePoints(rot_mat, pnts.begin(), pnts.end());
295}
296
297Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const& n)
298{
299 Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
300 // check if normal points already in the right direction
301 if (n[0] == 0 && n[1] == 0)
302 {
303 rot_mat(1, 1) = 1.0;
304
305 if (n[2] > 0)
306 {
307 // identity matrix
308 rot_mat(0, 0) = 1.0;
309 rot_mat(2, 2) = 1.0;
310 }
311 else
312 {
313 // rotate by pi about the y-axis
314 rot_mat(0, 0) = -1.0;
315 rot_mat(2, 2) = -1.0;
316 }
317
318 return rot_mat;
319 }
320
321 // sqrt (n_1^2 + n_3^2)
322 double const h0(std::sqrt(n[0] * n[0] + n[2] * n[2]));
323
324 // In case the x and z components of the normal are both zero the rotation
325 // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
326 // required. The angle is either pi/2 or 3/2*pi. Thus the components of
327 // rot_mat are as follows.
328 if (h0 < std::numeric_limits<double>::epsilon())
329 {
330 rot_mat(0, 0) = 1.0;
331 if (n[1] > 0)
332 {
333 rot_mat(1, 2) = -1.0;
334 rot_mat(2, 1) = 1.0;
335 }
336 else
337 {
338 rot_mat(1, 2) = 1.0;
339 rot_mat(2, 1) = -1.0;
340 }
341 return rot_mat;
342 }
343
344 double const h1(1 / n.norm());
345
346 // general case: calculate entries of rotation matrix
347 rot_mat(0, 0) = n[2] / h0;
348 rot_mat(0, 1) = 0;
349 rot_mat(0, 2) = -n[0] / h0;
350 rot_mat(1, 0) = -n[1] * n[0] / h0 * h1;
351 rot_mat(1, 1) = h0 * h1;
352 rot_mat(1, 2) = -n[1] * n[2] / h0 * h1;
353 rot_mat(2, 0) = n[0] * h1;
354 rot_mat(2, 1) = n[1] * h1;
355 rot_mat(2, 2) = n[2] * h1;
356
357 return rot_mat;
358}
359
360Eigen::Matrix3d rotatePointsToXY(std::vector<GeoLib::Point*>& pnts)
361{
362 return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
363}
364
365std::unique_ptr<GeoLib::Point> triangleLineIntersection(
366 MathLib::Point3d const& a, MathLib::Point3d const& b,
367 MathLib::Point3d const& c, MathLib::Point3d const& p,
368 MathLib::Point3d const& q)
369{
370 Eigen::Vector3d const pq = q.asEigenVector3d() - p.asEigenVector3d();
371 Eigen::Vector3d const pa = a.asEigenVector3d() - p.asEigenVector3d();
372 Eigen::Vector3d const pb = b.asEigenVector3d() - p.asEigenVector3d();
373 Eigen::Vector3d const pc = c.asEigenVector3d() - p.asEigenVector3d();
374
375 double u = pq.cross(pc).dot(pb);
376 if (u < 0)
377 {
378 return nullptr;
379 }
380 double v = pq.cross(pa).dot(pc);
381 if (v < 0)
382 {
383 return nullptr;
384 }
385 double w = pq.cross(pb).dot(pa);
386 if (w < 0)
387 {
388 return nullptr;
389 }
390
391 const double denom(1.0 / (u + v + w));
392 u *= denom;
393 v *= denom;
394 w *= denom;
395 return std::make_unique<GeoLib::Point>(u * a[0] + v * b[0] + w * c[0],
396 u * a[1] + v * b[1] + w * c[1],
397 u * a[2] + v * b[2] + w * c[2]);
398}
399
401 std::vector<GeoLib::Polyline*>& plys)
402{
403 auto computeSegmentIntersections =
404 [&pnt_vec](GeoLib::Polyline& poly0, GeoLib::Polyline& poly1)
405 {
406 for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
407 {
408 for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
409 {
410 GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
411 if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
412 {
413 std::size_t const id(
414 pnt_vec.push_back(new GeoLib::Point(s)));
415 poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
416 poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
417 }
418 }
419 }
420 };
421
422 for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
423 {
424 auto it1(it0);
425 ++it1;
426 for (; it1 != plys.end(); ++it1)
427 {
428 computeSegmentIntersections(*(*it0), *(*it1));
429 }
430 }
431}
432
433std::tuple<std::vector<GeoLib::Point*>, Eigen::Vector3d>
435{
436 // 1 copy all points
437 std::vector<GeoLib::Point*> polygon_points;
438 polygon_points.reserve(polygon_in.getNumberOfPoints());
439 for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
440 {
441 polygon_points.push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
442 }
443
444 // 2 rotate points
445 auto [plane_normal, d_polygon] = GeoLib::getNewellPlane(polygon_points);
446 Eigen::Matrix3d const rot_mat =
448 GeoLib::rotatePoints(rot_mat, polygon_points);
449
450 // 3 set z coord to zero
451 std::for_each(polygon_points.begin(), polygon_points.end(),
452 [](GeoLib::Point* p) { (*p)[2] = 0.0; });
453
454 return {polygon_points, plane_normal};
455}
456
457std::vector<MathLib::Point3d> lineSegmentIntersect2d(
458 GeoLib::LineSegment const& ab, GeoLib::LineSegment const& cd)
459{
460 GeoLib::Point const& a{ab.getBeginPoint()};
461 GeoLib::Point const& b{ab.getEndPoint()};
462 GeoLib::Point const& c{cd.getBeginPoint()};
463 GeoLib::Point const& d{cd.getEndPoint()};
464
465 double const orient_abc(getOrientation(a, b, c));
466 double const orient_abd(getOrientation(a, b, d));
467
468 // check if the segment (cd) lies on the left or on the right of (ab)
469 if ((orient_abc > 0 && orient_abd > 0) ||
470 (orient_abc < 0 && orient_abd < 0))
471 {
472 return std::vector<MathLib::Point3d>();
473 }
474
475 // check: (cd) and (ab) are on the same line
476 if (orient_abc == 0.0 && orient_abd == 0.0)
477 {
478 double const eps(std::numeric_limits<double>::epsilon());
479 if (MathLib::sqrDist2d(a, c) < eps && MathLib::sqrDist2d(b, d) < eps)
480 {
481 return {{a, b}};
482 }
483 if (MathLib::sqrDist2d(a, d) < eps && MathLib::sqrDist2d(b, c) < eps)
484 {
485 return {{a, b}};
486 }
487
488 // Since orient_ab and orient_abd vanish, a, b, c, d are on the same
489 // line and for this reason it is enough to check the x-component.
490 auto isPointOnSegment = [](double q, double p0, double p1)
491 {
492 double const t((q - p0) / (p1 - p0));
493 return 0 <= t && t <= 1;
494 };
495
496 // check if c in (ab)
497 if (isPointOnSegment(c[0], a[0], b[0]))
498 {
499 // check if a in (cd)
500 if (isPointOnSegment(a[0], c[0], d[0]))
501 {
502 return {{a, c}};
503 }
504 // check b == c
505 if (MathLib::sqrDist2d(b, c) < eps)
506 {
507 return {{b}};
508 }
509 // check if b in (cd)
510 if (isPointOnSegment(b[0], c[0], d[0]))
511 {
512 return {{b, c}};
513 }
514 // check d in (ab)
515 if (isPointOnSegment(d[0], a[0], b[0]))
516 {
517 return {{c, d}};
518 }
519 std::stringstream err;
520 err.precision(std::numeric_limits<double>::max_digits10);
521 err << ab << " x " << cd;
522 OGS_FATAL(
523 "The case of parallel line segments ({:s}) is not handled yet. "
524 "Aborting.",
525 err.str());
526 }
527
528 // check if d in (ab)
529 if (isPointOnSegment(d[0], a[0], b[0]))
530 {
531 // check if a in (cd)
532 if (isPointOnSegment(a[0], c[0], d[0]))
533 {
534 return {{a, d}};
535 }
536 // check if b==d
537 if (MathLib::sqrDist2d(b, d) < eps)
538 {
539 return {{b}};
540 }
541 // check if b in (cd)
542 if (isPointOnSegment(b[0], c[0], d[0]))
543 {
544 return {{b, d}};
545 }
546 // d in (ab), b not in (cd): check c in (ab)
547 if (isPointOnSegment(c[0], a[0], b[0]))
548 {
549 return {{c, d}};
550 }
551
552 std::stringstream err;
553 err.precision(std::numeric_limits<double>::max_digits10);
554 err << ab << " x " << cd;
555 OGS_FATAL(
556 "The case of parallel line segments ({:s}) is not handled yet. "
557 "Aborting.",
558 err.str());
559 }
560 return std::vector<MathLib::Point3d>();
561 }
562
563 // precondition: points a, b, c are collinear
564 // the function checks if the point c is onto the line segment (a,b)
565 auto isCollinearPointOntoLineSegment = [](MathLib::Point3d const& a,
566 MathLib::Point3d const& b,
567 MathLib::Point3d const& c)
568 {
569 if (b[0] - a[0] != 0)
570 {
571 double const t = (c[0] - a[0]) / (b[0] - a[0]);
572 return 0.0 <= t && t <= 1.0;
573 }
574 if (b[1] - a[1] != 0)
575 {
576 double const t = (c[1] - a[1]) / (b[1] - a[1]);
577 return 0.0 <= t && t <= 1.0;
578 }
579 if (b[2] - a[2] != 0)
580 {
581 double const t = (c[2] - a[2]) / (b[2] - a[2]);
582 return 0.0 <= t && t <= 1.0;
583 }
584 return false;
585 };
586
587 if (orient_abc == 0.0)
588 {
589 if (isCollinearPointOntoLineSegment(a, b, c))
590 {
591 return {{c}};
592 }
593 return std::vector<MathLib::Point3d>();
594 }
595
596 if (orient_abd == 0.0)
597 {
598 if (isCollinearPointOntoLineSegment(a, b, d))
599 {
600 return {{d}};
601 }
602 return std::vector<MathLib::Point3d>();
603 }
604
605 // check if the segment (ab) lies on the left or on the right of (cd)
606 double const orient_cda(getOrientation(c, d, a));
607 double const orient_cdb(getOrientation(c, d, b));
608 if ((orient_cda > 0 && orient_cdb > 0) ||
609 (orient_cda < 0 && orient_cdb < 0))
610 {
611 return std::vector<MathLib::Point3d>();
612 }
613
614 // at this point it is sure that there is an intersection and the system of
615 // linear equations will be invertible
616 // solve the two linear equations (b-a, c-d) (t, s)^T = (c-a) simultaneously
617 Eigen::Matrix2d mat;
618 mat(0, 0) = b[0] - a[0];
619 mat(0, 1) = c[0] - d[0];
620 mat(1, 0) = b[1] - a[1];
621 mat(1, 1) = c[1] - d[1];
622 Eigen::Vector2d rhs{c[0] - a[0], c[1] - a[1]};
623
624 rhs = mat.partialPivLu().solve(rhs);
625 if (0 <= rhs[1] && rhs[1] <= 1.0)
626 {
627 return {MathLib::Point3d{std::array<double, 3>{
628 {c[0] + rhs[1] * (d[0] - c[0]), c[1] + rhs[1] * (d[1] - c[1]),
629 c[2] + rhs[1] * (d[2] - c[2])}}}};
630 }
631 return std::vector<MathLib::Point3d>(); // parameter s not in the valid
632 // range
633}
634
635void sortSegments(MathLib::Point3d const& seg_beg_pnt,
636 std::vector<GeoLib::LineSegment>& sub_segments)
637{
638 double const eps(std::numeric_limits<double>::epsilon());
639
640 auto findNextSegment =
641 [&eps](MathLib::Point3d const& seg_beg_pnt,
642 std::vector<GeoLib::LineSegment>& sub_segments,
643 std::vector<GeoLib::LineSegment>::iterator& sub_seg_it)
644 {
645 if (sub_seg_it == sub_segments.end())
646 {
647 return;
648 }
649 // find appropriate segment for the given segment begin point
650 auto act_beg_seg_it = std::find_if(
651 sub_seg_it, sub_segments.end(),
652 [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg)
653 {
654 return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) <
655 eps ||
656 MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps;
657 });
658 if (act_beg_seg_it == sub_segments.end())
659 {
660 return;
661 }
662 // if necessary correct orientation of segment, i.e. swap beg and
663 // end
664 if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) <
665 MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint()))
666 {
667 std::swap(act_beg_seg_it->getBeginPoint(),
668 act_beg_seg_it->getEndPoint());
669 }
670 assert(sub_seg_it != sub_segments.end());
671 // exchange segments within the container
672 if (sub_seg_it != act_beg_seg_it)
673 {
674 std::swap(*sub_seg_it, *act_beg_seg_it);
675 }
676 };
677
678 // find start segment
679 auto seg_it = sub_segments.begin();
680 findNextSegment(seg_beg_pnt, sub_segments, seg_it);
681
682 while (seg_it != sub_segments.end())
683 {
684 MathLib::Point3d& new_seg_beg_pnt(seg_it->getEndPoint());
685 seg_it++;
686 if (seg_it != sub_segments.end())
687 {
688 findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
689 }
690 }
691}
692
693Eigen::Matrix3d compute2DRotationMatrixToX(Eigen::Vector3d const& v)
694{
695 Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
696 const double cos_theta = v[0];
697 const double sin_theta = v[1];
698 rot_mat(0, 0) = rot_mat(1, 1) = cos_theta;
699 rot_mat(0, 1) = sin_theta;
700 rot_mat(1, 0) = -sin_theta;
701 rot_mat(2, 2) = 1.0;
702 return rot_mat;
703}
704
705Eigen::Matrix3d compute3DRotationMatrixToX(Eigen::Vector3d const& v)
706{
707 // a vector on the plane
708 Eigen::Vector3d yy = Eigen::Vector3d::Zero();
709 auto const eps = std::numeric_limits<double>::epsilon();
710 if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
711 {
712 yy[2] = 1.0;
713 }
714 else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
715 {
716 yy[0] = 1.0;
717 }
718 else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
719 {
720 yy[1] = 1.0;
721 }
722 else
723 {
724 for (unsigned i = 0; i < 3; i++)
725 {
726 if (std::abs(v[i]) > 0.0)
727 {
728 yy[i] = -v[i];
729 break;
730 }
731 }
732 }
733 // z"_vec
734 Eigen::Vector3d const zz = v.cross(yy).normalized();
735 // y"_vec
736 yy = zz.cross(v).normalized();
737
738 Eigen::Matrix3d rot_mat;
739 rot_mat.row(0) = v;
740 rot_mat.row(1) = yy;
741 rot_mat.row(2) = zz;
742 return rot_mat;
743}
744
745} // end namespace GeoLib
Definition of analytical geometry functions.
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the PointVec class.
Definition of the PolyLine class.
Definition of string helper functions.
GeoLib::Point const & getBeginPoint() const
GeoLib::Point const & getEndPoint() const
This class manages pointers to Points in a std::vector along with a name. It also handles the deletio...
Definition PointVec.h:36
std::size_t push_back(Point *pnt)
Definition PointVec.cpp:133
std::size_t getSegmentNumber() const
Definition Polyline.cpp:387
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition Polyline.h:40
std::size_t getNumberOfSegments() const
Definition Polyline.cpp:114
std::size_t getNumberOfPoints() const
Definition Polyline.cpp:109
const Point * getPoint(std::size_t i) const
returns the i-th point contained in the polyline
Definition Polyline.cpp:179
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
Definition Polyline.cpp:55
SegmentIterator begin() const
Definition Polyline.h:188
SegmentIterator end() const
Definition Polyline.h:190
std::size_t size() const
Definition TemplateVec.h:99
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
const double * data() const
Definition Point3d.h:60
double getOrientation2d(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)
double getOrientation2dFast(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)
Orientation getOrientationFast(MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
Eigen::Matrix3d compute3DRotationMatrixToX(Eigen::Vector3d const &v)
bool lineSegmentsIntersect(const GeoLib::Polyline *ply, GeoLib::Polyline::SegmentIterator &seg_it0, GeoLib::Polyline::SegmentIterator &seg_it1, GeoLib::Point &intersection_pnt)
void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec, std::vector< GeoLib::Polyline * > &plys)
void rotatePoints(Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
bool parallel(Eigen::Vector3d v, Eigen::Vector3d w)
std::tuple< std::vector< GeoLib::Point * >, Eigen::Vector3d > rotatePolygonPointsToXY(GeoLib::Polygon const &polygon_in)
Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const &n)
Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)
void sortSegments(MathLib::Point3d const &seg_beg_pnt, std::vector< GeoLib::LineSegment > &sub_segments)
std::vector< MathLib::Point3d > lineSegmentIntersect2d(GeoLib::LineSegment const &ab, GeoLib::LineSegment const &cd)
std::unique_ptr< GeoLib::Point > triangleLineIntersection(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, MathLib::Point3d const &p, MathLib::Point3d const &q)
std::pair< Eigen::Vector3d, double > getNewellPlane(InputIterator pnts_begin, InputIterator pnts_end)
Eigen::Matrix3d compute2DRotationMatrixToX(Eigen::Vector3d const &v)
bool lineSegmentIntersect(GeoLib::LineSegment const &s0, GeoLib::LineSegment const &s1, GeoLib::Point &s)
Orientation getOrientation(MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
double sqrDist2d(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.h:123
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:26
Definitions of the predicates functions.
double orient2dfast(double *, double *, double *)
double orient2d(double *, double *, double *)