#include <cstddef>

#include <iostream>

#include <boost/geometry.hpp>

namespace bg = boost::geometry;

typedef bg::model::point<double, 2, bg::cs::cartesian> point;
typedef bg::model::polygon<point> polygon; // cw, closed
typedef bg::model::multi_polygon<polygon> multi_polygon;


// Display some info about a multi_polygon
void check(const multi_polygon& mp)
{
    std::cout << bg::wkt(mp) << std::endl;
    std::cout << "  Size: " << mp.size() << std::endl;
    for(std::size_t i = 0; i < mp.size(); ++i)
    {
        std::cout << "    Num points: "
                  << boost::geometry::num_points(mp[i])
                  << std::endl;
    }
    std::cout << "  Intersects: "
              << boost::geometry::intersects(mp)
              << std::endl;
}


int main()
{
   polygon green, blue;

   // Using polygons from example - details are unimportant, any overlapping
   //polygons would do
   boost::geometry::read_wkt(
      "POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2,3.7 1.6,3.4 2,4.1 3,5.3 2.6,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3),(4.0 2.0, 4.2 1.4, 4.8 1.9, 4.4 2.2, 4.0 2.0))", green);

   boost::geometry::read_wkt(
      "POLYGON((4.0 -0.5 , 3.5 1.0 , 2.0 1.5 , 3.5 2.0 , 4.0 3.5 , 4.5 2.0 ,6.0 1.5 , 4.5 1.0 , 4.0 -0.5))", blue);

   // Try 1 [WORKS]: Take union of two polygons and put result into
   //multi_polygon
   multi_polygon output1;
   boost::geometry::union_(green, blue, output1);
   std::cout << "output1 = union(green,blue)" << std::endl;
   check(output1);

   std::cout << std::endl;

   // Try 2 [WORKS]: Build up multi_polygon by repeatedly unioning in
   // each polygon
   multi_polygon output2;
   boost::geometry::union_(output2, green, output2);
   std::cout << "output2 = union(output2,green)" << std::endl;
   check(output2);
    
   std::cout << std::endl;

   multi_polygon output3;
   boost::geometry::union_(output2, blue, output3);
   std::cout << "output3 = union(output2,blue)" << std::endl;
   check(output3);

   return 0;
}