Update wagyu to 0.5.0

Closes https://github.com/postgis/postgis/pull/563
Closes #4620

LICENSE.TXT

@@ -3,7 +3,7 @@ PostGIS Licensing
 
 This file attempts to include all licenses that apply within the PostGIS
 source tree, in particular any that are supposed to be exposed to the end
-user for credit requirements for instance. 
+user for credit requirements for instance.
 
 PostGIS General
 ----------------
@@ -44,7 +44,8 @@ license. Copyright (c) 2005-2018, Troy D. Hanson.
 
 -- deps/wagyu (see LICENSE.README in folder) includes code released under the
 ISC license (copyrighted by Mapbox) and the Boost Software License (copyrighted
-by Angus Johnson and Mapbox).
+by Angus Johnson and Mapbox). Copyright for "almost_equal.hpp" file is held by
+Google Inc and its license is listed at the top of that file.
 
 The documentation for PostGIS is under a creative commons share-alike
 3.0 license.  http://creativecommons.org/licenses/by-sa/3.0/
@@ -58,4 +59,4 @@ Data used in documentation falls in one of the following categories
   utilizes ImageMagick) as described
   http://trac.osgeo.org/postgis/wiki/DevWikiDocNewFeature
 - PostGIS raster output functions such as the ST_AsPNG etc.
-- raster/vector data from MassGIS: http://www.mass.gov/mgis/laylist.htm 
+- raster/vector data from MassGIS: http://www.mass.gov/mgis/laylist.htm

NEWS

@@ -19,6 +19,7 @@ Only tickets not included in 3.1.0alpha1
   - #4672, Cache getSRSbySRID and getSRIDbySRS (Raúl Marín)
   - #4676, Avoid decompressing toasted geometries to read only the header (Raúl Marín)
            Optimize cast to Postgresql point type (Raúl Marín)
+  - #4620, Update internal wagyu to 0.5.0 (Raúl Marín)
 
 * Bug fixes *
   - #4652, Fix several memory related bugs in ST_GeomFromGML (Raúl Marín)

deps/wagyu/README.md

@@ -25,3 +25,4 @@ It is integrated in the project as an static library inside postgis.so, so it do
 
   - 2019-02-05 - [Wagyu] Library extraction from https://github.com/mapbox/wagyu
   - 2019-02-05 - [geometry.hpp] Library extraction from https://github.com/mapbox/geometry.hpp
+  - 2020-05-08 - [Wagyu] Update to 0.5.0

deps/wagyu/include/mapbox/geometry/wagyu/active_bound_list.hpp

