Test for overlay operations on random grids with integer coordinates.

test/robustness/overlay/areal_areal/Jamfile

@@ -25,4 +25,4 @@ exe recursive_polygons : recursive_polygons.cpp ;
 exe star_comb : star_comb.cpp ;
 
 exe ticket_9081 : ticket_9081.cpp ;
-
+exe random_integer_grids : random_integer_grids.cpp ;

test/robustness/overlay/areal_areal/random_integer_grids.cpp

@@ -0,0 +1,274 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Robustness Test
+
+// Copyright (c) 2025 Tinko Bartels, Shenzhen, China.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#define BOOST_GEOMETRY_NO_BOOST_TEST
+
+#include <bitset>
+#include <iostream>
+#include <random>
+#include <vector>
+
+#include <boost/program_options.hpp>
+
+#include <boost/geometry/algorithms/difference.hpp>
+#include <boost/geometry/algorithms/union.hpp>
+#include <boost/geometry/algorithms/sym_difference.hpp>
+#include <boost/geometry/algorithms/intersection.hpp>
+#include <boost/geometry/algorithms/is_valid.hpp>
+#include <boost/geometry/io/wkt/write.hpp>
+#include <boost/geometry/geometries/geometries.hpp>
+
+#define BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
+#define BOOST_GEOMETRY_DEFAULT_TEST_TYPE int
+#include <geometry_test_common.hpp>
+
+constexpr int chunk_size = 64;
+using bits = std::vector<std::bitset<chunk_size>>;
+using prng = std::mt19937_64;
+static_assert(sizeof(prng::result_type) >= chunk_size / 8, "PRNG output too narrow.");
+
+namespace bg = boost::geometry;
+
+using point = bg::model::d2::point_xy<int>;
+using box   = bg::model::box<point>;
+using mp    = bg::model::multi_polygon<bg::model::polygon<point>>;
+
+struct grid_settings
+{
+    int width              = 5;
+    int height             = 5;
+    int count              = 1;
+    prng::result_type seed = prng::default_seed;
+};
+
+constexpr int cell_width = 2;
+
+std::vector<box> grid_cells(grid_settings const& settings)
+{
+    std::vector<box> out;
+    for (int y = settings.height - 1; y >= 0; --y)
+    {
+        for (int x = 0; x < settings.width; ++x)
+        {
+            out.push_back(box{point{x * cell_width, y * cell_width},
+                              point{x * cell_width + cell_width, y * cell_width + cell_width}});
+        }
+    }
+    return out;
+}
+
+std::vector<point> test_points(grid_settings const& settings)
+{
+    std::vector<point> out;
+    for (int y = settings.height - 1; y >= 0; --y)
+    {
+        for (int x = 0; x < settings.width; ++x)
+        {
+            out.push_back(point{x * cell_width + cell_width / 2, y * cell_width + cell_width / 2});
+        }
+    }
+    return out;
+}
+
+std::pair<int, int> xalloc_index()
+{
+    static int width_index = std::ios_base::xalloc();
+    static int height_index = std::ios_base::xalloc();
+    return { width_index, height_index };
+} 
+
+std::ostream& operator<<(std::ostream& os, bits const& b)
+{
+    int width = os.iword(xalloc_index().first);
+    int height = os.iword(xalloc_index().second);
+    for (int y = 0; y < height; ++y)
+    {
+        for (int x = 0; x < width; ++x)
+        {
+            int index = y * width + x;
+            os << b[index / chunk_size][index % chunk_size];
+        }
+        if (y != height - 1) os << '\n';
+    }
+    return os;
+}
+
+bits geometry_to_bits(mp const& g, std::vector<point> const& tp)
+{
+    bits out(tp.size() / chunk_size + (tp.size() % chunk_size == 0 ? 0 : 1));
+    for (size_t i = 0; i < tp.size(); ++i)
+    {
+        out[i / chunk_size][i % chunk_size] = bg::within(tp[i], g);
+    }
+    return out;
+}
+
+mp bits_to_geometry(bits const& b, std::vector<box> const& grid, std::vector<point> const& points,
+                    grid_settings const& settings, bool& all_success)
+{
+    mp out({{{point{0, 0}, point{0, cell_width * settings.height},
+              point{cell_width * settings.width, cell_width * settings.height},
+              point{cell_width * settings.width, 0},
+              point{0, 0}}}});
+    for (size_t i = 0; i < grid.size(); ++i)
+    {
+        if( ! b[i / chunk_size][i % chunk_size] )
+        {
+            mp temp;
+            bg::difference(out, grid[i], temp);
+            out = temp;
+        }
+    }
+    if ( ! bg::is_valid(out) )
+    {
+        std::cout << "Generating grid from pattern\n\n" << b
+                  << "\n\nresults in invalid geometry: " << bg::wkt(out) << ".\n\n";
+        all_success = false;
+    }
+    if (geometry_to_bits(out, points) != b)
+    {
+        std::cout << "Generating grid from pattern\n\n" << b
+                  << "\n\nresults in mismatching geometry: "
+                  << bg::wkt(out) << ".\n\n";
+        all_success = false;
+    }
+    return out;
+}
+
+bits gen_bits(prng& gen, int bits_size)
+{
+    bits b(bits_size / chunk_size + (bits_size % chunk_size == 0 ? 0 : 1));
+    std::generate(b.begin(), b.end(), std::ref(gen));
+    if (bits_size % chunk_size != 0)
+    {
+        std::bitset<chunk_size> bm;
+        bm.set();
+        bm >>= chunk_size - bits_size % chunk_size;
+        b.back() &= bm;
+    }
+    return b;
+}
+
+template <typename BitOp>
+bits apply_for_each(bits a, bits const& b, BitOp bit_op)
+{
+    for (size_t i = 0; i < a.size(); ++i) a[i] = bit_op(a[i], b[i]);
+    return a;
+}
+
+template<typename BitOp, typename GeoOp>
+void test_op(bits const& bits1, bits const& bits2, mp const& geo1, mp const& geo2,
+             std::string op_label, BitOp bit_op, GeoOp geo_op, std::vector<point> const& points,
+             std::vector<box> const& grid, grid_settings const& settings, bool& success)
+{
+    auto test_geo = geo_op(geo1, geo2);
+    if ( ! bg::is_valid(test_geo) )
+    {
+        std::cout << op_label << "(\n " << bg::wkt(geo1) << ",\n " << bg::wkt(geo2)
+                  << "\n), generated from\n\n" << bits1 << "\n\nand\n\n" << bits2 << "\n\nis invalid.\n\n";
+        success = false;
+    }
+    bits expected = apply_for_each(bits1, bits2, bit_op);
+    bits obtained = geometry_to_bits(test_geo, points);
+    if (obtained != expected)
+    {
+        std::cout << op_label << "(\n " << bg::wkt(geo1) << ",\n " << bg::wkt(geo2)
+                  << "\n), generated from\n\n" << bits1 << "\n\nand\n\n" << bits2 << "\n\nis incorrect.\n"
+                  << "\nExpected: "
+                  << bg::wkt(bits_to_geometry(expected, grid, points, settings, success))
+                  << "\n\n" << expected << "\n\nObtained: " << bg::wkt(test_geo) << "\n\n" << obtained
+                  << "\n\n";
+        success = false;
+    }
+}
+
+bool test_all(grid_settings const& settings)
+{
+    prng gen(settings.seed);
+    const auto grid = grid_cells(settings);
+    const auto points = test_points(settings);
+    bool all_success = true;
+    for (int i = 0; i < settings.count || settings.count == -1; i++)
+    {
+        const bits bits1 = gen_bits(gen, settings.width * settings.height);
+        const bits bits2 = gen_bits(gen, settings.width * settings.height);
+        const auto geo1 = bits_to_geometry(bits1, grid, points, settings, all_success);
+        const auto geo2 = bits_to_geometry(bits2, grid, points, settings, all_success);
+        test_op(bits1, bits2, geo1, geo2, "union", std::bit_or<>{},
+                [](mp const& g1, mp const& g2) { mp g; bg::union_(g1, g2, g); return g; },
+                points, grid, settings, all_success);
+        test_op(bits1, bits2, geo1, geo2, "intersection", std::bit_and<>{},
+                [](mp const& g1, mp const& g2) { mp g; bg::intersection(g1, g2, g); return g; },
+                points, grid, settings, all_success);
+        test_op(bits1, bits2, geo1, geo2, "sym_difference", std::bit_xor<>{},
+                [](mp const& g1, mp const& g2) { mp g; bg::sym_difference(g1, g2, g); return g; },
+                points, grid, settings, all_success);
+        test_op(bits1, bits2, geo1, geo2, "difference g1 \\ g2", 
+                [](std::bitset<chunk_size> b1, std::bitset<chunk_size> b2) { return b1 & (~b2); },
+                [](mp const& g1, mp const& g2) { mp g; bg::difference(g1, g2, g); return g; },
+                points, grid, settings, all_success);
+        test_op(bits1, bits2, geo1, geo2, "difference g2 \\ g1", 
+                [](std::bitset<chunk_size> b1, std::bitset<chunk_size> b2) { return b2 & (~b1); },
+                [](mp const& g1, mp const& g2) { mp g; bg::difference(g2, g1, g); return g; },
+                points, grid, settings, all_success);
+    }
+    return all_success;
+}
+
+int main(int argc, char** argv)
+{
+    BoostGeometryWriteTestConfiguration();
+    try
+    {
+        namespace po = boost::program_options;
+        po::options_description description("=== random_integer_grids ===\nAllowed options");
+
+        grid_settings settings;
+
+        description.add_options()
+            ("help", "Help message")
+            ("seed",
+             po::value<decltype(settings.seed)>(&settings.seed)->default_value(settings.seed),
+             "Initialization seed for random generator")
+            ("count",
+             po::value<decltype(settings.count)>(&settings.count)->default_value(settings.count),
+             "Number of tests (-1 for infinite loop)")
+            ("width",
+             po::value<decltype(settings.width)>(&settings.width)->default_value(settings.width),
+             "Width of grid")
+            ("height",
+             po::value<decltype(settings.height)>(&settings.height)->default_value(settings.height),
+             "Height of grid")
+        ;
+
+        po::variables_map varmap;
+        po::store(po::parse_command_line(argc, argv, description), varmap);
+        po::notify(varmap);
+
+        if (varmap.count("help"))
+        {
+            std::cout << description << std::endl;
+            return 1;
+        }
+        std::cout.iword(xalloc_index().first) = settings.width;
+        std::cout.iword(xalloc_index().second) = settings.height;
+        if ( ! test_all(settings) )
+            return 1;
+    }
+    catch(std::exception const& e)
+    {
+        std::cout << "Exception " << e.what() << '\n';
+    }
+    catch(...)
+    {
+        std::cout << "Other exception" << '\n';
+    }
+
+    return 0;
+}