diff --git a/mapbox/geometry/wagyu/active_bound_list.hpp b/mapbox/geometry/wagyu/active_bound_list.hpp new file mode 100644 index 0000000..2f1c970 --- /dev/null +++ b/mapbox/geometry/wagyu/active_bound_list.hpp @@ -0,0 +1,526 @@ +#pragma once + +#ifdef DEBUG +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +using active_bound_list = std::list>; + +template +using active_bound_list_itr = typename active_bound_list::iterator; + +template +using active_bound_list_rev_itr = typename active_bound_list::reverse_iterator; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const active_bound_list& bnds) { + std::size_t c = 0; + for (auto const& bnd : bnds) { + out << "Index: " << c++ << std::endl; + out << *bnd; + } + return out; +} + +template +std::string output_edges(active_bound_list const& bnds) { + std::ostringstream out; + out << "["; + bool first = true; + for (auto const& bnd : bnds) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << bnd->current_edge->bot.x << "," << bnd->current_edge->bot.y << "],["; + out << bnd->current_edge->top.x << "," << bnd->current_edge->top.y << "]]"; + } + out << "]"; + return out.str(); +} + +#endif + +template +bool is_even_odd_fill_type(bound const& bound, + fill_type subject_fill_type, + fill_type clip_fill_type) { + if (bound.poly_type == polygon_type_subject) { + return subject_fill_type == fill_type_even_odd; + } else { + return clip_fill_type == fill_type_even_odd; + } +} + +template +bool is_even_odd_alt_fill_type(bound const& bound, + fill_type subject_fill_type, + fill_type clip_fill_type) { + if (bound.poly_type == polygon_type_subject) { + return clip_fill_type == fill_type_even_odd; + } else { + return subject_fill_type == fill_type_even_odd; + } +} + +template +inline bool bound2_inserts_before_bound1(bound const& bound1, bound const& bound2) { + if (values_are_equal(bound2.current_x, bound1.current_x)) { + if (bound2.current_edge->top.y > bound1.current_edge->top.y) { + return bound2.current_edge->top.x < + get_current_x(*(bound1.current_edge), bound2.current_edge->top.y); + } else { + return bound1.current_edge->top.x > + get_current_x(*(bound2.current_edge), bound1.current_edge->top.y); + } + } else { + return bound2.current_x < bound1.current_x; + } +} + +template +active_bound_list_itr insert_bound_into_ABL(bound& bnd, active_bound_list& active_bounds) { + auto itr = active_bounds.begin(); + while (itr != active_bounds.end() && !bound2_inserts_before_bound1(*(*itr), bnd)) { + ++itr; + } + return active_bounds.insert(itr, &bnd); +} + +template +active_bound_list_itr insert_bound_into_ABL(bound& bnd, + active_bound_list_itr itr, + active_bound_list& active_bounds) { + while (itr != active_bounds.end() && !bound2_inserts_before_bound1(*(*itr), bnd)) { + ++itr; + } + return active_bounds.insert(itr, &bnd); +} + +template +inline bool is_maxima(bound& bnd, T y) { + return bnd.next_edge == bnd.edges.end() && + bnd.current_edge->top.y == y; +} + +template +inline bool is_maxima(active_bound_list_itr& bnd, T y) { + return is_maxima(*(*bnd), y); +} + +template +inline bool is_intermediate(bound& bnd, T y) { + return bnd.next_edge != bnd.edges.end() && + bnd.current_edge->top.y == y; +} + +template +inline bool is_intermediate(active_bound_list_itr& bnd, T y) { + return is_intermediate(*(*bnd), y); +} + +template +inline bool current_edge_is_horizontal(active_bound_list_itr& bnd) { + return is_horizontal(*((*bnd)->current_edge)); +} + +template +inline bool next_edge_is_horizontal(active_bound_list_itr& bnd) { + return is_horizontal(*((*bnd)->next_edge)); +} + +template +inline void swap_positions_in_ABL(active_bound_list_itr& bnd1, + active_bound_list_itr& bnd2, + active_bound_list& active_bounds) { + if (std::next(bnd2) == bnd1) { + active_bounds.splice(bnd2, active_bounds, bnd1); + } else { + active_bounds.splice(bnd1, active_bounds, bnd2); + } +} + +template +void next_edge_in_bound(active_bound_list_itr& bnd, scanbeam_list& scanbeam) { + ++((*bnd)->current_edge); + if ((*bnd)->current_edge != (*bnd)->edges.end()) { + ++((*bnd)->next_edge); + (*bnd)->current_x = static_cast((*bnd)->current_edge->bot.x); + if (!current_edge_is_horizontal(bnd)) { + scanbeam.push((*bnd)->current_edge->top.y); + } + } +} + +template +active_bound_list_itr get_maxima_pair(active_bound_list_itr bnd, + active_bound_list& active_bounds) { + auto bnd_itr = active_bounds.begin(); + while (bnd_itr != active_bounds.end()) { + if (*bnd_itr == (*bnd)->maximum_bound) { + break; + } + ++bnd_itr; + } + return bnd_itr; +} + +template +void set_winding_count(active_bound_list_itr& bnd_itr, + active_bound_list& active_bounds, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + auto rev_bnd_itr = active_bound_list_rev_itr(bnd_itr); + if (rev_bnd_itr == active_bounds.rend()) { + if ((*bnd_itr)->winding_delta == 0) { + fill_type pft = ((*bnd_itr)->poly_type == polygon_type_subject) ? subject_fill_type + : clip_fill_type; + (*bnd_itr)->winding_count = (pft == fill_type_negative ? -1 : 1); + } else { + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + } + (*bnd_itr)->winding_count2 = 0; + return; + } + + // find the edge of the same polytype that immediately preceeds 'edge' in + // AEL + while (rev_bnd_itr != active_bounds.rend() && + ((*rev_bnd_itr)->poly_type != (*bnd_itr)->poly_type || + (*rev_bnd_itr)->winding_delta == 0)) { + ++rev_bnd_itr; + } + if (rev_bnd_itr == active_bounds.rend()) { + if ((*bnd_itr)->winding_delta == 0) { + fill_type pft = ((*bnd_itr)->poly_type == polygon_type_subject) ? subject_fill_type + : clip_fill_type; + (*bnd_itr)->winding_count = (pft == fill_type_negative ? -1 : 1); + } else { + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + } + (*bnd_itr)->winding_count2 = 0; + } else if ((*bnd_itr)->winding_delta == 0 && cliptype != clip_type_union) { + (*bnd_itr)->winding_count = 1; + (*bnd_itr)->winding_count2 = (*rev_bnd_itr)->winding_count2; + } else if (is_even_odd_fill_type(*(*bnd_itr), subject_fill_type, clip_fill_type)) { + // EvenOdd filling ... + if ((*bnd_itr)->winding_delta == 0) { + // are we inside a subj polygon ... + bool inside = true; + auto rev2 = std::next(rev_bnd_itr); + while (rev2 != active_bounds.rend()) { + if ((*rev2)->poly_type == (*rev_bnd_itr)->poly_type && + (*rev2)->winding_delta != 0) { + inside = !inside; + } + ++rev2; + } + (*bnd_itr)->winding_count = (inside ? 0 : 1); + } else { + (*bnd_itr)->winding_count = (*bnd_itr)->winding_delta; + } + (*bnd_itr)->winding_count2 = (*rev_bnd_itr)->winding_count2; + } else { + // nonZero, Positive or Negative filling ... + if ((*rev_bnd_itr)->winding_count * (*rev_bnd_itr)->winding_delta < 0) { + // prev edge is 'decreasing' WindCount (WC) toward zero + // so we're outside the previous polygon ... + if (std::abs((*rev_bnd_itr)->winding_count) > 1) { + // outside prev poly but still inside another. + // when reversing direction of prev poly use the same WC + if ((*rev_bnd_itr)->winding_delta * (*bnd_itr)->winding_delta < 0) { + (*bnd_itr)->winding_count = (*rev_bnd_itr)->winding_count; + } else { + // otherwise continue to 'decrease' WC ... + (*bnd_itr)->winding_count = + (*rev_bnd_itr)->winding_count + (*bnd_itr)->winding_delta; + } + } else { + // now outside all polys of same polytype so set own WC ... + (*bnd_itr)->winding_count = + ((*bnd_itr)->winding_delta == 0 ? 1 : (*bnd_itr)->winding_delta); + } + } else { + // prev edge is 'increasing' WindCount (WC) away from zero + // so we're inside the previous polygon ... + if ((*bnd_itr)->winding_delta == 0) { + (*bnd_itr)->winding_count = + ((*rev_bnd_itr)->winding_count < 0 ? (*rev_bnd_itr)->winding_count - 1 + : (*rev_bnd_itr)->winding_count + 1); + } else if ((*rev_bnd_itr)->winding_delta * (*bnd_itr)->winding_delta < 0) { + // if wind direction is reversing prev then use same WC + (*bnd_itr)->winding_count = (*rev_bnd_itr)->winding_count; + } else { + // otherwise add to WC ... + (*bnd_itr)->winding_count = + (*rev_bnd_itr)->winding_count + (*bnd_itr)->winding_delta; + } + } + (*bnd_itr)->winding_count2 = (*rev_bnd_itr)->winding_count2; + } + + // update winding_count2 ... + auto bnd_itr_forward = rev_bnd_itr.base(); + if (is_even_odd_alt_fill_type(*(*bnd_itr), subject_fill_type, clip_fill_type)) { + // EvenOdd filling ... + while (bnd_itr_forward != bnd_itr) { + if ((*bnd_itr_forward)->winding_delta != 0) { + (*bnd_itr)->winding_count2 = ((*bnd_itr)->winding_count2 == 0 ? 1 : 0); + } + ++bnd_itr_forward; + } + } else { + // nonZero, Positive or Negative filling ... + while (bnd_itr_forward != bnd_itr) { + (*bnd_itr)->winding_count2 += (*bnd_itr_forward)->winding_delta; + ++bnd_itr_forward; + } + } +} + +template +bool is_contributing(bound const& bnd, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + fill_type pft = subject_fill_type; + fill_type pft2 = clip_fill_type; + if (bnd.poly_type != polygon_type_subject) { + pft = clip_fill_type; + pft2 = subject_fill_type; + } + + switch (pft) { + case fill_type_even_odd: + // return false if a subj line has been flagged as inside a subj + // polygon + if (bnd.winding_delta == 0 && bnd.winding_count != 1) { + return false; + } + break; + case fill_type_non_zero: + if (std::abs(bnd.winding_count) != 1) { + return false; + } + break; + case fill_type_positive: + if (bnd.winding_count != 1) { + return false; + } + break; + case fill_type_negative: + default: + if (bnd.winding_count != -1) { + return false; + } + } + + switch (cliptype) { + case clip_type_intersection: + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 != 0); + case fill_type_positive: + return (bnd.winding_count2 > 0); + case fill_type_negative: + default: + return (bnd.winding_count2 < 0); + } + break; + case clip_type_union: + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 == 0); + case fill_type_positive: + return (bnd.winding_count2 <= 0); + case fill_type_negative: + default: + return (bnd.winding_count2 >= 0); + } + break; + case clip_type_difference: + if (bnd.poly_type == polygon_type_subject) { + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 == 0); + case fill_type_positive: + return (bnd.winding_count2 <= 0); + case fill_type_negative: + default: + return (bnd.winding_count2 >= 0); + } + } else { + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 != 0); + case fill_type_positive: + return (bnd.winding_count2 > 0); + case fill_type_negative: + default: + return (bnd.winding_count2 < 0); + } + } + break; + case clip_type_x_or: + if (bnd.winding_delta == 0) { + // XOr always contributing unless open + switch (pft2) { + case fill_type_even_odd: + case fill_type_non_zero: + return (bnd.winding_count2 == 0); + case fill_type_positive: + return (bnd.winding_count2 <= 0); + case fill_type_negative: + default: + return (bnd.winding_count2 >= 0); + } + } else { + return true; + } + break; + default: + return true; + } +} + +template +void insert_lm_only_one_bound(bound& bnd, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + auto abl_itr = insert_bound_into_ABL(bnd, active_bounds); + set_winding_count(abl_itr, active_bounds, cliptype, subject_fill_type, clip_fill_type); + if (is_contributing(bnd, cliptype, subject_fill_type, clip_fill_type)) { + add_first_point(abl_itr, active_bounds, (*abl_itr)->current_edge->bot, rings); + } + if (!current_edge_is_horizontal(abl_itr)) { + scanbeam.push((*abl_itr)->current_edge->top.y); + } +} + +template +void insert_lm_left_and_right_bound(bound& left_bound, + bound& right_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + // Both left and right bound + auto lb_abl_itr = insert_bound_into_ABL(left_bound, active_bounds); + auto rb_abl_itr = insert_bound_into_ABL(right_bound, lb_abl_itr, active_bounds); + set_winding_count(lb_abl_itr, active_bounds, cliptype, subject_fill_type, clip_fill_type); + (*rb_abl_itr)->winding_count = (*lb_abl_itr)->winding_count; + (*rb_abl_itr)->winding_count2 = (*lb_abl_itr)->winding_count2; + if (is_contributing(left_bound, cliptype, subject_fill_type, clip_fill_type)) { + add_local_minimum_point(lb_abl_itr, rb_abl_itr, active_bounds, + (*lb_abl_itr)->current_edge->bot, rings); + } + + // Add top of edges to scanbeam + scanbeam.push((*lb_abl_itr)->current_edge->top.y); + + if (!current_edge_is_horizontal(rb_abl_itr)) { + scanbeam.push((*rb_abl_itr)->current_edge->top.y); + } + auto abl_itr = std::next(lb_abl_itr); + while (abl_itr != rb_abl_itr && abl_itr != active_bounds.end()) { + // We call intersect_bounds here, but we do not swap positions in the ABL + // this is the logic that was copied from angus, but it might be correct + // to swap the positions in the ABL following this or at least move + // lb and rb to be next to each other in the ABL. + intersect_bounds(rb_abl_itr, abl_itr, (*lb_abl_itr)->current_edge->bot, cliptype, + subject_fill_type, clip_fill_type, rings, active_bounds); + ++abl_itr; + } +} + +template +void insert_local_minima_into_ABL(T const bot_y, + local_minimum_ptr_list const& minima_sorted, + local_minimum_ptr_list_itr& current_lm, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + while (current_lm != minima_sorted.end() && bot_y == (*current_lm)->y) { + initialize_lm(current_lm); + auto& left_bound = (*current_lm)->left_bound; + auto& right_bound = (*current_lm)->right_bound; + if (left_bound.edges.empty() && !right_bound.edges.empty()) { + insert_lm_only_one_bound(right_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } else if (right_bound.edges.empty() && !left_bound.edges.empty()) { + insert_lm_only_one_bound(left_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } else { + insert_lm_left_and_right_bound(left_bound, right_bound, active_bounds, rings, scanbeam, + cliptype, subject_fill_type, clip_fill_type); + } + ++current_lm; + } +} + +template +void insert_horizontal_local_minima_into_ABL(T const top_y, + local_minimum_ptr_list const& minima_sorted, + local_minimum_ptr_list_itr& current_lm, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + while (current_lm != minima_sorted.end() && top_y == (*current_lm)->y && + (*current_lm)->minimum_has_horizontal) { + initialize_lm(current_lm); + auto& left_bound = (*current_lm)->left_bound; + auto& right_bound = (*current_lm)->right_bound; + if (left_bound.edges.empty() && !right_bound.edges.empty()) { + insert_lm_only_one_bound(right_bound, active_bounds, rings, scanbeam, cliptype, + subject_fill_type, clip_fill_type); + } else if (right_bound.edges.empty() && !left_bound.edges.empty()) { + throw clipper_exception( + "There should only be horizontal local minimum on right bounds!"); + } else { + insert_lm_left_and_right_bound(left_bound, right_bound, active_bounds, rings, scanbeam, + cliptype, subject_fill_type, clip_fill_type); + } + ++current_lm; + } +} +} +} +} diff --git a/mapbox/geometry/wagyu/bound.hpp b/mapbox/geometry/wagyu/bound.hpp new file mode 100644 index 0000000..f0e3159 --- /dev/null +++ b/mapbox/geometry/wagyu/bound.hpp @@ -0,0 +1,95 @@ +#pragma once + +#include + +#include + +#include +#include +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct bound { + + edge_list edges; + edge_list_itr current_edge; + edge_list_itr next_edge; + mapbox::geometry::point last_point; + ring_ptr ring; + bound_ptr maximum_bound; // the bound who's maximum connects with this bound + double current_x; + std::size_t pos; + std::int32_t winding_count; + std::int32_t winding_count2; // winding count of the opposite polytype + std::int8_t winding_delta; // 1 or -1 depending on winding direction - 0 for linestrings + polygon_type poly_type; + edge_side side; // side only refers to current side of solution poly + + bound() noexcept + : edges(), + current_edge(edges.end()), + last_point({ 0, 0 }), + ring(nullptr), + maximum_bound(nullptr), + current_x(0.0), + pos(0), + winding_count(0), + winding_count2(0), + winding_delta(0), + poly_type(polygon_type_subject), + side(edge_left) { + } + + bound(bound && b) noexcept + : edges(std::move(b.edges)), + current_edge(std::move(b.current_edge)), + last_point(std::move(b.last_point)), + ring(std::move(b.ring)), + maximum_bound(std::move(b.maximum_bound)), + current_x(std::move(b.current_x)), + pos(std::move(b.pos)), + winding_count(std::move(b.winding_count)), + winding_count2(std::move(b.winding_count2)), + winding_delta(std::move(b.winding_delta)), + poly_type(std::move(b.poly_type)), + side(std::move(b.side)) { + } +}; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const bound& bnd) { + out << " Bound: " << &bnd << std::endl; + out << " current_x: " << bnd.current_x << std::endl; + out << " last_point: " << bnd.last_point.x << ", " << bnd.last_point.y << std::endl; + out << *(bnd.current_edge); + out << " winding count: " << bnd.winding_count << std::endl; + out << " winding_count2: " << bnd.winding_count2 << std::endl; + out << " winding_delta: " << static_cast(bnd.winding_delta) << std::endl; + out << " maximum_bound: " << bnd.maximum_bound << std::endl; + if (bnd.side == edge_left) { + out << " side: left" << std::endl; + } else { + out << " side: right" << std::endl; + } + out << " ring: " << bnd.ring << std::endl; + if (bnd.ring) { + out << " ring index: " << bnd.ring->ring_index << std::endl; + } + return out; +} + +#endif +} +} +} diff --git a/mapbox/geometry/wagyu/build_edges.hpp b/mapbox/geometry/wagyu/build_edges.hpp new file mode 100644 index 0000000..b7ce6d4 --- /dev/null +++ b/mapbox/geometry/wagyu/build_edges.hpp @@ -0,0 +1,216 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +bool build_edge_list(mapbox::geometry::line_string const& path_geometry, + edge_list& edges, + bool& is_flat) { + if (path_geometry.size() < 2) { + return false; + } + + auto itr_next = path_geometry.begin(); + ++itr_next; + auto itr = path_geometry.begin(); + while (itr_next != path_geometry.end()) { + if (*itr_next == *itr) { + // Duplicate point advance itr_next, but do not + // advance itr + ++itr_next; + continue; + } + + if (is_flat && itr_next->y != itr->y) { + is_flat = false; + } + edges.emplace_back(*itr, *itr_next); + itr = itr_next; + ++itr_next; + } + + if (edges.size() < 2) { + return false; + } + + return true; +} + +template +bool point_2_is_between_point_1_and_point_3(mapbox::geometry::point const& pt1, + mapbox::geometry::point const& pt2, + mapbox::geometry::point const& pt3) { + if ((pt1 == pt3) || (pt1 == pt2) || (pt3 == pt2)) { + return false; + } else if (pt1.x != pt3.x) { + return (pt2.x > pt1.x) == (pt2.x < pt3.x); + } else { + return (pt2.y > pt1.y) == (pt2.y < pt3.y); + } +} + +template +bool build_edge_list(mapbox::geometry::linear_ring const& path_geometry, edge_list& edges) { + using value_type = T; + + if (path_geometry.size() < 3) { + return false; + } + + // As this is a loop, we need to first go backwards from end to try and find + // the proper starting point for the iterators before the beginning + + auto itr_rev = path_geometry.rbegin(); + auto itr = path_geometry.begin(); + mapbox::geometry::point pt1 = *itr_rev; + mapbox::geometry::point pt2 = *itr; + + // Find next non repeated point going backwards from + // end for pt1 + while (pt1 == pt2) { + pt1 = *(++itr_rev); + if (itr_rev == path_geometry.rend()) { + return false; + } + } + ++itr; + mapbox::geometry::point pt3 = *itr; + auto itr_last = itr_rev.base(); + mapbox::geometry::point front_pt; + mapbox::geometry::point back_pt; + while (true) { + if (pt3 == pt2) { + // Duplicate point advance itr, but do not + // advance other points + if (itr == itr_last) { + break; + } + ++itr; + if (itr == itr_last) { + if (edges.empty()) { + break; + } + pt3 = front_pt; + } else { + pt3 = *itr; + } + continue; + } + + // Now check if slopes are equal between two segments - either + // a spike or a collinear point - if so drop point number 2. + if (slopes_equal(pt1, pt2, pt3)) { + // We need to reconsider previously added points + // because the point it was using was found to be collinear + // or a spike + pt2 = pt1; + if (!edges.empty()) { + edges.pop_back(); // remove previous edge (pt1) + } + if (!edges.empty()) { + if (back_pt == edges.back().top) { + pt1 = edges.back().bot; + } else { + pt1 = edges.back().top; + } + back_pt = pt1; + } else { + // If this occurs we must look to the back of the + // ring for new points. + do { + ++itr_rev; + if ((itr + 1) == itr_rev.base()) { + return false; + } + pt1 = *itr_rev; + } while (pt1 == pt2); + itr_last = itr_rev.base(); + } + continue; + } + + if (edges.empty()) { + front_pt = pt2; + } + edges.emplace_back(pt2, pt3); + back_pt = pt2; + if (itr == itr_last) { + break; + } + pt1 = pt2; + pt2 = pt3; + ++itr; + if (itr == itr_last) { + if (edges.empty()) { + break; + } + pt3 = front_pt; + } else { + pt3 = *itr; + } + } + + bool modified = false; + do { + modified = false; + if (edges.size() < 3) { + return false; + } + auto& f = edges.front(); + auto& b = edges.back(); + if (slopes_equal(f, b)) { + if (f.bot == b.top) { + if (f.top == b.bot) { + edges.pop_back(); + edges.erase(edges.begin()); + } else { + f.bot = b.bot; + edges.pop_back(); + } + modified = true; + } else if (f.top == b.bot) { + f.top = b.top; + edges.pop_back(); + modified = true; + } else if (f.top == b.top && f.bot == b.bot) { + edges.pop_back(); + edges.erase(edges.begin()); + modified = true; + } else if (f.top == b.top) { + if (point_2_is_between_point_1_and_point_3(f.top, f.bot, b.bot)) { + b.top = f.bot; + edges.erase(edges.begin()); + } else { + f.top = b.bot; + edges.pop_back(); + } + modified = true; + } else if (f.bot == b.bot) { + if (point_2_is_between_point_1_and_point_3(f.bot, f.top, b.top)) { + b.bot = f.top; + edges.erase(edges.begin()); + } else { + f.bot = b.top; + edges.pop_back(); + } + modified = true; + } + } + } while (modified); + + return true; +} +} +} +} diff --git a/mapbox/geometry/wagyu/build_local_minima_list.hpp b/mapbox/geometry/wagyu/build_local_minima_list.hpp new file mode 100644 index 0000000..595e7c4 --- /dev/null +++ b/mapbox/geometry/wagyu/build_local_minima_list.hpp @@ -0,0 +1,40 @@ +#pragma once + +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +bool add_line_string(mapbox::geometry::line_string const& path_geometry, + local_minimum_list& minima_list) { + bool is_flat = true; + edge_list new_edges; + new_edges.reserve(path_geometry.size()); + if (!build_edge_list(path_geometry, new_edges, is_flat) || new_edges.empty()) { + return false; + } + add_line_to_local_minima_list(new_edges, minima_list, polygon_type_subject); + return true; +} + +template +bool add_linear_ring(mapbox::geometry::linear_ring const& path_geometry, + local_minimum_list& minima_list, + polygon_type p_type) { + edge_list new_edges; + new_edges.reserve(path_geometry.size()); + if (!build_edge_list(path_geometry, new_edges) || new_edges.empty()) { + return false; + } + add_ring_to_local_minima_list(new_edges, minima_list, p_type); + return true; +} + +} +} +} diff --git a/mapbox/geometry/wagyu/build_result.hpp b/mapbox/geometry/wagyu/build_result.hpp new file mode 100644 index 0000000..5bc78ac --- /dev/null +++ b/mapbox/geometry/wagyu/build_result.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void push_ring_to_polygon(mapbox::geometry::polygon& poly, ring_ptr& r, bool reverse_output) { + mapbox::geometry::linear_ring lr; + lr.reserve(r->size + 1); + auto firstPt = r->points; + auto ptIt = r->points; + if (reverse_output) { + do { + lr.push_back({ ptIt->x, ptIt->y }); + ptIt = ptIt->next; + } while (ptIt != firstPt); + } else { + do { + lr.push_back({ ptIt->x, ptIt->y }); + ptIt = ptIt->prev; + } while (ptIt != firstPt); + } + lr.push_back({ firstPt->x, firstPt->y }); // close the ring + poly.push_back(lr); +} + +template +void build_result_polygons(std::vector>& solution, + ring_list& rings, + bool reverse_output) { + + for (auto& r : rings) { + assert(r->points); + std::size_t cnt = point_count(r->points); + if ((r->is_open && cnt < 2) || (!r->is_open && cnt < 3)) { + continue; + } + solution.emplace_back(); + push_ring_to_polygon(solution.back(), r, reverse_output); + for (auto& c : r->children) { + assert(c->points); + cnt = point_count(c->points); + if ((c->is_open && cnt < 2) || (!c->is_open && cnt < 3)) { + continue; + } + push_ring_to_polygon(solution.back(), c, reverse_output); + } + for (auto& c : r->children) { + if (!c->children.empty()) { + build_result_polygons(solution, c->children, reverse_output); + } + } + } +} + +template +void build_result(std::vector>& solution, + ring_manager& rings, + bool reverse_output) { + build_result_polygons(solution, rings.children, reverse_output); +} +} +} +} diff --git a/mapbox/geometry/wagyu/config.hpp b/mapbox/geometry/wagyu/config.hpp new file mode 100644 index 0000000..4dbdff7 --- /dev/null +++ b/mapbox/geometry/wagyu/config.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +enum clip_type : std::uint8_t { + clip_type_intersection = 0, + clip_type_union, + clip_type_difference, + clip_type_x_or +}; + +enum polygon_type : std::uint8_t { polygon_type_subject = 0, polygon_type_clip }; + +enum fill_type : std::uint8_t { + fill_type_even_odd = 0, + fill_type_non_zero, + fill_type_positive, + fill_type_negative +}; + +static double const def_arc_tolerance = 0.25; + +static int const EDGE_UNASSIGNED = -1; // edge not currently 'owning' a solution +static int const EDGE_SKIP = -2; // edge that would otherwise close a path +static std::int64_t const LOW_RANGE = 0x3FFFFFFF; +static std::int64_t const HIGH_RANGE = 0x3FFFFFFFFFFFFFFFLL; + +enum horizontal_direction : std::uint8_t { right_to_left = 0, left_to_right = 1 }; + +enum edge_side : std::uint8_t { edge_left = 0, edge_right }; + +enum join_type : std::uint8_t { join_type_square = 0, join_type_round, join_type_miter }; + +enum end_type { + end_type_closed_polygon = 0, + end_type_closed_line, + end_type_open_butt, + end_type_open_square, + end_type_open_round +}; + +template +using maxima_list = std::list; +} +} +} diff --git a/mapbox/geometry/wagyu/edge.hpp b/mapbox/geometry/wagyu/edge.hpp new file mode 100644 index 0000000..7f4f475 --- /dev/null +++ b/mapbox/geometry/wagyu/edge.hpp @@ -0,0 +1,122 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct bound; + +template +using bound_ptr = bound*; + +template +struct edge { + mapbox::geometry::point bot; + mapbox::geometry::point top; + double dx; + + edge(edge && e) noexcept + : bot(std::move(e.bot)), + top(std::move(e.top)), + dx(std::move(e.dx)) {} + + edge& operator=(edge && e) noexcept { + bot = std::move(e.bot); + top = std::move(e.top); + dx = std::move(e.dx); + return *this; + } + + edge(mapbox::geometry::point const& current, mapbox::geometry::point const& next_pt) noexcept + : bot(current), top(current), dx(0.0) { + if (current.y >= next_pt.y) { + top = next_pt; + } else { + bot = next_pt; + } + double dy = static_cast(top.y - bot.y); + if (value_is_zero(dy)) { + dx = std::numeric_limits::infinity(); + } else { + dx = static_cast(top.x - bot.x) / dy; + } + } + +}; + +template +using edge_ptr = edge*; + +template +using edge_list = std::vector>; + +template +using edge_list_itr = typename edge_list::iterator; + +template +bool slopes_equal(edge const& e1, edge const& e2) { + return (e1.top.y - e1.bot.y) * (e2.top.x - e2.bot.x) == + (e1.top.x - e1.bot.x) * (e2.top.y - e2.bot.y); +} + +template +inline bool is_horizontal(edge const& e) { + return std::isinf(e.dx); +} + +template +inline double get_current_x(edge const& edge, const T current_y) { + if (current_y == edge.top.y) { + return static_cast(edge.top.x); + } else { + return static_cast(edge.bot.x) + + edge.dx * static_cast(current_y - edge.bot.y); + } +} + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const edge& e) { + out << " Edge: " << std::endl; + out << " bot x: " << e.bot.x << " y: " << e.bot.y << std::endl; + out << " top x: " << e.top.x << " y: " << e.top.y << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + edge_list const& edges) { + out << "["; + bool first = true; + for (auto const& e : edges) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << e.bot.x << "," << e.bot.y << "],["; + out << e.top.x << "," << e.top.y << "]]"; + } + out << "]"; + return out; +} + +#endif +} +} +} diff --git a/mapbox/geometry/wagyu/exceptions.hpp b/mapbox/geometry/wagyu/exceptions.hpp new file mode 100644 index 0000000..f05b390 --- /dev/null +++ b/mapbox/geometry/wagyu/exceptions.hpp @@ -0,0 +1,23 @@ +#pragma once + +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { +class clipper_exception : public std::exception { +private: + std::string m_descr; + +public: + clipper_exception(const char* description) : m_descr(description) { + } + virtual ~clipper_exception() noexcept { + } + virtual const char* what() const noexcept { + return m_descr.c_str(); + } +}; +} +} +} diff --git a/mapbox/geometry/wagyu/intersect.hpp b/mapbox/geometry/wagyu/intersect.hpp new file mode 100644 index 0000000..cdb8571 --- /dev/null +++ b/mapbox/geometry/wagyu/intersect.hpp @@ -0,0 +1,73 @@ +#pragma once + +#include + +#include + +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct intersect_node { + + active_bound_list_itr bound1; + active_bound_list_itr bound2; + mapbox::geometry::point pt; + + intersect_node(intersect_node && n) + : bound1(std::move(n.bound1)), + bound2(std::move(n.bound2)), + pt(std::move(n.pt)) { } + + intersect_node& operator=(intersect_node && n) { + bound1 = std::move(n.bound1); + bound2 = std::move(n.bound2); + pt = std::move(n.pt); + return *this; + } + + intersect_node(active_bound_list_itr const& bound1_, + active_bound_list_itr const& bound2_, + mapbox::geometry::point const& pt_) + : bound1(bound1_), bound2(bound2_), pt(pt_) { + } +}; + +template +using intersect_list = std::vector>; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const intersect_node& e) { + out << " point x: " << e.pt.x << " y: " << e.pt.y << std::endl; + out << " bound 1: " << std::endl; + out << *(*e.bound1) << std::endl; + out << " bound 2: " << std::endl; + out << *(*e.bound2) << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const intersect_list& ints) { + std::size_t c = 0; + for (auto const& i : ints) { + out << "Intersection: " << c++ << std::endl; + out << i; + } + return out; +} + +#endif +} +} +} diff --git a/mapbox/geometry/wagyu/intersect_util.hpp b/mapbox/geometry/wagyu/intersect_util.hpp new file mode 100644 index 0000000..77be594 --- /dev/null +++ b/mapbox/geometry/wagyu/intersect_util.hpp @@ -0,0 +1,387 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct intersect_list_sorter { + inline bool operator()(intersect_node const& node1, intersect_node const& node2) { + if (!values_are_equal(node2.pt.y, node1.pt.y)) { + return node2.pt.y < node1.pt.y; + } else { + return ((*node2.bound1)->winding_count2 + (*node2.bound2)->winding_count2) > + ((*node1.bound1)->winding_count2 + (*node1.bound2)->winding_count2); + } + } +}; + +template +inline mapbox::geometry::point round_point(mapbox::geometry::point const& pt) { + return mapbox::geometry::point(round_towards_max(pt.x), round_towards_max(pt.y)); +} + +template +inline void swap_rings(bound& b1, bound& b2) { + ring_ptr ring = b1.ring; + b1.ring = b2.ring; + b2.ring = ring; +} + +template +inline void swap_sides(bound& b1, bound& b2) { + edge_side side = b1.side; + b1.side = b2.side; + b2.side = side; +} + +template +bool get_edge_intersection(edge const& e1, + edge const& e2, + mapbox::geometry::point& pt) { + T2 p0_x = static_cast(e1.bot.x); + T2 p0_y = static_cast(e1.bot.y); + T2 p1_x = static_cast(e1.top.x); + T2 p1_y = static_cast(e1.top.y); + T2 p2_x = static_cast(e2.bot.x); + T2 p2_y = static_cast(e2.bot.y); + T2 p3_x = static_cast(e2.top.x); + T2 p3_y = static_cast(e2.top.y); + T2 s1_x, s1_y, s2_x, s2_y; + s1_x = p1_x - p0_x; + s1_y = p1_y - p0_y; + s2_x = p3_x - p2_x; + s2_y = p3_y - p2_y; + + T2 s, t; + s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / (-s2_x * s1_y + s1_x * s2_y); + t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / (-s2_x * s1_y + s1_x * s2_y); + + if (s >= 0.0 && s <= 1.0 && t >= 0.0 && t <= 1.0) { + pt.x = p0_x + (t * s1_x); + pt.y = p0_y + (t * s1_y); + return true; + } + return false; +} + +template +void build_intersect_list(active_bound_list& active_bounds, + intersect_list& intersects) { + // bubblesort ... + bool isModified = false; + do { + isModified = false; + auto bnd = active_bounds.begin(); + auto bnd_next = std::next(bnd); + while (bnd_next != active_bounds.end()) { + if ((*bnd)->current_x > (*bnd_next)->current_x && + !slopes_equal(*((*bnd)->current_edge), *((*bnd_next)->current_edge))) { + mapbox::geometry::point pt; + if (!get_edge_intersection(*((*bnd)->current_edge), + *((*bnd_next)->current_edge), pt)) { + throw std::runtime_error( + "Trying to find intersection of lines that do not intersect"); + } + intersects.emplace_back(bnd, bnd_next, pt); + swap_positions_in_ABL(bnd, bnd_next, active_bounds); + bnd_next = std::next(bnd); + isModified = true; + } else { + bnd = bnd_next; + ++bnd_next; + } + } + } while (isModified); +} + +template +void intersect_bounds(active_bound_list_itr& b1, + active_bound_list_itr& b2, + mapbox::geometry::point const& pt, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings, + active_bound_list& active_bounds) { + bool b1Contributing = ((*b1)->ring != nullptr); + bool b2Contributing = ((*b2)->ring != nullptr); + // if either bound is on an OPEN path ... + if ((*b1)->winding_delta == 0 || (*b2)->winding_delta == 0) { + // ignore subject-subject open path intersections UNLESS they + // are both open paths, AND they are both 'contributing maximas' ... + if ((*b1)->winding_delta == 0 && (*b2)->winding_delta == 0) { + return; + } + + // if intersecting a subj line with a subj poly ... + else if ((*b1)->poly_type == (*b2)->poly_type && + (*b1)->winding_delta != (*b2)->winding_delta && cliptype == clip_type_union) { + if ((*b1)->winding_delta == 0) { + if (b2Contributing) { + add_point(b1, active_bounds, pt, rings); + if (b1Contributing) { + (*b1)->ring = nullptr; + } + } + } else { + if (b1Contributing) { + add_point(b2, active_bounds, pt, rings); + if (b2Contributing) { + (*b2)->ring = nullptr; + } + } + } + } else if ((*b1)->poly_type != (*b2)->poly_type) { + // toggle subj open path index on/off when std::abs(clip.WndCnt) == 1 + if (((*b1)->winding_delta == 0) && std::abs((*b2)->winding_count) == 1 && + (cliptype != clip_type_union || (*b2)->winding_count2 == 0)) { + add_point(b1, active_bounds, pt, rings); + if (b1Contributing) { + (*b1)->ring = nullptr; + } + } else if (((*b2)->winding_delta == 0) && (std::abs((*b1)->winding_count) == 1) && + (cliptype != clip_type_union || (*b1)->winding_count2 == 0)) { + add_point(b2, active_bounds, pt, rings); + if (b2Contributing) { + (*b2)->ring = nullptr; + } + } + } + return; + } + + // update winding counts... + // assumes that b1 will be to the Right of b2 ABOVE the intersection + if ((*b1)->poly_type == (*b2)->poly_type) { + if (is_even_odd_fill_type(*(*b1), subject_fill_type, clip_fill_type)) { + std::int32_t oldE1winding_count = (*b1)->winding_count; + (*b1)->winding_count = (*b2)->winding_count; + (*b2)->winding_count = oldE1winding_count; + } else { + if ((*b1)->winding_count + (*b2)->winding_delta == 0) { + (*b1)->winding_count = -(*b1)->winding_count; + } else { + (*b1)->winding_count += (*b2)->winding_delta; + } + if ((*b2)->winding_count - (*b1)->winding_delta == 0) { + (*b2)->winding_count = -(*b2)->winding_count; + } else { + (*b2)->winding_count -= (*b1)->winding_delta; + } + } + } else { + if (!is_even_odd_fill_type(*(*b2), subject_fill_type, clip_fill_type)) { + (*b1)->winding_count2 += (*b2)->winding_delta; + } else { + (*b1)->winding_count2 = ((*b1)->winding_count2 == 0) ? 1 : 0; + } + if (!is_even_odd_fill_type(*(*b1), subject_fill_type, clip_fill_type)) { + (*b2)->winding_count2 -= (*b1)->winding_delta; + } else { + (*b2)->winding_count2 = ((*b2)->winding_count2 == 0) ? 1 : 0; + } + } + + fill_type b1FillType, b2FillType, b1FillType2, b2FillType2; + if ((*b1)->poly_type == polygon_type_subject) { + b1FillType = subject_fill_type; + b1FillType2 = clip_fill_type; + } else { + b1FillType = clip_fill_type; + b1FillType2 = subject_fill_type; + } + if ((*b2)->poly_type == polygon_type_subject) { + b2FillType = subject_fill_type; + b2FillType2 = clip_fill_type; + } else { + b2FillType = clip_fill_type; + b2FillType2 = subject_fill_type; + } + + std::int32_t b1Wc, b2Wc; + switch (b1FillType) { + case fill_type_positive: + b1Wc = (*b1)->winding_count; + break; + case fill_type_negative: + b1Wc = -(*b1)->winding_count; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b1Wc = std::abs((*b1)->winding_count); + } + switch (b2FillType) { + case fill_type_positive: + b2Wc = (*b2)->winding_count; + break; + case fill_type_negative: + b2Wc = -(*b2)->winding_count; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b2Wc = std::abs((*b2)->winding_count); + } + if (b1Contributing && b2Contributing) { + if ((b1Wc != 0 && b1Wc != 1) || (b2Wc != 0 && b2Wc != 1) || + ((*b1)->poly_type != (*b2)->poly_type && cliptype != clip_type_x_or)) { + add_local_maximum_point(b1, b2, pt, rings, active_bounds); + } else { + add_point(b1, active_bounds, pt, rings); + add_point(b2, active_bounds, pt, rings); + swap_sides(*(*b1), *(*b2)); + swap_rings(*(*b1), *(*b2)); + } + } else if (b1Contributing) { + if (b2Wc == 0 || b2Wc == 1) { + add_point(b1, active_bounds, pt, rings); + (*b2)->last_point = pt; + swap_sides(*(*b1), *(*b2)); + swap_rings(*(*b1), *(*b2)); + } + } else if (b2Contributing) { + if (b1Wc == 0 || b1Wc == 1) { + (*b1)->last_point = pt; + add_point(b2, active_bounds, pt, rings); + swap_sides(*(*b1), *(*b2)); + swap_rings(*(*b1), *(*b2)); + } + } else if ((b1Wc == 0 || b1Wc == 1) && (b2Wc == 0 || b2Wc == 1)) { + // neither bound is currently contributing ... + + std::int32_t b1Wc2, b2Wc2; + switch (b1FillType2) { + case fill_type_positive: + b1Wc2 = (*b1)->winding_count2; + break; + case fill_type_negative: + b1Wc2 = -(*b1)->winding_count2; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b1Wc2 = std::abs((*b1)->winding_count2); + } + switch (b2FillType2) { + case fill_type_positive: + b2Wc2 = (*b2)->winding_count2; + break; + case fill_type_negative: + b2Wc2 = -(*b2)->winding_count2; + break; + case fill_type_even_odd: + case fill_type_non_zero: + default: + b2Wc2 = std::abs((*b2)->winding_count2); + } + + if ((*b1)->poly_type != (*b2)->poly_type) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } else if (b1Wc == 1 && b2Wc == 1) { + switch (cliptype) { + case clip_type_intersection: + if (b1Wc2 > 0 && b2Wc2 > 0) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + break; + default: + case clip_type_union: + if (b1Wc2 <= 0 && b2Wc2 <= 0) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + break; + case clip_type_difference: + if ((((*b1)->poly_type == polygon_type_clip) && (b1Wc2 > 0) && (b2Wc2 > 0)) || + (((*b1)->poly_type == polygon_type_subject) && (b1Wc2 <= 0) && (b2Wc2 <= 0))) { + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + break; + case clip_type_x_or: + add_local_minimum_point(b1, b2, active_bounds, pt, rings); + } + } else { + swap_sides(*(*b1), *(*b2)); + } + } +} + +template +bool bounds_adjacent(intersect_node const& inode) { + return (std::next(inode.bound1) == inode.bound2) || (std::next(inode.bound2) == inode.bound1); +} + +template +void process_intersect_list(intersect_list& intersects, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings, + active_bound_list& active_bounds) { + for (auto node_itr = intersects.begin(); node_itr != intersects.end(); ++node_itr) { + if (!bounds_adjacent(*node_itr)) { + auto next_itr = std::next(node_itr); + while (next_itr != intersects.end() && !bounds_adjacent(*next_itr)) { + ++next_itr; + } + if (next_itr == intersects.end()) { + throw std::runtime_error("Could not properly correct intersection order."); + } + std::iter_swap(node_itr, next_itr); + } + mapbox::geometry::point pt = round_point(node_itr->pt); + intersect_bounds(node_itr->bound1, node_itr->bound2, pt, cliptype, subject_fill_type, + clip_fill_type, rings, active_bounds); + swap_positions_in_ABL(node_itr->bound1, node_itr->bound2, active_bounds); + } +} + +template +void update_current_x(active_bound_list& active_bounds, T top_y) { + std::size_t pos = 0; + for (auto & bnd : active_bounds) { + bnd->pos = pos++; + bnd->current_x = get_current_x(*bnd->current_edge, top_y); + } +} + +template +void process_intersections(T top_y, + active_bound_list& active_bounds, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings) { + if (active_bounds.empty()) { + return; + } + update_current_x(active_bounds, top_y); + intersect_list intersects; + build_intersect_list(active_bounds, intersects); + + if (intersects.empty()) { + return; + } + + // Restore order of active bounds list + active_bounds.sort([] (bound_ptr const& b1, bound_ptr const& b2) { + return b1->pos < b2->pos; + }); + + // Sort the intersection list + std::stable_sort(intersects.begin(), intersects.end(), intersect_list_sorter()); + + process_intersect_list(intersects, cliptype, subject_fill_type, clip_fill_type, rings, + active_bounds); +} +} +} +} diff --git a/mapbox/geometry/wagyu/local_minimum.hpp b/mapbox/geometry/wagyu/local_minimum.hpp new file mode 100644 index 0000000..afd0a69 --- /dev/null +++ b/mapbox/geometry/wagyu/local_minimum.hpp @@ -0,0 +1,116 @@ +#pragma once + +#ifdef DEBUG +#include +#endif +#include + +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct local_minimum { + bound left_bound; + bound right_bound; + T y; + bool minimum_has_horizontal; + + local_minimum(bound&& left_bound_, bound&& right_bound_, T y_, bool has_horz_) + : left_bound(std::move(left_bound_)), + right_bound(std::move(right_bound_)), + y(y_), + minimum_has_horizontal(has_horz_) { + } +}; + +template +using local_minimum_list = std::deque>; + +template +using local_minimum_itr = typename local_minimum_list::iterator; + +template +using local_minimum_ptr = local_minimum*; + +template +using local_minimum_ptr_list = std::vector>; + +template +using local_minimum_ptr_list_itr = typename local_minimum_ptr_list::iterator; + +template +struct local_minimum_sorter { + inline bool operator()(local_minimum_ptr const& locMin1, + local_minimum_ptr const& locMin2) { + if (locMin2->y == locMin1->y) { + return locMin2->minimum_has_horizontal != locMin1->minimum_has_horizontal && + locMin1->minimum_has_horizontal; + } + return locMin2->y < locMin1->y; + } +}; + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const local_minimum& lm) { + out << " Local Minimum:" << std::endl; + out << " y: " << lm.y << std::endl; + if (lm.minimum_has_horizontal) { + out << " minimum_has_horizontal: true" << std::endl; + } else { + out << " minimum_has_horizontal: false" << std::endl; + } + out << " left_bound: " << std::endl; + out << lm.left_bound << std::endl; + out << " right_bound: " << std::endl; + out << lm.right_bound << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const local_minimum_ptr_list& lms) { + for (auto const& lm : lms) { + out << *lm; + } + return out; +} + +template +std::string output_all_edges(local_minimum_ptr_list const& lms) { + std::ostringstream out; + out << "["; + bool first = true; + for (auto const& lm : lms) { + for (auto const& e : lm->left_bound.edges) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << e.bot.x << "," << e.bot.y << "],["; + out << e.top.x << "," << e.top.y << "]]"; + } + for (auto const& e : lm->right_bound.edges) { + if (first) { + first = false; + } else { + out << ","; + } + out << "[[" << e.bot.x << "," << e.bot.y << "],["; + out << e.top.x << "," << e.top.y << "]]"; + } + } + out << "]"; + return out.str(); +} + +#endif +} +} +} diff --git a/mapbox/geometry/wagyu/local_minimum_util.hpp b/mapbox/geometry/wagyu/local_minimum_util.hpp new file mode 100644 index 0000000..b9863c8 --- /dev/null +++ b/mapbox/geometry/wagyu/local_minimum_util.hpp @@ -0,0 +1,419 @@ +#pragma once + +#include +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +inline void reverse_horizontal(edge& e) { + // swap horizontal edges' top and bottom x's so they follow the natural + // progression of the bounds - ie so their xbots will align with the + // adjoining lower edge. [Helpful in the process_horizontal() method.] + std::swap(e.top.x, e.bot.x); +} + +// Make a list start on a local maximum by +// shifting all the points not on a local maximum to the +template +void start_list_on_local_maximum(edge_list& edges) { + if (edges.size() <= 2) { + return; + } + // Find the first local maximum going forward in the list + auto prev_edge = edges.end(); + --prev_edge; + bool prev_edge_is_horizontal = is_horizontal(*prev_edge); + auto edge = edges.begin(); + bool edge_is_horizontal; + bool y_decreasing_before_last_horizontal = false; // assume false at start + + while (edge != edges.end()) { + edge_is_horizontal = is_horizontal(*edge); + if ((!prev_edge_is_horizontal && !edge_is_horizontal && edge->top == prev_edge->top)) { + break; + } + if (!edge_is_horizontal && prev_edge_is_horizontal) { + if (y_decreasing_before_last_horizontal && + (edge->top == prev_edge->bot || edge->top == prev_edge->top)) { + break; + } + } else if (!y_decreasing_before_last_horizontal && !prev_edge_is_horizontal && + edge_is_horizontal && + (prev_edge->top == edge->top || prev_edge->top == edge->bot)) { + y_decreasing_before_last_horizontal = true; + } + prev_edge_is_horizontal = edge_is_horizontal; + prev_edge = edge; + ++edge; + } + std::rotate(edges.begin(), edge, edges.end()); +} + +template +bound create_bound_towards_minimum(edge_list& edges) { + if (edges.size() == 1) { + if (is_horizontal(edges.front())) { + reverse_horizontal(edges.front()); + } + bound bnd; + std::swap(bnd.edges, edges); + return bnd; + } + auto next_edge = edges.begin(); + auto edge = next_edge; + ++next_edge; + bool edge_is_horizontal = is_horizontal(*edge); + if (edge_is_horizontal) { + reverse_horizontal(*edge); + } + bool next_edge_is_horizontal; + bool y_increasing_before_last_horizontal = false; // assume false at start + + while (next_edge != edges.end()) { + next_edge_is_horizontal = is_horizontal(*next_edge); + if ((!next_edge_is_horizontal && !edge_is_horizontal && edge->bot == next_edge->bot)) { + break; + } + if (!next_edge_is_horizontal && edge_is_horizontal) { + if (y_increasing_before_last_horizontal && + (next_edge->bot == edge->bot || next_edge->bot == edge->top)) { + break; + } + } else if (!y_increasing_before_last_horizontal && !edge_is_horizontal && + next_edge_is_horizontal && + (edge->bot == next_edge->top || edge->bot == next_edge->bot)) { + y_increasing_before_last_horizontal = true; + } + edge_is_horizontal = next_edge_is_horizontal; + edge = next_edge; + if (edge_is_horizontal) { + reverse_horizontal(*edge); + } + ++next_edge; + } + bound bnd; + if (next_edge == edges.end()) { + std::swap(edges, bnd.edges); + } else { + bnd.edges.reserve(std::distance(edges.begin(), next_edge)); + std::move(edges.begin(), next_edge, std::back_inserter(bnd.edges)); + edges.erase(edges.begin(), next_edge); + } + std::reverse(bnd.edges.begin(), bnd.edges.end()); + return bnd; +} + +template +bound create_bound_towards_maximum(edge_list& edges) { + if (edges.size() == 1) { + bound bnd; + std::swap(bnd.edges, edges); + return bnd; + } + auto next_edge = edges.begin(); + auto edge = next_edge; + ++next_edge; + bool edge_is_horizontal = is_horizontal(*edge); + bool next_edge_is_horizontal; + bool y_decreasing_before_last_horizontal = false; // assume false at start + + while (next_edge != edges.end()) { + next_edge_is_horizontal = is_horizontal(*next_edge); + if ((!next_edge_is_horizontal && !edge_is_horizontal && edge->top == next_edge->top)) { + break; + } + if (!next_edge_is_horizontal && edge_is_horizontal) { + if (y_decreasing_before_last_horizontal && + (next_edge->top == edge->bot || next_edge->top == edge->top)) { + break; + } + } else if (!y_decreasing_before_last_horizontal && !edge_is_horizontal && + next_edge_is_horizontal && + (edge->top == next_edge->top || edge->top == next_edge->bot)) { + y_decreasing_before_last_horizontal = true; + } + edge_is_horizontal = next_edge_is_horizontal; + edge = next_edge; + ++next_edge; + } + bound bnd; + if (next_edge == edges.end()) { + std::swap(bnd.edges, edges); + } else { + bnd.edges.reserve(std::distance(edges.begin(), next_edge)); + std::move(edges.begin(), next_edge, std::back_inserter(bnd.edges)); + edges.erase(edges.begin(), next_edge); + } + return bnd; +} + +template +void fix_horizontals(bound& bnd) { + + auto edge_itr = bnd.edges.begin(); + auto next_itr = std::next(edge_itr); + if (next_itr == bnd.edges.end()) { + return; + } + if (is_horizontal(*edge_itr) && next_itr->bot != edge_itr->top) { + reverse_horizontal(*edge_itr); + } + auto prev_itr = edge_itr++; + while (edge_itr != bnd.edges.end()) { + if (is_horizontal(*edge_itr) && prev_itr->top != edge_itr->bot) { + reverse_horizontal(*edge_itr); + } + prev_itr = edge_itr; + ++edge_itr; + } +} + +template +void move_horizontals_on_left_to_right(bound& left_bound, bound& right_bound) { + // We want all the horizontal segments that are at the same Y as the minimum to be on the right + // bound + auto edge_itr = left_bound.edges.begin(); + while (edge_itr != left_bound.edges.end()) { + if (!is_horizontal(*edge_itr)) { + break; + } + reverse_horizontal(*edge_itr); + ++edge_itr; + } + if (edge_itr == left_bound.edges.begin()) { + return; + } + std::reverse(left_bound.edges.begin(), edge_itr); + auto dist = std::distance(left_bound.edges.begin(), edge_itr); + std::move(left_bound.edges.begin(), edge_itr, std::back_inserter(right_bound.edges)); + left_bound.edges.erase(left_bound.edges.begin(), edge_itr); + std::rotate(right_bound.edges.begin(), std::prev(right_bound.edges.end(), dist), right_bound.edges.end()); +} + +template +void add_line_to_local_minima_list(edge_list& edges, local_minimum_list& minima_list) { + + if (edges.empty()) { + return; + } + // Adjust the order of the ring so we start on a local maximum + // therefore we start right away on a bound. + start_list_on_local_maximum(edges); + bound_ptr last_maximum = nullptr; + while (!edges.empty()) { + bool lm_minimum_has_horizontal = false; + auto to_minimum = create_bound_towards_minimum(edges); + assert(!to_minimum.edges.empty()); + fix_horizontals(to_minimum); + to_minimum.poly_type = polygon_type_subject; + to_minimum.maximum_bound = last_maximum; + to_minimum.winding_delta = 0; + auto to_min_first_non_horizontal = to_minimum.edges.begin(); + while (to_min_first_non_horizontal != to_minimum.edges.end() && + is_horizontal(*to_min_first_non_horizontal)) { + lm_minimum_has_horizontal = true; + ++to_min_first_non_horizontal; + } + if (edges.empty()) { + if (to_min_first_non_horizontal != to_minimum.edges.end() && + to_min_first_non_horizontal->dx > 0.0) { + to_minimum.side = edge_left; + bound right_bound; + right_bound.winding_delta = 0; + right_bound.side = edge_right; + right_bound.poly_type = polygon_type_subject; + move_horizontals_on_left_to_right(to_minimum, right_bound); + auto const& min_front = to_minimum.edges.front(); + minima_list.emplace_back(std::move(to_minimum), std::move(right_bound), min_front.y, + lm_minimum_has_horizontal); + if (last_maximum) { + last_maximum->maximum_bound = &(minima_list.back().left_bound); + last_maximum = nullptr; + } + } else { + to_minimum.side = edge_right; + bound left_bound; + left_bound.winding_delta = 0; + left_bound.side = edge_left; + left_bound.poly_type = polygon_type_subject; + auto const& min_front = to_minimum.edges.front(); + minima_list.emplace_back(std::move(left_bound), std::move(to_minimum), min_front.y); + if (last_maximum) { + last_maximum->maximum_bound = &(minima_list.back().right_bound); + last_maximum = nullptr; + } + } + break; + } + bool minimum_is_left = true; + auto to_maximum = create_bound_towards_maximum(edges); + assert(!to_maximum.edges.empty()); + fix_horizontals(to_maximum); + auto to_max_first_non_horizontal = to_minimum.edges.begin(); + while (to_max_first_non_horizontal != to_maximum.edges.end() && + is_horizontal(*to_max_first_non_horizontal)) { + lm_minimum_has_horizontal = true; + ++to_max_first_non_horizontal; + } + if (to_max_first_non_horizontal != to_maximum.edges.end() && + (to_min_first_non_horizontal == to_minimum.edges.end() || + to_max_first_non_horizontal->dx > to_min_first_non_horizontal->dx)) { + minimum_is_left = false; + move_horizontals_on_left_to_right(to_maximum, to_minimum); + } else { + minimum_is_left = true; + move_horizontals_on_left_to_right(to_minimum, to_maximum); + } + auto const& min_front = to_minimum.edges.front(); + to_maximum.poly_type = polygon_type_subject; + to_maximum.winding_delta = 0; + if (!minimum_is_left) { + to_minimum.side = edge_right; + to_maximum.side = edge_left; + minima_list.emplace_back(std::move(to_maximum), std::move(to_minimum), min_front.bot.y, + lm_minimum_has_horizontal); + if (last_maximum) { + last_maximum->maximum_bound = &(minima_list.back().right_bound); + } + last_maximum = &(minima_list.back().left_bound); + } else { + to_minimum.side = edge_left; + to_maximum.side = edge_right; + minima_list.emplace_back(std::move(to_minimum), std::move(to_maximum), min_front.bot.y, + lm_minimum_has_horizontal); + if (last_maximum) { + last_maximum->maximum_bound = &(minima_list.back().left_bound); + } + last_maximum = &(minima_list.back().right_bound); + } + } +} + +template +void add_ring_to_local_minima_list(edge_list& edges, + local_minimum_list& minima_list, + polygon_type poly_type) { + + if (edges.empty()) { + return; + } + // Adjust the order of the ring so we start on a local maximum + // therefore we start right away on a bound. + start_list_on_local_maximum(edges); + + bound_ptr first_minimum = nullptr; + bound_ptr last_maximum = nullptr; + while (!edges.empty()) { + bool lm_minimum_has_horizontal = false; + auto to_minimum = create_bound_towards_minimum(edges); + if (edges.empty()) { + throw std::runtime_error("Edges is empty after only creating a single bound."); + } + auto to_maximum = create_bound_towards_maximum(edges); + fix_horizontals(to_minimum); + fix_horizontals(to_maximum); + auto to_max_first_non_horizontal = to_maximum.edges.begin(); + auto to_min_first_non_horizontal = to_minimum.edges.begin(); + bool minimum_is_left = true; + while (to_max_first_non_horizontal != to_maximum.edges.end() && + is_horizontal(*to_max_first_non_horizontal)) { + lm_minimum_has_horizontal = true; + ++to_max_first_non_horizontal; + } + while (to_min_first_non_horizontal != to_minimum.edges.end() && + is_horizontal(*to_min_first_non_horizontal)) { + lm_minimum_has_horizontal = true; + ++to_min_first_non_horizontal; + } +#ifdef DEBUG + if (to_max_first_non_horizontal == to_maximum.edges.end() || + to_min_first_non_horizontal == to_minimum.edges.end()) { + throw clipper_exception("should not have a horizontal only bound for a ring"); + } +#endif + if (lm_minimum_has_horizontal) { + if (to_max_first_non_horizontal->bot.x > to_min_first_non_horizontal->bot.x) { + minimum_is_left = true; + move_horizontals_on_left_to_right(to_minimum, to_maximum); + } else { + minimum_is_left = false; + move_horizontals_on_left_to_right(to_maximum, to_minimum); + } + } else { + if (to_max_first_non_horizontal->dx > to_min_first_non_horizontal->dx) { + minimum_is_left = false; + } else { + minimum_is_left = true; + } + } + assert(!to_minimum.edges.empty()); + assert(!to_maximum.edges.empty()); + auto const& min_front = to_minimum.edges.front(); + if (last_maximum) { + to_minimum.maximum_bound = last_maximum; + } + to_minimum.poly_type = poly_type; + to_maximum.poly_type = poly_type; + if (!minimum_is_left) { + to_minimum.side = edge_right; + to_maximum.side = edge_left; + to_minimum.winding_delta = -1; + to_maximum.winding_delta = 1; + minima_list.emplace_back(std::move(to_maximum), std::move(to_minimum), min_front.bot.y, + lm_minimum_has_horizontal); + if (!last_maximum) { + first_minimum = &(minima_list.back().right_bound); + } else { + last_maximum->maximum_bound = &(minima_list.back().right_bound); + } + last_maximum = &(minima_list.back().left_bound); + } else { + to_minimum.side = edge_left; + to_maximum.side = edge_right; + to_minimum.winding_delta = -1; + to_maximum.winding_delta = 1; + minima_list.emplace_back(std::move(to_minimum), std::move(to_maximum), min_front.bot.y, + lm_minimum_has_horizontal); + if (!last_maximum) { + first_minimum = &(minima_list.back().left_bound); + } else { + last_maximum->maximum_bound = &(minima_list.back().left_bound); + } + last_maximum = &(minima_list.back().right_bound); + } + } + last_maximum->maximum_bound = first_minimum; + first_minimum->maximum_bound = last_maximum; +} + +template +void initialize_lm(local_minimum_ptr_list_itr& lm) { + if (!(*lm)->left_bound.edges.empty()) { + (*lm)->left_bound.current_edge = (*lm)->left_bound.edges.begin(); + (*lm)->left_bound.next_edge = std::next((*lm)->left_bound.current_edge); + (*lm)->left_bound.current_x = static_cast((*lm)->left_bound.current_edge->bot.x); + (*lm)->left_bound.winding_count = 0; + (*lm)->left_bound.winding_count2 = 0; + (*lm)->left_bound.side = edge_left; + (*lm)->left_bound.ring = nullptr; + } + if (!(*lm)->right_bound.edges.empty()) { + (*lm)->right_bound.current_edge = (*lm)->right_bound.edges.begin(); + (*lm)->right_bound.next_edge = std::next((*lm)->right_bound.current_edge); + (*lm)->right_bound.current_x = static_cast((*lm)->right_bound.current_edge->bot.x); + (*lm)->right_bound.winding_count = 0; + (*lm)->right_bound.winding_count2 = 0; + (*lm)->right_bound.side = edge_right; + (*lm)->right_bound.ring = nullptr; + } +} +} +} +} diff --git a/mapbox/geometry/wagyu/point.hpp b/mapbox/geometry/wagyu/point.hpp new file mode 100644 index 0000000..8811319 --- /dev/null +++ b/mapbox/geometry/wagyu/point.hpp @@ -0,0 +1,110 @@ +#pragma once + +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct point; + +template +using point_ptr = point*; + +template +using const_point_ptr = point* const; + +template +struct ring; + +template +using ring_ptr = ring*; + +template +using const_ring_ptr = ring* const; + +template +struct point { + using coordinate_type = T; + ring_ptr ring; + T x; + T y; + point_ptr next; + point_ptr prev; + + point(point && p) + : ring(std::move(p.ring)), + x(std::move(p.x)), + y(std::move(p.y)), + next(std::move(p.next)), + prev(std::move(p.prev)) { } + + point() : ring(nullptr), x(0), y(0), prev(this), next(this) { + } + point(T x_, T y_) : ring(nullptr), x(x_), y(y_), next(this), prev(this) { + } + point(ring_ptr ring_, mapbox::geometry::point const& pt) + : ring(ring_), x(pt.x), y(pt.y), next(this), prev(this) { + } + + point(ring_ptr ring_, mapbox::geometry::point const& pt, point_ptr before_this_point) + : ring(ring_), x(pt.x), y(pt.y), next(before_this_point), prev(before_this_point->prev) { + before_this_point->prev = this; + prev->next = this; + } +}; + +template +bool operator==(point const& lhs, point const& rhs) { + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +template +bool operator==(mapbox::geometry::point const& lhs, point const& rhs) { + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +template +bool operator==(point const& lhs, mapbox::geometry::point const& rhs) { + return lhs.x == rhs.x && lhs.y == rhs.y; +} + +template +bool operator!=(point const& lhs, point const& rhs) { + return lhs.x != rhs.x || lhs.y != rhs.y; +} + +template +bool operator!=(mapbox::geometry::point const& lhs, point const& rhs) { + return lhs.x != rhs.x || lhs.y != rhs.y; +} + +template +bool operator!=(point const& lhs, mapbox::geometry::point const& rhs) { + return lhs.x != rhs.x || lhs.y != rhs.y; +} + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const point& p) { + out << " point at: " << p.x << ", " << p.y; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const mapbox::geometry::point& p) { + out << " point at: " << p.x << ", " << p.y; + return out; +} +#endif +} +} +} diff --git a/mapbox/geometry/wagyu/process_horizontal.hpp b/mapbox/geometry/wagyu/process_horizontal.hpp new file mode 100644 index 0000000..ad8774f --- /dev/null +++ b/mapbox/geometry/wagyu/process_horizontal.hpp @@ -0,0 +1,285 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +active_bound_list_itr process_horizontal_left_to_right(T scanline_y, + active_bound_list_itr& horz_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + auto horizontal_itr_behind = horz_bound; + bool is_open = (*horz_bound)->winding_delta == 0; + bool is_maxima_edge = is_maxima(horz_bound, scanline_y); + auto bound_max_pair = active_bounds.end(); + if (is_maxima_edge) { + bound_max_pair = get_maxima_pair(horz_bound, active_bounds); + } + + auto hp_itr = rings.current_hp_itr; + while (hp_itr != rings.hot_pixels.end() && (hp_itr->y > scanline_y || (hp_itr->y == scanline_y && hp_itr->x < (*horz_bound)->current_edge->bot.x))) { + ++hp_itr; + } + + auto bnd = std::next(horz_bound); + + while (bnd != active_bounds.end()) { + // this code block inserts extra coords into horizontal edges (in output + // polygons) wherever hot pixels touch these horizontal edges. This helps + //'simplifying' polygons (ie if the Simplify property is set). + while (hp_itr != rings.hot_pixels.end() && hp_itr->y == scanline_y && hp_itr->x < std::llround((*bnd)->current_x) && + hp_itr->x < (*horz_bound)->current_edge->top.x) { + if ((*horz_bound)->ring && !is_open) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + } + ++hp_itr; + } + + if ((*bnd)->current_x > static_cast((*horz_bound)->current_edge->top.x)) { + break; + } + + // Also break if we've got to the end of an intermediate horizontal edge ... + // nb: Smaller Dx's are to the right of larger Dx's ABOVE the horizontal. + if (std::llround((*bnd)->current_x) == (*horz_bound)->current_edge->top.x && + (*horz_bound)->next_edge != (*horz_bound)->edges.end() && + (*horz_bound)->current_edge->dx < (*horz_bound)->next_edge->dx) { + break; + } + + // note: may be done multiple times + if ((*horz_bound)->ring && !is_open) { + add_point_to_ring(*(*horz_bound), + mapbox::geometry::point(std::llround((*bnd)->current_x), + scanline_y), + rings); + } + + // OK, so far we're still in range of the horizontal Edge but make sure + // we're at the last of consec. horizontals when matching with eMaxPair + if (is_maxima_edge && bnd == bound_max_pair) { + if ((*horz_bound)->ring) { + add_local_maximum_point(horz_bound, bound_max_pair, + (*horz_bound)->current_edge->top, rings, active_bounds); + } + active_bounds.erase(bound_max_pair); + auto after_horz = active_bounds.erase(horz_bound); + if (horizontal_itr_behind != horz_bound) { + return horizontal_itr_behind; + } else { + return after_horz; + } + } + + intersect_bounds(horz_bound, bnd, + mapbox::geometry::point(std::llround((*bnd)->current_x), + scanline_y), + cliptype, subject_fill_type, clip_fill_type, rings, active_bounds); + auto next_bnd = std::next(bnd); + swap_positions_in_ABL(horz_bound, bnd, active_bounds); + if (current_edge_is_horizontal(bnd) && horizontal_itr_behind == horz_bound) { + horizontal_itr_behind = bnd; + } + bnd = next_bnd; + } // end while (bnd != active_bounds.end()) + + if ((*horz_bound)->ring && !is_open) { + while (hp_itr != rings.hot_pixels.end() && hp_itr->y == scanline_y && + hp_itr->x < std::llround((*horz_bound)->current_edge->top.x)) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + ++hp_itr; + } + } + + if ((*horz_bound)->next_edge != (*horz_bound)->edges.end()) { + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), (*horz_bound)->current_edge->top, rings); + next_edge_in_bound(horz_bound, scanbeam); + + if ((*horz_bound)->winding_delta == 0) { + if (horizontal_itr_behind != horz_bound) { + return horizontal_itr_behind; + } else { + return std::next(horz_bound); + } + } + } else { + next_edge_in_bound(horz_bound, scanbeam); + } + if (horizontal_itr_behind != horz_bound) { + return horizontal_itr_behind; + } else { + return std::next(horz_bound); + } + } else { + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), (*horz_bound)->current_edge->top, rings); + } + auto after_horz = active_bounds.erase(horz_bound); + if (horizontal_itr_behind != horz_bound) { + return horizontal_itr_behind; + } else { + return after_horz; + } + } +} + +template +active_bound_list_itr process_horizontal_right_to_left(T scanline_y, + active_bound_list_itr& horz_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + bool is_open = (*horz_bound)->winding_delta == 0; + bool is_maxima_edge = is_maxima(horz_bound, scanline_y); + auto bound_max_pair = active_bounds.end(); + if (is_maxima_edge) { + bound_max_pair = get_maxima_pair(horz_bound, active_bounds); + } + auto hp_itr_fwd = rings.current_hp_itr; + while (hp_itr_fwd != rings.hot_pixels.end() && (hp_itr_fwd->y < scanline_y || (hp_itr_fwd->y == scanline_y && hp_itr_fwd->x < (*horz_bound)->current_edge->top.x))) { + ++hp_itr_fwd; + } + auto hp_itr = hot_pixel_rev_itr(hp_itr_fwd); + + auto bnd = active_bound_list_rev_itr(horz_bound); + while (bnd != active_bounds.rend()) { + // this code block inserts extra coords into horizontal edges (in output + // polygons) wherever hot pixels touch these horizontal edges. + while (hp_itr != rings.hot_pixels.rend() && hp_itr->y == scanline_y && hp_itr->x > std::llround((*bnd)->current_x) && + hp_itr->x > (*horz_bound)->current_edge->top.x) { + if ((*horz_bound)->ring && !is_open) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + } + ++hp_itr; + } + + if ((*bnd)->current_x < static_cast((*horz_bound)->current_edge->top.x)) { + break; + } + + // Also break if we've got to the end of an intermediate horizontal edge ... + // nb: Smaller Dx's are to the right of larger Dx's ABOVE the horizontal. + if (std::llround((*bnd)->current_x) == (*horz_bound)->current_edge->top.x && + (*horz_bound)->next_edge != (*horz_bound)->edges.end() && + (*horz_bound)->current_edge->dx < (*horz_bound)->next_edge->dx) { + break; + } + + // note: may be done multiple times + if ((*horz_bound)->ring && !is_open) { + add_point_to_ring(*(*horz_bound), + mapbox::geometry::point(std::llround((*bnd)->current_x), + scanline_y), + rings); + } + auto bnd_forward = --(bnd.base()); + + // OK, so far we're still in range of the horizontal Edge but make sure + // we're at the last of consec. horizontals when matching with eMaxPair + if (is_maxima_edge && bnd_forward == bound_max_pair) { + if ((*horz_bound)->ring) { + add_local_maximum_point(horz_bound, bound_max_pair, + (*horz_bound)->current_edge->top, rings, active_bounds); + } + active_bounds.erase(bound_max_pair); + return active_bounds.erase(horz_bound); + } + + intersect_bounds(bnd_forward, horz_bound, + mapbox::geometry::point(std::llround((*bnd)->current_x), + scanline_y), + cliptype, subject_fill_type, clip_fill_type, rings, active_bounds); + swap_positions_in_ABL(horz_bound, bnd_forward, active_bounds); + // Why are we not incrementing the bnd iterator here: + // It is because reverse iterators point to a `base()` iterator that is a forward + // iterator that is one ahead of the reverse bound. This will always be the horizontal + // bound, + // so what the reverse bound points to will have changed. + } // end while (bnd != active_bounds.rend()) + + if ((*horz_bound)->ring && !is_open) { + while (hp_itr != rings.hot_pixels.rend() && hp_itr->y == scanline_y && hp_itr->x > (*horz_bound)->current_edge->top.x) { + add_point_to_ring(*(*horz_bound), *hp_itr, rings); + ++hp_itr; + } + } + + if ((*horz_bound)->next_edge != (*horz_bound)->edges.end()) { + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), (*horz_bound)->current_edge->top, rings); + next_edge_in_bound(horz_bound, scanbeam); + + if ((*horz_bound)->winding_delta == 0) { + return std::next(horz_bound); + } + } else { + next_edge_in_bound(horz_bound, scanbeam); + } + return std::next(horz_bound); + } else { + if ((*horz_bound)->ring) { + add_point_to_ring(*(*horz_bound), (*horz_bound)->current_edge->top, rings); + } + return active_bounds.erase(horz_bound); + } +} + +template +active_bound_list_itr process_horizontal(T scanline_y, + active_bound_list_itr& horz_bound, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + if ((*horz_bound)->current_edge->bot.x < (*horz_bound)->current_edge->top.x) { + return process_horizontal_left_to_right(scanline_y, horz_bound, active_bounds, rings, + scanbeam, cliptype, subject_fill_type, + clip_fill_type); + } else { + return process_horizontal_right_to_left(scanline_y, horz_bound, active_bounds, rings, + scanbeam, cliptype, subject_fill_type, + clip_fill_type); + } +} + +template +void process_horizontals(T scanline_y, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + for (auto bnd_itr = active_bounds.begin(); bnd_itr != active_bounds.end();) { + if (current_edge_is_horizontal(bnd_itr)) { + bnd_itr = process_horizontal(scanline_y, bnd_itr, active_bounds, rings, scanbeam, + cliptype, subject_fill_type, clip_fill_type); + } else { + ++bnd_itr; + } + } +} +} +} +} diff --git a/mapbox/geometry/wagyu/process_maxima.hpp b/mapbox/geometry/wagyu/process_maxima.hpp new file mode 100644 index 0000000..edac1c3 --- /dev/null +++ b/mapbox/geometry/wagyu/process_maxima.hpp @@ -0,0 +1,133 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +active_bound_list_itr do_maxima(active_bound_list_itr& bnd, + active_bound_list_itr& bndMaxPair, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type, + ring_manager& rings, + active_bound_list& active_bounds) { + if (bndMaxPair == active_bounds.end()) { + if ((*bnd)->ring) { + add_point_to_ring(*(*bnd), (*bnd)->current_edge->top, rings); + } + return active_bounds.erase(bnd); + } + auto bnd_next = std::next(bnd); + auto return_bnd = bnd_next; + bool skipped = false; + while (bnd_next != active_bounds.end() && bnd_next != bndMaxPair) { + skipped = true; + intersect_bounds(bnd, bnd_next, (*bnd)->current_edge->top, cliptype, subject_fill_type, + clip_fill_type, rings, active_bounds); + swap_positions_in_ABL(bnd, bnd_next, active_bounds); + bnd_next = std::next(bnd); + } + + if (!(*bnd)->ring && !(*bndMaxPair)->ring) { + active_bounds.erase(bndMaxPair); + } else if ((*bnd)->ring && (*bndMaxPair)->ring) { + add_local_maximum_point(bnd, bndMaxPair, (*bnd)->current_edge->top, rings, active_bounds); + active_bounds.erase(bndMaxPair); + } else if ((*bnd)->winding_delta == 0 && (*bnd)->ring) { + add_point_to_ring(*(*bnd), (*bnd)->current_edge->top, rings); + active_bounds.erase(bndMaxPair); + } else if ((*bnd)->winding_delta == 0 && (*bndMaxPair)->ring) { + add_point_to_ring(*(*bndMaxPair), (*bnd)->current_edge->top, rings); + active_bounds.erase(bndMaxPair); + } else { + throw clipper_exception("DoMaxima error"); + } + auto prev_itr = active_bounds.erase(bnd); + if (skipped) { + return return_bnd; + } else { + return prev_itr; + } +} + +template +void process_edges_at_top_of_scanbeam(T top_y, + active_bound_list& active_bounds, + scanbeam_list& scanbeam, + local_minimum_ptr_list const& minima_sorted, + local_minimum_ptr_list_itr& current_lm, + ring_manager& rings, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + for (auto bnd = active_bounds.begin(); bnd != active_bounds.end();) { + // 1. Process maxima, treating them as if they are "bent" horizontal edges, + // but exclude maxima with horizontal edges. + + bool is_maxima_edge = is_maxima(bnd, top_y); + + if (is_maxima_edge) { + auto bnd_max_pair = get_maxima_pair(bnd, active_bounds); + is_maxima_edge = ((bnd_max_pair == active_bounds.end() || + !current_edge_is_horizontal(bnd_max_pair)) && + is_maxima(bnd_max_pair, top_y)); + if (is_maxima_edge) { + bnd = do_maxima(bnd, bnd_max_pair, cliptype, subject_fill_type, clip_fill_type, rings, + active_bounds); + continue; + } + } + + // 2. Promote horizontal edges. + if (is_intermediate(bnd, top_y) && next_edge_is_horizontal(bnd)) { + if ((*bnd)->ring) { + insert_hot_pixels_in_path(*(*bnd), (*bnd)->current_edge->top, rings, false); + } + next_edge_in_bound(bnd, scanbeam); + if ((*bnd)->ring) { + add_point_to_ring(*(*bnd), (*bnd)->current_edge->bot, rings); + } + } else { + (*bnd)->current_x = get_current_x(*((*bnd)->current_edge), top_y); + } + + ++bnd; + } + + insert_horizontal_local_minima_into_ABL(top_y, minima_sorted, current_lm, active_bounds, rings, + scanbeam, cliptype, subject_fill_type, clip_fill_type); + + process_horizontals(top_y, active_bounds, rings, scanbeam, cliptype, subject_fill_type, + clip_fill_type); + + // 4. Promote intermediate vertices + + for (auto bnd = active_bounds.begin(); bnd != active_bounds.end(); ++bnd) { + if (is_intermediate(bnd, top_y)) { + if ((*bnd)->ring) { + add_point_to_ring(*(*bnd), (*bnd)->current_edge->top, rings); + insert_hot_pixels_in_path(*(*bnd), (*bnd)->current_edge->top, rings, false); + } + next_edge_in_bound(bnd, scanbeam); + } + } +} + +} +} +} diff --git a/mapbox/geometry/wagyu/ring.hpp b/mapbox/geometry/wagyu/ring.hpp new file mode 100644 index 0000000..b4eda27 --- /dev/null +++ b/mapbox/geometry/wagyu/ring.hpp @@ -0,0 +1,479 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEBUG +#include +#include +#include +// +// void* callstack[128]; +// int i, frames = backtrace(callstack, 128); +// char** strs = backtrace_symbols(callstack, frames); +// for (i = 0; i < frames; ++i) { +// printf("%s\n", strs[i]); +// } +// free(strs); +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +// NOTE: ring and ring_ptr are forward declared in wagyu/point.hpp + +template +using ring_vector = std::vector>; + +template +using ring_list = std::list>; + +template +struct ring { + std::size_t ring_index; // To support unset 0 is undefined and indexes offset by 1 + std::size_t size; + double area; + ring_ptr parent; + ring_list children; + point_ptr points; + point_ptr bottom_point; + bool is_open; + + ring( ring const& ) = delete; + ring& operator=(ring const& ) = delete; + + ring() + : ring_index(0), + size(0), + area(std::numeric_limits::quiet_NaN()), + parent(nullptr), + children(), + points(nullptr), + bottom_point(nullptr), + is_open(false) { + } +}; + +template +using hot_pixel_vector = std::vector>; + +template +using hot_pixel_itr = typename hot_pixel_vector::iterator; + +template +using hot_pixel_rev_itr = typename hot_pixel_vector::reverse_iterator; + +template +struct ring_manager { + + ring_list children; + std::vector> all_points; + hot_pixel_vector hot_pixels; + hot_pixel_itr current_hp_itr; + std::deque> points; + std::deque> rings; + std::vector> storage; + std::size_t index; + + ring_manager( ring_manager const& ) = delete; + ring_manager& operator=(ring_manager const& ) = delete; + + ring_manager() + : children(), + all_points(), + hot_pixels(), + current_hp_itr(hot_pixels.end()), + points(), + rings(), + storage(), + index(0) { + } +}; + +template +void preallocate_point_memory(ring_manager& rings, std::size_t size) { + rings.storage.reserve(size); + rings.all_points.reserve(size); +} + +template +ring_ptr create_new_ring(ring_manager& rings) { + rings.rings.emplace_back(); + ring_ptr result = &rings.rings.back(); + result->ring_index = rings.index++; + return result; +} + +template +point_ptr +create_new_point(ring_ptr r, mapbox::geometry::point const& pt, ring_manager& rings) { + point_ptr point; + if (rings.storage.size() < rings.storage.capacity()) { + rings.storage.emplace_back(r, pt); + point = &rings.storage.back(); + } else { + rings.points.emplace_back(r, pt); + point = &rings.points.back(); + } + rings.all_points.push_back(point); + return point; +} + +template +point_ptr create_new_point(ring_ptr r, + mapbox::geometry::point const& pt, + point_ptr before_this_point, + ring_manager& rings) { + point_ptr point; + if (rings.storage.size() < rings.storage.capacity()) { + rings.storage.emplace_back(r, pt, before_this_point); + point = &rings.storage.back(); + } else { + rings.points.emplace_back(r, pt, before_this_point); + point = &rings.points.back(); + } + rings.all_points.push_back(point); + return point; +} + +template +void ring1_child_of_ring2(ring_ptr ring1, ring_ptr ring2, ring_manager& manager) { + assert(ring1 != ring2); + if (ring1->parent == ring2) { + return; + } + if (ring1->parent == nullptr) { + manager.children.remove(ring1); + } else { + ring1->parent->children.remove(ring1); + } + if (ring2 == nullptr) { + ring1->parent = nullptr; + manager.children.push_back(ring1); + } else { + ring1->parent = ring2; + ring2->children.push_back(ring1); + } +} + +template +void ring1_sibling_of_ring2(ring_ptr ring1, ring_ptr ring2, ring_manager& manager) { + assert(ring1 != ring2); + if (ring1->parent == ring2->parent) { + return; + } + if (ring1->parent == nullptr) { + manager.children.remove(ring1); + } else { + ring1->parent->children.remove(ring1); + } + if (ring2->parent == nullptr) { + manager.children.push_back(ring1); + } else { + ring2->parent->children.push_back(ring1); + } + ring1->parent = ring2->parent; +} + +template +void ring1_replaces_ring2(ring_ptr ring1, ring_ptr ring2, ring_manager& manager) { + assert(ring1 != ring2); + if (ring2->parent == nullptr) { + manager.children.remove(ring2); + } else { + ring2->parent->children.remove(ring2); + } + for (auto& c : ring2->children) { + c->parent = ring1; + } + if (ring1 == nullptr) { + manager.children.splice(manager.children.end(), ring2->children); + } else { + ring1->children.splice(ring1->children.end(), ring2->children); + } + ring2->parent = nullptr; +} + +template +void remove_ring(ring_ptr r, ring_manager& manager) { + if (r->parent == nullptr) { + manager.children.remove(r); + for (auto& c : r->children) { + c->parent = nullptr; + } + manager.children.splice(manager.children.end(), r->children); + } else { + r->parent->children.remove(r); + for (auto& c : r->children) { + c->parent = r->parent; + } + r->parent->children.splice(r->parent->children.end(), r->children); + r->parent = nullptr; + } +} + +template +inline std::size_t ring_depth(ring_ptr r) { + std::size_t depth = 0; + if (!r) { + return depth; + } + while (r->parent) { + depth++; + r = r->parent; + } + return depth; +} + +template +inline bool ring_is_hole(ring_ptr r) { + return ring_depth(r) & 1; +} + +template +void set_next(const_point_ptr& node, const const_point_ptr& next_node) { + node->next = next_node; +} + +template +point_ptr get_next(const_point_ptr& node) { + return node->next; +} + +template +point_ptr get_prev(const_point_ptr& node) { + return node->prev; +} + +template +void set_prev(const_point_ptr& node, const const_point_ptr& prev_node) { + node->prev = prev_node; +} + +template +void init(const_point_ptr& node) { + set_next(node, node); + set_prev(node, node); +} + +template +std::size_t point_count(const const_point_ptr& orig_node) { + std::size_t size = 0; + point_ptr n = orig_node; + do { + n = get_next(n); + ++size; + } while (n != orig_node); + return size; +} + +template +void link_before(point_ptr& node, point_ptr& new_node) { + point_ptr prev_node = get_prev(node); + set_prev(new_node, prev_node); + set_next(new_node, node); + set_prev(node, new_node); + set_next(prev_node, new_node); +} + +template +void link_after(point_ptr& node, point_ptr& new_node) { + point_ptr next_node = get_next(node); + set_prev(new_node, node); + set_next(new_node, next_node); + set_next(node, new_node); + set_prev(next_node, new_node); +} + +template +void transfer_point(point_ptr& p, point_ptr& b, point_ptr& e) { + if (b != e) { + point_ptr prev_p = get_prev(p); + point_ptr prev_b = get_prev(b); + point_ptr prev_e = get_prev(e); + set_next(prev_e, p); + set_prev(p, prev_e); + set_next(prev_b, e); + set_prev(e, prev_b); + set_next(prev_p, b); + set_prev(b, prev_p); + } else { + link_before(p, b); + } +} + +template +void reverse_ring(point_ptr pp) { + if (!pp) { + return; + } + point_ptr pp1; + point_ptr pp2; + pp1 = pp; + do { + pp2 = pp1->next; + pp1->next = pp1->prev; + pp1->prev = pp2; + pp1 = pp2; + } while (pp1 != pp); +} + +template +double area_from_point(point_ptr op, std::size_t & size) { + point_ptr startOp = op; + size = 1; + if (!op) { + return 0.0; + } + double a = 0.0; + do { + ++size; + a += static_cast(op->prev->x + op->x) * static_cast(op->prev->y - op->y); + op = op->next; + } while (op != startOp); + return a * 0.5; +} + +template +double area(ring_ptr r) { + assert(r != nullptr); + if (std::isnan(r->area)) { + r->area = area_from_point(r->points, r->size); + } + return r->area; +} + +#ifdef DEBUG + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const ring& r) { + out << " ring_index: " << r.ring_index << std::endl; + if (!r.parent) { + // out << " parent_ring ptr: nullptr" << std::endl; + out << " parent_index: -----" << std::endl; + } else { + // out << " parent_ring ptr: " << r.parent << std::endl; + out << " parent_ring idx: " << r.parent->ring_index << std::endl; + } + ring_ptr n = const_cast>(&r); + if (ring_is_hole(n)) { + out << " is_hole: true " << std::endl; + } else { + out << " is_hole: false " << std::endl; + } + auto pt_itr = r.points; + if (pt_itr) { + out << " area: " << r.area << std::endl; + out << " points:" << std::endl; + out << " [[[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + while (pt_itr != r.points) { + out << "[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + } + out << "[" << pt_itr->x << "," << pt_itr->y << "]]]" << std::endl; + } else { + out << " area: NONE" << std::endl; + out << " points: NONE" << std::endl; + } + return out; +} + +template +std::string output_as_polygon(ring_ptr r) { + std::ostringstream out; + + auto pt_itr = r->points; + if (pt_itr) { + out << "["; + out << "[[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + while (pt_itr != r->points) { + out << "[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + } + out << "[" << pt_itr->x << "," << pt_itr->y << "]]"; + for (auto const& c : r->children) { + pt_itr = c->points; + if (pt_itr) { + out << ",[[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + while (pt_itr != c->points) { + out << "[" << pt_itr->x << "," << pt_itr->y << "],"; + pt_itr = pt_itr->next; + } + out << "[" << pt_itr->x << "," << pt_itr->y << "]]"; + } + } + out << "]" << std::endl; + } else { + out << "[]" << std::endl; + } + + return out.str(); +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const ring_list& rings) { + out << "START RING LIST" << std::endl; + for (auto& r : rings) { + out << " ring: " << r->ring_index << " - " << r << std::endl; + out << *r; + } + out << "END RING LIST" << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const ring_vector& rings) { + out << "START RING VECTOR" << std::endl; + for (auto& r : rings) { + if (!r->points) { + continue; + } + out << " ring: " << r->ring_index << " - " << r << std::endl; + out << *r; + } + out << "END RING VECTOR" << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const std::deque>& rings) { + out << "START RING VECTOR" << std::endl; + for (auto& r : rings) { + if (!r.points) { + continue; + } + out << " ring: " << r.ring_index << std::endl; + out << r; + } + out << "END RING VECTOR" << std::endl; + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const hot_pixel_vector& hp_vec) { + out << "Hot Pixels: " << std::endl; + for (auto& hp : hp_vec) { + out << hp << std::endl; + } + return out; +} +#endif +} +} +} diff --git a/mapbox/geometry/wagyu/ring_util.hpp b/mapbox/geometry/wagyu/ring_util.hpp new file mode 100644 index 0000000..0988a77 --- /dev/null +++ b/mapbox/geometry/wagyu/ring_util.hpp @@ -0,0 +1,915 @@ +#pragma once + +#ifdef DEBUG +#include +// Example debug print for backtrace - only works on IOS +#include +#include +// +// void* callstack[128]; +// int i, frames = backtrace(callstack, 128); +// char** strs = backtrace_symbols(callstack, frames); +// for (i = 0; i < frames; ++i) { +// printf("%s\n", strs[i]); +// } +// free(strs); +#endif + +#include + +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void set_hole_state(active_bound_list_itr& bnd, + active_bound_list& active_bounds, + ring_manager& rings) { + auto bnd2 = active_bound_list_rev_itr(bnd); + bound_ptr bndTmp = nullptr; + // Find first non line ring to the left of current bound. + while (bnd2 != active_bounds.rend()) { + if ((*bnd2)->ring && (*bnd2)->winding_delta != 0) { + if (!bndTmp) { + bndTmp = (*bnd2); + } else if (bndTmp->ring == (*bnd2)->ring) { + bndTmp = nullptr; + } + } + ++bnd2; + } + if (!bndTmp) { + (*bnd)->ring->parent = nullptr; + rings.children.push_back((*bnd)->ring); + } else { + (*bnd)->ring->parent = bndTmp->ring; + bndTmp->ring->children.push_back((*bnd)->ring); + } +} + +template +void set_hole_state(active_bound_list_rev_itr& bnd, + active_bound_list& active_bounds, + ring_manager& rings) { + auto bnd2 = std::next(bnd); + bound_ptr bndTmp = nullptr; + // Find first non line ring to the left of current bound. + while (bnd2 != active_bounds.rend()) { + if ((*bnd2)->ring && (*bnd2)->winding_delta != 0) { + if (!bndTmp) { + bndTmp = (*bnd2); + } else if (bndTmp->ring == (*bnd2)->ring) { + bndTmp = nullptr; + } + } + ++bnd2; + } + + if (!bndTmp) { + (*bnd)->ring->parent = nullptr; + rings.children.push_back((*bnd)->ring); + } else { + (*bnd)->ring->parent = bndTmp->ring; + bndTmp->ring->children.push_back((*bnd)->ring); + } +} + +template +void update_current_hp_itr(T scanline_y, ring_manager& rings) { + while (rings.current_hp_itr->y > scanline_y) { + ++rings.current_hp_itr; + } +} + +template +struct hot_pixel_sorter { + inline bool operator()(mapbox::geometry::point const& pt1, mapbox::geometry::point const& pt2) { + if (pt1.y == pt2.y) { + return pt1.x < pt2.x; + } else { + return pt1.y > pt2.y; + } + } +}; + +// Due to the nature of floating point calculations +// and the high likely hood of values around X.5, we +// need to fudge what is X.5 some for our rounding. +const double rounding_offset = 1e-12; +const double rounding_offset_y = 5e-13; + +template +T round_towards_min(double val) { + // 0.5 rounds to 0 + // 0.0 rounds to 0 + // -0.5 rounds to -1 + return static_cast(std::ceil(val - 0.5 + rounding_offset)); +} + +template +T round_towards_max(double val) { + // 0.5 rounds to 1 + // 0.0 rounds to 0 + // -0.5 rounds to 0 + return static_cast(std::floor(val + 0.5 + rounding_offset)); +} + +template +inline T get_edge_min_x(edge const& edge, const T current_y) { + if (is_horizontal(edge)) { + if (edge.bot.x < edge.top.x) { + return edge.bot.x; + } else { + return edge.top.x; + } + } else if (edge.dx > 0.0) { + if (current_y == edge.top.y) { + return edge.top.x; + } else { + double lower_range_y = static_cast(current_y - edge.bot.y) - 0.5; + double return_val = static_cast(edge.bot.x) + edge.dx * lower_range_y; + T value = round_towards_min(return_val); + return value; + } + } else { + if (current_y == edge.bot.y) { + return edge.bot.x; + } else { + double return_val = static_cast(edge.bot.x) + + edge.dx * (static_cast(current_y - edge.bot.y) + 0.5 - rounding_offset_y); + T value = round_towards_min(return_val); + return value; + } + } +} + +template +inline T get_edge_max_x(edge const& edge, const T current_y) { + if (is_horizontal(edge)) { + if (edge.bot.x > edge.top.x) { + return edge.bot.x; + } else { + return edge.top.x; + } + } else if (edge.dx < 0.0) { + if (current_y == edge.top.y) { + return edge.top.x; + } else { + double lower_range_y = static_cast(current_y - edge.bot.y) - 0.5; + double return_val = static_cast(edge.bot.x) + edge.dx * lower_range_y; + T value = round_towards_max(return_val); + return value; + } + } else { + if (current_y == edge.bot.y) { + return edge.bot.x; + } else { + double return_val = static_cast(edge.bot.x) + + edge.dx * (static_cast(current_y - edge.bot.y) + 0.5 - rounding_offset_y); + T value = round_towards_max(return_val); + return value; + } + } +} + +template +void hot_pixel_set_left_to_right(T y, + T start_x, + T end_x, + bound& bnd, + ring_manager& rings, + hot_pixel_itr & itr, + hot_pixel_itr & end, + bool add_end_point) { + T x_min = get_edge_min_x(*(bnd.current_edge), y); + x_min = std::max(x_min, start_x); + T x_max = get_edge_max_x(*(bnd.current_edge), y); + x_max = std::min(x_max, end_x); + for (;itr != end; ++itr) { + if (itr->x < x_min) { + continue; + } + if (itr->x > x_max) { + break; + } + if (!add_end_point && itr->x == end_x) { + continue; + } + point_ptr op = bnd.ring->points; + bool to_front = (bnd.side == edge_left); + if (to_front && (*itr == *op)) { + continue; + } else if (!to_front && (*itr == *op->prev)) { + continue; + } + point_ptr new_point = create_new_point(bnd.ring, *itr, op, rings); + if (to_front) { + bnd.ring->points = new_point; + } + } +} + +template +void hot_pixel_set_right_to_left(T y, + T start_x, + T end_x, + bound& bnd, + ring_manager& rings, + hot_pixel_rev_itr & itr, + hot_pixel_rev_itr & end, + bool add_end_point) { + T x_min = get_edge_min_x(*(bnd.current_edge), y); + x_min = std::max(x_min, end_x); + T x_max = get_edge_max_x(*(bnd.current_edge), y); + x_max = std::min(x_max, start_x); + for (;itr != end; ++itr) { + if (itr->x > x_max) { + continue; + } + if (itr->x < x_min) { + break; + } + if (!add_end_point && itr->x == end_x) { + continue; + } + point_ptr op = bnd.ring->points; + bool to_front = (bnd.side == edge_left); + if (to_front && (*itr == *op)) { + continue; + } else if (!to_front && (*itr == *op->prev)) { + continue; + } + point_ptr new_point = create_new_point(bnd.ring, *itr, op, rings); + if (to_front) { + bnd.ring->points = new_point; + } + } +} + +template +void sort_hot_pixels(ring_manager& rings) { + std::sort(rings.hot_pixels.begin(), rings.hot_pixels.end(), hot_pixel_sorter()); + auto last = std::unique(rings.hot_pixels.begin(), rings.hot_pixels.end()); + rings.hot_pixels.erase(last, rings.hot_pixels.end()); +} + +template +void insert_hot_pixels_in_path(bound& bnd, + mapbox::geometry::point const& end_pt, + ring_manager& rings, + bool add_end_point) { + if (end_pt == bnd.last_point) { + return; + } + if (!bnd.ring) { + bnd.last_point = end_pt; + return; + } + + T start_y = bnd.last_point.y; + T start_x = bnd.last_point.x; + T end_y = end_pt.y; + T end_x = end_pt.x; + + auto itr = rings.current_hp_itr; + while (itr->y <= start_y && itr != rings.hot_pixels.begin()) { + --itr; + } + if (start_x > end_x) { + for (; itr != rings.hot_pixels.end();) { + if (itr->y > start_y) { + ++itr; + continue; + } + if (itr->y < end_y) { + break; + } + T y = itr->y; + auto last_itr = hot_pixel_rev_itr(itr); + while (itr != rings.hot_pixels.end() && itr->y == y) { + ++itr; + } + auto first_itr = hot_pixel_rev_itr(itr); + bool add_end_point_itr = (y != end_pt.y || add_end_point); + hot_pixel_set_right_to_left(y, start_x, end_x, bnd, rings, first_itr, last_itr, + add_end_point_itr); + } + } else { + for (; itr != rings.hot_pixels.end();) { + if (itr->y > start_y) { + ++itr; + continue; + } + if (itr->y < end_y) { + break; + } + T y = itr->y; + auto first_itr = itr; + while (itr != rings.hot_pixels.end() && itr->y == y) { + ++itr; + } + auto last_itr = itr; + bool add_end_point_itr = (y != end_pt.y || add_end_point); + hot_pixel_set_left_to_right(y, start_x, end_x, bnd, rings, first_itr, last_itr, + add_end_point_itr); + } + } + bnd.last_point = end_pt; +} + +template +void add_to_hot_pixels(mapbox::geometry::point const& pt, ring_manager& rings) { + rings.hot_pixels.push_back(pt); +} + +template +void add_first_point(active_bound_list_itr& bnd, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + + ring_ptr r = create_new_ring(rings); + (*bnd)->ring = r; + r->is_open = ((*bnd)->winding_delta == 0); + r->points = create_new_point(r, pt, rings); + if (!r->is_open) { + set_hole_state(bnd, active_bounds, rings); + } + (*bnd)->last_point = pt; +} + +template +void add_first_point(active_bound_list_rev_itr& bnd, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + ring_ptr r = create_new_ring(rings); + // no ring currently set! + (*bnd)->ring = r; + r->is_open = ((*bnd)->winding_delta == 0); + r->points = create_new_point(r, pt, rings); + if (!r->is_open) { + set_hole_state(bnd, active_bounds, rings); + } + (*bnd)->last_point = pt; +} + +template +void add_point_to_ring(bound& bnd, + mapbox::geometry::point const& pt, + ring_manager& rings) { + assert(bnd.ring); + // Handle hot pixels + insert_hot_pixels_in_path(bnd, pt, rings, false); + + // bnd.ring->points is the 'Left-most' point & bnd.ring->points->prev is the + // 'Right-most' + point_ptr op = bnd.ring->points; + bool to_front = (bnd.side == edge_left); + if (to_front && (pt == *op)) { + return; + } else if (!to_front && (pt == *op->prev)) { + return; + } + point_ptr new_point = create_new_point(bnd.ring, pt, bnd.ring->points, rings); + if (to_front) { + bnd.ring->points = new_point; + } +} + +template +void add_point(active_bound_list_itr& bnd, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + if (!(*bnd)->ring) { + add_first_point(bnd, active_bounds, pt, rings); + } else { + add_point_to_ring(*(*bnd), pt, rings); + } +} + +template +void add_point(active_bound_list_rev_itr& bnd, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + if (!(*bnd)->ring) { + add_first_point(bnd, active_bounds, pt, rings); + } else { + add_point_to_ring(*(*bnd), pt, rings); + } +} + +template +void add_local_minimum_point(active_bound_list_itr b1, + active_bound_list_itr b2, + active_bound_list& active_bounds, + mapbox::geometry::point const& pt, + ring_manager& rings) { + active_bound_list_itr b; + active_bound_list_rev_itr prev_bound; + active_bound_list_rev_itr prev_b1(b1); + active_bound_list_rev_itr prev_b2(b2); + if (is_horizontal(*((*b2)->current_edge)) || + ((*b1)->current_edge->dx > (*b2)->current_edge->dx)) { + add_point(b1, active_bounds, pt, rings); + (*b2)->last_point = pt; + (*b2)->ring = (*b1)->ring; + (*b1)->side = edge_left; + (*b2)->side = edge_right; + b = b1; + if (prev_b1 != active_bounds.rend() && std::prev(b) == b2) { + prev_bound = prev_b2; + } else { + prev_bound = prev_b1; + } + } else { + add_point(b2, active_bounds, pt, rings); + (*b1)->last_point = pt; + (*b1)->ring = (*b2)->ring; + (*b1)->side = edge_right; + (*b2)->side = edge_left; + b = b2; + if (prev_b2 != active_bounds.rend() && std::prev(b) == b1) { + prev_bound = prev_b1; + } else { + prev_bound = prev_b2; + } + } +} + +template +inline double get_dx(point const& pt1, point const& pt2) { + if (pt1.y == pt2.y) { + return std::numeric_limits::infinity(); + } else { + return static_cast(pt2.x - pt2.x) / static_cast(pt2.y - pt1.y); + } +} + +template +bool first_is_bottom_point(const_point_ptr btmPt1, const_point_ptr btmPt2) { + point_ptr p = btmPt1->prev; + while ((*p == *btmPt1) && (p != btmPt1)) { + p = p->prev; + } + double dx1p = std::fabs(get_dx(*btmPt1, *p)); + + p = btmPt1->next; + while ((*p == *btmPt1) && (p != btmPt1)) { + p = p->next; + } + double dx1n = std::fabs(get_dx(*btmPt1, *p)); + + p = btmPt2->prev; + while ((*p == *btmPt2) && (p != btmPt2)) { + p = p->prev; + } + double dx2p = std::fabs(get_dx(*btmPt2, *p)); + + p = btmPt2->next; + while ((*p == *btmPt2) && (p != btmPt2)) { + p = p->next; + } + double dx2n = std::fabs(get_dx(*btmPt2, *p)); + + if (values_are_equal(std::max(dx1p, dx1n), std::max(dx2p, dx2n)) && + values_are_equal(std::min(dx1p, dx1n), std::min(dx2p, dx2n))) { + std::size_t s = 0; + return area_from_point(btmPt1, s) > 0.0; // if otherwise identical use orientation + } else { + return (greater_than_or_equal(dx1p, dx2p) && greater_than_or_equal(dx1p, dx2n)) || + (greater_than_or_equal(dx1n, dx2p) && greater_than_or_equal(dx1n, dx2n)); + } +} + +template +point_ptr get_bottom_point(point_ptr pp) { + point_ptr dups = 0; + point_ptr p = pp->next; + while (p != pp) { + if (p->y > pp->y) { + pp = p; + dups = 0; + } else if (p->y == pp->y && p->x <= pp->x) { + if (p->x < pp->x) { + dups = 0; + pp = p; + } else { + if (p->next != pp && p->prev != pp) { + dups = p; + } + } + } + p = p->next; + } + if (dups) { + // there appears to be at least 2 vertices at bottom_point so ... + while (dups != p) { + if (!first_is_bottom_point(p, dups)) { + pp = dups; + } + dups = dups->next; + while (*dups != *pp) { + dups = dups->next; + } + } + } + return pp; +} + +template +ring_ptr get_lower_most_ring(ring_ptr outRec1, ring_ptr outRec2) { + // work out which polygon fragment has the correct hole state ... + if (!outRec1->bottom_point) { + outRec1->bottom_point = get_bottom_point(outRec1->points); + } + if (!outRec2->bottom_point) { + outRec2->bottom_point = get_bottom_point(outRec2->points); + } + point_ptr OutPt1 = outRec1->bottom_point; + point_ptr OutPt2 = outRec2->bottom_point; + if (OutPt1->y > OutPt2->y) { + return outRec1; + } else if (OutPt1->y < OutPt2->y) { + return outRec2; + } else if (OutPt1->x < OutPt2->x) { + return outRec1; + } else if (OutPt1->x > OutPt2->x) { + return outRec2; + } else if (OutPt1->next == OutPt1) { + return outRec2; + } else if (OutPt2->next == OutPt2) { + return outRec1; + } else if (first_is_bottom_point(OutPt1, OutPt2)) { + return outRec1; + } else { + return outRec2; + } +} + +template +bool ring1_right_of_ring2(ring_ptr ring1, ring_ptr ring2) { + do { + ring1 = ring1->parent; + if (ring1 == ring2) { + return true; + } + } while (ring1); + return false; +} + +template +void update_points_ring(ring_ptr ring) { + point_ptr op = ring->points; + do { + op->ring = ring; + op = op->prev; + } while (op != ring->points); +} + +template +void append_ring(active_bound_list_itr& b1, + active_bound_list_itr& b2, + active_bound_list& active_bounds, + ring_manager& manager) { + // get the start and ends of both output polygons ... + ring_ptr outRec1 = (*b1)->ring; + ring_ptr outRec2 = (*b2)->ring; + + ring_ptr keep_ring; + bound_ptr keep_bound; + ring_ptr remove_ring; + bound_ptr remove_bound; + if (ring1_right_of_ring2(outRec1, outRec2)) { + keep_ring = outRec2; + keep_bound = *b2; + remove_ring = outRec1; + remove_bound = *b1; + } else if (ring1_right_of_ring2(outRec2, outRec1)) { + keep_ring = outRec1; + keep_bound = *b1; + remove_ring = outRec2; + remove_bound = *b2; + } else if (outRec1 == get_lower_most_ring(outRec1, outRec2)) { + keep_ring = outRec1; + keep_bound = *b1; + remove_ring = outRec2; + remove_bound = *b2; + } else { + keep_ring = outRec2; + keep_bound = *b2; + remove_ring = outRec1; + remove_bound = *b1; + } + + // get the start and ends of both output polygons and + // join b2 poly onto b1 poly and delete pointers to b2 ... + + point_ptr p1_lft = keep_ring->points; + point_ptr p1_rt = p1_lft->prev; + point_ptr p2_lft = remove_ring->points; + point_ptr p2_rt = p2_lft->prev; + + // join b2 poly onto b1 poly and delete pointers to b2 ... + if (keep_bound->side == edge_left) { + if (remove_bound->side == edge_left) { + // z y x a b c + reverse_ring(p2_lft); + p2_lft->next = p1_lft; + p1_lft->prev = p2_lft; + p1_rt->next = p2_rt; + p2_rt->prev = p1_rt; + keep_ring->points = p2_rt; + } else { + // x y z a b c + p2_rt->next = p1_lft; + p1_lft->prev = p2_rt; + p2_lft->prev = p1_rt; + p1_rt->next = p2_lft; + keep_ring->points = p2_lft; + } + } else { + if (remove_bound->side == edge_right) { + // a b c z y x + reverse_ring(p2_lft); + p1_rt->next = p2_rt; + p2_rt->prev = p1_rt; + p2_lft->next = p1_lft; + p1_lft->prev = p2_lft; + } else { + // a b c x y z + p1_rt->next = p2_lft; + p2_lft->prev = p1_rt; + p1_lft->prev = p2_rt; + p2_rt->next = p1_lft; + } + } + + keep_ring->bottom_point = nullptr; + bool keep_is_hole = ring_is_hole(keep_ring); + bool remove_is_hole = ring_is_hole(remove_ring); + + remove_ring->points = nullptr; + remove_ring->bottom_point = nullptr; + if (keep_is_hole != remove_is_hole) { + ring1_replaces_ring2(keep_ring->parent, remove_ring, manager); + } else { + ring1_replaces_ring2(keep_ring, remove_ring, manager); + } + + update_points_ring(keep_ring); + + // nb: safe because we only get here via AddLocalMaxPoly + keep_bound->ring = nullptr; + remove_bound->ring = nullptr; + + for (auto& b : active_bounds) { + if (b->ring == remove_ring) { + b->ring = keep_ring; + b->side = keep_bound->side; + break; // Not sure why there is a break here but was transfered logic from angus + } + } +} + +template +void add_local_maximum_point(active_bound_list_itr& b1, + active_bound_list_itr& b2, + mapbox::geometry::point const& pt, + ring_manager& rings, + active_bound_list& active_bounds) { + insert_hot_pixels_in_path(*(*b2), pt, rings, false); + add_point(b1, active_bounds, pt, rings); + if ((*b2)->winding_delta == 0) { + add_point(b2, active_bounds, pt, rings); + } + if ((*b1)->ring == (*b2)->ring) { + (*b1)->ring = nullptr; + (*b2)->ring = nullptr; + // I am not certain that order is important here? + } else if ((*b1)->ring->ring_index < (*b2)->ring->ring_index) { + append_ring(b1, b2, active_bounds, rings); + } else { + append_ring(b2, b1, active_bounds, rings); + } +} + +enum point_in_polygon_result : std::int8_t { + point_on_polygon = -1, + point_inside_polygon = 0, + point_outside_polygon = 1 +}; + +template +point_in_polygon_result point_in_polygon(point const& pt, point_ptr op) { + // returns 0 if false, +1 if true, -1 if pt ON polygon boundary + point_in_polygon_result result = point_outside_polygon; + point_ptr startOp = op; + do { + if (op->next->y == pt.y) { + if ((op->next->x == pt.x) || + (op->y == pt.y && ((op->next->x > pt.x) == (op->x < pt.x)))) { + return point_on_polygon; + } + } + if ((op->y < pt.y) != (op->next->y < pt.y)) { + if (op->x >= pt.x) { + if (op->next->x > pt.x) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } else { + double d = + static_cast(op->x - pt.x) * + static_cast(op->next->y - pt.y) - + static_cast(op->next->x - pt.x) * static_cast(op->y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0) == (op->next->y > op->y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } else { + if (op->next->x > pt.x) { + double d = + static_cast(op->x - pt.x) * + static_cast(op->next->y - pt.y) - + static_cast(op->next->x - pt.x) * static_cast(op->y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0) == (op->next->y > op->y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } + } + op = op->next; + } while (startOp != op); + return result; +} + +template +point_in_polygon_result point_in_polygon(mapbox::geometry::point const& pt, + point_ptr op) { + // returns 0 if false, +1 if true, -1 if pt ON polygon boundary + point_in_polygon_result result = point_outside_polygon; + point_ptr startOp = op; + do { + double op_x = static_cast(op->x); + double op_y = static_cast(op->y); + double op_next_x = static_cast(op->next->x); + double op_next_y = static_cast(op->next->y); + if (values_are_equal(op_next_y, pt.y)) { + if (values_are_equal(op_next_x, pt.x) || + (values_are_equal(op_y, pt.y) && ((op_next_x > pt.x) == (op_x < pt.x)))) { + return point_on_polygon; + } + } + if ((op_y < pt.y) != (op_next_y < pt.y)) { + if (greater_than_or_equal(op_x, pt.x)) { + if (op_next_x > pt.x) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } else { + double d = + (op_x - pt.x) * (op_next_y - pt.y) - (op_next_x - pt.x) * (op_y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0.0) == (op_next_y > op->y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } else { + if (op_next_x > pt.x) { + double d = + (op_x - pt.x) * (op_next_y - pt.y) - (op_next_x - pt.x) * (op_y - pt.y); + if (value_is_zero(d)) { + return point_on_polygon; + } + if ((d > 0.0) == (op_next_y > op->y)) { + // Switch between point outside polygon and point inside + // polygon + if (result == point_outside_polygon) { + result = point_inside_polygon; + } else { + result = point_outside_polygon; + } + } + } + } + } + op = op->next; + } while (startOp != op); + return result; +} + +template +point_in_polygon_result inside_or_outside_special(point_ptr first_pt, point_ptr other_poly) { + + if (value_is_zero(area(first_pt->ring))) { + return point_inside_polygon; + } + if (value_is_zero(area(other_poly->ring))) { + return point_outside_polygon; + } + point_ptr pt = first_pt; + do { + if (*pt == *(pt->prev) || *pt == *(pt->next) || *(pt->next) == *(pt->prev) || + slopes_equal(*(pt->prev), *pt, *(pt->next))) { + pt = pt->next; + continue; + } + double dx = ((pt->prev->x - pt->x) / 3.0) + ((pt->next->x - pt->x) / 3.0); + double dy = ((pt->prev->y - pt->y) / 3.0) + ((pt->next->y - pt->y) / 3.0); + mapbox::geometry::point offset_pt(pt->x + dx, pt->y + dy); + point_in_polygon_result res = point_in_polygon(offset_pt, pt); + if (res != point_inside_polygon) { + offset_pt.x = pt->x - dx; + offset_pt.y = pt->y - dy; + res = point_in_polygon(offset_pt, pt); + if (res != point_inside_polygon) { + pt = pt->next; + continue; + } + } + res = point_in_polygon(offset_pt, other_poly); + if (res == point_on_polygon) { + pt = pt->next; + continue; + } + return res; + } while (pt != first_pt); + return point_inside_polygon; +} + +template +bool poly2_contains_poly1(point_ptr outpt1, point_ptr outpt2) { + point_ptr op = outpt1; + do { + // nb: PointInPolygon returns 0 if false, +1 if true, -1 if pt on polygon + point_in_polygon_result res = point_in_polygon(*op, outpt2); + if (res != point_on_polygon) { + return res == point_inside_polygon; + } + op = op->next; + } while (op != outpt1); + point_in_polygon_result res = inside_or_outside_special(outpt1, outpt2); + return res == point_inside_polygon; +} + +template +void dispose_out_points(point_ptr& pp) { + if (pp == nullptr) { + return; + } + pp->prev->next = nullptr; + while (pp) { + point_ptr tmpPp = pp; + pp = pp->next; + tmpPp->next = tmpPp; + tmpPp->prev = tmpPp; + tmpPp->ring = nullptr; + // delete tmpPp; + } +} +} +} +} diff --git a/mapbox/geometry/wagyu/scanbeam.hpp b/mapbox/geometry/wagyu/scanbeam.hpp new file mode 100644 index 0000000..a3bf8f7 --- /dev/null +++ b/mapbox/geometry/wagyu/scanbeam.hpp @@ -0,0 +1,37 @@ +#pragma once + +#include + +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +using scanbeam_list = std::priority_queue; + +template +bool pop_from_scanbeam(T& Y, scanbeam_list& scanbeam) { + if (scanbeam.empty()) { + return false; + } + Y = scanbeam.top(); + scanbeam.pop(); + while (!scanbeam.empty() && Y == scanbeam.top()) { + scanbeam.pop(); + } // Pop duplicates. + return true; +} + +template +void setup_scanbeam(local_minimum_list& minima_list, scanbeam_list& scanbeam) { + + for (auto lm = minima_list.begin(); lm != minima_list.end(); ++lm) { + scanbeam.push(lm->y); + } +} +} +} +} diff --git a/mapbox/geometry/wagyu/snap_rounding.hpp b/mapbox/geometry/wagyu/snap_rounding.hpp new file mode 100644 index 0000000..fc65d27 --- /dev/null +++ b/mapbox/geometry/wagyu/snap_rounding.hpp @@ -0,0 +1,174 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +void process_hot_pixel_intersections(T top_y, + active_bound_list& active_bounds, + ring_manager& rings) { + if (active_bounds.empty()) { + return; + } + + update_current_x(active_bounds, top_y); + // bubblesort ... + bool isModified; + do { + isModified = false; + auto bnd = active_bounds.begin(); + auto bnd_next = std::next(bnd); + while (bnd_next != active_bounds.end()) { + if ((*bnd)->current_x > (*bnd_next)->current_x && !slopes_equal(*(*bnd)->current_edge, *(*bnd_next)->current_edge)) { + mapbox::geometry::point pt; + if (!get_edge_intersection(*((*bnd)->current_edge), + *((*bnd_next)->current_edge), pt)) { + throw std::runtime_error("Edges do not intersect!"); + } + add_to_hot_pixels(round_point(pt), rings); + swap_positions_in_ABL(bnd, bnd_next, active_bounds); + bnd_next = std::next(bnd); + isModified = true; + } else { + bnd = bnd_next; + ++bnd_next; + } + } + } while (isModified); +} + +template +void process_hot_pixel_edges_at_top_of_scanbeam(T top_y, + scanbeam_list& scanbeam, + active_bound_list& active_bounds, + ring_manager& rings) { + for (auto bnd = active_bounds.begin(); bnd != active_bounds.end();) { + auto bnd_2 = std::next(bnd); + while ((*bnd)->current_edge != (*bnd)->edges.end() && (*bnd)->current_edge->top.y == top_y) { + add_to_hot_pixels((*bnd)->current_edge->top, rings); + if (current_edge_is_horizontal(bnd)) { + (*bnd)->current_x = static_cast((*bnd)->current_edge->top.x); + if ((*bnd)->current_edge->bot.x < (*bnd)->current_edge->top.x) { + // left to right + auto bnd_next = std::next(bnd); + while (bnd_next != active_bounds.end() && (*bnd_next)->current_x < (*bnd)->current_x) { + if (std::llround((*bnd_next)->current_edge->top.y) != top_y && std::llround((*bnd_next)->current_edge->bot.y) != top_y) { + mapbox::geometry::point pt(std::llround((*bnd_next)->current_x), top_y); + add_to_hot_pixels(pt, rings); + } + swap_positions_in_ABL(bnd, bnd_next, active_bounds); + bnd_next = std::next(bnd); + } + } else { + // right to left + if (bnd != active_bounds.begin()) { + auto bnd_prev = std::prev(bnd); + while (bnd != active_bounds.begin() && (*bnd_prev)->current_x > (*bnd)->current_x) { + if (std::llround((*bnd_prev)->current_edge->top.y) != top_y && std::llround((*bnd_prev)->current_edge->bot.y) != top_y) { + mapbox::geometry::point pt(std::llround((*bnd_prev)->current_x), top_y); + add_to_hot_pixels(pt, rings); + } + swap_positions_in_ABL(bnd, bnd_prev, active_bounds); + bnd_prev = std::prev(bnd); + } + } + } + } + next_edge_in_bound(bnd, scanbeam); + } + if ((*bnd)->current_edge == (*bnd)->edges.end()) { + active_bounds.erase(bnd); + } + bnd = bnd_2; + } +} + +template +void insert_local_minima_into_ABL_hot_pixel(T top_y, + local_minimum_ptr_list & minima_sorted, + local_minimum_ptr_list_itr & lm, + active_bound_list& active_bounds, + ring_manager& rings, + scanbeam_list& scanbeam) { + while (lm != minima_sorted.end() && (*lm)->y == top_y) { + if ((*lm)->left_bound.edges.empty() || (*lm)->right_bound.edges.empty()) { + ++lm; + continue; + } + + add_to_hot_pixels((*lm)->left_bound.edges.front().bot, rings); + auto& left_bound = (*lm)->left_bound; + left_bound.current_edge = left_bound.edges.begin(); + left_bound.current_x = static_cast(left_bound.current_edge->bot.x); + auto lb_abl_itr = insert_bound_into_ABL(left_bound, active_bounds); + if (!current_edge_is_horizontal(lb_abl_itr)) { + scanbeam.push((*lb_abl_itr)->current_edge->top.y); + } + auto& right_bound = (*lm)->right_bound; + right_bound.current_edge = right_bound.edges.begin(); + right_bound.current_x = static_cast(right_bound.current_edge->bot.x); + auto rb_abl_itr = insert_bound_into_ABL(right_bound, lb_abl_itr, active_bounds); + if (!current_edge_is_horizontal(rb_abl_itr)) { + scanbeam.push((*rb_abl_itr)->current_edge->top.y); + } + ++lm; + } +} + +template +void build_hot_pixels(local_minimum_list& minima_list, + ring_manager& rings) { + if (minima_list.empty()) { + return; + } + + active_bound_list active_bounds; + scanbeam_list scanbeam; + T scanline_y = std::numeric_limits::max(); + + local_minimum_ptr_list minima_sorted; + minima_sorted.reserve(minima_list.size()); + for (auto& lm : minima_list) { + minima_sorted.push_back(&lm); + } + std::stable_sort(minima_sorted.begin(), minima_sorted.end(), local_minimum_sorter()); + local_minimum_ptr_list_itr current_lm = minima_sorted.begin(); + + setup_scanbeam(minima_list, scanbeam); + + // Estimate size for reserving hot pixels + std::size_t reserve = 0; + for (auto & lm : minima_list) { + reserve += lm.left_bound.edges.size() + 2; + reserve += lm.right_bound.edges.size() + 2; + } + rings.hot_pixels.reserve(reserve); + + while (pop_from_scanbeam(scanline_y, scanbeam) || current_lm != minima_sorted.end()) { + + process_hot_pixel_intersections(scanline_y, active_bounds, rings); + + insert_local_minima_into_ABL_hot_pixel(scanline_y, minima_sorted, current_lm, active_bounds, + rings, scanbeam); + + process_hot_pixel_edges_at_top_of_scanbeam(scanline_y, scanbeam, active_bounds, rings); + + } + preallocate_point_memory(rings, rings.hot_pixels.size()); + sort_hot_pixels(rings); +} + +} +} +} diff --git a/mapbox/geometry/wagyu/topology_correction.hpp b/mapbox/geometry/wagyu/topology_correction.hpp new file mode 100644 index 0000000..d65e1ff --- /dev/null +++ b/mapbox/geometry/wagyu/topology_correction.hpp @@ -0,0 +1,1932 @@ +#pragma once + +#define _USE_MATH_DEFINES +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#ifdef DEBUG +#include +#endif + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +struct point_ptr_pair { + point_ptr op1; + point_ptr op2; + + constexpr point_ptr_pair(point_ptr o1, point_ptr o2) + : op1(o1), + op2(o2) {} + + point_ptr_pair(point_ptr_pair const& p) = default; + + point_ptr_pair(point_ptr_pair && p) + : op1(std::move(p.op1)), + op2(std::move(p.op2)) {} + + point_ptr_pair& operator=(point_ptr_pair && p) { + op1 = std::move(p.op1); + op2 = std::move(p.op2); + return *this; + } + +}; + +#ifdef DEBUG + +template +inline std::basic_ostream& +operator<<(std::basic_ostream& out, + const std::unordered_multimap, point_ptr_pair>& dupe_ring) { + + out << " BEGIN CONNECTIONS: " << std::endl; + for (auto& r : dupe_ring) { + out << " Ring: "; + if (r.second.op1->ring) { + out << r.second.op1->ring->ring_index; + } else { + out << "---"; + } + out << " to "; + if (r.second.op2->ring) { + out << r.second.op2->ring->ring_index; + } else { + out << "---"; + } + out << " ( at " << r.second.op1->x << ", " << r.second.op1->y << " )"; + out << " Ring1 ( "; + if (r.second.op1->ring) { + out << "area: " << r.second.op1->ring->area << " parent: "; + if (r.second.op1->ring->parent) { + out << r.second.op1->ring->parent->ring_index; + } else { + out << "---"; + } + } else { + out << "---"; + } + out << " )"; + out << " Ring2 ( "; + if (r.second.op2->ring) { + out << "area: " << r.second.op2->ring->area << " parent: "; + if (r.second.op2->ring->parent) { + out << r.second.op2->ring->parent->ring_index; + } else { + out << "---"; + } + } else { + out << "---"; + } + out << " )"; + out << std::endl; + } + out << " END CONNECTIONS: " << std::endl; + return out; +} + +#endif + +template +bool find_intersect_loop(std::unordered_multimap, point_ptr_pair>& dupe_ring, + std::list, point_ptr_pair>>& iList, + ring_ptr ring_parent, + ring_ptr ring_origin, + ring_ptr ring_search, + std::set>& visited, + point_ptr orig_pt, + point_ptr prev_pt, + ring_manager& rings) { + { + auto range = dupe_ring.equal_range(ring_search); + // Check for direct connection + for (auto& it = range.first; it != range.second;) { + ring_ptr it_ring1 = it->second.op1->ring; + ring_ptr it_ring2 = it->second.op2->ring; + if (!it_ring1 || !it_ring2 || it_ring1 != ring_search || + (!ring_is_hole(it_ring1) && !ring_is_hole(it_ring2))) { + it = dupe_ring.erase(it); + continue; + } + if (it_ring2 == ring_origin && + (ring_parent == it_ring2 || ring_parent == it_ring2->parent) && + *prev_pt != *it->second.op2 && *orig_pt != *it->second.op2) { + iList.emplace_front(ring_search, it->second); + return true; + } + ++it; + } + } + auto range = dupe_ring.equal_range(ring_search); + visited.insert(ring_search); + // Check for connection through chain of other intersections + for (auto& it = range.first; + it != range.second && it != dupe_ring.end() && it->first == ring_search; ++it) { + ring_ptr it_ring = it->second.op2->ring; + if (visited.count(it_ring) > 0 || it_ring == nullptr || + (ring_parent != it_ring && ring_parent != it_ring->parent) || + value_is_zero(area(it_ring)) || *prev_pt == *it->second.op2) { + continue; + } + if (find_intersect_loop(dupe_ring, iList, ring_parent, ring_origin, it_ring, visited, + orig_pt, it->second.op2, rings)) { + iList.emplace_front(ring_search, it->second); + return true; + } + } + return false; +} + +template +void remove_spikes(point_ptr& pt) { + ring_ptr r = pt->ring; + while (true) { + if (pt->next == pt) { + r->points = nullptr; + r->area = std::numeric_limits::quiet_NaN(); + pt->ring = nullptr; + pt = nullptr; + break; + } else if (*(pt) == *(pt->next)) { + point_ptr old_next = pt->next; + old_next->next->prev = pt; + pt->next = old_next->next; + old_next->next = old_next; + old_next->prev = old_next; + if (r->points == old_next) { + r->points = pt; + } + r->area = std::numeric_limits::quiet_NaN(); + old_next->ring = nullptr; + } else if (*(pt) == *(pt->prev)) { + point_ptr old_prev = pt->prev; + old_prev->prev->next = pt; + pt->prev = old_prev->prev; + old_prev->next = old_prev; + old_prev->prev = old_prev; + if (r->points == old_prev) { + r->points = pt; + } + r->area = std::numeric_limits::quiet_NaN(); + old_prev->ring = nullptr; + } else if (*(pt->next) == *(pt->prev)) { + point_ptr next = pt->next; + point_ptr prev = pt->prev; + next->prev = prev; + prev->next = next; + if (r->points == pt) { + r->points = prev; + } + r->area = std::numeric_limits::quiet_NaN(); + pt->ring = nullptr; + pt->next = pt; + pt->prev = pt; + pt = next; + } else { + break; + } + } +} + +template +void fixup_children(ring_ptr old_ring, ring_ptr new_ring) { + // Tests if any of the children from the old ring are now children of the new ring + assert(old_ring != new_ring); + for (auto r = old_ring->children.begin(); r != old_ring->children.end();) { + assert((*r)->points); + assert((*r) != old_ring); + if ((*r) != new_ring && !ring1_right_of_ring2(new_ring, (*r)) && + poly2_contains_poly1((*r)->points, new_ring->points)) { + (*r)->parent = new_ring; + new_ring->children.push_back((*r)); + r = old_ring->children.erase(r); + } else { + ++r; + } + } +} + +template +bool fix_intersects(std::unordered_multimap, point_ptr_pair>& dupe_ring, + point_ptr op_j, + point_ptr op_k, + ring_manager& rings, + mapbox::geometry::point& rewind_point) { + ring_ptr ring_j = op_j->ring; + ring_ptr ring_k = op_k->ring; + if (ring_j == ring_k) { + return false; + } + + if (!ring_is_hole(ring_j) && !ring_is_hole(ring_k)) { + // Both are not holes, return nothing to do. + return false; + } + + ring_ptr ring_origin; + ring_ptr ring_search; + ring_ptr ring_parent; + point_ptr op_origin_1; + point_ptr op_origin_2; + if (!ring_is_hole(ring_j)) { + ring_origin = ring_j; + ring_parent = ring_origin; + ring_search = ring_k; + op_origin_1 = op_j; + op_origin_2 = op_k; + } else if (!ring_is_hole(ring_k)) { + ring_origin = ring_k; + ring_parent = ring_origin; + ring_search = ring_j; + op_origin_1 = op_k; + op_origin_2 = op_j; + + } else { + // both are holes + // Order doesn't matter + ring_origin = ring_j; + ring_parent = ring_origin->parent; + ring_search = ring_k; + op_origin_1 = op_j; + op_origin_2 = op_k; + } + if (ring_parent != ring_search->parent) { + // The two holes do not have the same parent, do not add them + // simply return! + if (ring_parent->parent != ring_search && + poly2_contains_poly1(ring_search->points, ring_parent->points) && + !ring1_right_of_ring2(ring_search, ring_parent)) { + ring_ptr old_parent = ring_search->parent; + ring_search->parent = ring_parent; + old_parent->children.remove(ring_search); + ring_parent->children.push_back(ring_search); + } else { + return false; + } + } + bool found = false; + std::list, point_ptr_pair>> iList; + { + auto range = dupe_ring.equal_range(ring_search); + // Check for direct connection + for (auto& it = range.first; it != range.second;) { + if (!it->second.op1->ring) { + it = dupe_ring.erase(it); + continue; + } + if (!it->second.op2->ring) { + it = dupe_ring.erase(it); + continue; + } + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring2 == ring_origin) { + found = true; + if (*op_origin_1 != *(it->second.op2)) { + iList.emplace_back(ring_search, it->second); + break; + } + } + ++it; + } + } + if (iList.empty()) { + auto range = dupe_ring.equal_range(ring_search); + std::set> visited; + visited.insert(ring_search); + // Check for connection through chain of other intersections + for (auto& it = range.first; + it != range.second && it != dupe_ring.end() && it->first == ring_search; ++it) { + ring_ptr it_ring = it->second.op2->ring; + if (it_ring != ring_search && *op_origin_2 != *it->second.op2 && it_ring != nullptr && + (ring_parent == it_ring || ring_parent == it_ring->parent) && + !value_is_zero(area(it_ring)) && + find_intersect_loop(dupe_ring, iList, ring_parent, ring_origin, it_ring, visited, + op_origin_2, it->second.op2, rings)) { + found = true; + iList.emplace_front(ring_search, it->second); + break; + } + } + } + if (!found) { + point_ptr_pair intPt_origin = { op_origin_1, op_origin_2 }; + point_ptr_pair intPt_search = { op_origin_2, op_origin_1 }; + dupe_ring.emplace(ring_origin, std::move(intPt_origin)); + dupe_ring.emplace(ring_search, std::move(intPt_search)); + return false; + } + + if (iList.empty()) { + // The situation where both origin and search are holes might have a missing + // search condition, we must check if a new pair must be added. + bool missing = true; + auto rng = dupe_ring.equal_range(ring_origin); + // Check for direct connection + for (auto& it = rng.first; it != rng.second; ++it) { + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring2 == ring_search) { + missing = false; + } + } + if (missing) { + point_ptr_pair intPt_origin = { op_origin_1, op_origin_2 }; + dupe_ring.emplace(ring_origin, std::move(intPt_origin)); + } + return false; + } + + if (ring_is_hole(ring_origin)) { + for (auto& iRing : iList) { + ring_ptr ring_itr = iRing.first; + if (!ring_is_hole(ring_itr)) { + // Make the hole the origin! + point_ptr op1 = op_origin_1; + op_origin_1 = iRing.second.op1; + iRing.second.op1 = op1; + point_ptr op2 = op_origin_2; + op_origin_2 = iRing.second.op2; + iRing.second.op2 = op2; + iRing.first = ring_origin; + ring_origin = ring_itr; + ring_parent = ring_origin; + break; + } + } + } + + // Switch + point_ptr op_origin_1_next = op_origin_1->next; + point_ptr op_origin_2_next = op_origin_2->next; + op_origin_1->next = op_origin_2_next; + op_origin_2->next = op_origin_1_next; + op_origin_1_next->prev = op_origin_2; + op_origin_2_next->prev = op_origin_1; + + for (auto& iRing : iList) { + mapbox::geometry::point possible_rewind_point = find_rewind_point(iRing.second.op2); + if (possible_rewind_point.y > rewind_point.y || + (possible_rewind_point.y == rewind_point.y && + possible_rewind_point.x < rewind_point.x)) { + rewind_point.x = possible_rewind_point.x; + rewind_point.y = possible_rewind_point.y; + } + } + + for (auto& iRing : iList) { + point_ptr op_search_1 = iRing.second.op1; + point_ptr op_search_2 = iRing.second.op2; + point_ptr op_search_1_next = op_search_1->next; + point_ptr op_search_2_next = op_search_2->next; + op_search_1->next = op_search_2_next; + op_search_2->next = op_search_1_next; + op_search_1_next->prev = op_search_2; + op_search_2_next->prev = op_search_1; + } + + remove_spikes(op_origin_1); + remove_spikes(op_origin_2); + + if (op_origin_1 == nullptr || op_origin_2 == nullptr) { + if (op_origin_1 == nullptr && op_origin_2 == nullptr) { + // Self destruction! + ring_origin->points = nullptr; + ring_origin->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring_origin, rings); + for (auto& iRing : iList) { + ring_ptr ring_itr = iRing.first; + ring_itr->points = nullptr; + ring_itr->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring_itr, rings); + } + } else { + if (op_origin_1 == nullptr) { + ring_origin->points = op_origin_2; + } else { + //(op_origin_2 == nullptr) + ring_origin->points = op_origin_1; + } + ring_origin->area = std::numeric_limits::quiet_NaN(); + update_points_ring(ring_origin); + for (auto& iRing : iList) { + ring_ptr ring_itr = iRing.first; + ring_itr->points = nullptr; + ring_itr->area = std::numeric_limits::quiet_NaN(); + ring_itr->bottom_point = nullptr; + if (ring_is_hole(ring_origin)) { + ring1_replaces_ring2(ring_origin, ring_itr, rings); + } else { + ring1_replaces_ring2(ring_origin->parent, ring_itr, rings); + } + } + } + } else { + ring_ptr ring_new = create_new_ring(rings); + std::size_t size_1 = 0; + std::size_t size_2 = 0; + double area_1 = area_from_point(op_origin_1, size_1); + double area_2 = area_from_point(op_origin_2, size_2); + if (ring_is_hole(ring_origin) && ((area_1 < 0.0))) { + ring_origin->points = op_origin_1; + ring_origin->area = area_1; + ring_origin->size = size_1; + ring_new->points = op_origin_2; + ring_new->area = area_2; + ring_new->size = size_2; + } else { + ring_origin->points = op_origin_2; + ring_origin->area = area_2; + ring_origin->size = size_2; + ring_new->points = op_origin_1; + ring_new->area = area_1; + ring_new->size = size_1; + } + + update_points_ring(ring_origin); + update_points_ring(ring_new); + + ring_origin->bottom_point = nullptr; + + for (auto& iRing : iList) { + ring_ptr ring_itr = iRing.first; + ring_itr->points = nullptr; + ring_itr->area = std::numeric_limits::quiet_NaN(); + ring_itr->bottom_point = nullptr; + if (ring_is_hole(ring_origin)) { + ring1_replaces_ring2(ring_origin, ring_itr, rings); + } else { + ring1_replaces_ring2(ring_origin->parent, ring_itr, rings); + } + } + if (ring_is_hole(ring_origin)) { + ring_new->parent = ring_origin; + ring_new->parent->children.push_back(ring_new); + fixup_children(ring_origin, ring_new); + fixup_children(ring_parent, ring_new); + } else { + ring_new->parent = ring_origin->parent; + if (ring_new->parent == nullptr) { + rings.children.push_back(ring_new); + } else { + ring_new->parent->children.push_back(ring_new); + } + fixup_children(ring_origin, ring_new); + } + } + + std::list, point_ptr_pair>> move_list; + + for (auto& iRing : iList) { + auto range_itr = dupe_ring.equal_range(iRing.first); + if (range_itr.first != range_itr.second) { + for (auto& it = range_itr.first; it != range_itr.second; ++it) { + ring_ptr it_ring = it->second.op1->ring; + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring == nullptr || it_ring2 == nullptr || it_ring == it_ring2) { + continue; + } + if ((ring_is_hole(it_ring) || ring_is_hole(it_ring2))) { + move_list.emplace_back(it_ring, it->second); + } + } + dupe_ring.erase(iRing.first); + } + } + + auto range_itr = dupe_ring.equal_range(ring_origin); + for (auto& it = range_itr.first; it != range_itr.second;) { + ring_ptr it_ring = it->second.op1->ring; + ring_ptr it_ring2 = it->second.op2->ring; + if (it_ring == nullptr || it_ring2 == nullptr || it_ring == it_ring2) { + it = dupe_ring.erase(it); + continue; + } + if (it_ring != ring_origin) { + if ((ring_is_hole(it_ring) || ring_is_hole(it_ring2))) { + move_list.emplace_back(it_ring, it->second); + } + it = dupe_ring.erase(it); + } else { + if ((ring_is_hole(it_ring) || ring_is_hole(it_ring2))) { + ++it; + } else { + it = dupe_ring.erase(it); + } + } + } + + if (!move_list.empty()) { + dupe_ring.insert(move_list.begin(), move_list.end()); + } + + return true; +} + +template +struct point_ptr_cmp { + inline bool operator()(point_ptr op1, point_ptr op2) { + if (op1->y != op2->y) { + return (op1->y > op2->y); + } else if (op1->x != op2->x) { + return (op1->x < op2->x); + } else { + std::size_t depth_1 = ring_depth(op1->ring); + std::size_t depth_2 = ring_depth(op2->ring); + return depth_1 > depth_2; + } + } +}; + +template +struct point_ptr_depth_cmp { + inline bool operator()(point_ptr op1, point_ptr op2) { + std::size_t depth_1 = ring_depth(op1->ring); + std::size_t depth_2 = ring_depth(op2->ring); + return depth_1 > depth_2; + } +}; + +template +void update_duplicate_point_entries( + ring_ptr ring, std::unordered_multimap, point_ptr_pair>& dupe_ring) { + auto range = dupe_ring.equal_range(ring); + std::list, point_ptr_pair>> move_list; + for (auto& it = range.first; it != range.second;) { + ring_ptr it_ring = it->second.op1->ring; + ring_ptr it_ring_2 = it->second.op2->ring; + if (it_ring == nullptr || it_ring_2 == nullptr) { + it = dupe_ring.erase(it); + continue; + } + if (it_ring != ring) { + if ((ring_is_hole(it_ring) || ring_is_hole(it_ring_2))) { + move_list.emplace_back(it_ring, it->second); + } + it = dupe_ring.erase(it); + } else { + if ((ring_is_hole(it_ring) || ring_is_hole(it_ring_2))) { + ++it; + } else { + it = dupe_ring.erase(it); + } + } + } + if (!move_list.empty()) { + dupe_ring.insert(move_list.begin(), move_list.end()); + } +} + +template +bool parent_in_tree(ring_ptr r, ring_ptr possible_parent) { + + ring_ptr current_ring = r->parent; + while (current_ring != nullptr) { + if (current_ring == possible_parent) { + return true; + } + current_ring = current_ring->parent; + } + return false; +} + +template +void fixup_children_new_interior_ring(ring_ptr old_ring, + ring_ptr new_ring, + ring_manager& rings) { + bool old_ring_area_is_positive = area(old_ring) > 0.0; + // Now we must search the siblings of the old ring + // This is slow, but there is no other way I know of currently + // to solve these problems. + if (old_ring->parent == nullptr) { + for (auto r = rings.children.begin(); r != rings.children.end();) { + assert((*r)->points); + bool ring_area_is_positive = area((*r)) > 0.0; + if ((*r) != new_ring && ring_area_is_positive == old_ring_area_is_positive && + poly2_contains_poly1((*r)->points, new_ring->points)) { + (*r)->parent = new_ring; + new_ring->children.push_back((*r)); + r = rings.children.erase(r); + } else { + ++r; + } + } + } else { + ring_ptr parent = old_ring->parent; + for (auto r = parent->children.begin(); r != parent->children.end();) { + assert((*r)->points); + assert((*r) != parent); + bool ring_area_is_positive = area((*r)) > 0.0; + if ((*r) != new_ring && ring_area_is_positive == old_ring_area_is_positive && + poly2_contains_poly1((*r)->points, new_ring->points)) { + (*r)->parent = new_ring; + new_ring->children.push_back((*r)); + r = parent->children.erase(r); + } else { + ++r; + } + } + } +} + +#ifdef DEBUG + +template +void check_if_intersections_cross(point_ptr p1, point_ptr p2) { + // LCOV_EXCL_START + point_ptr p1_next = p1->next; + point_ptr p2_next = p2->next; + point_ptr p1_prev = p1->prev; + point_ptr p2_prev = p2->prev; + while (*p1_next == *p1) { + if (p1_next == p1) { + return; + } + p1_next = p1_next->next; + } + while (*p2_next == *p2) { + if (p2_next == p2) { + return; + } + p2_next = p2_next->next; + } + while (*p1_prev == *p1) { + if (p1_prev == p1) { + return; + } + p1_prev = p1_prev->prev; + } + while (*p2_prev == *p2) { + if (p2_prev == p2) { + return; + } + p2_prev = p2_prev->prev; + } + double a1_p1 = std::atan2(static_cast(p1_prev->y - p1->y), + static_cast(p1_prev->x - p1->x)); + double a2_p1 = std::atan2(static_cast(p1_next->y - p1->y), + static_cast(p1_next->x - p1->x)); + double a1_p2 = std::atan2(static_cast(p2_prev->y - p2->y), + static_cast(p2_prev->x - p2->x)); + double a2_p2 = std::atan2(static_cast(p2_next->y - p2->y), + static_cast(p2_next->x - p2->x)); + double min_p1 = std::min(a1_p1, a2_p1); + double max_p1 = std::max(a1_p1, a2_p1); + double min_p2 = std::min(a1_p2, a2_p2); + double max_p2 = std::max(a1_p2, a2_p2); + if ((min_p1 < max_p2 && min_p1 > min_p2 && max_p1 > max_p2) || + (min_p2 < max_p1 && min_p2 > min_p1 && max_p2 > max_p1)) { + throw std::runtime_error("Paths are found to be crossing"); + } + // LCOV_EXCL_END +} + +#endif + +template +void handle_self_intersections(point_ptr op, + point_ptr op2, + std::unordered_multimap, point_ptr_pair>& dupe_ring, + ring_manager& rings) { + // Check that are same ring + assert(op->ring == op2->ring); + ring_ptr ring = op->ring; + double original_area = area(ring); + bool original_is_positive = (original_area > 0.0); +#ifdef DEBUG + check_if_intersections_cross(op, op2); +#endif + + // split the polygon into two ... + point_ptr op3 = op->prev; + point_ptr op4 = op2->prev; + op->prev = op4; + op4->next = op; + op2->prev = op3; + op3->next = op2; + + remove_spikes(op); + remove_spikes(op2); + + if (op == nullptr && op2 == nullptr) { + // Self destruction! + // I am not positive that it could ever reach this point -- but leaving + // the logic in here for now. + ring->points = nullptr; + ring->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring, rings); + update_duplicate_point_entries(ring, dupe_ring); + return; + } else if (op == nullptr) { + ring->points = op2; + ring->area = std::numeric_limits::quiet_NaN(); + update_duplicate_point_entries(ring, dupe_ring); + return; + } else if (op2 == nullptr) { + ring->points = op; + ring->area = std::numeric_limits::quiet_NaN(); + update_duplicate_point_entries(ring, dupe_ring); + return; + } + + ring_ptr new_ring = create_new_ring(rings); + std::size_t size_1 = 0; + std::size_t size_2 = 0; + double area_1 = area_from_point(op, size_1); + double area_2 = area_from_point(op2, size_2); + bool area_1_is_positive = (area_1 > 0.0); + bool area_2_is_positive = (area_2 > 0.0); + bool area_1_is_zero = value_is_zero(area_1); + bool area_2_is_zero = value_is_zero(area_2); + + // Situation # 1 - Orientations are NOT the same: + // - One ring contains the other and MUST be a child of that ring + // - The one that changed orientation is the child of the other ring + // + // Situation # 2 - Orientations are the same + // - The rings are now split, such a new ring of the same orientation + // must be created. + // - If the new ring is WITHIN the old ring: + // * It WILL be the child of a hole of that ring (this ring may not yet be created) + // or possible the child of a child of a child of the ring (an so on)... + // - If the new ring is OUTSIDE the old ring: + // * It may contain any of the children of the old ring. + if (area_2_is_zero || area_1_is_zero || area_1_is_positive != area_2_is_positive) { + // Situation #1 - new_ring is contained by ring ... + if (area_2_is_zero || (!area_1_is_zero && area_1_is_positive == original_is_positive)) { + ring->points = op; + ring->area = area_1; + ring->size = size_1; + new_ring->points = op2; + new_ring->area = area_2; + new_ring->size = size_2; + } else { + ring->points = op2; + ring->area = area_2; + ring->size = size_2; + new_ring->points = op; + new_ring->area = area_1; + new_ring->size = size_1; + } + update_points_ring(ring); + update_points_ring(new_ring); + new_ring->parent = ring; + new_ring->parent->children.push_back(new_ring); + fixup_children_new_interior_ring(ring, new_ring, rings); + } else { + // Situation #2 - create new ring + // The largest absolute area is the parent + if (std::fabs(area_1) > std::fabs(area_2)) { + ring->points = op; + ring->area = area_1; + ring->size = size_1; + new_ring->points = op2; + new_ring->area = area_2; + new_ring->size = size_2; + } else { + ring->points = op2; + ring->area = area_2; + ring->size = size_2; + new_ring->points = op; + new_ring->area = area_1; + new_ring->size = size_1; + } + update_points_ring(ring); + update_points_ring(new_ring); + if (poly2_contains_poly1(new_ring->points, ring->points)) { + // This is the situation where there is the new ring is + // created inside the ring. Later on this should be inherited + // as child of a newly created hole. However, we should check existing + // holes of this polygon to see if they might belong inside this polygon. + new_ring->parent = ring; + new_ring->parent->children.push_back(new_ring); + fixup_children(ring, new_ring); + } else { + // Polygons are completely seperate + new_ring->parent = ring->parent; + if (new_ring->parent == nullptr) { + rings.children.push_back(new_ring); + } else { + new_ring->parent->children.push_back(new_ring); + } + fixup_children(ring, new_ring); + } + } + update_duplicate_point_entries(ring, dupe_ring); +} + +template +mapbox::geometry::point find_rewind_point(point_ptr pt) { + mapbox::geometry::point rewind; + rewind.x = pt->x; + rewind.y = pt->y; + point_ptr itr = pt->next; + while (pt != itr) { + if (itr->y > rewind.y || (itr->y == rewind.y && itr->x < rewind.x)) { + rewind.x = itr->x; + rewind.y = itr->y; + } + itr = itr->next; + } + return rewind; +} + +template +bool handle_collinear_edges(point_ptr pt1, + point_ptr pt2, + std::unordered_multimap, point_ptr_pair>& dupe_ring, + ring_manager& rings, + mapbox::geometry::point& rewind_point) { + ring_ptr ring1 = pt1->ring; + ring_ptr ring2 = pt2->ring; + if (ring1 == ring2) { + return false; + } + + bool valid = (ring1 != ring2 && (ring1->parent == ring2->parent || ring2->parent == ring1 || + ring1->parent == ring2)); + if (!valid) { + return false; + } + + if (*(pt1->next) != *(pt2->prev) && *(pt2->next) != *(pt1->prev)) { + return false; + } + + if (ring1->parent == ring2) { + // switch ring1 and ring2 + std::swap(pt1, pt2); + std::swap(ring1, ring2); + } + + mapbox::geometry::point rewind_1 = find_rewind_point(pt1); + mapbox::geometry::point rewind_2 = find_rewind_point(pt2); + + // The lower right of the two points is the rewind point. + mapbox::geometry::point possible_rewind; + if (rewind_1.y > rewind_2.y) { + possible_rewind = rewind_2; + } else if (rewind_1.y < rewind_2.y) { + possible_rewind = rewind_1; + } else if (rewind_1.x > rewind_2.x) { + possible_rewind = rewind_1; + } else { + possible_rewind = rewind_2; + } + if (possible_rewind.y > rewind_point.y || + (possible_rewind.y == rewind_point.y && possible_rewind.x < rewind_point.x)) { + rewind_point.x = possible_rewind.x; + rewind_point.y = possible_rewind.y; + } + + // swap points + point_ptr pt3 = pt1->prev; + point_ptr pt4 = pt2->prev; + pt1->prev = pt4; + pt4->next = pt1; + pt2->prev = pt3; + pt3->next = pt2; + + // remove spikes + remove_spikes(pt1); + if (!pt1) { + // rings self destructed + ring1->points = nullptr; + ring1->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring1, rings); + ring2->points = nullptr; + ring2->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring2, rings); + return false; + } + if (pt2->ring) { + remove_spikes(pt2); + if (!pt2) { + // rings self destructed + // Not sure this logic is reachable, but leaving it in for now. + ring1->points = nullptr; + ring1->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring1, rings); + ring2->points = nullptr; + ring2->area = std::numeric_limits::quiet_NaN(); + remove_ring(ring2, rings); + return false; + } + // We could have removed pt1 during this process.. + // if we did we need to reassign pt2 to pt1 + if (!pt1->ring) { + pt1 = pt2; + } + } + ring1->points = pt1; + ring2->points = nullptr; + ring1->area = std::numeric_limits::quiet_NaN(); + ring2->area = std::numeric_limits::quiet_NaN(); + if (ring2->parent == ring1) { + ring1_replaces_ring2(ring1->parent, ring2, rings); + } else { + ring1_replaces_ring2(ring1, ring2, rings); + } + update_points_ring(ring1); + update_duplicate_point_entries(ring2, dupe_ring); + + return true; +} + +template +bool point_2_is_between_point_1_and_point_3(point_ptr pt1, point_ptr pt2, point_ptr pt3) { + if ((*pt1 == *pt3) || (*pt1 == *pt2) || (*pt3 == *pt2)) { + return false; + } else if (pt1->x != pt3->x) { + return (pt2->x > pt1->x) == (pt2->x < pt3->x); + } else { + return (pt2->y > pt1->y) == (pt2->y < pt3->y); + } +} + +enum orientation_type : std::uint8_t { + orientation_collinear_spike = 0, + orientation_clockwise, + orientation_collinear_line, + orientation_counter_clockwise +}; + +// To find orientation of ordered triplet (p, q, r) +// Orientation between q and r with respect to p. +template +inline orientation_type orientation_of_points(point_ptr p, point_ptr q, point_ptr r) { + T val = (q->y - p->y) * (r->x - q->x) - (q->x - p->x) * (r->y - q->y); + if (val == 0) { + if (point_2_is_between_point_1_and_point_3(q, p, r)) { + return orientation_collinear_line; + } else { + return orientation_collinear_spike; + } + } + return (val > 0) ? orientation_clockwise : orientation_counter_clockwise; +} + +// Self intersection point vector +template +using si_point_vector = std::vector, bool>>; + +#ifdef DEBUG + +template +std::string output_si_angles(point_ptr pt) { + // LCOV_EXCL_START + std::ostringstream out; + double prev_angle = std::atan2(static_cast(pt->prev->y - pt->y), + static_cast(pt->prev->x - pt->x)); + prev_angle = (180.0 / M_PI) * prev_angle; + if (prev_angle < 0.0) { + prev_angle += 360.0; + } + double next_angle = std::atan2(static_cast(pt->next->y - pt->y), + static_cast(pt->next->x - pt->x)); + next_angle = (180.0 / M_PI) * next_angle; + if (next_angle < 0.0) { + next_angle += 360.0; + } + out << " angles: " << prev_angle << " , " << next_angle; + out << " - [[" << pt->prev->x << "," << pt->prev->y << "],["; + out << pt->x << "," << pt->y << "],["; + out << pt->next->x << "," << pt->next->y << "]]"; + return out.str(); + // LCOV_EXCL_END +} + +template +std::string output_compare_si_angles(point_ptr pt, point_ptr compare) { + std::ostringstream out; + double cmp_prev_angle = std::atan2(static_cast(compare->prev->y - compare->y), + static_cast(compare->prev->x - compare->x)); + double prev_angle = std::atan2(static_cast(pt->prev->y - pt->y), + static_cast(pt->prev->x - pt->x)); + prev_angle = (180.0 / M_PI) * prev_angle; + if (prev_angle < 0.0) { + prev_angle += 360.0; + } + cmp_prev_angle = (180.0 / M_PI) * cmp_prev_angle; + if (cmp_prev_angle < 0.0) { + cmp_prev_angle += 360.0; + } + double cmp_next_angle = std::atan2(static_cast(compare->next->y - compare->y), + static_cast(compare->next->x - compare->x)); + double next_angle = std::atan2(static_cast(pt->next->y - pt->y), + static_cast(pt->next->x - pt->x)); + next_angle = (180.0 / M_PI) * next_angle; + if (next_angle < 0.0) { + next_angle += 360.0; + } + cmp_next_angle = (180.0 / M_PI) * cmp_next_angle; + if (cmp_next_angle < 0.0) { + cmp_next_angle += 360.0; + } + out << " compared to prev: " << prev_angle - cmp_prev_angle << ", " + << next_angle - cmp_prev_angle << std::endl; + out << " compared to next: " << prev_angle - cmp_next_angle << ", " + << next_angle - cmp_next_angle << std::endl; + return out.str(); +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const si_point_vector& pts) { + out << "Self Intersection Point Vector:" << std::endl; + for (auto& pt : pts) { + out << output_si_angles(pt.first); + if (pt.second) { + out << " - clockwise" << std::endl; + } else { + out << " - counter clockwise" << std::endl; + } + } + return out; +} + +template +inline std::basic_ostream& operator<<(std::basic_ostream& out, + const orientation_type& ot) { + switch (ot) { + case orientation_collinear_spike: + out << "collinear - spike"; + break; + case orientation_clockwise: + out << "clockwise"; + break; + case orientation_collinear_line: + out << "collinear - line"; + break; + case orientation_counter_clockwise: + out << "counter clockwise"; + break; + }; + return out; +} + +#endif + +template +bool clockwise_of_next(point_ptr const& origin, point_ptr pt) { + + // Determine if the prev and next of pt is either clockwise of next from + // origin (and therefore ccw of previous) or if it is ccw of next and clockwise + // of previous + + // First inspect orientation of points, keep in mind this is the orientations of + // the three points + orientation_type ot_origin = orientation_of_points(origin, origin->next, origin->prev); + if (ot_origin == orientation_collinear_spike) { + return true; + } else if (ot_origin == orientation_clockwise) { + orientation_type ot_prev_next = orientation_of_points(origin, origin->next, pt->prev); + if (ot_prev_next == orientation_collinear_spike) { + orientation_type ot_next_next = orientation_of_points(origin, origin->next, pt->next); + if (ot_next_next == orientation_collinear_spike) { + // Pt forms a spike on origin next + // so lets assume it is CW. + // I am not sure that we would ever encounter a spike so this will + // be missing from coverage, but left in code. + return true; + } else if (ot_next_next == orientation_clockwise) { + // We need to check which is after "origin->next" traveling + // clockwise -- origin->prev or pt->next + orientation_type ot_next_prev = + orientation_of_points(origin, origin->prev, pt->next); + if (ot_next_prev == orientation_collinear_spike) { + // Both origin and pt follow the same paths. + // so we will call this clockwise + return true; + } else if (ot_next_prev == orientation_clockwise) { + // pt->next is clockwise of prev + return false; + } else if (ot_next_prev == orientation_collinear_line) { + // This shouldn't happen + // LCOV_EXCL_START + throw std::runtime_error("Impossible situation reached in clockwise_of_next"); + // LCOV_EXCL_END + } else { + // Pt next is counter clockwise of origin prev + return true; + } + } else { + // ot_next_next == orientation_collinear_line + // ot_next_next == orientation_counter_clockwise + return false; + } + } else if (ot_prev_next == orientation_clockwise) { + // We need to check which is after "origin->next" traveling + // clockwise -- origin->prev or pt->prev + orientation_type ot_prev_prev = orientation_of_points(origin, origin->prev, pt->prev); + if (ot_prev_prev == orientation_clockwise || + ot_prev_prev == orientation_collinear_spike) { + // pt->prev is before this, so.. between the two. + return false; + } else { + // ot_prev_prev == orientation_clockwise + // ot_prev_prev == orientation_collinear_line + return true; + } + } else { + // ot_prev_next == orientation_collinear_line + // ot_prev_next == orientation_counter_clockwise + return false; + } + } else if (ot_origin == orientation_collinear_line) { + orientation_type ot_prev_next = orientation_of_points(origin, origin->next, pt->prev); + if (ot_prev_next == orientation_collinear_spike || + ot_prev_next == orientation_collinear_line) { + // prev and next on top of each other + orientation_type ot_next_next = orientation_of_points(origin, origin->next, pt->next); + if (ot_next_next == orientation_counter_clockwise) { + return false; + } else { + // ot_next_next == orientation_collinear_spike + // ot_next_next == orientation_clockwise + // ot_next_next == orientation_collinear_line + return true; + } + } else if (ot_prev_next == orientation_clockwise) { + return true; + } else { + // ot_prev_next == orientation_counter_clockwise + return false; + } + } else { + // orientation_counter_clockwise + orientation_type ot_prev_next = orientation_of_points(origin, origin->next, pt->prev); + if (ot_prev_next == orientation_collinear_spike) { + orientation_type ot_next_next = orientation_of_points(origin, origin->next, pt->next); + if (ot_next_next == orientation_collinear_spike) { + // Pt forms a spike on origin next + // so lets assume it is CW. + // I am not sure that we would ever encounter a spike so this will + // be missing from coverage, but left in code. + return true; + } else if (ot_next_next == orientation_counter_clockwise) { + // We need to check which is after "origin->next" traveling + // counter clockwise -- origin->prev or pt->next + orientation_type ot_next_prev = + orientation_of_points(origin, origin->prev, pt->next); + if (ot_next_prev == orientation_collinear_spike) { + // Both origin and pt follow the same paths. + // so we will call this clockwise + return true; + } else if (ot_next_prev == orientation_clockwise) { + // pt->next is clockwise of prev + return false; + } else if (ot_next_prev == orientation_collinear_line) { + // This shouldn't happen? + // LCOV_EXCL_START + throw std::runtime_error("Impossible situation reached in clockwise_of_next - 2"); + // LCOV_EXCL_END + } else { + // Pt next is counter clockwise of origin prev + return true; + } + } else { + // ot_next_next == orientation_clockwise + // ot_next_next == orientation_collinear_line + return true; + } + } else if (ot_prev_next == orientation_counter_clockwise) { + // We need to check which is after "origin->next" traveling + // counter clockwise -- origin->prev or pt->prev + orientation_type ot_prev_prev = orientation_of_points(origin, origin->prev, pt->prev); + if (ot_prev_prev == orientation_clockwise || + ot_prev_prev == orientation_collinear_spike) { + // pt->prev is before this, so.. between the two. + return false; + } else { + // ot_prev_prev == orientation_counter_clockwise + // ot_prev_prev == orientation_collinear_line + return true; + } + } else { + // ot_prev_next == orientation_collinear_line + // ot_prev_next == orientation_clockwise + return true; + } + } +} + +template +inline bool cw_p1p2_prev_collinear_spike(point_ptr const& origin, + point_ptr const& next, + point_ptr const& p1, + point_ptr const& p2) { + + // we must compare the nexts to determine the order between the two. + orientation_type ot_p1_next = orientation_of_points(origin, next, p1->next); + orientation_type ot_p2_next = orientation_of_points(origin, next, p2->next); + if (ot_p1_next == orientation_collinear_spike) { + // p1 is a spike on origin next + if (ot_p2_next == orientation_collinear_spike) { + // I am not sure that a spike will ever occur at this point in the processing! + return true; + } else { + return false; + } + } else if (ot_p1_next == orientation_clockwise) { + if (ot_p2_next == orientation_collinear_spike) { + return true; + } else if (ot_p2_next == orientation_clockwise) { + // Both are clockwise we have to compare. + orientation_type ot = orientation_of_points(origin, p1->next, p2->next); + if (ot == orientation_collinear_spike) { + return true; + } else if (ot == orientation_clockwise) { + return false; + } else if (ot == orientation_collinear_line) { + return false; + } else { + return true; + } + } else { + // ot_p2_next == orienation_collinear_line + // ot_p2_next == orienation_counter_clockwise + return false; + } + } else if (ot_p1_next == orientation_collinear_line) { + if (ot_p2_next == orientation_counter_clockwise) { + return false; + } else { + // ot_p2_next == orienation_collinear_spike + // ot_p2_next == orienation_clockwise + // ot_p2_next == orienation_collinear_line + return true; + } + } else { + // ot_p1_next == orientation_counter_clockwise + if (ot_p2_next == orientation_counter_clockwise) { + // Both are counter clockwise we have to compare. + orientation_type ot = orientation_of_points(origin, p1->next, p2->next); + if (ot == orientation_collinear_spike) { + return true; + } else if (ot == orientation_clockwise) { + return false; + } else if (ot == orientation_collinear_line) { + return false; + } else { + return true; + } + } else { + // ot_p2_next == orienation_collinear_spike + // ot_p2_next == orienation_clockwise + // ot_p2_next == orienation_collinear_line + return true; + } + } +} + +template +inline bool ccw_p1p2_prev_collinear_spike(point_ptr const& origin, + point_ptr const& next, + point_ptr const& p1, + point_ptr const& p2) { + + // we must compare the nexts to determine the order between the two. + orientation_type ot_p1_next = orientation_of_points(origin, next, p1->next); + orientation_type ot_p2_next = orientation_of_points(origin, next, p2->next); + if (ot_p1_next == orientation_collinear_spike) { + // p1 is a spike on origin next + if (ot_p2_next == orientation_collinear_spike) { + return true; + } else { + return false; + } + } else if (ot_p1_next == orientation_clockwise) { + if (ot_p2_next == orientation_collinear_spike) { + return false; + } else if (ot_p2_next == orientation_clockwise) { + // Both are clockwise we have to compare. + orientation_type ot = orientation_of_points(origin, p1->next, p2->next); + if (ot == orientation_collinear_spike) { + return true; + } else if (ot == orientation_clockwise) { + return true; + } else if (ot == orientation_collinear_line) { + return true; + } else { + return false; + } + } else { + // ot_p2_next == orienation_collinear_line + // ot_p2_next == orienation_counter_clockwise + return true; + } + } else if (ot_p1_next == orientation_collinear_line) { + if (ot_p2_next == orientation_clockwise) { + return false; + } else { + // ot_p2_next == orienation_collinear_spike + // ot_p2_next == orienation_counter_clockwise + // ot_p2_next == orienation_collinear_line + return true; + } + } else { + // ot_p1_next == orientation_counter_clockwise + if (ot_p2_next == orientation_counter_clockwise) { + // Both are counter clockwise we have to compare. + orientation_type ot = orientation_of_points(origin, p1->next, p2->next); + if (ot == orientation_collinear_spike) { + return true; + } else if (ot == orientation_clockwise) { + return true; + } else if (ot == orientation_collinear_line) { + return true; + } else { + return false; + } + } else { + // ot_p2_next == orienation_collinear_spike + // ot_p2_next == orienation_clockwise + // ot_p2_next == orienation_collinear_line + return false; + } + } +} + +template +struct si_point_sorter { + + point_ptr origin; + point_ptr next; + + si_point_sorter(point_ptr origin_) : origin(origin_), next(origin_->next) { + } + + // Sorting order + // + // Primary Sort: + // * Left of next, right of previous + // * Right of next, left of previous + // + // Secondary Sort: + // * Magnitude of angle (direction based on primary sort) + // between item's previous and origin's next + + inline bool operator()(std::pair, bool> const& pp1, + std::pair, bool> const& pp2) { + // Because a next must be paired with a previous, we are only + // caring first about previous segments for ordering + point_ptr p1 = pp1.first; + point_ptr p2 = pp2.first; + if (pp1.second != pp2.second) { + return pp1.second; + } + + orientation_type ot_p1 = orientation_of_points(origin, next, p1->prev); + orientation_type ot_p2 = orientation_of_points(origin, next, p2->prev); + if (pp1.second) { + if (ot_p1 == orientation_collinear_spike) { + // p1 prev lines up with origin next + if (ot_p2 != orientation_collinear_spike) { + // ot_p2 == orienation_clockwise + // ot_p2 == orienation_collinear_line + // ot_p2 == orienation_counter_clockwise + return true; + } else { + // ot_p2 == orientation_collinear_spike + return cw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } + } else if (ot_p1 == orientation_clockwise) { + if (ot_p2 == orientation_collinear_spike) { + return false; + } else if (ot_p2 == orientation_clockwise) { + orientation_type ot_p1p2 = orientation_of_points(origin, p1->prev, p2->prev); + if (ot_p1p2 == orientation_collinear_spike) { + // Both p1 prev and p2 prev are on top of each other + return cw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } else if (ot_p1p2 == orientation_clockwise) { + return true; + } else { + // ot_p1p2 == orientation_counter_clockwise + // ot_p1p2 == orientation_collinear_line + return false; + } + } else { + // ot_p2 == orientation_collinear_line + // ot_p2 == orientation_counter_clockwise + return true; + } + } else if (ot_p1 == orientation_collinear_line) { + if (ot_p2 == orientation_collinear_spike || ot_p2 == orientation_clockwise) { + return false; + } else if (ot_p2 == orientation_collinear_line) { + return cw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } else { + // ot_p2 == orientation_counter_clockwise + return true; + } + } else { + // ot_p1 == orientation_counter_clockwise + if (ot_p2 == orientation_counter_clockwise) { + orientation_type ot_p1p2 = orientation_of_points(origin, p1->prev, p2->prev); + if (ot_p1p2 == orientation_collinear_spike) { + // Both p1 prev and p2 prev are on top of each other + return cw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } else if (ot_p1p2 == orientation_clockwise) { + return true; + } else { + // ot_p1p2 == orientation_counter_clockwise + // ot_p1p2 == orientation_collinear_line + return false; + } + } else { + // ot_p2 == orientation_collinear_spike + // ot_p2 == orientation_counter_clockwise + // ot_p2 == orientation_collinear_line + return false; + } + } + } else { + if (ot_p1 == orientation_collinear_spike) { + // p1 prev lines up with origin next + if (ot_p2 != orientation_collinear_spike) { + // ot_p2 == orienation_clockwise + // ot_p2 == orienation_collinear_line + // ot_p2 == orienation_counter_clockwise + return true; + } else { + // ot_p2 == orientation_collinear_spike + return ccw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } + } else if (ot_p1 == orientation_clockwise) { + if (ot_p2 == orientation_collinear_spike) { + return false; + } else if (ot_p2 == orientation_clockwise) { + orientation_type ot_p1p2 = orientation_of_points(origin, p1->prev, p2->prev); + if (ot_p1p2 == orientation_collinear_spike) { + // Both p1 prev and p2 prev are on top of each other + return ccw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } else if (ot_p1p2 == orientation_counter_clockwise) { + return true; + } else { + // ot_p1p2 == orientation_clockwise + // ot_p1p2 == orientation_collinear_line + return false; + } + } else { + // ot_p2 == orientation_collinear_line + // ot_p2 == orientation_counter_clockwise + return false; + } + } else if (ot_p1 == orientation_collinear_line) { + if (ot_p2 == orientation_collinear_spike || + ot_p2 == orientation_counter_clockwise) { + return false; + } else if (ot_p2 == orientation_collinear_line) { + return ccw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } else { + // ot_p2 == orientation_clockwise + return true; + } + } else { + // ot_p1 == orientation_counter_clockwise + if (ot_p2 == orientation_collinear_spike) { + return false; + } else if (ot_p2 == orientation_counter_clockwise) { + orientation_type ot_p1p2 = orientation_of_points(origin, p1->prev, p2->prev); + if (ot_p1p2 == orientation_collinear_spike) { + // Both p1 prev and p2 prev are on top of each other + return ccw_p1p2_prev_collinear_spike(origin, next, p1, p2); + } else if (ot_p1p2 == orientation_counter_clockwise) { + return true; + } else { + // ot_p1p2 == orientation_clockwise + // ot_p1p2 == orientation_collinear_line + return false; + } + } else { + // ot_p2 == orientation_clockwise + // ot_p2 == orientation_collinear_line + return true; + } + } + } + } +}; + +template +si_point_vector build_si_point_vector(std::size_t first_index, + std::size_t last_index, + std::size_t current_index, + ring_ptr match_ring, + ring_manager& rings) { + si_point_vector point_vec; + point_ptr origin = rings.all_points[current_index]; + for (std::size_t j = first_index; j <= last_index; ++j) { + if (j == current_index) { + continue; + } + point_ptr op_j = rings.all_points[j]; + if (op_j->ring == match_ring) { + bool clockwise = clockwise_of_next(origin, op_j); + point_vec.emplace_back(op_j, clockwise); + } + } + return point_vec; +} + +template +bool process_repeated_point_set(std::size_t first_index, + std::size_t last_index, + std::size_t current_index, + std::unordered_multimap, point_ptr_pair>& dupe_ring, + ring_manager& rings) { + point_ptr point_1 = rings.all_points[current_index]; + + if (point_1->ring == nullptr) { + return false; + } + + // Build point vector + si_point_vector vec = + build_si_point_vector(first_index, last_index, current_index, point_1->ring, rings); + + if (vec.empty()) { + return false; + } + + // Sort points in vector + std::stable_sort(vec.begin(), vec.end(), si_point_sorter(point_1)); + + auto vec_itr = vec.begin(); + point_ptr point_2 = vec_itr->first; + + // If there are collinear sets of lines, we might not be able to just pick + // the first point in the sorted list (and we will have to do special processing). + if (vec.size() > 2) { + ++vec_itr; + point_ptr point_3 = vec_itr->first; + orientation_type ot_next = orientation_of_points(point_2, point_2->next, point_3->next); + if (ot_next == orientation_collinear_spike) { + orientation_type ot_prev = orientation_of_points(point_2, point_2->prev, point_3->prev); + if (ot_prev == orientation_collinear_spike) { + // Start at point_1 and we will travel around the circle until we find another point at + // the same location. We will measure its area, and then continue till the next point at + // the same location. The smallest absolute value of area is the one we will use. + point_ptr point_a = point_1; + point_ptr min_a = nullptr; + point_ptr min_b = nullptr; + point_ptr pt = point_1->next; + double area = 0.0; + double min_area = std::numeric_limits::max(); + while (pt != point_1) { + area += static_cast(pt->prev->x + pt->x) * static_cast(pt->prev->y - pt->y); + if (*pt == *point_1) { + if (std::fabs(area) < min_area) { + min_area = std::fabs(area); + min_a = point_a; + min_b = pt; + } + point_a = pt; + area = 0.0; + } + pt = pt->next; + } + if (point_a == point_1) { + // LCOV_EXCL_START + throw std::runtime_error("No other point was between point_1 on the path"); + // LCOV_EXCL_END + } + area += static_cast(pt->prev->x + pt->x) * static_cast(pt->prev->y - pt->y); + if (std::fabs(area) < min_area) { + min_area = std::fabs(area); + min_a = point_a; + min_b = pt; + } + handle_self_intersections(min_a, min_b, dupe_ring, rings); + return true; + } + } + } + handle_self_intersections(point_1, point_2, dupe_ring, rings); + return true; +} + +template +void process_repeated_points(std::size_t first_index, + std::size_t last_index, + std::unordered_multimap, point_ptr_pair>& dupe_ring, + ring_manager& rings) { + for (std::size_t j = first_index; j <= last_index; ++j) { + while (process_repeated_point_set(first_index, last_index, j, dupe_ring, rings)) + ; + } + +#ifdef DEBUG + // LCOV_EXCL_START + for (std::size_t j = first_index; j <= last_index; ++j) { + point_ptr op_j = rings.all_points[j]; + if (!op_j->ring) { + continue; + } + double ring_area = area(op_j->ring); + bool ring_is_positive = ring_area > 0.0; + bool ring_is_zero = value_is_zero(ring_area); + if (!ring_is_zero) { + if (op_j->ring->parent) { + double parent_area = area(op_j->ring->parent); + bool parent_is_positive = parent_area > 0.0; + bool parent_is_zero = value_is_zero(parent_area); + if (!parent_is_zero && ring_is_positive == parent_is_positive) { + throw std::runtime_error( + "Created a ring with a parent having the same orientation (sign of area)"); + } + } + for (auto& c : op_j->ring->children) { + double c_area = area(c); + bool c_is_positive = c_area > 0.0; + bool c_is_zero = value_is_zero(c_area); + if (!c_is_zero && ring_is_positive == c_is_positive) { + throw std::runtime_error( + "Created a ring with a child having the same orientation (sign of area)"); + } + } + } + } + // LCOV_EXCL_STOP +#endif +} + +template +bool process_chains(std::size_t first_index, + std::size_t last_index, + std::unordered_multimap, point_ptr_pair>& dupe_ring, + ring_manager& rings, + mapbox::geometry::point& rewind_point) { + bool rewind = false; + for (std::size_t j = first_index; j <= last_index; ++j) { + point_ptr op_j = rings.all_points[j]; + if (!op_j->ring) { + continue; + } + for (std::size_t k = j + 1; k <= last_index; ++k) { + point_ptr op_k = rings.all_points[k]; + if (!op_k->ring || !op_j->ring) { + continue; + } + if (fix_intersects(dupe_ring, op_j, op_k, rings, rewind_point)) { + rewind = true; + } + } + } + return rewind; +} + +template +bool process_collinear_edges(std::size_t first_index, + std::size_t last_index, + std::unordered_multimap, point_ptr_pair>& dupe_ring, + ring_manager& rings, + mapbox::geometry::point& rewind_point) { + bool rewind = false; + for (std::size_t j = first_index; j <= last_index; ++j) { + point_ptr op_j = rings.all_points[j]; + if (!op_j->ring) { + continue; + } + for (std::size_t k = j + 1; k <= last_index; ++k) { + point_ptr op_k = rings.all_points[k]; + if (!op_k->ring || !op_j->ring) { + continue; + } + if (handle_collinear_edges(op_j, op_k, dupe_ring, rings, rewind_point)) { + rewind = true; + } + } + } + return rewind; +} + +template +bool index_is_after_point(std::size_t const& i, + mapbox::geometry::point const& pt, + ring_manager const& rings) { + if (i == 0) { + return false; + } + + if (rings.all_points[i]->y < pt.y) { + return true; + } else if (rings.all_points[i]->y > pt.y) { + return false; + } else if (rings.all_points[i]->x >= pt.x) { + return true; + } else { + return false; + } +} + +template +void rewind_to_point(std::size_t& i, + mapbox::geometry::point const& pt, + ring_manager const& rings) { + if (i >= rings.rings.size()) { + i = rings.rings.size() - 1; + } + while (index_is_after_point(i, pt, rings)) { + --i; + } +} + +template +void remove_spikes_in_polygons(ring_ptr r, ring_manager& rings) { + + point_ptr first_point = r->points; + remove_spikes(first_point); + if (!first_point) { + r->points = nullptr; + r->area = std::numeric_limits::quiet_NaN(); + remove_ring(r, rings); + return; + } + point_ptr p = first_point->next; + while (p != first_point) { + remove_spikes(p); + if (!p) { + r->points = nullptr; + r->area = std::numeric_limits::quiet_NaN(); + remove_ring(r, rings); + return; + } + if (p->ring && !first_point->ring) { + first_point = p; + } + p = p->next; + } +} + +template +void fixup_out_polyline(ring& ring, ring_manager& rings) { + point_ptr pp = ring.points; + point_ptr lastPP = pp->prev; + while (pp != lastPP) { + pp = pp->next; + if (*pp == *pp->prev) { + if (pp == lastPP) + lastPP = pp->prev; + point_ptr tmpPP = pp->prev; + tmpPP->next = pp->next; + pp->next->prev = tmpPP; + // delete pp; + pp->next = pp; + pp->prev = pp; + pp->ring = nullptr; + pp = tmpPP; + } + } + + if (pp == pp->prev) { + remove_ring(&ring, rings); + dispose_out_points(pp); + ring.points = nullptr; + return; + } +} + + +template +void fixup_out_polygon(ring& ring, ring_manager& rings) { + // FixupOutPolygon() - removes duplicate points and simplifies consecutive + // parallel edges by removing the middle vertex. + point_ptr lastOK = nullptr; + ring.bottom_point = nullptr; + point_ptr pp = ring.points; + + for (;;) { + if (pp->prev == pp || pp->prev == pp->next) { + // We now need to make sure any children rings to this are promoted and their hole + // status is changed + // promote_children_of_removed_ring(&ring, rings); + remove_ring(&ring, rings); + dispose_out_points(pp); + ring.points = nullptr; + return; + } + + // test for duplicate points and collinear edges ... + if ((*pp == *pp->next) || (*pp == *pp->prev) || + (slopes_equal(*pp->prev, *pp, *pp->next))) { + lastOK = nullptr; + point_ptr tmp = pp; + pp->prev->next = pp->next; + pp->next->prev = pp->prev; + pp = pp->prev; + tmp->ring = nullptr; + tmp->next = tmp; + tmp->prev = tmp; + } else if (pp == lastOK) { + break; + } else { + if (!lastOK) { + lastOK = pp; + } + pp = pp->next; + } + } + ring.points = pp; +} + + +template +void do_simple_polygons(ring_manager& rings) { + + // fix orientations ... + for (auto& r : rings.rings) { + if (!r.points || r.is_open) { + continue; + } + std::size_t s = 0; + if (ring_is_hole(&r) == (area_from_point(r.points, s) > 0)) { + reverse_ring(r.points); + } + remove_spikes_in_polygons(&r, rings); + r.area = std::numeric_limits::quiet_NaN(); + } + + std::stable_sort(rings.all_points.begin(), rings.all_points.end(), point_ptr_cmp()); + std::unordered_multimap, point_ptr_pair> dupe_ring; + dupe_ring.reserve(rings.rings.size()); + + // Find sets of repeated points and process them + std::size_t count = 0; + for (std::size_t i = 1; i < rings.all_points.size(); ++i) { + if (*rings.all_points[i] == *rings.all_points[i - 1]) { + ++count; + if (i < (rings.all_points.size() - 1)) { + continue; + } else { + ++i; + } + } + + if (count == 0) { + continue; + } + std::size_t first_index = i - count - 1; + std::size_t last_index = i - 1; + auto sort_begin = rings.all_points.begin(); + std::advance(sort_begin, first_index); + auto sort_end = rings.all_points.begin(); + std::advance(sort_end, i); + std::stable_sort(sort_begin, sort_end, point_ptr_depth_cmp()); + process_repeated_points(first_index, last_index, dupe_ring, rings); + mapbox::geometry::point rewind_point(std::numeric_limits::min(), + std::numeric_limits::min()); + bool do_rewind = false; + if (process_chains(first_index, last_index, dupe_ring, rings, rewind_point)) { + do_rewind = true; + } + if (process_collinear_edges(first_index, last_index, dupe_ring, rings, rewind_point)) { + do_rewind = true; + } + if (do_rewind) { + rewind_to_point(i, rewind_point, rings); + } + count = 0; + } + +#if DEBUG + // LCOV_EXCL_START + for (auto& r : rings.rings) { + if (!r.points || r.is_open) { + continue; + } + double stored_area = area(&r); + std::size_t s = 0; + double calculated_area = area_from_point(r.points, s); + if (!values_near_equal(stored_area, calculated_area)) { + throw std::runtime_error("Difference in stored area vs calculated area!"); + } + } + // LCOV_EXCL_END +#endif + + for (auto& r : rings.rings) { + if (!r.points || r.is_open) { + continue; + } + fixup_out_polygon(r, rings); + if (ring_is_hole(&r) == (area(&r) > 0.0)) { + reverse_ring(r.points); + r.area = std::numeric_limits::quiet_NaN(); + } + } +} +} +} +} diff --git a/mapbox/geometry/wagyu/util.hpp b/mapbox/geometry/wagyu/util.hpp new file mode 100644 index 0000000..45938f6 --- /dev/null +++ b/mapbox/geometry/wagyu/util.hpp @@ -0,0 +1,83 @@ +#pragma once + +#include + +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +double area(mapbox::geometry::linear_ring const& poly) { + std::size_t size = poly.size(); + if (size < 3) { + return 0.0; + } + + double a = 0.0; + auto itr = poly.begin(); + auto itr_prev = poly.end(); + --itr_prev; + a += static_cast(itr_prev->x + itr->x) * static_cast(itr_prev->y - itr->y); + ++itr; + itr_prev = poly.begin(); + for (; itr != poly.end(); ++itr, ++itr_prev) { + a += static_cast(itr_prev->x + itr->x) * static_cast(itr_prev->y - itr->y); + } + return -a * 0.5; +} + +inline bool value_is_zero(double val) { + return std::fabs(val) < std::numeric_limits::epsilon(); +} + +inline bool values_are_equal(double x, double y) { + return value_is_zero(x - y); +} + +inline bool values_near_equal(double x, double y) { + return std::fabs(x - y) < (5.0 * std::numeric_limits::epsilon()); +} + +inline bool greater_than_or_equal(double x, double y) { + return x > y || values_are_equal(x, y); +} + +inline bool less_than_or_equal(double x, double y) { + return x < y || values_are_equal(x, y); +} + +template +bool slopes_equal(mapbox::geometry::point const& pt1, + mapbox::geometry::point const& pt2, + mapbox::geometry::point const& pt3) { + return (pt1.y - pt2.y) * (pt2.x - pt3.x) == (pt1.x - pt2.x) * (pt2.y - pt3.y); +} + +template +bool slopes_equal(mapbox::geometry::wagyu::point const& pt1, + mapbox::geometry::wagyu::point const& pt2, + mapbox::geometry::point const& pt3) { + return (pt1.y - pt2.y) * (pt2.x - pt3.x) == (pt1.x - pt2.x) * (pt2.y - pt3.y); +} + +template +bool slopes_equal(mapbox::geometry::wagyu::point const& pt1, + mapbox::geometry::wagyu::point const& pt2, + mapbox::geometry::wagyu::point const& pt3) { + return (pt1.y - pt2.y) * (pt2.x - pt3.x) == (pt1.x - pt2.x) * (pt2.y - pt3.y); +} + +template +bool slopes_equal(mapbox::geometry::point const& pt1, + mapbox::geometry::point const& pt2, + mapbox::geometry::point const& pt3, + mapbox::geometry::point const& pt4) { + return (pt1.y - pt2.y) * (pt3.x - pt4.x) == (pt1.x - pt2.x) * (pt3.y - pt4.y); +} +} +} +} diff --git a/mapbox/geometry/wagyu/vatti.hpp b/mapbox/geometry/wagyu/vatti.hpp new file mode 100644 index 0000000..5e32fa0 --- /dev/null +++ b/mapbox/geometry/wagyu/vatti.hpp @@ -0,0 +1,77 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +bool execute_vatti(local_minimum_list& minima_list, + ring_manager& rings, + clip_type cliptype, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + if (minima_list.empty()) { + return false; + } + + active_bound_list active_bounds; + scanbeam_list scanbeam; + T scanline_y = std::numeric_limits::max(); + + local_minimum_ptr_list minima_sorted; + minima_sorted.reserve(minima_list.size()); + for (auto& lm : minima_list) { + minima_sorted.push_back(&lm); + } + std::stable_sort(minima_sorted.begin(), minima_sorted.end(), local_minimum_sorter()); + local_minimum_ptr_list_itr current_lm = minima_sorted.begin(); + // std::clog << output_all_edges(minima_sorted) << std::endl; + + setup_scanbeam(minima_list, scanbeam); + rings.current_hp_itr = rings.hot_pixels.begin(); + + while (pop_from_scanbeam(scanline_y, scanbeam) || current_lm != minima_sorted.end()) { + + process_intersections(scanline_y, active_bounds, cliptype, + subject_fill_type, clip_fill_type, rings); + + update_current_hp_itr(scanline_y, rings); + + // First we process bounds that has already been added to the active bound list -- + // if the active bound list is empty local minima that are at this scanline_y and + // have a horizontal edge at the local minima will be processed + process_edges_at_top_of_scanbeam(scanline_y, active_bounds, scanbeam, minima_sorted, + current_lm, rings, cliptype, subject_fill_type, + clip_fill_type); + + // Next we will add local minima bounds to the active bounds list that are on the local + // minima queue at + // this current scanline_y + insert_local_minima_into_ABL(scanline_y, minima_sorted, current_lm, active_bounds, + rings, scanbeam, cliptype, subject_fill_type, + clip_fill_type); + + } + // std::clog << rings.rings << std::endl; + // std::clog << output_as_polygon(rings.all_rings[0]); + return true; +} +} +} +} diff --git a/mapbox/geometry/wagyu/wagyu.hpp b/mapbox/geometry/wagyu/wagyu.hpp new file mode 100644 index 0000000..730a565 --- /dev/null +++ b/mapbox/geometry/wagyu/wagyu.hpp @@ -0,0 +1,141 @@ +#pragma once + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace mapbox { +namespace geometry { +namespace wagyu { + +template +class wagyu { +private: + using value_type = T; + + local_minimum_list minima_list; + bool has_open_paths; + bool reverse_output; + + wagyu( wagyu const& ) = delete; + wagyu& operator=(wagyu const& ) = delete; + +public: + wagyu() : minima_list(), has_open_paths(false), reverse_output(false) { + } + + ~wagyu() { + clear(); + } + + bool add_line(mapbox::geometry::line_string const& pg) { + bool success = add_line_string(pg, minima_list); + if (success) { + has_open_paths = true; + } + return success; + } + + bool add_ring(mapbox::geometry::linear_ring const& pg, + polygon_type p_type = polygon_type_subject) { + return add_linear_ring(pg, minima_list, p_type); + } + + bool add_polygon(mapbox::geometry::polygon const& ppg, + polygon_type p_type = polygon_type_subject) { + bool result = false; + for (auto const& r : ppg) { + if (add_ring(r, p_type)) { + result = true; + } + } + return result; + } + + void reverse_rings(bool value) { + reverse_output = value; + } + + void clear() { + minima_list.clear(); + has_open_paths = false; + } + + mapbox::geometry::box get_bounds() { + mapbox::geometry::point min = { 0, 0 }; + mapbox::geometry::point max = { 0, 0 }; + if (minima_list.empty()) { + return mapbox::geometry::box(min, max); + } + bool first_set = false; + for (auto const& lm : minima_list) { + if (!lm.left_bound.edges.empty()) { + if (!first_set) { + min = lm.left_bound.edges.front().top; + max = lm.left_bound.edges.back().bot; + first_set = true; + } else { + min.y = std::min(min.y, lm.left_bound.edges.front().top.y); + max.y = std::max(max.y, lm.left_bound.edges.back().bot.y); + max.x = std::max(max.x, lm.left_bound.edges.back().top.x); + min.x = std::min(min.x, lm.left_bound.edges.back().top.x); + } + for (auto const& e : lm.left_bound.edges) { + max.x = std::max(max.x, e.bot.x); + min.x = std::min(min.x, e.bot.x); + } + } + if (!lm.right_bound.edges.empty()) { + if (!first_set) { + min = lm.right_bound.edges.front().top; + max = lm.right_bound.edges.back().bot; + first_set = true; + } else { + min.y = std::min(min.y, lm.right_bound.edges.front().top.y); + max.y = std::max(max.y, lm.right_bound.edges.back().bot.y); + max.x = std::max(max.x, lm.right_bound.edges.back().top.x); + min.x = std::min(min.x, lm.right_bound.edges.back().top.x); + } + for (auto const& e : lm.right_bound.edges) { + max.x = std::max(max.x, e.bot.x); + min.x = std::min(min.x, e.bot.x); + } + } + } + return mapbox::geometry::box(min, max); + } + + bool execute(clip_type cliptype, + mapbox::geometry::multi_polygon& solution, + fill_type subject_fill_type, + fill_type clip_fill_type) { + + ring_manager rings; + + build_hot_pixels(minima_list, rings); + + if (!execute_vatti(minima_list, rings, cliptype, subject_fill_type, clip_fill_type)) { + return false; + } + + do_simple_polygons(rings); + + build_result(solution, rings, reverse_output); + + return true; + } +}; +} +} +}