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