Quick test of bg set ops vs corresponding logic operations

test_bg_set_ops.cpp

#include <iostream>
#include <bitset>
#include <fstream>
#include <random>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>

namespace bg = boost::geometry;

using point = bg::model::point<int, 2, bg::cs::cartesian>;
using polygon = bg::model::polygon<point>;
using multipoly = bg::model::multi_polygon<polygon>;


template <int width, int height, bool verbose = false>
multipoly bits_to_geometry(std::bitset<width * height> const& bits)
{
    const polygon full({{point{0, 0}, point{0, 2 * height}, point{2 * width, 2 * height}, point{2 * width, 0}, point{0, 0}}});
    if( !bg::is_valid(full) ) std::cout << "Full invalid.\n";
    multipoly out({full});
    for (int i = 0; i < width; ++i)
        for (int j = 0; j < height; ++j)
        {
            if( ! bits[j * width + i] )
            {
                multipoly temp;
                polygon hole({{{2 * i, 2 * j}, {2 * i, 2 * j + 2}, {2 * i + 2, 2* j + 2}, {2 * i + 2, 2 * j}, {2 * i, 2 * j}}});
                bg::difference(out, hole, temp);
                if constexpr(verbose)
                    if( !bg::is_valid(temp) )
                    {
                        std::cout << "out: " << bg::wkt(out) << "\n";
                        std::cout << "hole: " << bg::wkt(hole) << "\n";
                        std::cout << "out \\ hole: " << bg::wkt(temp) << "\n";
                    }
                out = temp;
            }
        }
    return out;
}

template <int width, int height>
std::bitset<width * height> geometry_to_bits(multipoly const& geo)
{
    std::bitset<width * height> out;
    for (int i = 0; i < width; ++i)
        for (int j = 0; j < height; ++j)
        {
            point test({2 * i + 1, 2 * j + 1});
            out[j * width + i] = bg::within(test, geo);
        }
    return out;
}

template <int width, int height, bool svg = false>
void print_failure(auto bits1, auto bits2, auto resbits, auto geo1, auto geo2, auto resgeo, auto op)
{
    if constexpr(svg)
    {
        std::ofstream svg1("my_map1.svg");
        std::ofstream svg2("my_map2.svg");
        std::ofstream svg3("my_map3.svg");
        bg::svg_mapper<point> mapper1(svg1, 400, 400);
        bg::svg_mapper<point> mapper2(svg2, 400, 400);
        bg::svg_mapper<point> mapper3(svg3, 400, 400);
        mapper1.add(point{0,0});
        mapper1.add(point{2 * width, 2 * height});
        mapper2.add(point{0,0});
        mapper2.add(point{2 * width, 2 * height});
        mapper3.add(point{0,0});
        mapper3.add(point{2 * width, 2 * height});
        mapper1.map(geo1, "fill: red");
        mapper2.map(geo2, "fill: red");
        mapper3.map(resgeo, "fill: red");
    }
    std::cout << "geo1: " << bg::wkt(geo1) << "\n";
    std::cout << "geo2: " << bg::wkt(geo2) << "\n";
    std::cout << op << ": " << bg::wkt(resgeo) << "\n";
    std::cout << "expected: " << bg::wkt(bits_to_geometry<width, height>(resbits)) << "\n";
    std::cout << bits1 << "\n";
    std::cout << bits2 << "\n";
    std::cout << resbits << "\n";
}

template <int width, int height, bool validity = true>
bool test(unsigned long long val)
{
    bool correct = true;
    auto first = val;
    auto second = val >> (width * height);

    std::bitset<width * height> bits1(first), bits2(second);
    auto geo1 = bits_to_geometry<width, height/*, true*/>(bits1);
    auto geo2 = bits_to_geometry<width, height/*, true*/>(bits2);
    if constexpr (validity)
    {
        correct &= bg::is_valid(geo1);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits1, geo1, geo2, geo1, "geo1 invalid"); return false; }
        correct &= bg::is_valid(geo2);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits2, geo1, geo2, geo2, "geo2 invalid"); return false; }
    }

    multipoly testgeo;
    bg::union_(geo1, geo2, testgeo);
    auto geo_union = geometry_to_bits<width, height>(testgeo);
    auto bits_or = bits1 | bits2;
    correct &= (bits_or == geo_union);
    if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_or, geo1, geo2, testgeo, "geo1 union_ geo2 wrong"); return false; }
    if constexpr (validity)
    {
        correct &= bg::is_valid(testgeo);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_or, geo1, geo2, testgeo, "geo1 union_ geo2 invalid"); return false; }
    }

    bg::clear(testgeo);
    bg::intersection(geo1, geo2, testgeo);
    auto geo_intersection = geometry_to_bits<width, height>(testgeo);
    auto bits_and = bits1 & bits2;
    correct &= (bits_and == geo_intersection);
    if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_and, geo1, geo2, testgeo, "geo1 intersection geo2 wrong"); return false; }
    if constexpr (validity)
    {
        correct &= bg::is_valid(testgeo);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_and, geo1, geo2, testgeo, "geo1 intersection geo2 invalid"); return false; }
    }

    bg::clear(testgeo);
    bg::sym_difference(geo1, geo2, testgeo);
    auto geo_symdiff = geometry_to_bits<width, height>(testgeo);
    auto bits_xor = bits1 ^ bits2;
    correct &= (bits_xor == geo_symdiff);
    if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_xor, geo1, geo2, testgeo, "geo1 sym_difference geo2 wrong"); return false; }
    if constexpr (validity)
    {
        correct &= bg::is_valid(testgeo);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_xor, geo1, geo2, testgeo, "geo1 sym_difference geo2 invalid"); return false; }
    }

    bg::clear(testgeo);
    bg::difference(geo1, geo2, testgeo);
    auto geo_d12 = geometry_to_bits<width, height>(testgeo);
    auto bits_d12 = bits1 & (~bits2);
    correct &= (bits_d12 == geo_d12);
    if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_d12, geo1, geo2, testgeo, "geo1 \\ geo2 wrong"); return false; }
    if constexpr (validity)
    {
        correct &= bg::is_valid(testgeo);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_d12, geo1, geo2, testgeo, "geo1 \\ geo2 invalid"); return false; }
    }

    bg::clear(testgeo);
    bg::difference(geo2, geo1, testgeo);
    auto geo_d21 = geometry_to_bits<width, height>(testgeo);
    auto bits_d21 = bits2 & (~bits1);
    correct &= (bits_d21 == geo_d21);
    if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_d21, geo1, geo2, testgeo, "geo2 \\ geo1 wrong"); return false; }
    if constexpr (validity)
    {
        correct &= bg::is_valid(testgeo);
        if ( !correct ) { print_failure<width, height>(bits1, bits2, bits_d21, geo1, geo2, testgeo, "geo2 \\ geo1 invalid"); return false; }
    }

    return correct;
}

template <int dim1, int dim2> auto get_prng() {
    if constexpr (dim1 * dim2 * 2 <= 32) return std::mt19937{}; 
    else return std::mt19937_64{};
}

int main() {
    constexpr int width = 3, height = 5;

    auto prng = get_prng<width, height>();
    for ( /*unsigned long long i = 0ULL; i < (1ULL << (width * height * 2); ++i*/ ;;)
    {
        auto correct = test<width, height, false>( /*i*/ prng());
    }
    return 0;
}