OGS 6.2.1-97-g73d1aeda3
AnalyticalGeometry.cpp
Go to the documentation of this file.
1 
15 #include "AnalyticalGeometry.h"
16 
17 #include <algorithm>
18 #include <cmath>
19 #include <limits>
20 
21 #include <logog/include/logog.hpp>
22 
23 #include <Eigen/Dense>
24 
25 #include "BaseLib/StringTools.h"
26 
27 #include "Polyline.h"
28 #include "PointVec.h"
29 
31 
32 extern double orient2d(double *, double *, double *);
33 extern double orient2dfast(double*, double*, double*);
34 
35 namespace ExactPredicates
36 {
38  MathLib::Point3d const& b, MathLib::Point3d const& c)
39 {
40  return orient2d(const_cast<double*>(a.getCoords()),
41  const_cast<double*>(b.getCoords()),
42  const_cast<double*>(c.getCoords()));
43 }
44 
46  MathLib::Point3d const& b,
47  MathLib::Point3d const& c)
48 {
49  return orient2dfast(const_cast<double*>(a.getCoords()),
50  const_cast<double*>(b.getCoords()),
51  const_cast<double*>(c.getCoords()));
52 }
53 } // namespace ExactPredicates
54 
55 namespace GeoLib
56 {
58  MathLib::Point3d const& p1,
59  MathLib::Point3d const& p2)
60 {
61  double const orientation = ExactPredicates::getOrientation2d(p0, p1, p2);
62  if (orientation > 0)
63  {
64  return CCW;
65  }
66  if (orientation < 0)
67  {
68  return CW;
69  }
70  return COLLINEAR;
71 }
72 
74  MathLib::Point3d const& p1,
75  MathLib::Point3d const& p2)
76 {
77  double const orientation =
79  if (orientation > 0)
80  {
81  return CCW;
82  }
83  if (orientation < 0)
84  {
85  return CW;
86  }
87  return COLLINEAR;
88 }
89 
91 {
92  const double eps(std::numeric_limits<double>::epsilon());
93 
94  // check degenerated cases
95  if (v.getLength() < eps)
96  {
97  return false;
98  }
99 
100  if (w.getLength() < eps)
101  {
102  return false;
103  }
104 
105  v.normalize();
106  w.normalize();
107 
108  bool parallel(true);
109  if (std::abs(v[0] - w[0]) > eps)
110  {
111  parallel = false;
112  }
113  if (std::abs(v[1] - w[1]) > eps)
114  {
115  parallel = false;
116  }
117  if (std::abs(v[2] - w[2]) > eps)
118  {
119  parallel = false;
120  }
121 
122  if (! parallel) {
123  parallel = true;
124  // change sense of direction of v_normalised
125  v *= -1.0;
126  // check again
127  if (std::abs(v[0] - w[0]) > eps)
128  {
129  parallel = false;
130  }
131  if (std::abs(v[1] - w[1]) > eps)
132  {
133  parallel = false;
134  }
135  if (std::abs(v[2] - w[2]) > eps)
136  {
137  parallel = false;
138  }
139  }
140 
141  return parallel;
142 }
143 
145  GeoLib::LineSegment const& s1,
146  GeoLib::Point& s)
147 {
148  GeoLib::Point const& a{s0.getBeginPoint()};
149  GeoLib::Point const& b{s0.getEndPoint()};
150  GeoLib::Point const& c{s1.getBeginPoint()};
151  GeoLib::Point const& d{s1.getEndPoint()};
152 
153  if (!isCoplanar(a, b, c, d))
154  {
155  return false;
156  }
157 
158  // handle special cases here to avoid computing intersection numerical
159  if (MathLib::sqrDist(a, c) < std::numeric_limits<double>::epsilon() ||
160  MathLib::sqrDist(a, d) < std::numeric_limits<double>::epsilon()) {
161  s = a;
162  return true;
163  }
164  if (MathLib::sqrDist(b, c) < std::numeric_limits<double>::epsilon() ||
165  MathLib::sqrDist(b, d) < std::numeric_limits<double>::epsilon()) {
166  s = b;
167  return true;
168  }
169 
170  MathLib::Vector3 const v(a, b);
171  MathLib::Vector3 const w(c, d);
172  MathLib::Vector3 const qp(a, c);
173  MathLib::Vector3 const pq(c, a);
174 
175  auto isLineSegmentIntersectingAB = [&v](MathLib::Vector3 const& ap,
176  std::size_t i)
177  {
178  // check if p is located at v=(a,b): (ap = t*v, t in [0,1])
179  return 0.0 <= ap[i] / v[i] && ap[i] / v[i] <= 1.0;
180  };
181 
182  if (parallel(v,w)) { // original line segments (a,b) and (c,d) are parallel
183  if (parallel(pq,v)) { // line segment (a,b) and (a,c) are also parallel
184  // Here it is already checked that the line segments (a,b) and (c,d)
185  // are parallel. At this point it is also known that the line
186  // segment (a,c) is also parallel to (a,b). In that case it is
187  // possible to express c as c(t) = a + t * (b-a) (analog for the
188  // point d). Since the evaluation of all three coordinate equations
189  // (x,y,z) have to lead to the same solution for the parameter t it
190  // is sufficient to evaluate t only once.
191 
192  // Search id of coordinate with largest absolute value which is will
193  // be used in the subsequent computations. This prevents division by
194  // zero in case the line segments are parallel to one of the
195  // coordinate axis.
196  std::size_t i_max(std::abs(v[0]) <= std::abs(v[1]) ? 1 : 0);
197  i_max = std::abs(v[i_max]) <= std::abs(v[2]) ? 2 : i_max;
198  if (isLineSegmentIntersectingAB(qp, i_max)) {
199  s = c;
200  return true;
201  }
202  MathLib::Vector3 const ad(a, d);
203  if (isLineSegmentIntersectingAB(ad, i_max)) {
204  s = d;
205  return true;
206  }
207  return false;
208  }
209  return false;
210  }
211 
212  // general case
213  const double sqr_len_v(v.getSqrLength());
214  const double sqr_len_w(w.getSqrLength());
215 
216  Eigen::Matrix2d mat;
217  mat(0,0) = sqr_len_v;
218  mat(0,1) = -1.0 * MathLib::scalarProduct(v,w);
219  mat(1,1) = sqr_len_w;
220  mat(1,0) = mat(0,1);
221 
222  Eigen::Vector2d rhs;
223  rhs << MathLib::scalarProduct(v, qp), MathLib::scalarProduct(w, pq);
224 
225  rhs = mat.partialPivLu().solve(rhs);
226 
227  // no theory for the following tolerances, determined by testing
228  // lower tolerance: little bit smaller than zero
229  const double l(-1.0*std::numeric_limits<float>::epsilon());
230  // upper tolerance a little bit greater than one
231  const double u(1.0+std::numeric_limits<float>::epsilon());
232  if (rhs[0] < l || u < rhs[0] || rhs[1] < l || u < rhs[1]) {
233  return false;
234  }
235 
236  // compute points along line segments with minimal distance
237  GeoLib::Point const p0(a[0]+rhs[0]*v[0], a[1]+rhs[0]*v[1], a[2]+rhs[0]*v[2]);
238  GeoLib::Point const p1(c[0]+rhs[1]*w[0], c[1]+rhs[1]*w[1], c[2]+rhs[1]*w[2]);
239 
240  double const min_dist(sqrt(MathLib::sqrDist(p0, p1)));
241  double const min_seg_len(std::min(sqrt(sqr_len_v), sqrt(sqr_len_w)));
242  if (min_dist < min_seg_len * 1e-6) {
243  s[0] = 0.5 * (p0[0] + p1[0]);
244  s[1] = 0.5 * (p0[1] + p1[1]);
245  s[2] = 0.5 * (p0[2] + p1[2]);
246  return true;
247  }
248 
249  return false;
250 }
251 
255  GeoLib::Point& intersection_pnt)
256 {
257  std::size_t const n_segs(ply->getNumberOfSegments());
258  // Neighbouring segments always intersects at a common vertex. The algorithm
259  // checks for intersections of non-neighbouring segments.
260  for (seg_it0 = ply->begin(); seg_it0 != ply->end() - 2; ++seg_it0)
261  {
262  seg_it1 = seg_it0+2;
263  std::size_t const seg_num_0 = seg_it0.getSegmentNumber();
264  for ( ; seg_it1 != ply->end(); ++seg_it1) {
265  // Do not check first and last segment, because they are
266  // neighboured.
267  if (!(seg_num_0 == 0 && seg_it1.getSegmentNumber() == n_segs - 1)) {
268  if (lineSegmentIntersect(*seg_it0, *seg_it1, intersection_pnt)) {
269  return true;
270  }
271  }
272  }
273  }
274  return false;
275 }
276 
278 {
279  // *** some frequently used terms ***
280  // n_1^2 + n_2^2
281  const double h0(plane_normal[0] * plane_normal[0] + plane_normal[1] * plane_normal[1]);
282  // 1 / sqrt (n_1^2 + n_2^2)
283  const double h1(1.0 / sqrt(h0));
284  // 1 / sqrt (n_1^2 + n_2^2 + n_3^2)
285  const double h2(1.0 / sqrt(h0 + plane_normal[2] * plane_normal[2]));
286 
287  // calc rotation matrix
288  rot_mat(0, 0) = plane_normal[1] * h1;
289  rot_mat(0, 1) = -plane_normal[0] * h1;
290  rot_mat(0, 2) = 0.0;
291  rot_mat(1, 0) = plane_normal[0] * h2;
292  rot_mat(1, 1) = plane_normal[1] * h2;
293  rot_mat(1, 2) = plane_normal[2] * h2;
294  rot_mat(2, 0) = plane_normal[0] * plane_normal[2] * h1 * h2;
295  rot_mat(2, 1) = plane_normal[1] * plane_normal[2] * h1 * h2;
296  rot_mat(2, 2) = -sqrt(h0) * h2;
297 }
298 
299 void rotatePoints(MathLib::DenseMatrix<double> const& rot_mat, std::vector<GeoLib::Point*> &pnts)
300 {
301  rotatePoints(rot_mat, pnts.begin(), pnts.end());
302 }
303 
304 MathLib::DenseMatrix<double> rotatePointsToXY(std::vector<GeoLib::Point*>& pnts)
305 {
306  return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
307 }
308 
309 std::unique_ptr<GeoLib::Point> triangleLineIntersection(
310  MathLib::Point3d const& a, MathLib::Point3d const& b,
311  MathLib::Point3d const& c, MathLib::Point3d const& p,
312  MathLib::Point3d const& q)
313 {
314  const MathLib::Vector3 pq(p, q);
315  const MathLib::Vector3 pa(p, a);
316  const MathLib::Vector3 pb(p, b);
317  const MathLib::Vector3 pc(p, c);
318 
319  double u (MathLib::scalarTriple(pq, pc, pb));
320  if (u < 0)
321  {
322  return nullptr;
323  }
324  double v (MathLib::scalarTriple(pq, pa, pc));
325  if (v < 0)
326  {
327  return nullptr;
328  }
329  double w (MathLib::scalarTriple(pq, pb, pa));
330  if (w < 0)
331  {
332  return nullptr;
333  }
334 
335  const double denom (1.0/(u+v+w));
336  u*=denom;
337  v*=denom;
338  w*=denom;
339  return std::make_unique<GeoLib::Point>(u * a[0] + v * b[0] + w * c[0],
340  u * a[1] + v * b[1] + w * c[1],
341  u * a[2] + v * b[2] + w * c[2]);
342 }
343 
345  std::vector<GeoLib::Polyline*> & plys)
346 {
347  auto computeSegmentIntersections = [&pnt_vec](GeoLib::Polyline& poly0,
348  GeoLib::Polyline& poly1)
349  {
350  for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it)
351  {
352  for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it)
353  {
354  GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
355  if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
356  {
357  std::size_t const id(
358  pnt_vec.push_back(new GeoLib::Point(s)));
359  poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
360  poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
361  }
362  }
363  }
364  };
365 
366  for (auto it0(plys.begin()); it0 != plys.end(); ++it0) {
367  auto it1(it0);
368  ++it1;
369  for (; it1 != plys.end(); ++it1) {
370  computeSegmentIntersections(*(*it0), *(*it1));
371  }
372  }
373 }
374 
376  MathLib::Vector3 & plane_normal)
377 {
378  // 1 copy all points
379  auto* polygon_pnts(new std::vector<GeoLib::Point*>);
380  for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
381  {
382  polygon_pnts->push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
383  }
384 
385  // 2 rotate points
386  double d_polygon (0.0);
387  GeoLib::getNewellPlane (*polygon_pnts, plane_normal, d_polygon);
388  MathLib::DenseMatrix<double> rot_mat(3,3);
389  GeoLib::computeRotationMatrixToXY(plane_normal, rot_mat);
390  GeoLib::rotatePoints(rot_mat, *polygon_pnts);
391 
392  // 3 set z coord to zero
393  std::for_each(polygon_pnts->begin(), polygon_pnts->end(),
394  [] (GeoLib::Point* p) { (*p)[2] = 0.0; }
395  );
396 
397  // 4 create new polygon
398  GeoLib::Polyline rot_polyline(*polygon_pnts);
399  for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
400  {
401  rot_polyline.addPoint(k);
402  }
403  rot_polyline.addPoint(0);
404  return GeoLib::Polygon(rot_polyline);
405 }
406 
407 std::vector<MathLib::Point3d> lineSegmentIntersect2d(
408  GeoLib::LineSegment const& ab, GeoLib::LineSegment const& cd)
409 {
410  GeoLib::Point const& a{ab.getBeginPoint()};
411  GeoLib::Point const& b{ab.getEndPoint()};
412  GeoLib::Point const& c{cd.getBeginPoint()};
413  GeoLib::Point const& d{cd.getEndPoint()};
414 
415  double const orient_abc(getOrientation(a, b, c));
416  double const orient_abd(getOrientation(a, b, d));
417 
418  // check if the segment (cd) lies on the left or on the right of (ab)
419  if ((orient_abc > 0 && orient_abd > 0) || (orient_abc < 0 && orient_abd < 0)) {
420  return std::vector<MathLib::Point3d>();
421  }
422 
423  // check: (cd) and (ab) are on the same line
424  if (orient_abc == 0.0 && orient_abd == 0.0) {
425  double const eps(std::numeric_limits<double>::epsilon());
426  if (MathLib::sqrDist2d(a, c) < eps && MathLib::sqrDist2d(b, d) < eps)
427  {
428  return {{a, b}};
429  }
430  if (MathLib::sqrDist2d(a, d) < eps && MathLib::sqrDist2d(b, c) < eps)
431  {
432  return {{a, b}};
433  }
434 
435  // Since orient_ab and orient_abd vanish, a, b, c, d are on the same
436  // line and for this reason it is enough to check the x-component.
437  auto isPointOnSegment = [](double q, double p0, double p1)
438  {
439  double const t((q - p0) / (p1 - p0));
440  return 0 <= t && t <= 1;
441  };
442 
443  // check if c in (ab)
444  if (isPointOnSegment(c[0], a[0], b[0])) {
445  // check if a in (cd)
446  if (isPointOnSegment(a[0], c[0], d[0])) {
447  return {{a, c}};
448  }
449  // check b == c
450  if (MathLib::sqrDist2d(b,c) < eps) {
451  return {{b}};
452  }
453  // check if b in (cd)
454  if (isPointOnSegment(b[0], c[0], d[0])) {
455  return {{b, c}};
456  }
457  // check d in (ab)
458  if (isPointOnSegment(d[0], a[0], b[0])) {
459  return {{c, d}};
460  }
461  std::stringstream err;
462  err.precision(std::numeric_limits<double>::digits10);
463  err << ab << " x " << cd;
464  OGS_FATAL(
465  "The case of parallel line segments (%s) is not handled yet. "
466  "Aborting.",
467  err.str().c_str());
468  }
469 
470  // check if d in (ab)
471  if (isPointOnSegment(d[0], a[0], b[0])) {
472  // check if a in (cd)
473  if (isPointOnSegment(a[0], c[0], d[0])) {
474  return {{a, d}};
475  }
476  // check if b==d
477  if (MathLib::sqrDist2d(b, d) < eps) {
478  return {{b}};
479  }
480  // check if b in (cd)
481  if (isPointOnSegment(b[0], c[0], d[0])) {
482  return {{b, d}};
483  }
484  // d in (ab), b not in (cd): check c in (ab)
485  if (isPointOnSegment(c[0], a[0], b[0])) {
486  return {{c, d}};
487  }
488 
489  std::stringstream err;
490  err.precision(std::numeric_limits<double>::digits10);
491  err << ab << " x " << cd;
492  OGS_FATAL(
493  "The case of parallel line segments (%s) "
494  "is not handled yet. Aborting.",
495  err.str().c_str());
496  }
497  return std::vector<MathLib::Point3d>();
498  }
499 
500  // precondition: points a, b, c are collinear
501  // the function checks if the point c is onto the line segment (a,b)
502  auto isCollinearPointOntoLineSegment = [](MathLib::Point3d const& a,
503  MathLib::Point3d const& b,
504  MathLib::Point3d const& c) {
505  if (b[0] - a[0] != 0)
506  {
507  double const t = (c[0] - a[0]) / (b[0] - a[0]);
508  return 0.0 <= t && t <= 1.0;
509  }
510  if (b[1] - a[1] != 0)
511  {
512  double const t = (c[1] - a[1]) / (b[1] - a[1]);
513  return 0.0 <= t && t <= 1.0;
514  }
515  if (b[2] - a[2] != 0)
516  {
517  double const t = (c[2] - a[2]) / (b[2] - a[2]);
518  return 0.0 <= t && t <= 1.0;
519  }
520  return false;
521  };
522 
523  if (orient_abc == 0.0) {
524  if (isCollinearPointOntoLineSegment(a, b, c))
525  {
526  return {{c}};
527  }
528  return std::vector<MathLib::Point3d>();
529  }
530 
531  if (orient_abd == 0.0) {
532  if (isCollinearPointOntoLineSegment(a, b, d))
533  {
534  return {{d}};
535  }
536  return std::vector<MathLib::Point3d>();
537  }
538 
539  // check if the segment (ab) lies on the left or on the right of (cd)
540  double const orient_cda(getOrientation(c, d, a));
541  double const orient_cdb(getOrientation(c, d, b));
542  if ((orient_cda > 0 && orient_cdb > 0) || (orient_cda < 0 && orient_cdb < 0)) {
543  return std::vector<MathLib::Point3d>();
544  }
545 
546  // at this point it is sure that there is an intersection and the system of
547  // linear equations will be invertible
548  // solve the two linear equations (b-a, c-d) (t, s)^T = (c-a) simultaneously
549  Eigen::Matrix2d mat;
550  mat(0,0) = b[0]-a[0];
551  mat(0,1) = c[0]-d[0];
552  mat(1,0) = b[1]-a[1];
553  mat(1,1) = c[1]-d[1];
554  Eigen::Vector2d rhs;
555  rhs << c[0] - a[0], c[1] - a[1];
556 
557  rhs = mat.partialPivLu().solve(rhs);
558  if (0 <= rhs[1] && rhs[1] <= 1.0) {
559  return { MathLib::Point3d{std::array<double,3>{{
560  c[0]+rhs[1]*(d[0]-c[0]), c[1]+rhs[1]*(d[1]-c[1]),
561  c[2]+rhs[1]*(d[2]-c[2])}} } };
562  }
563  return std::vector<MathLib::Point3d>(); // parameter s not in the valid
564  // range
565 }
566 
568  MathLib::Point3d const& seg_beg_pnt,
569  std::vector<GeoLib::LineSegment>& sub_segments)
570 {
571  double const eps(std::numeric_limits<double>::epsilon());
572 
573  auto findNextSegment = [&eps](
574  MathLib::Point3d const& seg_beg_pnt,
575  std::vector<GeoLib::LineSegment>& sub_segments,
576  std::vector<GeoLib::LineSegment>::iterator&
577  sub_seg_it) {
578  if (sub_seg_it == sub_segments.end())
579  {
580  return;
581  }
582  // find appropriate segment for the given segment begin point
583  auto act_beg_seg_it = std::find_if(
584  sub_seg_it, sub_segments.end(),
585  [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg)
586  {
587  return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) < eps ||
588  MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps;
589  });
590  if (act_beg_seg_it == sub_segments.end())
591  {
592  return;
593  }
594  // if necessary correct orientation of segment, i.e. swap beg and end
595  if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) <
596  MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint()))
597  {
598  std::swap(act_beg_seg_it->getBeginPoint(),
599  act_beg_seg_it->getEndPoint());
600  }
601  assert(sub_seg_it != sub_segments.end());
602  // exchange segments within the container
603  if (sub_seg_it != act_beg_seg_it)
604  {
605  std::swap(*sub_seg_it, *act_beg_seg_it);
606  }
607  };
608 
609  // find start segment
610  auto seg_it = sub_segments.begin();
611  findNextSegment(seg_beg_pnt, sub_segments, seg_it);
612 
613  while (seg_it != sub_segments.end())
614  {
615  MathLib::Point3d & new_seg_beg_pnt(seg_it->getEndPoint());
616  seg_it++;
617  if (seg_it != sub_segments.end())
618  {
619  findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
620  }
621  }
622 }
623 
624 } // end namespace GeoLib
double getOrientation2d(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)
double getLength() const
Returns the length.
Definition: Vector3.h:139
Definition of the PolyLine class.
void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec, std::vector< GeoLib::Polyline *> &plys)
MathLib::DenseMatrix< double > rotatePointsToXY(std::vector< GeoLib::Point *> &pnts)
double sqrDist(const double *p0, const double *p1)
Definition: MathTools.h:83
void rotatePoints(MathLib::DenseMatrix< double > const &rot_mat, std::vector< GeoLib::Point *> &pnts)
Definition of string helper functions.
std::vector< MathLib::Point3d > lineSegmentIntersect2d(GeoLib::LineSegment const &ab, GeoLib::LineSegment const &cd)
double getOrientation2dFast(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)
GeoLib::Point const & getBeginPoint() const
Definition: LineSegment.cpp:62
std::size_t size() const
Definition: TemplateVec.h:107
std::size_t getNumberOfSegments() const
Definition: Polyline.cpp:203
void rotatePoints(MathLib::DenseMatrix< double > const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
double orient2d(double *, double *, double *)
double scalarTriple(MathLib::Vector3 const &u, MathLib::Vector3 const &v, MathLib::Vector3 const &w)
Calculates the scalar triple (u x v) . w.
Definition: Vector3.cpp:16
double sqrDist2d(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition: Point3d.h:58
bool lineSegmentsIntersect(const GeoLib::Polyline *ply, GeoLib::Polyline::SegmentIterator &seg_it0, GeoLib::Polyline::SegmentIterator &seg_it1, GeoLib::Point &intersection_pnt)
Orientation getOrientation(MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
const Point * getPoint(std::size_t i) const
returns the i-th point contained in the polyline
Definition: Polyline.cpp:271
Definition of analytical geometry functions.
bool lineSegmentIntersect(GeoLib::LineSegment const &s0, GeoLib::LineSegment const &s1, GeoLib::Point &s)
double getSqrLength() const
Returns the squared length.
Definition: Vector3.h:133
void computeRotationMatrixToXZ(MathLib::Vector3 const &plane_normal, MathLib::DenseMatrix< double > &rot_mat)
std::size_t getNumberOfPoints() const
Definition: Polyline.cpp:198
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)
Definition of the GEOObjects class.
Definition: BaseItem.h:20
GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const &polygon_in, MathLib::Vector3 &plane_normal)
void sortSegments(MathLib::Point3d const &seg_beg_pnt, std::vector< GeoLib::LineSegment > &sub_segments)
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.
Orientation getOrientationFast(MathLib::Point3d const &p0, MathLib::Point3d const &p1, MathLib::Point3d const &p2)
static const double q
T scalarProduct(T const *const v0, T const *const v1)
Definition: MathTools.h:28
GeoLib::Point const & getEndPoint() const
Definition: LineSegment.cpp:72
const T * getCoords() const
Definition: TemplatePoint.h:77
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
Definition: Polyline.cpp:73
void computeRotationMatrixToXY(MathLib::Vector3 const &n, T_MATRIX &rot_mat)
void getNewellPlane(InputIterator pnts_begin, InputIterator pnts_end, MathLib::Vector3 &plane_normal, double &d)
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition: Polyline.h:50
SegmentIterator begin() const
Definition: Polyline.h:191
This class manages pointers to Points in a std::vector along with a name. It also handles the deletin...
Definition: PointVec.h:39
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
std::size_t push_back(Point *pnt)
Definition: PointVec.cpp:131
std::size_t getSegmentNumber() const
Definition: Polyline.cpp:501
Definition of the PointVec class.
bool parallel(MathLib::Vector3 v, MathLib::Vector3 w)
double orient2dfast(double *, double *, double *)
SegmentIterator end() const
Definition: Polyline.h:196