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