Copy of wagyu-0.2.0 (b3fa127)

mapbox/geometry/wagyu/active_bound_list.hpp

@@ -0,0 +1,526 @@
+#pragma once
+
+#ifdef DEBUG
+#include <iostream>
+#endif
+
+#include <mapbox/geometry/wagyu/bound.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/exceptions.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+#include <mapbox/geometry/wagyu/local_minimum_util.hpp>
+#include <mapbox/geometry/wagyu/ring.hpp>
+#include <mapbox/geometry/wagyu/scanbeam.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+using active_bound_list = std::list<bound_ptr<T>>;
+
+template <typename T>
+using active_bound_list_itr = typename active_bound_list<T>::iterator;
+
+template <typename T>
+using active_bound_list_rev_itr = typename active_bound_list<T>::reverse_iterator;
+
+#ifdef DEBUG
+
+template <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const active_bound_list<T>& bnds) {
+    std::size_t c = 0;
+    for (auto const& bnd : bnds) {
+        out << "Index: " << c++ << std::endl;
+        out << *bnd;
+    }
+    return out;
+}
+
+template <typename T>
+std::string output_edges(active_bound_list<T> 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 <typename T>
+bool is_even_odd_fill_type(bound<T> 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 <typename T>
+bool is_even_odd_alt_fill_type(bound<T> 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 <typename T>
+inline bool bound2_inserts_before_bound1(bound<T> const& bound1, bound<T> 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 <typename T>
+active_bound_list_itr<T> insert_bound_into_ABL(bound<T>& bnd, active_bound_list<T>& 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 <typename T>
+active_bound_list_itr<T> insert_bound_into_ABL(bound<T>& bnd,
+                                               active_bound_list_itr<T> itr,
+                                               active_bound_list<T>& active_bounds) {
+    while (itr != active_bounds.end() && !bound2_inserts_before_bound1(*(*itr), bnd)) {
+        ++itr;
+    }
+    return active_bounds.insert(itr, &bnd);
+}
+
+template <typename T>
+inline bool is_maxima(bound<T>& bnd, T y) {
+    return bnd.next_edge == bnd.edges.end() &&
+           bnd.current_edge->top.y == y;
+}
+
+template <typename T>
+inline bool is_maxima(active_bound_list_itr<T>& bnd, T y) {
+    return is_maxima(*(*bnd), y);
+}
+
+template <typename T>
+inline bool is_intermediate(bound<T>& bnd, T y) {
+    return bnd.next_edge != bnd.edges.end() &&
+           bnd.current_edge->top.y == y;
+}
+
+template <typename T>
+inline bool is_intermediate(active_bound_list_itr<T>& bnd, T y) {
+    return is_intermediate(*(*bnd), y);
+}
+
+template <typename T>
+inline bool current_edge_is_horizontal(active_bound_list_itr<T>& bnd) {
+    return is_horizontal(*((*bnd)->current_edge));
+}
+
+template <typename T>
+inline bool next_edge_is_horizontal(active_bound_list_itr<T>& bnd) {
+    return is_horizontal(*((*bnd)->next_edge));
+}
+
+template <typename T>
+inline void swap_positions_in_ABL(active_bound_list_itr<T>& bnd1,
+                                  active_bound_list_itr<T>& bnd2,
+                                  active_bound_list<T>& active_bounds) {
+    if (std::next(bnd2) == bnd1) {
+        active_bounds.splice(bnd2, active_bounds, bnd1);
+    } else { 
+        active_bounds.splice(bnd1, active_bounds, bnd2);
+    } 
+}
+
+template <typename T>
+void next_edge_in_bound(active_bound_list_itr<T>& bnd, scanbeam_list<T>& scanbeam) {
+    ++((*bnd)->current_edge);
+    if ((*bnd)->current_edge != (*bnd)->edges.end()) {
+        ++((*bnd)->next_edge);
+        (*bnd)->current_x = static_cast<double>((*bnd)->current_edge->bot.x);
+        if (!current_edge_is_horizontal<T>(bnd)) {
+            scanbeam.push((*bnd)->current_edge->top.y);
+        }
+    }
+}
+
+template <typename T>
+active_bound_list_itr<T> get_maxima_pair(active_bound_list_itr<T> bnd,
+                                         active_bound_list<T>& 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 <typename T>
+void set_winding_count(active_bound_list_itr<T>& bnd_itr,
+                       active_bound_list<T>& active_bounds,
+                       clip_type cliptype,
+                       fill_type subject_fill_type,
+                       fill_type clip_fill_type) {
+
+    auto rev_bnd_itr = active_bound_list_rev_itr<T>(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 <typename T>
+bool is_contributing(bound<T> 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 <typename T>
+void insert_lm_only_one_bound(bound<T>& bnd,
+                              active_bound_list<T>& active_bounds,
+                              ring_manager<T>& rings,
+                              scanbeam_list<T>& 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<T>(abl_itr)) {
+        scanbeam.push((*abl_itr)->current_edge->top.y);
+    }
+}
+
+template <typename T>
+void insert_lm_left_and_right_bound(bound<T>& left_bound,
+                                    bound<T>& right_bound,
+                                    active_bound_list<T>& active_bounds,
+                                    ring_manager<T>& rings,
+                                    scanbeam_list<T>& 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<T>(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 <typename T>
+void insert_local_minima_into_ABL(T const bot_y,
+                                  local_minimum_ptr_list<T> const& minima_sorted,
+                                  local_minimum_ptr_list_itr<T>& current_lm,
+                                  active_bound_list<T>& active_bounds,
+                                  ring_manager<T>& rings,
+                                  scanbeam_list<T>& 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<T>(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 <typename T>
+void insert_horizontal_local_minima_into_ABL(T const top_y,
+                                             local_minimum_ptr_list<T> const& minima_sorted,
+                                             local_minimum_ptr_list_itr<T>& current_lm,
+                                             active_bound_list<T>& active_bounds,
+                                             ring_manager<T>& rings,
+                                             scanbeam_list<T>& 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<T>(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;
+    }
+}
+}
+}
+}

mapbox/geometry/wagyu/bound.hpp

@@ -0,0 +1,95 @@
+#pragma once
+
+#include <list>
+
+#include <mapbox/geometry/point.hpp>
+
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/ring.hpp>
+
+#ifdef DEBUG
+#include <iostream>
+#endif
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+struct bound {
+
+    edge_list<T> edges;
+    edge_list_itr<T> current_edge;
+    edge_list_itr<T> next_edge;
+    mapbox::geometry::point<T> last_point;
+    ring_ptr<T> ring;
+    bound_ptr<T> 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<T> && 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const bound<T>& 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<int>(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
+}
+}
+}

mapbox/geometry/wagyu/build_edges.hpp

@@ -0,0 +1,216 @@
+#pragma once
+
+#include <mapbox/geometry/line_string.hpp>
+#include <mapbox/geometry/point.hpp>
+#include <mapbox/geometry/polygon.hpp>
+
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/exceptions.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+bool build_edge_list(mapbox::geometry::line_string<T> const& path_geometry,
+                     edge_list<T>& 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 <typename T>
+bool point_2_is_between_point_1_and_point_3(mapbox::geometry::point<T> const& pt1, 
+                                            mapbox::geometry::point<T> const& pt2, 
+                                            mapbox::geometry::point<T> 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 <typename T>
+bool build_edge_list(mapbox::geometry::linear_ring<T> const& path_geometry, edge_list<T>& 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<value_type> pt1 = *itr_rev;
+    mapbox::geometry::point<value_type> 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<value_type> pt3 = *itr;
+    auto itr_last = itr_rev.base();
+    mapbox::geometry::point<value_type> front_pt;
+    mapbox::geometry::point<value_type> 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;
+}
+}
+}
+}

mapbox/geometry/wagyu/build_local_minima_list.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+#include <mapbox/geometry/wagyu/build_edges.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+#include <mapbox/geometry/wagyu/local_minimum_util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+bool add_line_string(mapbox::geometry::line_string<T> const& path_geometry,
+                     local_minimum_list<T>& minima_list) {
+    bool is_flat = true;
+    edge_list<T> 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 <typename T>
+bool add_linear_ring(mapbox::geometry::linear_ring<T> const& path_geometry,
+                     local_minimum_list<T>& minima_list,
+                     polygon_type p_type) {
+    edge_list<T> 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;
+}
+
+}
+}
+}

mapbox/geometry/wagyu/build_result.hpp

@@ -0,0 +1,68 @@
+#pragma once
+
+#include <mapbox/geometry/wagyu/ring.hpp>
+#include <mapbox/geometry/wagyu/ring_util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+void push_ring_to_polygon(mapbox::geometry::polygon<T>& poly, ring_ptr<T>& r, bool reverse_output) {
+    mapbox::geometry::linear_ring<T> 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 <typename T>
+void build_result_polygons(std::vector<mapbox::geometry::polygon<T>>& solution,
+                           ring_list<T>& 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 <typename T>
+void build_result(std::vector<mapbox::geometry::polygon<T>>& solution,
+                  ring_manager<T>& rings,
+                  bool reverse_output) {
+    build_result_polygons(solution, rings.children, reverse_output);
+}
+}
+}
+}

mapbox/geometry/wagyu/config.hpp

@@ -0,0 +1,52 @@
+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <list>
+
+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 <typename T>
+using maxima_list = std::list<T>;
+}
+}
+}