@@ -90,11 +90,11 @@ struct bound_insert_location {
         auto const& bound1 = *b;
         if (values_are_equal(bound2.current_x, bound1.current_x)) {
             if (bound2.current_edge->top.y > bound1.current_edge->top.y) {
-                return static_cast<double>(bound2.current_edge->top.x) <
-                       get_current_x(*(bound1.current_edge), bound2.current_edge->top.y);
+                return less_than(static_cast<double>(bound2.current_edge->top.x),
+                                 get_current_x(*(bound1.current_edge), bound2.current_edge->top.y));
             } else {
-                return static_cast<double>(bound1.current_edge->top.x) >
-                       get_current_x(*(bound2.current_edge), bound1.current_edge->top.y);
+                return greater_than(static_cast<double>(bound1.current_edge->top.x),
+                                    get_current_x(*(bound2.current_edge), bound1.current_edge->top.y));
             }
         } else {
             return bound2.current_x < bound1.current_x;

deps/wagyu/include/mapbox/geometry/wagyu/almost_equal.hpp

@@ -0,0 +1,277 @@
+// Copyright 2005, Google Inc.
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+//     * Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+//     * Redistributions in binary form must reproduce the above
+// copyright notice, this list of conditions and the following disclaimer
+// in the documentation and/or other materials provided with the
+// distribution.
+//     * Neither the name of Google Inc. nor the names of its
+// contributors may be used to endorse or promote products derived from
+// this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Authors: wan@google.com (Zhanyong Wan), eefacm@gmail.com (Sean Mcafee)
+//
+// The Google C++ Testing Framework (Google Test)
+
+// This template class serves as a compile-time function from size to
+// type.  It maps a size in bytes to a primitive type with that
+// size. e.g.
+//
+//   TypeWithSize<4>::UInt
+//
+// is typedef-ed to be unsigned int (unsigned integer made up of 4
+// bytes).
+//
+// Such functionality should belong to STL, but I cannot find it
+// there.
+//
+// Google Test uses this class in the implementation of floating-point
+// comparison.
+//
+// For now it only handles UInt (unsigned int) as that's all Google Test
+// needs.  Other types can be easily added in the future if need
+// arises.
+namespace mapbox {
+namespace geometry {
+namespace wagyu {
+namespace util {
+
+template <size_t size>
+class TypeWithSize {
+public:
+    // This prevents the user from using TypeWithSize<N> with incorrect
+    // values of N.
+    typedef void UInt;
+};
+
+// The specialization for size 4.
+template <>
+class TypeWithSize<4> {
+public:
+    // unsigned int has size 4 in both gcc and MSVC.
+    //
+    // As base/basictypes.h doesn't compile on Windows, we cannot use
+    // uint32, uint64, and etc here.
+    typedef int Int;
+    typedef unsigned int UInt;
+};
+
+// The specialization for size 8.
+template <>
+class TypeWithSize<8> {
+public:
+#if GTEST_OS_WINDOWS
+    typedef __int64 Int;
+    typedef unsigned __int64 UInt;
+#else
+    typedef long long Int;           // NOLINT
+    typedef unsigned long long UInt; // NOLINT
+#endif // GTEST_OS_WINDOWS
+};
+
+// This template class represents an IEEE floating-point number
+// (either single-precision or double-precision, depending on the
+// template parameters).
+//
+// The purpose of this class is to do more sophisticated number
+// comparison.  (Due to round-off error, etc, it's very unlikely that
+// two floating-points will be equal exactly.  Hence a naive
+// comparison by the == operation often doesn't work.)
+//
+// Format of IEEE floating-point:
+//
+//   The most-significant bit being the leftmost, an IEEE
+//   floating-point looks like
+//
+//     sign_bit exponent_bits fraction_bits
+//
+//   Here, sign_bit is a single bit that designates the sign of the
+//   number.
+//
+//   For float, there are 8 exponent bits and 23 fraction bits.
+//
+//   For double, there are 11 exponent bits and 52 fraction bits.
+//
+//   More details can be found at
+//   http://en.wikipedia.org/wiki/IEEE_floating-point_standard.
+//
+// Template parameter:
+//
+//   RawType: the raw floating-point type (either float or double)
+template <typename RawType>
+class FloatingPoint {
+public:
+    // Defines the unsigned integer type that has the same size as the
+    // floating point number.
+    typedef typename TypeWithSize<sizeof(RawType)>::UInt Bits;
+
+    // Constants.
+
+    // # of bits in a number.
+    static const size_t kBitCount = 8 * sizeof(RawType);
+
+    // # of fraction bits in a number.
+    static const size_t kFractionBitCount = std::numeric_limits<RawType>::digits - 1;
+
+    // # of exponent bits in a number.
+    static const size_t kExponentBitCount = kBitCount - 1 - kFractionBitCount;
+
+    // The mask for the sign bit.
+    static const Bits kSignBitMask = static_cast<Bits>(1) << (kBitCount - 1);
+
+    // The mask for the fraction bits.
+    static const Bits kFractionBitMask = ~static_cast<Bits>(0) >> (kExponentBitCount + 1);
+
+    // The mask for the exponent bits.
+    static const Bits kExponentBitMask = ~(kSignBitMask | kFractionBitMask);
+
+    // How many ULP's (Units in the Last Place) we want to tolerate when
+    // comparing two numbers.  The larger the value, the more error we
+    // allow.  A 0 value means that two numbers must be exactly the same
+    // to be considered equal.
+    //
+    // The maximum error of a single floating-point operation is 0.5
+    // units in the last place.  On Intel CPU's, all floating-point
+    // calculations are done with 80-bit precision, while double has 64
+    // bits.  Therefore, 4 should be enough for ordinary use.
+    //
+    // See the following article for more details on ULP:
+    // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm.
+    static const size_t kMaxUlps = 4;
+
+    // Constructs a FloatingPoint from a raw floating-point number.
+    //
+    // On an Intel CPU, passing a non-normalized NAN (Not a Number)
+    // around may change its bits, although the new value is guaranteed
+    // to be also a NAN.  Therefore, don't expect this constructor to
+    // preserve the bits in x when x is a NAN.
+    explicit FloatingPoint(const RawType& x) : u_(x) {
+    }
+
+    // Static methods
+
+    // Reinterprets a bit pattern as a floating-point number.
+    //
+    // This function is needed to test the AlmostEquals() method.
+    static RawType ReinterpretBits(const Bits bits) {
+        FloatingPoint fp(0);
+        fp.u_.bits_ = bits;
+        return fp.u_.value_;
+    }
+
+    // Returns the floating-point number that represent positive infinity.
+    static RawType Infinity() {
+        return ReinterpretBits(kExponentBitMask);
+    }
+
+    // Non-static methods
+
+    // Returns the bits that represents this number.
+    const Bits& bits() const {
+        return u_.bits_;
+    }
+
+    // Returns the exponent bits of this number.
+    Bits exponent_bits() const {
+        return kExponentBitMask & u_.bits_;
+    }
+
+    // Returns the fraction bits of this number.
+    Bits fraction_bits() const {
+        return kFractionBitMask & u_.bits_;
+    }
+
+    // Returns the sign bit of this number.
+    Bits sign_bit() const {
+        return kSignBitMask & u_.bits_;
+    }
+
+    // Returns true iff this is NAN (not a number).
+    bool is_nan() const {
+        // It's a NAN if the exponent bits are all ones and the fraction
+        // bits are not entirely zeros.
+        return (exponent_bits() == kExponentBitMask) && (fraction_bits() != 0);
+    }
+
+    // Returns true iff this number is at most kMaxUlps ULP's away from
+    // rhs.  In particular, this function:
+    //
+    //   - returns false if either number is (or both are) NAN.
+    //   - treats really large numbers as almost equal to infinity.
+    //   - thinks +0.0 and -0.0 are 0 DLP's apart.
+    bool AlmostEquals(const FloatingPoint& rhs) const {
+        // The IEEE standard says that any comparison operation involving
+        // a NAN must return false.
+        if (is_nan() || rhs.is_nan())
+            return false;
+
+        return DistanceBetweenSignAndMagnitudeNumbers(u_.bits_, rhs.u_.bits_) <= kMaxUlps;
+    }
+
+private:
+    // The data type used to store the actual floating-point number.
+    union FloatingPointUnion {
+        explicit FloatingPointUnion(RawType val) : value_(val) {
+        }
+        RawType value_; // The raw floating-point number.
+        Bits bits_;     // The bits that represent the number.
+    };
+
+    // Converts an integer from the sign-and-magnitude representation to
+    // the biased representation.  More precisely, let N be 2 to the
+    // power of (kBitCount - 1), an integer x is represented by the
+    // unsigned number x + N.
+    //
+    // For instance,
+    //
+    //   -N + 1 (the most negative number representable using
+    //          sign-and-magnitude) is represented by 1;
+    //   0      is represented by N; and
+    //   N - 1  (the biggest number representable using
+    //          sign-and-magnitude) is represented by 2N - 1.
+    //
+    // Read http://en.wikipedia.org/wiki/Signed_number_representations
+    // for more details on signed number representations.
+    static Bits SignAndMagnitudeToBiased(const Bits& sam) {
+        if (kSignBitMask & sam) {
+            // sam represents a negative number.
+            return ~sam + 1;
+        } else {
+            // sam represents a positive number.
+            return kSignBitMask | sam;
+        }
+    }
+
+    // Given two numbers in the sign-and-magnitude representation,
+    // returns the distance between them as an unsigned number.
+    static Bits DistanceBetweenSignAndMagnitudeNumbers(const Bits& sam1, const Bits& sam2) {
+        const Bits biased1 = SignAndMagnitudeToBiased(sam1);
+        const Bits biased2 = SignAndMagnitudeToBiased(sam2);
+        return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1);
+    }
+
+    FloatingPointUnion u_;
+};
+
+} // namespace util
+} // namespace wagyu
+} // namespace geometry
+} // namespace mapbox

deps/wagyu/include/mapbox/geometry/wagyu/interrupt.hpp

@@ -9,7 +9,7 @@
 #ifdef USE_WAGYU_INTERRUPT
 
 namespace {
-bool WAGYU_INTERRUPT_REQUESTED = false;
+thread_local bool WAGYU_INTERRUPT_REQUESTED = false;
 }
 
 namespace mapbox {

deps/wagyu/include/mapbox/geometry/wagyu/process_horizontal.hpp

@@ -54,7 +54,7 @@ active_bound_list_itr<T> process_horizontal_left_to_right(T scanline_y,
             ++hp_itr;
         }
 
-        if ((*bnd)->current_x > static_cast<double>((*horz_bound)->current_edge->top.x)) {
+        if (greater_than((*bnd)->current_x, static_cast<double>((*horz_bound)->current_edge->top.x))) {
             break;
         }
 
@@ -75,7 +75,7 @@ active_bound_list_itr<T> process_horizontal_left_to_right(T scanline_y,
         // 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) {
+            if ((*horz_bound)->ring && (*bound_max_pair)->ring) {
                 add_local_maximum_point(*(*horz_bound), *(*bound_max_pair), (*horz_bound)->current_edge->top, rings,
                                         active_bounds);
             }
@@ -159,7 +159,7 @@ active_bound_list_itr<T> process_horizontal_right_to_left(T scanline_y,
             ++hp_itr;
         }
 
-        if ((*bnd)->current_x < static_cast<double>((*horz_bound)->current_edge->top.x)) {
+        if (less_than((*bnd)->current_x, static_cast<double>((*horz_bound)->current_edge->top.x))) {
             break;
         }
 
@@ -180,7 +180,7 @@ active_bound_list_itr<T> process_horizontal_right_to_left(T scanline_y,
         // 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) {
+            if ((*horz_bound)->ring && (*bound_max_pair)->ring) {
                 add_local_maximum_point(*(*horz_bound), *(*bound_max_pair), (*horz_bound)->current_edge->top, rings,
                                         active_bounds);
             }

deps/wagyu/include/mapbox/geometry/wagyu/ring_util.hpp

@@ -74,18 +74,17 @@ struct hot_pixel_sorter {
     }
 };
 
-// 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));
+    double half = std::floor(val) + 0.5;
+    if (values_are_equal(val, half)) {
+        return static_cast<T>(std::floor(val));
+    } else {
+        return static_cast<T>(std::llround(val));
+    }
 }
 
 template <typename T>
@@ -93,7 +92,12 @@ 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));
+    double half = std::floor(val) + 0.5;
+    if (values_are_equal(val, half)) {
+        return static_cast<T>(std::ceil(val));
+    } else {
+        return static_cast<T>(std::llround(val));
+    }
 }
 
 template <typename T>
@@ -117,8 +121,8 @@ inline T get_edge_min_x(edge<T> const& edge, const T current_y) {
         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);
+            double return_val =
+                static_cast<double>(edge.bot.x) + edge.dx * (static_cast<double>(current_y - edge.bot.y) + 0.5);
             T value = round_towards_min<T>(return_val);
             return value;
         }
