Add optional fallback polygon clipping using `wagyu`

[artemp]

Fallback on wagyu if boost::geometry::intersection fails. Could be important if dealing with VTs < v2.0

[springmeyer]

We are lucky that #84 and #82 were able to handle all the v1 edge cases we encountered, so closing this!

[springmeyer]

Adding @artemp's diff for the record:

diff --git a/src/feature_builder.hpp b/src/feature_builder.hpp
index 8196f8f..7ffc443 100644
--- a/src/feature_builder.hpp
+++ b/src/feature_builder.hpp
@@ -10,11 +10,37 @@
 #include <boost/core/ignore_unused.hpp>
 #include <boost/geometry/algorithms/intersects.hpp>
 #include <boost/geometry/algorithms/intersection.hpp>
+// wagyu
+#include <mapbox/geometry/wagyu/quick_clip.hpp>
+#include <mapbox/geometry/wagyu/wagyu.hpp>
 
 //BOOST_GEOMETRY_REGISTER_POINT_2D(mapbox::geometry::point<int>, int, boost::geometry::cs::cartesian, x, y)
 // ^ Uncomment to enable coordinate_type = int ^
 
-namespace vtile {
+
+#define USE_WAGYU
+
+namespace vtile { namespace detail {
+
+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;
+}
+}
 
 template <typename CoordinateType>
 struct feature_builder_visitor
@@ -111,7 +137,67 @@ struct feature_builder_visitor
         }
         if (valid) finalize(feature_builder);
     }
+#if defined(USE_WAGYU)
+
+    void operator()(mapbox::geometry::polygon<coordinate_type> & poly)
+    {
+        if (poly.empty() || (poly.front().size() < 3))
+        {
+            return;
+        }
+
+        mapbox::geometry::wagyu::wagyu<coordinate_type> clipper;
+        bool first = true;
+        for (auto & ring : poly)
+        {
+            if (ring.size() < 3) continue;
+            double area = detail::area(ring);
+            if (first)
+            {
+                first = false;
+                if (area < 0) std::reverse(ring.begin(), ring.end());
+
+                auto new_ring = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox_);
+                if (new_ring.empty()) return;
+                clipper.add_ring(new_ring);
+            }
+            else
+            {
+                if (area > 0) std::reverse(ring.begin(), ring.end());
+                auto new_ring = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox_);
+                if (new_ring.empty()) continue;
+                clipper.add_ring(new_ring);
+            }
+        }
 
+        mapbox::geometry::multi_polygon<std::int64_t> result;
+
+        clipper.execute(mapbox::geometry::wagyu::clip_type_union,
+                        result,
+                        mapbox::geometry::wagyu::fill_type_even_odd,
+                        mapbox::geometry::wagyu::fill_type_even_odd);
+
+        vtzero::polygon_feature_builder feature_builder{layer_builder_};
+        if (feature_.has_id()) feature_builder.set_id(feature_.id());
+        bool valid = false;
+        for (auto const& p : result)
+        {
+            for (auto const& ring : p)
+            {
+                if (ring.size() > 3)
+                {
+                    valid = true;
+                    feature_builder.add_ring(static_cast<unsigned>(ring.size()));
+                    for (auto const& pt : ring)
+                    {
+                        feature_builder.set_point(static_cast<int>(pt.x), static_cast<int>(pt.y));
+                    }
+                }
+            }
+        }
+        if (valid) finalize(feature_builder);
+    }
+#else
     void operator()(mapbox::geometry::polygon<coordinate_type> const& poly)
     {
         auto extent = mapbox::geometry::envelope(poly);
@@ -138,7 +224,65 @@ struct feature_builder_visitor
         }
         if (valid) finalize(feature_builder);
     }
+#endif
 
+#if defined(USE_WAGYU)
+
+    void operator()(mapbox::geometry::multi_polygon<coordinate_type> & mpoly)
+    {
+        if (mpoly.empty()) return;
+        //mapbox::geometry::multi_polygon<coordinate_type>> result;
+        vtzero::polygon_feature_builder feature_builder{layer_builder_};
+        if (feature_.has_id()) feature_builder.set_id(feature_.id());
+        bool valid = false;
+        for (auto & poly : mpoly)
+        {
+            mapbox::geometry::wagyu::wagyu<std::int64_t> clipper;
+            bool first = true;
+            for (auto & ring : poly)
+            {
+                if (ring.size() < 3) continue;
+                double area = detail::area(ring);
+                if (first)
+                {
+                    first = false;
+                    if (area < 0) std::reverse(ring.begin(), ring.end());
+                    auto new_ring = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox_);
+                    if (new_ring.empty()) return;
+                    clipper.add_ring(new_ring);
+                }
+                else
+                {
+                    if (area > 0) std::reverse(ring.begin(), ring.end());
+                    auto new_ring = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox_);
+                    if (new_ring.empty()) continue;
+                    clipper.add_ring(new_ring);
+                }
+            }
+            mapbox::geometry::multi_polygon<std::int64_t> result;
+            clipper.execute(mapbox::geometry::wagyu::clip_type_union,
+                            result,
+                            mapbox::geometry::wagyu::fill_type_even_odd,
+                            mapbox::geometry::wagyu::fill_type_even_odd);
+            for (auto const& p : result)
+            {
+                for (auto const& ring : p)
+                {
+                    if (ring.size() > 3)
+                    {
+                        valid = true;
+                        feature_builder.add_ring(static_cast<unsigned>(ring.size()));
+                        for (auto const& pt : ring)
+                        {
+                            feature_builder.set_point(static_cast<int>(pt.x), static_cast<int>(pt.y));
+                        }
+                    }
+                }
+            }
+        }
+        if (valid) finalize(feature_builder);
+    }
+#else
     void operator()(mapbox::geometry::multi_polygon<coordinate_type> const& mpoly)
     {
         std::vector<mapbox::geometry::polygon<coordinate_type>> result;
@@ -169,6 +313,7 @@ struct feature_builder_visitor
         }
         if (valid) finalize(feature_builder);
     }
+#endif
 
     template <typename T>
     void operator()(T const& geom)