OGS
Polygon.cpp
Go to the documentation of this file.
1 
15 #include "Polygon.h"
16 
17 #include <boost/math/constants/constants.hpp>
18 
19 #include "AnalyticalGeometry.h"
20 #include "BaseLib/quicksort.h"
21 
22 namespace GeoLib
23 {
24 Polygon::Polygon(const Polyline& ply, bool init)
25  : Polyline(ply), _aabb(ply.getPointsVec(), ply._ply_pnt_ids)
26 {
27  if (init)
28  {
29  initialise();
30  }
31  _simple_polygon_list.push_back(this);
32 }
33 
34 Polygon::Polygon(Polygon const& other) : Polyline(other), _aabb(other._aabb)
35 {
36  _simple_polygon_list.push_back(this);
37  auto sub_polygon_it(other._simple_polygon_list.begin());
38  for (sub_polygon_it++; // the first entry is the polygon itself, skip the
39  // entry
40  sub_polygon_it != other._simple_polygon_list.end();
41  ++sub_polygon_it)
42  {
43  _simple_polygon_list.emplace_back(new Polygon(*(*sub_polygon_it)));
44  }
45 }
46 
48 {
49  // remove polygons from list
50  for (auto& polygon : _simple_polygon_list)
51  {
52  // the first entry of the list can be a pointer the object itself
53  if (polygon != this)
54  {
55  delete polygon;
56  }
57  }
58 }
59 
61 {
62  if (this->isClosed())
63  {
65  return true;
66  }
67  WARN("Polygon::initialise(): base polyline is not closed.");
68  return false;
69 }
70 
72 {
73  auto const& min_aabb_pnt(_aabb.getMinPoint());
74  auto const& max_aabb_pnt(_aabb.getMaxPoint());
75 
76  if (pnt[0] < min_aabb_pnt[0] || max_aabb_pnt[0] < pnt[0] ||
77  pnt[1] < min_aabb_pnt[1] || max_aabb_pnt[1] < pnt[1])
78  {
79  return false;
80  }
81 
82  if (_simple_polygon_list.size() == 1)
83  {
84  std::size_t n_intersections(0);
85  const std::size_t n_nodes(getNumberOfPoints() - 1);
86  for (std::size_t k(0); k < n_nodes; k++)
87  {
88  if (((*(getPoint(k)))[1] <= pnt[1] &&
89  pnt[1] <= (*(getPoint(k + 1)))[1]) ||
90  ((*(getPoint(k + 1)))[1] <= pnt[1] &&
91  pnt[1] <= (*(getPoint(k)))[1]))
92  {
93  switch (getEdgeType(k, pnt))
94  {
95  case EdgeType::TOUCHING:
96  return true;
97  case EdgeType::CROSSING:
98  n_intersections++;
99  break;
101  break;
102  default:
103  // do nothing
104  ;
105  }
106  }
107  }
108  if (n_intersections % 2 == 1)
109  {
110  return true;
111  }
112  }
113  else
114  {
115  for (auto it(_simple_polygon_list.begin()++);
116  it != _simple_polygon_list.end();
117  ++it)
118  {
119  if ((*it)->isPntInPolygon(pnt))
120  {
121  return true;
122  }
123  }
124  }
125  return false;
126 }
127 
128 std::vector<GeoLib::Point> Polygon::getAllIntersectionPoints(
129  GeoLib::LineSegment const& segment) const
130 {
131  std::vector<GeoLib::Point> intersections;
132  GeoLib::Point s;
133  for (auto&& seg_it : *this)
134  {
135  if (GeoLib::lineSegmentIntersect(seg_it, segment, s))
136  {
137  intersections.push_back(s);
138  }
139  }
140 
141  return intersections;
142 }
143 
145 {
146  std::vector<GeoLib::Point> s(getAllIntersectionPoints(segment));
147 
148  GeoLib::Point const& a{segment.getBeginPoint()};
149  GeoLib::Point const& b{segment.getEndPoint()};
150  // no intersections -> check if at least one point of segment is in polygon
151  if (s.empty())
152  {
153  return (isPntInPolygon(a));
154  }
155 
156  const double tol(std::numeric_limits<float>::epsilon());
157 
158  // one intersection, intersection in line segment end point
159  if (s.size() == 1)
160  {
161  const double sqr_dist_as(MathLib::sqrDist(a, s[0]));
162  if (sqr_dist_as < tol)
163  {
164  return (isPntInPolygon(b));
165  }
166 
167  const double sqr_dist_bs(MathLib::sqrDist(b, s[0]));
168  if (sqr_dist_bs < tol)
169  {
170  return (isPntInPolygon(a));
171  }
172  }
173 
174  // Sorting the intersection with respect to the distance to the point a.
175  // This induces a partition of the line segment into sub segments.
176  std::sort(s.begin(), s.end(),
177  [&a](GeoLib::Point const& p0, GeoLib::Point const& p1)
178  { return MathLib::sqrDist(a, p0) < MathLib::sqrDist(a, p1); });
179 
180  // remove sub segments with almost zero length
181  for (std::size_t k(0); k < s.size() - 1;)
182  {
183  if (MathLib::sqrDist(s[k], s[k + 1]) < tol)
184  {
185  s.erase(s.begin() + k + 1);
186  }
187  else
188  {
189  k++;
190  }
191  }
192 
193  // Check if all sub segments are within the polygon.
194  if (!isPntInPolygon(GeoLib::Point(0.5 * (a[0] + s[0][0]),
195  0.5 * (a[1] + s[0][1]),
196  0.5 * (a[2] + s[0][2]))))
197  {
198  return false;
199  }
200  const std::size_t n_sub_segs(s.size() - 1);
201  for (std::size_t k(0); k < n_sub_segs; k++)
202  {
203  if (!isPntInPolygon(GeoLib::Point(0.5 * (s[k][0] + s[k + 1][0]),
204  0.5 * (s[k][1] + s[k + 1][1]),
205  0.5 * (s[k][2] + s[k + 1][2]))))
206  {
207  return false;
208  }
209  }
210  return isPntInPolygon(GeoLib::Point(0.5 * (s[0][0] + b[0]),
211  0.5 * (s[0][1] + b[1]),
212  0.5 * (s[0][2] + b[2])));
213 }
214 
216 {
217  return std::all_of(ply.begin(), ply.end(),
218  [this](auto const& segment)
219  { return containsSegment(segment); });
220 }
221 
223 {
224  const std::size_t ply_size(ply.getNumberOfPoints());
225  // check points
226  for (std::size_t k(0); k < ply_size; k++)
227  {
228  if (isPntInPolygon(*(ply.getPoint(k))))
229  {
230  return true;
231  }
232  }
233 
234  auto polygon_segment_intersects_line = [&](auto const& polygon_seg)
235  {
236  GeoLib::Point s;
237  return std::any_of(ply.begin(), ply.end(),
238  [&polygon_seg, &s](auto const& polyline_seg) {
239  return GeoLib::lineSegmentIntersect(
240  polyline_seg, polygon_seg, s);
241  });
242  };
243 
244  return any_of(std::cbegin(*this), std::cend(*this),
245  polygon_segment_intersects_line);
246 }
247 
249  GeoLib::LineSegment const& seg, GeoLib::Point& intersection_pnt,
250  std::size_t& seg_num) const
251 {
252  if (_simple_polygon_list.size() == 1)
253  {
254  for (auto seg_it(begin() + seg_num); seg_it != end(); ++seg_it)
255  {
256  if (GeoLib::lineSegmentIntersect(*seg_it, seg, intersection_pnt))
257  {
258  seg_num = seg_it.getSegmentNumber();
259  return true;
260  }
261  }
262  }
263  else
264  {
265  for (auto polygon : _simple_polygon_list)
266  {
267  for (auto seg_it(polygon->begin()); seg_it != polygon->end();
268  ++seg_it)
269  {
270  if (GeoLib::lineSegmentIntersect(*seg_it, seg,
271  intersection_pnt))
272  {
273  seg_num = seg_it.getSegmentNumber();
274  return true;
275  }
276  }
277  }
278  }
279  return false;
280 }
281 
282 EdgeType Polygon::getEdgeType(std::size_t k, MathLib::Point3d const& pnt) const
283 {
284  switch (getLocationOfPoint(k, pnt))
285  {
286  case Location::LEFT:
287  {
288  const GeoLib::Point& v(*(getPoint(k)));
289  const GeoLib::Point& w(*(getPoint(k + 1)));
290  if (v[1] < pnt[1] && pnt[1] <= w[1])
291  {
292  return EdgeType::CROSSING;
293  }
294 
295  return EdgeType::INESSENTIAL;
296  }
297  case Location::RIGHT:
298  {
299  const GeoLib::Point& v(*(getPoint(k)));
300  const GeoLib::Point& w(*(getPoint(k + 1)));
301  if (w[1] < pnt[1] && pnt[1] <= v[1])
302  {
303  return EdgeType::CROSSING;
304  }
305 
306  return EdgeType::INESSENTIAL;
307  }
308  case Location::BETWEEN:
309  case Location::SOURCE:
311  return EdgeType::TOUCHING;
312  default:
313  return EdgeType::INESSENTIAL;
314  }
315 }
316 
318 {
319  // *** pre processing: rotate points to xy-plan
320  // *** copy points to vector - last point is identical to the first
321  std::size_t n_pnts(this->getNumberOfPoints() - 1);
322  std::vector<GeoLib::Point*> tmp_polygon_pnts;
323  for (std::size_t k(0); k < n_pnts; k++)
324  {
325  tmp_polygon_pnts.push_back(new GeoLib::Point(*(this->getPoint(k))));
326  }
327 
328  // rotate copied points into x-y-plane
329  GeoLib::rotatePointsToXY(tmp_polygon_pnts);
330 
331  for (auto& tmp_polygon_pnt : tmp_polygon_pnts)
332  {
333  (*tmp_polygon_pnt)[2] =
334  0.0; // should be -= d but there are numerical errors
335  }
336 
337  // *** get the left most upper point
338  std::size_t min_x_max_y_idx(0); // for orientation check
339  for (std::size_t k(0); k < n_pnts; k++)
340  {
341  if ((*(tmp_polygon_pnts[k]))[0] <=
342  (*(tmp_polygon_pnts[min_x_max_y_idx]))[0])
343  {
344  if ((*(tmp_polygon_pnts[k]))[0] <
345  (*(tmp_polygon_pnts[min_x_max_y_idx]))[0])
346  {
347  min_x_max_y_idx = k;
348  }
349  else if ((*(tmp_polygon_pnts[k]))[1] >
350  (*(tmp_polygon_pnts[min_x_max_y_idx]))[1])
351  {
352  min_x_max_y_idx = k;
353  }
354  }
355  }
356  // *** determine orientation
357  GeoLib::Orientation orient;
358  if (0 < min_x_max_y_idx && min_x_max_y_idx < n_pnts - 2)
359  {
360  orient = GeoLib::getOrientation(*tmp_polygon_pnts[min_x_max_y_idx - 1],
361  *tmp_polygon_pnts[min_x_max_y_idx],
362  *tmp_polygon_pnts[min_x_max_y_idx + 1]);
363  }
364  else
365  {
366  if (0 == min_x_max_y_idx)
367  {
368  orient = GeoLib::getOrientation(*tmp_polygon_pnts[n_pnts - 1],
369  *tmp_polygon_pnts[0],
370  *tmp_polygon_pnts[1]);
371  }
372  else
373  {
374  orient = GeoLib::getOrientation(*tmp_polygon_pnts[n_pnts - 2],
375  *tmp_polygon_pnts[n_pnts - 1],
376  *tmp_polygon_pnts[0]);
377  }
378  }
379 
380  if (orient != GeoLib::CCW)
381  {
382  // switch orientation
383  std::size_t tmp_n_pnts(n_pnts);
384  tmp_n_pnts++; // include last point of polygon (which is identical to
385  // the first)
386  for (std::size_t k(0); k < tmp_n_pnts / 2; k++)
387  {
388  std::swap(_ply_pnt_ids[k], _ply_pnt_ids[tmp_n_pnts - 1 - k]);
389  }
390  }
391 
392  for (std::size_t k(0); k < n_pnts; k++)
393  {
394  delete tmp_polygon_pnts[k];
395  }
396 }
397 
399  const std::list<Polygon*>::const_iterator& polygon_it)
400 {
401  GeoLib::Polyline::SegmentIterator seg_it0((*polygon_it)->begin());
402  GeoLib::Polyline::SegmentIterator seg_it1((*polygon_it)->begin());
403  GeoLib::Point intersection_pnt;
404  if (!GeoLib::lineSegmentsIntersect(*polygon_it, seg_it0, seg_it1,
405  intersection_pnt))
406  {
407  return;
408  }
409 
410  std::size_t idx0(seg_it0.getSegmentNumber());
411  std::size_t idx1(seg_it1.getSegmentNumber());
412  // adding intersection point to pnt_vec
413  std::size_t const intersection_pnt_id(_ply_pnts.size());
414  const_cast<std::vector<Point*>&>(_ply_pnts).push_back(
415  new GeoLib::Point(intersection_pnt));
416 
417  // split Polygon
418  if (idx0 > idx1)
419  {
420  std::swap(idx0, idx1);
421  }
422 
423  GeoLib::Polyline polyline0{(*polygon_it)->getPointsVec()};
424  for (std::size_t k(0); k <= idx0; k++)
425  {
426  polyline0.addPoint((*polygon_it)->getPointID(k));
427  }
428  polyline0.addPoint(intersection_pnt_id);
429  for (std::size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
430  {
431  polyline0.addPoint((*polygon_it)->getPointID(k));
432  }
433 
434  GeoLib::Polyline polyline1{(*polygon_it)->getPointsVec()};
435  polyline1.addPoint(intersection_pnt_id);
436  for (std::size_t k(idx0 + 1); k <= idx1; k++)
437  {
438  polyline1.addPoint((*polygon_it)->getPointID(k));
439  }
440  polyline1.addPoint(intersection_pnt_id);
441 
442  // remove the polygon except the first
443  if (*polygon_it != this)
444  {
445  delete *polygon_it;
446  }
447  // erase polygon_it and add two new polylines
448  auto polygon1_it = _simple_polygon_list.insert(
449  _simple_polygon_list.erase(polygon_it), new GeoLib::Polygon(polyline1));
450  auto polygon0_it = _simple_polygon_list.insert(
451  polygon1_it, new GeoLib::Polygon(polyline0));
452 
453  splitPolygonAtIntersection(polygon0_it);
454  splitPolygonAtIntersection(polygon1_it);
455 }
456 
458  const std::list<GeoLib::Polygon*>::iterator& polygon_it)
459 {
460  std::size_t const n((*polygon_it)->getNumberOfPoints() - 1);
461  std::vector<std::size_t> id_vec(n);
462  std::vector<std::size_t> perm(n);
463  for (std::size_t k(0); k < n; k++)
464  {
465  id_vec[k] = (*polygon_it)->getPointID(k);
466  perm[k] = k;
467  }
468 
469  BaseLib::quicksort(id_vec, 0, n, perm);
470 
471  for (std::size_t k(0); k < n - 1; k++)
472  {
473  if (id_vec[k] == id_vec[k + 1])
474  {
475  std::size_t idx0 = perm[k];
476  std::size_t idx1 = perm[k + 1];
477 
478  if (idx0 > idx1)
479  {
480  std::swap(idx0, idx1);
481  }
482 
483  // create two closed polylines
484  GeoLib::Polyline polyline0{*(*polygon_it)};
485  for (std::size_t j(0); j <= idx0; j++)
486  {
487  polyline0.addPoint((*polygon_it)->getPointID(j));
488  }
489  for (std::size_t j(idx1 + 1);
490  j < (*polygon_it)->getNumberOfPoints();
491  j++)
492  {
493  polyline0.addPoint((*polygon_it)->getPointID(j));
494  }
495 
496  GeoLib::Polyline polyline1{*(*polygon_it)};
497  for (std::size_t j(idx0); j <= idx1; j++)
498  {
499  polyline1.addPoint((*polygon_it)->getPointID(j));
500  }
501 
502  // remove the polygon except the first
503  if (*polygon_it != this)
504  {
505  delete *polygon_it;
506  }
507  // erase polygon_it and add two new polygons
508  auto polygon1_it = _simple_polygon_list.insert(
509  _simple_polygon_list.erase(polygon_it), new Polygon(polyline1));
510  auto polygon0_it = _simple_polygon_list.insert(
511  polygon1_it, new Polygon(polyline0));
512 
513  splitPolygonAtPoint(polygon0_it);
514  splitPolygonAtPoint(polygon1_it);
515 
516  return;
517  }
518  }
519 }
520 
521 bool operator==(Polygon const& lhs, Polygon const& rhs)
522 {
523  if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
524  {
525  return false;
526  }
527 
528  const std::size_t n(lhs.getNumberOfPoints());
529  const std::size_t start_pnt(lhs.getPointID(0));
530 
531  // search start point of first polygon in second polygon
532  bool nfound(true);
533  std::size_t k(0);
534  for (; k < n - 1 && nfound; k++)
535  {
536  if (start_pnt == rhs.getPointID(k))
537  {
538  nfound = false;
539  break;
540  }
541  }
542 
543  // case: start point not found in second polygon
544  if (nfound)
545  {
546  return false;
547  }
548 
549  // *** determine direction
550  // opposite direction
551  if (k == n - 2)
552  {
553  for (k = 1; k < n - 1; k++)
554  {
555  if (lhs.getPointID(k) != rhs.getPointID(n - 1 - k))
556  {
557  return false;
558  }
559  }
560  return true;
561  }
562 
563  // same direction - start point of first polygon at arbitrary position in
564  // second polygon
565  if (lhs.getPointID(1) == rhs.getPointID(k + 1))
566  {
567  std::size_t j(k + 2);
568  for (; j < n - 1; j++)
569  {
570  if (lhs.getPointID(j - k) != rhs.getPointID(j))
571  {
572  return false;
573  }
574  }
575  j = 0; // new start point at second polygon
576  for (; j < k + 1; j++)
577  {
578  if (lhs.getPointID(n - (k + 2) + j + 1) != rhs.getPointID(j))
579  {
580  return false;
581  }
582  }
583  return true;
584  }
585  // opposite direction with start point of first polygon at arbitrary
586  // position
587  // *** ATTENTION
588  WARN(
589  "operator==(Polygon const& lhs, Polygon const& rhs) - not tested case "
590  "(implementation is probably buggy) - please contact "
591  "thomas.fischer@ufz.de mentioning the problem.");
592  // in second polygon
593  if (lhs.getPointID(1) == rhs.getPointID(k - 1))
594  {
595  std::size_t j(k - 2);
596  for (; j > 0; j--)
597  {
598  if (lhs.getPointID(k - 2 - j) != rhs.getPointID(j))
599  {
600  return false;
601  }
602  }
603  // new start point at second polygon - the point n-1 of a polygon is
604  // equal to the first point of the polygon (for this reason: n-2)
605  j = n - 2;
606  for (; j > k - 1; j--)
607  {
608  if (lhs.getPointID(n - 2 + j + k - 2) != rhs.getPointID(j))
609  {
610  return false;
611  }
612  }
613  return true;
614  }
615  return false;
616 }
617 
618 std::list<Polygon*> const& Polygon::computeListOfSimplePolygons()
619 {
622 
623  for (auto polygon : _simple_polygon_list)
624  {
625  polygon->initialise();
626  }
627  return _simple_polygon_list;
628 }
629 
630 } // end namespace GeoLib
Definition of analytical geometry functions.
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of the Polygon class.
Eigen::Vector3d const & getMinPoint() const
Definition: AABB.h:174
Eigen::Vector3d const & getMaxPoint() const
Definition: AABB.h:181
GeoLib::Point const & getBeginPoint() const
Definition: LineSegment.cpp:66
GeoLib::Point const & getEndPoint() const
Definition: LineSegment.cpp:76
std::list< Polygon * > const & computeListOfSimplePolygons()
Definition: Polygon.cpp:618
bool isPolylineInPolygon(const Polyline &ply) const
Definition: Polygon.cpp:215
void splitPolygonAtPoint(const std::list< Polygon * >::iterator &polygon_it)
Definition: Polygon.cpp:457
std::vector< GeoLib::Point > getAllIntersectionPoints(GeoLib::LineSegment const &segment) const
Definition: Polygon.cpp:128
void splitPolygonAtIntersection(const std::list< Polygon * >::const_iterator &polygon_it)
Definition: Polygon.cpp:398
bool isPartOfPolylineInPolygon(const Polyline &ply) const
Definition: Polygon.cpp:222
bool isPntInPolygon(const MathLib::Point3d &pnt) const
Definition: Polygon.cpp:71
bool containsSegment(GeoLib::LineSegment const &segment) const
Definition: Polygon.cpp:144
~Polygon() override
Definition: Polygon.cpp:47
EdgeType getEdgeType(std::size_t k, MathLib::Point3d const &pnt) const
Definition: Polygon.cpp:282
std::list< Polygon * > _simple_polygon_list
Definition: Polygon.h:132
bool getNextIntersectionPointPolygonLine(GeoLib::LineSegment const &seg, GeoLib::Point &intersection_pnt, std::size_t &seg_num) const
Definition: Polygon.cpp:248
bool initialise()
Definition: Polygon.cpp:60
void ensureCCWOrientation()
Definition: Polygon.cpp:317
std::size_t getSegmentNumber() const
Definition: Polyline.cpp:422
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition: Polyline.h:51
friend class Polygon
Definition: Polyline.h:92
std::size_t getPointID(std::size_t i) const
Definition: Polyline.cpp:150
std::size_t getNumberOfPoints() const
Definition: Polyline.cpp:99
const Point * getPoint(std::size_t i) const
returns the i-th point contained in the polyline
Definition: Polyline.cpp:176
std::vector< std::size_t > _ply_pnt_ids
Definition: Polyline.h:227
SegmentIterator begin() const
Definition: Polyline.h:189
bool isClosed() const
Definition: Polyline.cpp:109
virtual bool addPoint(std::size_t pnt_id)
Definition: Polyline.cpp:34
const std::vector< Point * > & _ply_pnts
Definition: Polyline.h:225
Location getLocationOfPoint(std::size_t k, MathLib::Point3d const &pnt) const
Definition: Polyline.cpp:312
SegmentIterator end() const
Definition: Polyline.h:194
std::vector< Point * > const & getPointsVec() const
Definition: Polyline.cpp:182
EdgeType
Definition: Polygon.h:34
@ INESSENTIAL
INESSENTIAL.
@ CROSSING
CROSSING.
@ TOUCHING
TOUCHING.
void quicksort(It1 first1, It1 last1, It2 first2, Comparator compare)
Definition: quicksort.h:26
bool lineSegmentsIntersect(const GeoLib::Polyline *ply, GeoLib::Polyline::SegmentIterator &seg_it0, GeoLib::Polyline::SegmentIterator &seg_it1, GeoLib::Point &intersection_pnt)
Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)
bool operator==(LineSegment const &s0, LineSegment const &s1)
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 sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition: Point3d.h:48
Definition of the quicksort function.