mapbox/geometry/wagyu/edge.hpp

@@ -0,0 +1,122 @@
+#pragma once
+
+#include <cmath>
+#include <limits>
+#include <list>
+
+#include <mapbox/geometry/point.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+#ifdef DEBUG
+#include <iostream>
+#endif
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+struct bound;
+
+template <typename T>
+using bound_ptr = bound<T>*;
+
+template <typename T>
+struct edge {
+    mapbox::geometry::point<T> bot;
+    mapbox::geometry::point<T> top;
+    double dx;
+
+    edge(edge<T> && e) noexcept
+        : bot(std::move(e.bot)),
+          top(std::move(e.top)),
+          dx(std::move(e.dx)) {}
+
+    edge& operator=(edge<T> && e) noexcept {
+        bot = std::move(e.bot);
+        top = std::move(e.top);
+        dx = std::move(e.dx);
+        return *this;
+    }
+
+    edge(mapbox::geometry::point<T> const& current, mapbox::geometry::point<T> 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<double>(top.y - bot.y);
+        if (value_is_zero(dy)) {
+            dx = std::numeric_limits<double>::infinity();
+        } else {
+            dx = static_cast<double>(top.x - bot.x) / dy;
+        }
+    }
+
+};
+
+template <typename T>
+using edge_ptr = edge<T>*;
+
+template <typename T>
+using edge_list = std::vector<edge<T>>;
+
+template <typename T>
+using edge_list_itr = typename edge_list<T>::iterator;
+
+template <typename T>
+bool slopes_equal(edge<T> const& e1, edge<T> 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 <typename T>
+inline bool is_horizontal(edge<T> const& e) {
+    return std::isinf(e.dx);
+}
+
+template <typename T>
+inline double get_current_x(edge<T> const& edge, const T current_y) {
+    if (current_y == edge.top.y) {
+        return static_cast<double>(edge.top.x);
+    } else {
+        return static_cast<double>(edge.bot.x) +
+               edge.dx * static_cast<double>(current_y - edge.bot.y);
+    }
+}
+
+#ifdef DEBUG
+
+template <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const edge<T>& 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     edge_list<T> 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
+}
+}
+}

mapbox/geometry/wagyu/exceptions.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <stdexcept>
+
+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();
+    }
+};
+}
+}
+}

mapbox/geometry/wagyu/intersect.hpp

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <set>
+
+#include <mapbox/geometry/point.hpp>
+
+#include <mapbox/geometry/wagyu/active_bound_list.hpp>
+
+#ifdef DEBUG
+#include <iostream>
+#endif
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+struct intersect_node {
+
+    active_bound_list_itr<T> bound1;
+    active_bound_list_itr<T> bound2;
+    mapbox::geometry::point<double> pt;
+
+    intersect_node(intersect_node<T> && n)
+        : bound1(std::move(n.bound1)),
+          bound2(std::move(n.bound2)),
+          pt(std::move(n.pt)) { }
+
+    intersect_node& operator=(intersect_node<T> && n) {
+        bound1 = std::move(n.bound1);
+        bound2 = std::move(n.bound2);
+        pt = std::move(n.pt);
+        return *this;
+    }
+
+    intersect_node(active_bound_list_itr<T> const& bound1_,
+                   active_bound_list_itr<T> const& bound2_,
+                   mapbox::geometry::point<double> const& pt_)
+        : bound1(bound1_), bound2(bound2_), pt(pt_) {
+    }
+};
+
+template <typename T>
+using intersect_list = std::vector<intersect_node<T>>;
+
+#ifdef DEBUG
+
+template <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const intersect_node<T>& 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const intersect_list<T>& ints) {
+    std::size_t c = 0;
+    for (auto const& i : ints) {
+        out << "Intersection: " << c++ << std::endl;
+        out << i;
+    }
+    return out;
+}
+
+#endif
+}
+}
+}

mapbox/geometry/wagyu/intersect_util.hpp

