18#include <boost/math/constants/constants.hpp>
60 long double const a[2] = {destination[0] - source[0],
61 destination[1] - source[1]};
62 long double const b[2] = {pnt[0] - source[0],
65 long double const det_2x2(a[0] * b[1] - a[1] * b[0]);
66 constexpr double eps = std::numeric_limits<double>::epsilon();
72 if (eps < std::abs(det_2x2))
76 if (a[0] * b[0] < 0.0 || a[1] * b[1] < 0.0)
80 if (a[0] * a[0] + a[1] * a[1] < b[0] * b[0] + b[1] * b[1])
111 if (a[1] < pnt[1] && pnt[1] <= b[1])
120 if (b[1] < pnt[1] && pnt[1] <= a[1])
137 :
Polyline(ply), _aabb(ply.getPointsVec(), ply.getPolylinePointIDs())
150 for (sub_polygon_it++;
179 WARN(
"Polygon::initialise(): base polyline is not closed.");
193 std::vector<GeoLib::Point> intersections;
195 for (
auto&& seg_it : polygon)
199 intersections.push_back(s);
203 return intersections;
210 if (pnt[0] < min_aabb_pnt[0] || max_aabb_pnt[0] < pnt[0] ||
211 pnt[1] < min_aabb_pnt[1] || max_aabb_pnt[1] < pnt[1])
218 std::size_t n_intersections(0);
220 for (std::size_t k(0); k < n_nodes; k++)
222 if (((*(
getPoint(k)))[1] <= pnt[1] &&
223 pnt[1] <= (*(
getPoint(k + 1)))[1]) ||
224 ((*(
getPoint(k + 1)))[1] <= pnt[1] &&
242 if (n_intersections % 2 == 1)
253 if ((*it)->isPntInPolygon(pnt))
274 const double tol(std::numeric_limits<float>::epsilon());
280 if (sqr_dist_as < tol)
286 if (sqr_dist_bs < tol)
294 std::sort(s.begin(), s.end(),
296 { return MathLib::sqrDist(a, p0) < MathLib::sqrDist(a, p1); });
299 for (std::size_t k(0); k < s.size() - 1;)
303 s.erase(s.begin() + k + 1);
313 0.5 * (a[1] + s[0][1]),
314 0.5 * (a[2] + s[0][2]))))
318 const std::size_t n_sub_segs(s.size() - 1);
319 for (std::size_t k(0); k < n_sub_segs; k++)
322 0.5 * (s[k][1] + s[k + 1][1]),
323 0.5 * (s[k][2] + s[k + 1][2]))))
329 0.5 * (s[0][1] + b[1]),
330 0.5 * (s[0][2] + b[2])));
335 return std::all_of(ply.
begin(), ply.
end(),
336 [
this](
auto const& segment)
337 { return containsSegment(segment); });
344 for (std::size_t k(0); k < ply_size; k++)
352 auto polygon_segment_intersects_line = [&](
auto const& polygon_seg)
355 return std::any_of(ply.
begin(), ply.
end(),
356 [&polygon_seg, &s](
auto const& polyline_seg) {
357 return GeoLib::lineSegmentIntersect(
358 polyline_seg, polygon_seg, s);
362 return std::any_of(std::cbegin(*
this), std::cend(*
this),
363 polygon_segment_intersects_line);
368 std::size_t& seg_num)
const
372 for (
auto seg_it(
begin() + seg_num); seg_it !=
end(); ++seg_it)
376 seg_num = seg_it.getSegmentNumber();
385 for (
auto seg_it(polygon->begin()); seg_it != polygon->end();
391 seg_num = seg_it.getSegmentNumber();
405 std::vector<GeoLib::Point*> tmp_polygon_pnts;
406 for (std::size_t k(0); k < n_pnts; k++)
414 for (
auto& tmp_polygon_pnt : tmp_polygon_pnts)
416 (*tmp_polygon_pnt)[2] =
421 std::size_t min_x_max_y_idx(0);
422 for (std::size_t k(0); k < n_pnts; k++)
424 if ((*(tmp_polygon_pnts[k]))[0] <=
425 (*(tmp_polygon_pnts[min_x_max_y_idx]))[0])
427 if ((*(tmp_polygon_pnts[k]))[0] <
428 (*(tmp_polygon_pnts[min_x_max_y_idx]))[0])
432 else if ((*(tmp_polygon_pnts[k]))[1] >
433 (*(tmp_polygon_pnts[min_x_max_y_idx]))[1])
441 if (0 < min_x_max_y_idx && min_x_max_y_idx < n_pnts - 2)
444 *tmp_polygon_pnts[min_x_max_y_idx],
445 *tmp_polygon_pnts[min_x_max_y_idx + 1]);
449 if (0 == min_x_max_y_idx)
452 *tmp_polygon_pnts[0],
453 *tmp_polygon_pnts[1]);
458 *tmp_polygon_pnts[n_pnts - 1],
459 *tmp_polygon_pnts[0]);
468 for (std::size_t k(0); k < n_pnts; k++)
470 delete tmp_polygon_pnts[k];
475 const std::list<Polygon*>::const_iterator& polygon_it)
489 std::size_t
const intersection_pnt_id(
_ply_pnts.size());
490 const_cast<std::vector<Point*>&
>(
_ply_pnts).push_back(
496 std::swap(idx0, idx1);
500 for (std::size_t k(0); k <= idx0; k++)
502 polyline0.addPoint((*polygon_it)->getPointID(k));
504 polyline0.addPoint(intersection_pnt_id);
505 for (std::size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
507 polyline0.addPoint((*polygon_it)->getPointID(k));
511 polyline1.addPoint(intersection_pnt_id);
512 for (std::size_t k(idx0 + 1); k <= idx1; k++)
514 polyline1.addPoint((*polygon_it)->getPointID(k));
516 polyline1.addPoint(intersection_pnt_id);
519 if (*polygon_it !=
this)
534 const std::list<GeoLib::Polygon*>::iterator& polygon_it)
536 std::size_t
const n((*polygon_it)->getNumberOfPoints() - 1);
537 std::vector<std::size_t> id_vec(n);
538 std::vector<std::size_t> perm(n);
539 for (std::size_t k(0); k < n; k++)
541 id_vec[k] = (*polygon_it)->getPointID(k);
547 for (std::size_t k(0); k < n - 1; k++)
549 if (id_vec[k] == id_vec[k + 1])
551 std::size_t idx0 = perm[k];
552 std::size_t idx1 = perm[k + 1];
556 std::swap(idx0, idx1);
561 for (std::size_t j(0); j <= idx0; j++)
563 polyline0.
addPoint((*polygon_it)->getPointID(j));
565 for (std::size_t j(idx1 + 1);
566 j < (*polygon_it)->getNumberOfPoints();
569 polyline0.addPoint((*polygon_it)->getPointID(j));
573 for (std::size_t j(idx0); j <= idx1; j++)
575 polyline1.
addPoint((*polygon_it)->getPointID(j));
579 if (*polygon_it !=
this)
587 polygon1_it,
new Polygon(polyline0));
605 const std::size_t start_pnt(lhs.
getPointID(0));
610 for (; k < n - 1 && nfound; k++)
629 for (k = 1; k < n - 1; k++)
643 std::size_t j(k + 2);
644 for (; j < n - 1; j++)
652 for (; j < k + 1; j++)
665 "operator==(Polygon const& lhs, Polygon const& rhs) - not tested case "
666 "(implementation is probably buggy) - please contact "
667 "thomas.fischer@ufz.de mentioning the problem.");
671 std::size_t j(k - 2);
682 for (; j > k - 1; j--)
701 polygon->initialise();
Definition of analytical geometry functions.
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the Polygon class.
MinMaxPoints getMinMaxPoints() const
GeoLib::Point const & getBeginPoint() const
GeoLib::Point const & getEndPoint() const
std::list< Polygon * > const & computeListOfSimplePolygons()
bool isPolylineInPolygon(const Polyline &ply) const
void splitPolygonAtPoint(const std::list< Polygon * >::iterator &polygon_it)
bool isPntInPolygon(MathLib::Point3d const &pnt) const
void splitPolygonAtIntersection(const std::list< Polygon * >::const_iterator &polygon_it)
bool isPartOfPolylineInPolygon(const Polyline &ply) const
bool containsSegment(GeoLib::LineSegment const &segment) const
std::list< Polygon * > _simple_polygon_list
bool getNextIntersectionPointPolygonLine(GeoLib::LineSegment const &seg, GeoLib::Point &intersection_pnt, std::size_t &seg_num) const
void ensureCCWOrientation()
std::size_t getSegmentNumber() const
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
std::size_t getPointID(std::size_t const i) const
std::size_t getNumberOfPoints() const
const Point * getPoint(std::size_t i) const
returns the i-th point contained in the polyline
void reverseOrientation()
SegmentIterator begin() const
virtual bool addPoint(std::size_t pnt_id)
const std::vector< Point * > & _ply_pnts
SegmentIterator end() const
std::vector< Point * > const & getPointsVec() const
void quicksort(It1 first1, It1 last1, It2 first2, Comparator compare)
bool lineSegmentsIntersect(const GeoLib::Polyline *ply, GeoLib::Polyline::SegmentIterator &seg_it0, GeoLib::Polyline::SegmentIterator &seg_it1, GeoLib::Point &intersection_pnt)
std::vector< GeoLib::Point > getAllIntersectionPoints(Polygon const &polygon, GeoLib::LineSegment const &segment)
Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)
Location getLocationOfPoint(MathLib::Point3d const &source, MathLib::Point3d const &destination, MathLib::Point3d const &pnt)
bool operator==(LineSegment const &s0, LineSegment const &s1)
EdgeType getEdgeType(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &pnt)
@ INESSENTIAL
INESSENTIAL.
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 of the quicksort function.