@@ -146,8 +150,8 @@ inline T get_edge_max_x(edge<T> const& edge, const T current_y) {
         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);
+            double return_val =
+                static_cast<double>(edge.bot.x) + edge.dx * (static_cast<double>(current_y - edge.bot.y) + 0.5);
             T value = round_towards_max<T>(return_val);
             return value;
         }

deps/wagyu/include/mapbox/geometry/wagyu/util.hpp

@@ -4,6 +4,7 @@
 
 #include <mapbox/geometry/point.hpp>
 #include <mapbox/geometry/polygon.hpp>
+#include <mapbox/geometry/wagyu/almost_equal.hpp>
 #include <mapbox/geometry/wagyu/point.hpp>
 
 namespace mapbox {
@@ -30,18 +31,26 @@ double area(mapbox::geometry::linear_ring<T> const& poly) {
     return -a * 0.5;
 }
 
-inline bool value_is_zero(double val) {
-    return std::fabs(val) < (5.0 * std::numeric_limits<double>::epsilon());
+inline bool values_are_equal(double x, double y) {
+    return util::FloatingPoint<double>(x).AlmostEquals(util::FloatingPoint<double>(y));
 }
 
-inline bool values_are_equal(double x, double y) {
-    return value_is_zero(x - y);
+inline bool value_is_zero(double val) {
+    return values_are_equal(val, static_cast<double>(0.0));
 }
 
 inline bool greater_than_or_equal(double x, double y) {
     return x > y || values_are_equal(x, y);
 }
 
+inline bool greater_than(double x, double y) {
+    return (!values_are_equal(x, y) && x > y);
+}
+
+inline bool less_than(double x, double y) {
+    return (!values_are_equal(x, y) && x < y);
+}
+
 template <typename T>
 bool slopes_equal(mapbox::geometry::point<T> const& pt1,
                   mapbox::geometry::point<T> const& pt2,

deps/wagyu/include/mapbox/geometry/wagyu/wagyu.hpp

@@ -17,8 +17,8 @@
 #include <mapbox/geometry/wagyu/vatti.hpp>
 
 #define WAGYU_MAJOR_VERSION 0
-#define WAGYU_MINOR_VERSION 4
-#define WAGYU_PATCH_VERSION 3
+#define WAGYU_MINOR_VERSION 5
+#define WAGYU_PATCH_VERSION 0
 
 #define WAGYU_VERSION (WAGYU_MAJOR_VERSION * 100000) + (WAGYU_MINOR_VERSION * 100) + (WAGYU_PATCH_VERSION)