@@ -0,0 +1,387 @@
+#pragma once
+
+#include <mapbox/geometry/wagyu/active_bound_list.hpp>
+#include <mapbox/geometry/wagyu/bound.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/intersect.hpp>
+#include <mapbox/geometry/wagyu/ring_util.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+struct intersect_list_sorter {
+    inline bool operator()(intersect_node<T> const& node1, intersect_node<T> 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 <typename T>
+inline mapbox::geometry::point<T> round_point(mapbox::geometry::point<double> const& pt) {
+    return mapbox::geometry::point<T>(round_towards_max<T>(pt.x), round_towards_max<T>(pt.y));
+}
+
+template <typename T>
+inline void swap_rings(bound<T>& b1, bound<T>& b2) {
+    ring_ptr<T> ring = b1.ring;
+    b1.ring = b2.ring;
+    b2.ring = ring;
+}
+
+template <typename T>
+inline void swap_sides(bound<T>& b1, bound<T>& b2) {
+    edge_side side = b1.side;
+    b1.side = b2.side;
+    b2.side = side;
+}
+
+template <typename T1, typename T2>
+bool get_edge_intersection(edge<T1> const& e1,
+                           edge<T1> const& e2,
+                           mapbox::geometry::point<T2>& pt) {
+    T2 p0_x = static_cast<T2>(e1.bot.x);
+    T2 p0_y = static_cast<T2>(e1.bot.y);
+    T2 p1_x = static_cast<T2>(e1.top.x);
+    T2 p1_y = static_cast<T2>(e1.top.y);
+    T2 p2_x = static_cast<T2>(e2.bot.x);
+    T2 p2_y = static_cast<T2>(e2.bot.y);
+    T2 p3_x = static_cast<T2>(e2.top.x);
+    T2 p3_y = static_cast<T2>(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 <typename T>
+void build_intersect_list(active_bound_list<T>& active_bounds,
+                          intersect_list<T>& 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<double> pt;
+                if (!get_edge_intersection<T, double>(*((*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 <typename T>
+void intersect_bounds(active_bound_list_itr<T>& b1,
+                      active_bound_list_itr<T>& b2,
+                      mapbox::geometry::point<T> const& pt,
+                      clip_type cliptype,
+                      fill_type subject_fill_type,
+                      fill_type clip_fill_type,
+                      ring_manager<T>& rings,
+                      active_bound_list<T>& 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 <typename T>
+bool bounds_adjacent(intersect_node<T> const& inode) {
+    return (std::next(inode.bound1) == inode.bound2) || (std::next(inode.bound2) == inode.bound1);
+}
+
+template <typename T>
+void process_intersect_list(intersect_list<T>& intersects,
+                            clip_type cliptype,
+                            fill_type subject_fill_type,
+                            fill_type clip_fill_type,
+                            ring_manager<T>& rings,
+                            active_bound_list<T>& 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<T> pt = round_point<T>(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 <typename T>
+void update_current_x(active_bound_list<T>& 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 <typename T>
+void process_intersections(T top_y,
+                           active_bound_list<T>& active_bounds,
+                           clip_type cliptype,
+                           fill_type subject_fill_type,
+                           fill_type clip_fill_type,
+                           ring_manager<T>& rings) {
+    if (active_bounds.empty()) {
+        return;
+    }
+    update_current_x(active_bounds, top_y);
+    intersect_list<T> intersects;
+    build_intersect_list(active_bounds, intersects);
+
+    if (intersects.empty()) {
+        return;
+    }
+
+    // Restore order of active bounds list
+    active_bounds.sort([] (bound_ptr<T> const& b1, bound_ptr<T> const& b2) {
+        return b1->pos < b2->pos;
+    });
+
+    // Sort the intersection list
+    std::stable_sort(intersects.begin(), intersects.end(), intersect_list_sorter<T>());
+
+    process_intersect_list(intersects, cliptype, subject_fill_type, clip_fill_type, rings,
+                           active_bounds);
+}
+}
+}
+}

mapbox/geometry/wagyu/local_minimum.hpp

@@ -0,0 +1,116 @@
+#pragma once
+
+#ifdef DEBUG
+#include <iostream>
+#endif
+#include <queue>
+
+#include <mapbox/geometry/wagyu/bound.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+struct local_minimum {
+    bound<T> left_bound;
+    bound<T> right_bound;
+    T y;
+    bool minimum_has_horizontal;
+
+    local_minimum(bound<T>&& left_bound_, bound<T>&& 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 <typename T>
+using local_minimum_list = std::deque<local_minimum<T>>;
+
+template <typename T>
+using local_minimum_itr = typename local_minimum_list<T>::iterator;
+
+template <typename T>
+using local_minimum_ptr = local_minimum<T>*;
+
+template <typename T>
+using local_minimum_ptr_list = std::vector<local_minimum_ptr<T>>;
+
+template <typename T>
+using local_minimum_ptr_list_itr = typename local_minimum_ptr_list<T>::iterator;
+
+template <typename T>
+struct local_minimum_sorter {
+    inline bool operator()(local_minimum_ptr<T> const& locMin1,
+                           local_minimum_ptr<T> 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const local_minimum<T>& 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const local_minimum_ptr_list<T>& lms) {
+    for (auto const& lm : lms) {
+        out << *lm;
+    }
+    return out;
+}
+
+template <typename T>
+std::string output_all_edges(local_minimum_ptr_list<T> 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
+}
+}
+}

mapbox/geometry/wagyu/local_minimum_util.hpp

@@ -0,0 +1,419 @@
+#pragma once
+
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+
+#ifdef DEBUG
+#include <mapbox/geometry/wagyu/exceptions.hpp>
+#endif
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+inline void reverse_horizontal(edge<T>& 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 <typename T>
+void start_list_on_local_maximum(edge_list<T>& 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 <typename T>
+bound<T> create_bound_towards_minimum(edge_list<T>& edges) {
+    if (edges.size() == 1) {
+        if (is_horizontal(edges.front())) {
+            reverse_horizontal(edges.front());
+        }
+        bound<T> 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<T> 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 <typename T>
+bound<T> create_bound_towards_maximum(edge_list<T>& edges) {
+    if (edges.size() == 1) {
+        bound<T> 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<T> 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 <typename T>
+void fix_horizontals(bound<T>& 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 <typename T>
+void move_horizontals_on_left_to_right(bound<T>& left_bound, bound<T>& 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 <typename T>
+void add_line_to_local_minima_list(edge_list<T>& edges, local_minimum_list<T>& 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<T> 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<T> 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<T> 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 <typename T>
+void add_ring_to_local_minima_list(edge_list<T>& edges,
+                                   local_minimum_list<T>& 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<T> first_minimum = nullptr;
+    bound_ptr<T> 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 <typename T>
+void initialize_lm(local_minimum_ptr_list_itr<T>& 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<double>((*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<double>((*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;
+    }
+}
+}
+}
+}

mapbox/geometry/wagyu/point.hpp

@@ -0,0 +1,110 @@
+#pragma once
+
+#include <mapbox/geometry/point.hpp>
+
+#ifdef DEBUG
+#include <iostream>
+#endif
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+struct point;
+
+template <typename T>
+using point_ptr = point<T>*;
+
+template <typename T>
+using const_point_ptr = point<T>* const;
+
+template <typename T>
+struct ring;
+
+template <typename T>
+using ring_ptr = ring<T>*;
+
+template <typename T>
+using const_ring_ptr = ring<T>* const;
+
+template <typename T>
+struct point {
+    using coordinate_type = T;
+    ring_ptr<T> ring;
+    T x;
+    T y;
+    point_ptr<T> next;
+    point_ptr<T> prev;
+
+    point(point<T> && 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<T> ring_, mapbox::geometry::point<T> const& pt)
+        : ring(ring_), x(pt.x), y(pt.y), next(this), prev(this) {
+    }
+
+    point(ring_ptr<T> ring_, mapbox::geometry::point<T> const& pt, point_ptr<T> 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 <typename T>
+bool operator==(point<T> const& lhs, point<T> const& rhs) {
+    return lhs.x == rhs.x && lhs.y == rhs.y;
+}
+
+template <typename T>
+bool operator==(mapbox::geometry::point<T> const& lhs, point<T> const& rhs) {
+    return lhs.x == rhs.x && lhs.y == rhs.y;
+}
+
+template <typename T>
+bool operator==(point<T> const& lhs, mapbox::geometry::point<T> const& rhs) {
+    return lhs.x == rhs.x && lhs.y == rhs.y;
+}
+
+template <typename T>
+bool operator!=(point<T> const& lhs, point<T> const& rhs) {
+    return lhs.x != rhs.x || lhs.y != rhs.y;
+}
+
+template <typename T>
+bool operator!=(mapbox::geometry::point<T> const& lhs, point<T> const& rhs) {
+    return lhs.x != rhs.x || lhs.y != rhs.y;
+}
+
+template <typename T>
+bool operator!=(point<T> const& lhs, mapbox::geometry::point<T> const& rhs) {
+    return lhs.x != rhs.x || lhs.y != rhs.y;
+}
+
+#ifdef DEBUG
+
+template <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const point<T>& p) {
+    out << "  point at: " << p.x << ", " << p.y;
+    return out;
+}
+
+template <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const mapbox::geometry::point<T>& p) {
+    out << "  point at: " << p.x << ", " << p.y;
+    return out;
+}
+#endif
+}
+}
+}

mapbox/geometry/wagyu/process_horizontal.hpp

@@ -0,0 +1,285 @@
+#pragma once
+
+#include <mapbox/geometry/line_string.hpp>
+#include <mapbox/geometry/point.hpp>
+#include <mapbox/geometry/polygon.hpp>
+
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/exceptions.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+active_bound_list_itr<T> process_horizontal_left_to_right(T scanline_y,
+                                                          active_bound_list_itr<T>& horz_bound,
+                                                          active_bound_list<T>& active_bounds,
+                                                          ring_manager<T>& rings,
+                                                          scanbeam_list<T>& 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<T>(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<double>((*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<T>(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<T>(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<T>(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 <typename T>
+active_bound_list_itr<T> process_horizontal_right_to_left(T scanline_y,
+                                                          active_bound_list_itr<T>& horz_bound,
+                                                          active_bound_list<T>& active_bounds,
+                                                          ring_manager<T>& rings,
+                                                          scanbeam_list<T>& 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<T>(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<T>(hp_itr_fwd);
+
+    auto bnd = active_bound_list_rev_itr<T>(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<double>((*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<T>(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<T>(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 <typename T>
+active_bound_list_itr<T> process_horizontal(T scanline_y,
+                                            active_bound_list_itr<T>& horz_bound,
+                                            active_bound_list<T>& active_bounds,
+                                            ring_manager<T>& rings,
+                                            scanbeam_list<T>& 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 <typename T>
+void process_horizontals(T scanline_y,
+                         active_bound_list<T>& active_bounds,
+                         ring_manager<T>& rings,
+                         scanbeam_list<T>& 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<T>(bnd_itr)) {
+            bnd_itr = process_horizontal(scanline_y, bnd_itr, active_bounds, rings, scanbeam,
+                                         cliptype, subject_fill_type, clip_fill_type);
+        } else {
+            ++bnd_itr;
+        }
+    }
+}
+}
+}
+}

mapbox/geometry/wagyu/process_maxima.hpp

@@ -0,0 +1,133 @@
+#pragma once
+
+#include <mapbox/geometry/wagyu/active_bound_list.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/exceptions.hpp>
+#include <mapbox/geometry/wagyu/intersect_util.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+#include <mapbox/geometry/wagyu/local_minimum_util.hpp>
+#include <mapbox/geometry/wagyu/process_horizontal.hpp>
+#include <mapbox/geometry/wagyu/ring.hpp>
+#include <mapbox/geometry/wagyu/ring_util.hpp>
+#include <mapbox/geometry/wagyu/topology_correction.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+active_bound_list_itr<T> do_maxima(active_bound_list_itr<T>& bnd,
+                                   active_bound_list_itr<T>& bndMaxPair,
+                                   clip_type cliptype,
+                                   fill_type subject_fill_type,
+                                   fill_type clip_fill_type,
+                                   ring_manager<T>& rings,
+                                   active_bound_list<T>& 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 <typename T>
+void process_edges_at_top_of_scanbeam(T top_y,
+                                      active_bound_list<T>& active_bounds,
+                                      scanbeam_list<T>& scanbeam,
+                                      local_minimum_ptr_list<T> const& minima_sorted,
+                                      local_minimum_ptr_list_itr<T>& current_lm,
+                                      ring_manager<T>& 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<T>(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<T>(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);
+        }
+    }
+}
+
+}
+}
+}

mapbox/geometry/wagyu/ring.hpp

@@ -0,0 +1,479 @@
+#pragma once
+
+#include <assert.h>
+#include <cmath>
+#include <list>
+#include <map>
+#include <set>
+#include <sstream>
+#include <vector>
+#include <deque>
+#include <mapbox/geometry/wagyu/point.hpp>
+
+#ifdef DEBUG
+#include <execinfo.h>
+#include <iostream>
+#include <stdio.h>
+//
+// 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 <typename T>
+using ring_vector = std::vector<ring_ptr<T>>;
+
+template <typename T>
+using ring_list = std::list<ring_ptr<T>>;
+
+template <typename T>
+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<T> parent;
+    ring_list<T> children;
+    point_ptr<T> points;
+    point_ptr<T> bottom_point;
+    bool is_open;
+
+    ring( ring const& ) = delete;
+    ring& operator=(ring const& ) = delete;
+
+    ring()
+        : ring_index(0),
+          size(0),
+          area(std::numeric_limits<double>::quiet_NaN()),
+          parent(nullptr),
+          children(),
+          points(nullptr),
+          bottom_point(nullptr),
+          is_open(false) {
+    }
+};
+
+template <typename T>
+using hot_pixel_vector = std::vector<mapbox::geometry::point<T>>;
+
+template <typename T>
+using hot_pixel_itr = typename hot_pixel_vector<T>::iterator;
+
+template <typename T>
+using hot_pixel_rev_itr = typename hot_pixel_vector<T>::reverse_iterator;
+
+template <typename T>
+struct ring_manager {
+    
+    ring_list<T> children;
+    std::vector<point_ptr<T>> all_points;
+    hot_pixel_vector<T> hot_pixels;
+    hot_pixel_itr<T> current_hp_itr;
+    std::deque<point<T>> points;
+    std::deque<ring<T>> rings;
+    std::vector<point<T>> 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 <typename T>
+void preallocate_point_memory(ring_manager<T>& rings, std::size_t size) {
+    rings.storage.reserve(size);
+    rings.all_points.reserve(size);
+}
+
+template <typename T>
+ring_ptr<T> create_new_ring(ring_manager<T>& rings) {
+    rings.rings.emplace_back();
+    ring_ptr<T> result = &rings.rings.back();
+    result->ring_index = rings.index++;
+    return result;
+}
+
+template <typename T>
+point_ptr<T>
+create_new_point(ring_ptr<T> r, mapbox::geometry::point<T> const& pt, ring_manager<T>& rings) {
+    point_ptr<T> 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 <typename T>
+point_ptr<T> create_new_point(ring_ptr<T> r,
+                              mapbox::geometry::point<T> const& pt,
+                              point_ptr<T> before_this_point,
+                              ring_manager<T>& rings) {
+    point_ptr<T> 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 <typename T>
+void ring1_child_of_ring2(ring_ptr<T> ring1, ring_ptr<T> ring2, ring_manager<T>& 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 <typename T>
+void ring1_sibling_of_ring2(ring_ptr<T> ring1, ring_ptr<T> ring2, ring_manager<T>& 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 <typename T>
+void ring1_replaces_ring2(ring_ptr<T> ring1, ring_ptr<T> ring2, ring_manager<T>& 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 <typename T>
+void remove_ring(ring_ptr<T> r, ring_manager<T>& 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 <typename T>
+inline std::size_t ring_depth(ring_ptr<T> r) {
+    std::size_t depth = 0;
+    if (!r) {
+        return depth;
+    }
+    while (r->parent) {
+        depth++;
+        r = r->parent;
+    }
+    return depth;
+}
+
+template <typename T>
+inline bool ring_is_hole(ring_ptr<T> r) {
+    return ring_depth(r) & 1;
+}
+
+template <typename T>
+void set_next(const_point_ptr<T>& node, const const_point_ptr<T>& next_node) {
+    node->next = next_node;
+}
+
+template <typename T>
+point_ptr<T> get_next(const_point_ptr<T>& node) {
+    return node->next;
+}
+
+template <typename T>
+point_ptr<T> get_prev(const_point_ptr<T>& node) {
+    return node->prev;
+}
+
+template <typename T>
+void set_prev(const_point_ptr<T>& node, const const_point_ptr<T>& prev_node) {
+    node->prev = prev_node;
+}
+
+template <typename T>
+void init(const_point_ptr<T>& node) {
+    set_next(node, node);
+    set_prev(node, node);
+}
+
+template <typename T>
+std::size_t point_count(const const_point_ptr<T>& orig_node) {
+    std::size_t size = 0;
+    point_ptr<T> n = orig_node;
+    do {
+        n = get_next(n);
+        ++size;
+    } while (n != orig_node);
+    return size;
+}
+
+template <typename T>
+void link_before(point_ptr<T>& node, point_ptr<T>& new_node) {
+    point_ptr<T> 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 <typename T>
+void link_after(point_ptr<T>& node, point_ptr<T>& new_node) {
+    point_ptr<T> 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 <typename T>
+void transfer_point(point_ptr<T>& p, point_ptr<T>& b, point_ptr<T>& e) {
+    if (b != e) {
+        point_ptr<T> prev_p = get_prev(p);
+        point_ptr<T> prev_b = get_prev(b);
+        point_ptr<T> 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 <typename T>
+void reverse_ring(point_ptr<T> pp) {
+    if (!pp) {
+        return;
+    }
+    point_ptr<T> pp1;
+    point_ptr<T> pp2;
+    pp1 = pp;
+    do {
+        pp2 = pp1->next;
+        pp1->next = pp1->prev;
+        pp1->prev = pp2;
+        pp1 = pp2;
+    } while (pp1 != pp);
+}
+
+template <typename T>
+double area_from_point(point_ptr<T> op, std::size_t & size) {
+    point_ptr<T> startOp = op;
+    size = 1;
+    if (!op) {
+        return 0.0;
+    }
+    double a = 0.0;
+    do {
+        ++size;
+        a += static_cast<double>(op->prev->x + op->x) * static_cast<double>(op->prev->y - op->y);
+        op = op->next;
+    } while (op != startOp);
+    return a * 0.5;
+}
+
+template <typename T>
+double area(ring_ptr<T> r) {
+    assert(r != nullptr);
+    if (std::isnan(r->area)) {
+        r->area = area_from_point(r->points, r->size);
+    }
+    return r->area;
+}
+
+#ifdef DEBUG
+
+template <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const ring<T>& 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<T> n = const_cast<ring_ptr<T>>(&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 <typename T>
+std::string output_as_polygon(ring_ptr<T> 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const ring_list<T>& 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const ring_vector<T>& 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const std::deque<ring<T>>& 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 <class charT, class traits, typename T>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& out,
+                                                     const hot_pixel_vector<T>& hp_vec) {
+    out << "Hot Pixels: " << std::endl;
+    for (auto& hp : hp_vec) {
+        out << hp << std::endl;
+    }
+    return out;
+}
+#endif
+}
+}
+}

mapbox/geometry/wagyu/ring_util.hpp

@@ -0,0 +1,915 @@
+#pragma once
+
+#ifdef DEBUG
+#include <iostream>
+// Example debug print for backtrace - only works on IOS
+#include <execinfo.h>
+#include <stdio.h>
+//
+// 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 <queue>
+
+#include <mapbox/geometry/wagyu/active_bound_list.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/ring.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+void set_hole_state(active_bound_list_itr<T>& bnd,
+                    active_bound_list<T>& active_bounds,
+                    ring_manager<T>& rings) {
+    auto bnd2 = active_bound_list_rev_itr<T>(bnd);
+    bound_ptr<T> 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 <typename T>
+void set_hole_state(active_bound_list_rev_itr<T>& bnd,
+                    active_bound_list<T>& active_bounds,
+                    ring_manager<T>& rings) {
+    auto bnd2 = std::next(bnd);
+    bound_ptr<T> 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 <typename T>
+void update_current_hp_itr(T scanline_y, ring_manager<T>& rings) {
+    while (rings.current_hp_itr->y > scanline_y) {
+        ++rings.current_hp_itr;       
+    }
+}
+
+template <typename T>
+struct hot_pixel_sorter {
+    inline bool operator()(mapbox::geometry::point<T> const& pt1, mapbox::geometry::point<T> 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 <typename T>
+T round_towards_min(double val) {
+    // 0.5 rounds to 0
+    // 0.0 rounds to 0
+    // -0.5 rounds to -1
+    return static_cast<T>(std::ceil(val - 0.5 + rounding_offset));
+}
+
+template <typename T>
+T round_towards_max(double val) {
+    // 0.5 rounds to 1
+    // 0.0 rounds to 0
+    // -0.5 rounds to 0
+    return static_cast<T>(std::floor(val + 0.5 + rounding_offset));
+}
+
+template <typename T>
+inline T get_edge_min_x(edge<T> 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<double>(current_y - edge.bot.y) - 0.5;
+            double return_val = static_cast<double>(edge.bot.x) + edge.dx * lower_range_y;
+            T value = round_towards_min<T>(return_val);
+            return value;
+        }
+    } else {
+        if (current_y == edge.bot.y) {
+            return edge.bot.x;
+        } else {
+            double return_val = static_cast<double>(edge.bot.x) +
+                                edge.dx * (static_cast<double>(current_y - edge.bot.y) + 0.5 - rounding_offset_y);
+            T value = round_towards_min<T>(return_val);
+            return value;
+        }
+    }
+}
+
+template <typename T>
+inline T get_edge_max_x(edge<T> 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<double>(current_y - edge.bot.y) - 0.5;
+            double return_val = static_cast<double>(edge.bot.x) + edge.dx * lower_range_y;
+            T value = round_towards_max<T>(return_val);
+            return value;
+        }
+    } else {
+        if (current_y == edge.bot.y) {
+            return edge.bot.x;
+        } else {
+            double return_val = static_cast<double>(edge.bot.x) +
+                                edge.dx * (static_cast<double>(current_y - edge.bot.y) + 0.5 - rounding_offset_y);
+            T value = round_towards_max<T>(return_val);
+            return value;
+        }
+    }
+}
+
+template <typename T>
+void hot_pixel_set_left_to_right(T y,
+                                 T start_x,
+                                 T end_x,
+                                 bound<T>& bnd,
+                                 ring_manager<T>& rings,
+                                 hot_pixel_itr<T> & itr,
+                                 hot_pixel_itr<T> & 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<T> 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<T> new_point = create_new_point(bnd.ring, *itr, op, rings);
+        if (to_front) {
+            bnd.ring->points = new_point;
+        }
+    }
+}
+
+template <typename T>
+void hot_pixel_set_right_to_left(T y,
+                                 T start_x,
+                                 T end_x,
+                                 bound<T>& bnd,
+                                 ring_manager<T>& rings,
+                                 hot_pixel_rev_itr<T> & itr,
+                                 hot_pixel_rev_itr<T> & 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<T> 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<T> new_point = create_new_point(bnd.ring, *itr, op, rings);
+        if (to_front) {
+            bnd.ring->points = new_point;
+        }
+    }
+}
+
+template <typename T>
+void sort_hot_pixels(ring_manager<T>& rings) {
+    std::sort(rings.hot_pixels.begin(), rings.hot_pixels.end(), hot_pixel_sorter<T>());
+    auto last = std::unique(rings.hot_pixels.begin(), rings.hot_pixels.end());
+    rings.hot_pixels.erase(last, rings.hot_pixels.end());
+}
+
+template <typename T>
+void insert_hot_pixels_in_path(bound<T>& bnd,
+                               mapbox::geometry::point<T> const& end_pt,
+                               ring_manager<T>& 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<T>(itr);
+            while (itr != rings.hot_pixels.end() && itr->y == y) {
+                ++itr;
+            }
+            auto first_itr = hot_pixel_rev_itr<T>(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 <typename T>
+void add_to_hot_pixels(mapbox::geometry::point<T> const& pt, ring_manager<T>& rings) {
+    rings.hot_pixels.push_back(pt);
+}
+
+template <typename T>
+void add_first_point(active_bound_list_itr<T>& bnd,
+                     active_bound_list<T>& active_bounds,
+                     mapbox::geometry::point<T> const& pt,
+                     ring_manager<T>& rings) {
+
+    ring_ptr<T> 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 <typename T>
+void add_first_point(active_bound_list_rev_itr<T>& bnd,
+                     active_bound_list<T>& active_bounds,
+                     mapbox::geometry::point<T> const& pt,
+                     ring_manager<T>& rings) {
+    ring_ptr<T> 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 <typename T>
+void add_point_to_ring(bound<T>& bnd,
+                       mapbox::geometry::point<T> const& pt,
+                       ring_manager<T>& 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<T> 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<T> new_point = create_new_point(bnd.ring, pt, bnd.ring->points, rings);
+    if (to_front) {
+        bnd.ring->points = new_point;
+    }
+}
+
+template <typename T>
+void add_point(active_bound_list_itr<T>& bnd,
+               active_bound_list<T>& active_bounds,
+               mapbox::geometry::point<T> const& pt,
+               ring_manager<T>& rings) {
+    if (!(*bnd)->ring) {
+        add_first_point(bnd, active_bounds, pt, rings);
+    } else {
+        add_point_to_ring(*(*bnd), pt, rings);
+    }
+}
+
+template <typename T>
+void add_point(active_bound_list_rev_itr<T>& bnd,
+               active_bound_list<T>& active_bounds,
+               mapbox::geometry::point<T> const& pt,
+               ring_manager<T>& rings) {
+    if (!(*bnd)->ring) {
+        add_first_point(bnd, active_bounds, pt, rings);
+    } else {
+        add_point_to_ring(*(*bnd), pt, rings);
+    }
+}
+
+template <typename T>
+void add_local_minimum_point(active_bound_list_itr<T> b1,
+                             active_bound_list_itr<T> b2,
+                             active_bound_list<T>& active_bounds,
+                             mapbox::geometry::point<T> const& pt,
+                             ring_manager<T>& rings) {
+    active_bound_list_itr<T> b;
+    active_bound_list_rev_itr<T> prev_bound;
+    active_bound_list_rev_itr<T> prev_b1(b1);
+    active_bound_list_rev_itr<T> 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 <typename T>
+inline double get_dx(point<T> const& pt1, point<T> const& pt2) {
+    if (pt1.y == pt2.y) {
+        return std::numeric_limits<double>::infinity();
+    } else {
+        return static_cast<double>(pt2.x - pt2.x) / static_cast<double>(pt2.y - pt1.y);
+    }
+}
+
+template <typename T>
+bool first_is_bottom_point(const_point_ptr<T> btmPt1, const_point_ptr<T> btmPt2) {
+    point_ptr<T> 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 <typename T>
+point_ptr<T> get_bottom_point(point_ptr<T> pp) {
+    point_ptr<T> dups = 0;
+    point_ptr<T> 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 <typename T>
+ring_ptr<T> get_lower_most_ring(ring_ptr<T> outRec1, ring_ptr<T> 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<T> OutPt1 = outRec1->bottom_point;
+    point_ptr<T> 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 <typename T>
+bool ring1_right_of_ring2(ring_ptr<T> ring1, ring_ptr<T> ring2) {
+    do {
+        ring1 = ring1->parent;
+        if (ring1 == ring2) {
+            return true;
+        }
+    } while (ring1);
+    return false;
+}
+
+template <typename T>
+void update_points_ring(ring_ptr<T> ring) {
+    point_ptr<T> op = ring->points;
+    do {
+        op->ring = ring;
+        op = op->prev;
+    } while (op != ring->points);
+}
+
+template <typename T>
+void append_ring(active_bound_list_itr<T>& b1,
+                 active_bound_list_itr<T>& b2,
+                 active_bound_list<T>& active_bounds,
+                 ring_manager<T>& manager) {
+    // get the start and ends of both output polygons ...
+    ring_ptr<T> outRec1 = (*b1)->ring;
+    ring_ptr<T> outRec2 = (*b2)->ring;
+
+    ring_ptr<T> keep_ring;
+    bound_ptr<T> keep_bound;
+    ring_ptr<T> remove_ring;
+    bound_ptr<T> 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<T> p1_lft = keep_ring->points;
+    point_ptr<T> p1_rt = p1_lft->prev;
+    point_ptr<T> p2_lft = remove_ring->points;
+    point_ptr<T> 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 <typename T>
+void add_local_maximum_point(active_bound_list_itr<T>& b1,
+                             active_bound_list_itr<T>& b2,
+                             mapbox::geometry::point<T> const& pt,
+                             ring_manager<T>& rings,
+                             active_bound_list<T>& 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 <typename T>
+point_in_polygon_result point_in_polygon(point<T> const& pt, point_ptr<T> op) {
+    // returns 0 if false, +1 if true, -1 if pt ON polygon boundary
+    point_in_polygon_result result = point_outside_polygon;
+    point_ptr<T> 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<double>(op->x - pt.x) *
+                            static_cast<double>(op->next->y - pt.y) -
+                        static_cast<double>(op->next->x - pt.x) * static_cast<double>(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<double>(op->x - pt.x) *
+                            static_cast<double>(op->next->y - pt.y) -
+                        static_cast<double>(op->next->x - pt.x) * static_cast<double>(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 <typename T>
+point_in_polygon_result point_in_polygon(mapbox::geometry::point<double> const& pt,
+                                         point_ptr<T> op) {
+    // returns 0 if false, +1 if true, -1 if pt ON polygon boundary
+    point_in_polygon_result result = point_outside_polygon;
+    point_ptr<T> startOp = op;
+    do {
+        double op_x = static_cast<double>(op->x);
+        double op_y = static_cast<double>(op->y);
+        double op_next_x = static_cast<double>(op->next->x);
+        double op_next_y = static_cast<double>(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 <typename T>
+point_in_polygon_result inside_or_outside_special(point_ptr<T> first_pt, point_ptr<T> 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<T> 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<double> 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 <typename T>
+bool poly2_contains_poly1(point_ptr<T> outpt1, point_ptr<T> outpt2) {
+    point_ptr<T> 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 <typename T>
+void dispose_out_points(point_ptr<T>& pp) {
+    if (pp == nullptr) {
+        return;
+    }
+    pp->prev->next = nullptr;
+    while (pp) {
+        point_ptr<T> tmpPp = pp;
+        pp = pp->next;
+        tmpPp->next = tmpPp;
+        tmpPp->prev = tmpPp;
+        tmpPp->ring = nullptr;
+        // delete tmpPp;
+    }
+}
+}
+}
+}

mapbox/geometry/wagyu/scanbeam.hpp

@@ -0,0 +1,37 @@
+#pragma once
+
+#include <queue>
+
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+using scanbeam_list = std::priority_queue<T>;
+
+template <typename T>
+bool pop_from_scanbeam(T& Y, scanbeam_list<T>& scanbeam) {
+    if (scanbeam.empty()) {
+        return false;
+    }
+    Y = scanbeam.top();
+    scanbeam.pop();
+    while (!scanbeam.empty() && Y == scanbeam.top()) {
+        scanbeam.pop();
+    } // Pop duplicates.
+    return true;
+}
+
+template <typename T>
+void setup_scanbeam(local_minimum_list<T>& minima_list, scanbeam_list<T>& scanbeam) {
+
+    for (auto lm = minima_list.begin(); lm != minima_list.end(); ++lm) {
+        scanbeam.push(lm->y);
+    }
+}
+}
+}
+}

mapbox/geometry/wagyu/snap_rounding.hpp

@@ -0,0 +1,174 @@
+#pragma once
+
+#include <mapbox/geometry/wagyu/active_bound_list.hpp>
+#include <mapbox/geometry/wagyu/edge.hpp>
+#include <mapbox/geometry/wagyu/bound.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/intersect.hpp>
+#include <mapbox/geometry/wagyu/intersect_util.hpp>
+#include <mapbox/geometry/wagyu/ring.hpp>
+#include <mapbox/geometry/wagyu/ring_util.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+void process_hot_pixel_intersections(T top_y,
+                                     active_bound_list<T>& active_bounds,
+                                     ring_manager<T>& 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<double> pt;
+                if (!get_edge_intersection<T, double>(*((*bnd)->current_edge),
+                                                     *((*bnd_next)->current_edge), pt)) {
+                    throw std::runtime_error("Edges do not intersect!");
+                }
+                add_to_hot_pixels(round_point<T>(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 <typename T>
+void process_hot_pixel_edges_at_top_of_scanbeam(T top_y,
+                                                scanbeam_list<T>& scanbeam,
+                                                active_bound_list<T>& active_bounds,
+                                                ring_manager<T>& 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<T>(bnd)) {
+                (*bnd)->current_x = static_cast<double>((*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<T> 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<T> 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 <typename T>
+void insert_local_minima_into_ABL_hot_pixel(T top_y,
+                                            local_minimum_ptr_list<T> & minima_sorted,
+                                            local_minimum_ptr_list_itr<T> & lm,
+                                            active_bound_list<T>& active_bounds,
+                                            ring_manager<T>& rings,
+                                            scanbeam_list<T>& 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<double>(left_bound.current_edge->bot.x);
+        auto lb_abl_itr = insert_bound_into_ABL(left_bound, active_bounds);
+        if (!current_edge_is_horizontal<T>(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<double>(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<T>(rb_abl_itr)) {
+            scanbeam.push((*rb_abl_itr)->current_edge->top.y);
+        }
+        ++lm;
+    }
+}
+
+template <typename T>
+void build_hot_pixels(local_minimum_list<T>& minima_list,
+                      ring_manager<T>& rings) {
+    if (minima_list.empty()) {
+        return;
+    }
+    
+    active_bound_list<T> active_bounds;
+    scanbeam_list<T> scanbeam;
+    T scanline_y = std::numeric_limits<T>::max();
+
+    local_minimum_ptr_list<T> 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<T>());
+    local_minimum_ptr_list_itr<T> 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);
+}
+
+}
+}
+}

mapbox/geometry/wagyu/util.hpp

@@ -0,0 +1,83 @@
+#pragma once
+
+#include <cmath>
+
+#include <mapbox/geometry/point.hpp>
+#include <mapbox/geometry/polygon.hpp>
+#include <mapbox/geometry/wagyu/point.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+double area(mapbox::geometry::linear_ring<T> 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<double>(itr_prev->x + itr->x) * static_cast<double>(itr_prev->y - itr->y);
+    ++itr;
+    itr_prev = poly.begin();
+    for (; itr != poly.end(); ++itr, ++itr_prev) {
+        a += static_cast<double>(itr_prev->x + itr->x) * static_cast<double>(itr_prev->y - itr->y);
+    }
+    return -a * 0.5;
+}
+
+inline bool value_is_zero(double val) {
+    return std::fabs(val) < std::numeric_limits<double>::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<double>::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 <typename T>
+bool slopes_equal(mapbox::geometry::point<T> const& pt1,
+                  mapbox::geometry::point<T> const& pt2,
+                  mapbox::geometry::point<T> const& pt3) {
+    return (pt1.y - pt2.y) * (pt2.x - pt3.x) == (pt1.x - pt2.x) * (pt2.y - pt3.y);
+}
+
+template <typename T>
+bool slopes_equal(mapbox::geometry::wagyu::point<T> const& pt1,
+                  mapbox::geometry::wagyu::point<T> const& pt2,
+                  mapbox::geometry::point<T> const& pt3) {
+    return (pt1.y - pt2.y) * (pt2.x - pt3.x) == (pt1.x - pt2.x) * (pt2.y - pt3.y);
+}
+
+template <typename T>
+bool slopes_equal(mapbox::geometry::wagyu::point<T> const& pt1,
+                  mapbox::geometry::wagyu::point<T> const& pt2,
+                  mapbox::geometry::wagyu::point<T> const& pt3) {
+    return (pt1.y - pt2.y) * (pt2.x - pt3.x) == (pt1.x - pt2.x) * (pt2.y - pt3.y);
+}
+
+template <typename T>
+bool slopes_equal(mapbox::geometry::point<T> const& pt1,
+                  mapbox::geometry::point<T> const& pt2,
+                  mapbox::geometry::point<T> const& pt3,
+                  mapbox::geometry::point<T> const& pt4) {
+    return (pt1.y - pt2.y) * (pt3.x - pt4.x) == (pt1.x - pt2.x) * (pt3.y - pt4.y);
+}
+}
+}
+}

mapbox/geometry/wagyu/vatti.hpp

@@ -0,0 +1,77 @@
+#pragma once
+
+#include <algorithm>
+#include <set>
+
+#include <mapbox/geometry/wagyu/active_bound_list.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/exceptions.hpp>
+#include <mapbox/geometry/wagyu/intersect_util.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+#include <mapbox/geometry/wagyu/local_minimum_util.hpp>
+#include <mapbox/geometry/wagyu/process_horizontal.hpp>
+#include <mapbox/geometry/wagyu/process_maxima.hpp>
+#include <mapbox/geometry/wagyu/ring.hpp>
+#include <mapbox/geometry/wagyu/ring_util.hpp>
+#include <mapbox/geometry/wagyu/util.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+bool execute_vatti(local_minimum_list<T>& minima_list,
+                   ring_manager<T>& rings,
+                   clip_type cliptype,
+                   fill_type subject_fill_type,
+                   fill_type clip_fill_type) {
+
+    if (minima_list.empty()) {
+        return false;
+    }
+    
+    active_bound_list<T> active_bounds;
+    scanbeam_list<T> scanbeam;
+    T scanline_y = std::numeric_limits<T>::max();
+
+    local_minimum_ptr_list<T> 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<T>());
+    local_minimum_ptr_list_itr<T> 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;
+}
+}
+}
+}

mapbox/geometry/wagyu/wagyu.hpp

@@ -0,0 +1,141 @@
+#pragma once
+
+#include <list>
+
+#include <mapbox/geometry/box.hpp>
+#include <mapbox/geometry/line_string.hpp>
+#include <mapbox/geometry/multi_polygon.hpp>
+#include <mapbox/geometry/polygon.hpp>
+
+#include <mapbox/geometry/wagyu/build_result.hpp>
+#include <mapbox/geometry/wagyu/build_local_minima_list.hpp>
+#include <mapbox/geometry/wagyu/config.hpp>
+#include <mapbox/geometry/wagyu/local_minimum.hpp>
+#include <mapbox/geometry/wagyu/snap_rounding.hpp>
+#include <mapbox/geometry/wagyu/topology_correction.hpp>
+#include <mapbox/geometry/wagyu/vatti.hpp>
+
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+
+template <typename T>
+class wagyu {
+private:
+    using value_type = T;
+
+    local_minimum_list<value_type> 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<value_type> 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<value_type> const& pg,
+                  polygon_type p_type = polygon_type_subject) {
+        return add_linear_ring(pg, minima_list, p_type);
+    }
+
+    bool add_polygon(mapbox::geometry::polygon<value_type> 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<value_type> get_bounds() {
+        mapbox::geometry::point<value_type> min = { 0, 0 };
+        mapbox::geometry::point<value_type> max = { 0, 0 };
+        if (minima_list.empty()) {
+            return mapbox::geometry::box<value_type>(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<value_type>(min, max);
+    }
+
+    bool execute(clip_type cliptype,
+                 mapbox::geometry::multi_polygon<value_type>& solution,
+                 fill_type subject_fill_type,
+                 fill_type clip_fill_type) {
+
+        ring_manager<T> 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;
+    }
+};
+}
+}
+}