Use a cascaded union for adding many shapes together

Makefile.am

@@ -21,6 +21,8 @@ pcb2gcode_SOURCES = \
     eulerian_paths.hpp \
     eulerian_paths.cpp \
     flatten.hpp \
+    geos_helpers.hpp \
+    geos_helpers.cpp \
     geometry.hpp \
     geometry_int.hpp \
     gerberimporter.hpp \
@@ -79,20 +81,20 @@ check_PROGRAMS = voronoi_tests eulerian_paths_tests segmentize_tests tsp_solver_
                  available_drills_tests gerberimporter_tests options_tests path_finding_tests \
                  autoleveller_tests common_tests backtrack_tests trim_paths_tests outline_bridges_tests
 voronoi_tests_SOURCES = voronoi.hpp voronoi.cpp voronoi_tests.cpp boost_unit_test.cpp
-eulerian_paths_tests_SOURCES = eulerian_paths_tests.cpp eulerian_paths.hpp geometry_int.hpp boost_unit_test.cpp  bg_operators.hpp bg_operators.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.cpp segmentize.cpp merge_near_points.cpp
-segmentize_tests_SOURCES = segmentize_tests.cpp segmentize.cpp segmentize.hpp merge_near_points.cpp merge_near_points.hpp boost_unit_test.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.cpp bg_operators.hpp bg_operators.cpp
-path_finding_tests_SOURCES = path_finding_tests.cpp path_finding.cpp path_finding.hpp boost_unit_test.cpp bg_helpers.cpp bg_helpers.hpp eulerian_paths.cpp eulerian_paths.hpp segmentize.hpp segmentize.cpp merge_near_points.cpp merge_near_points.hpp bg_operators.hpp bg_operators.cpp
+eulerian_paths_tests_SOURCES = eulerian_paths_tests.cpp eulerian_paths.hpp geometry_int.hpp boost_unit_test.cpp  bg_operators.hpp bg_operators.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.cpp segmentize.cpp merge_near_points.cpp geos_helpers.hpp geos_helpers.cpp
+segmentize_tests_SOURCES = segmentize_tests.cpp segmentize.cpp segmentize.hpp merge_near_points.cpp merge_near_points.hpp boost_unit_test.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.cpp bg_operators.hpp bg_operators.cpp geos_helpers.hpp geos_helpers.cpp
+path_finding_tests_SOURCES = path_finding_tests.cpp path_finding.cpp path_finding.hpp boost_unit_test.cpp bg_helpers.cpp bg_helpers.hpp eulerian_paths.cpp eulerian_paths.hpp segmentize.hpp segmentize.cpp merge_near_points.cpp merge_near_points.hpp bg_operators.hpp bg_operators.cpp geos_helpers.hpp geos_helpers.cpp
 tsp_solver_tests_SOURCES = tsp_solver_tests.cpp tsp_solver.hpp boost_unit_test.cpp
 units_tests_SOURCES = units_tests.cpp units.hpp boost_unit_test.cpp
 available_drills_tests_SOURCES = available_drills_tests.cpp available_drills.hpp boost_unit_test.cpp
-gerberimporter_tests_SOURCES = gerberimporter.hpp gerberimporter.cpp gerberimporter_tests.cpp merge_near_points.hpp merge_near_points.cpp eulerian_paths.cpp eulerian_paths.hpp segmentize.cpp segmentize.hpp boost_unit_test.cpp bg_helpers.cpp bg_helpers.hpp bg_operators.hpp bg_operators.cpp
+gerberimporter_tests_SOURCES = gerberimporter.hpp gerberimporter.cpp gerberimporter_tests.cpp merge_near_points.hpp merge_near_points.cpp eulerian_paths.cpp eulerian_paths.hpp segmentize.cpp segmentize.hpp boost_unit_test.cpp bg_helpers.cpp bg_helpers.hpp bg_operators.hpp bg_operators.cpp geos_helpers.hpp geos_helpers.cpp
 gerberimporter_tests_LDFLAGS = $(glibmm_LIBS) $(gdkmm_LIBS) $(rsvg_LIBS)
 options_tests_SOURCES = options_tests.cpp options.hpp options.cpp boost_unit_test.cpp
-autoleveller_tests_SOURCES = autoleveller_tests.cpp autoleveller.hpp autoleveller.cpp options.cpp options.hpp boost_unit_test.cpp bg_operators.hpp bg_operators.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp merge_near_points.hpp merge_near_points.cpp
+autoleveller_tests_SOURCES = autoleveller_tests.cpp autoleveller.hpp autoleveller.cpp options.cpp options.hpp boost_unit_test.cpp bg_operators.hpp bg_operators.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp merge_near_points.hpp merge_near_points.cpp geos_helpers.hpp geos_helpers.cpp
 common_tests_SOURCES = common.hpp common.cpp common_tests.cpp boost_unit_test.cpp
-backtrack_tests_SOURCES = backtrack.hpp backtrack.cpp backtrack_tests.cpp boost_unit_test.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp merge_near_points.hpp merge_near_points.cpp bg_operators.hpp bg_operators.cpp
-trim_paths_tests_SOURCES = trim_paths.hpp trim_paths.cpp trim_paths_tests.cpp boost_unit_test.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp merge_near_points.hpp merge_near_points.cpp bg_operators.hpp bg_operators.cpp
-outline_bridges_tests_SOURCES = outline_bridges_tests.cpp outline_bridges.hpp outline_bridges.cpp bg_operators.hpp bg_operators.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp boost_unit_test.cpp merge_near_points.hpp merge_near_points.cpp
+backtrack_tests_SOURCES = backtrack.hpp backtrack.cpp backtrack_tests.cpp boost_unit_test.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp merge_near_points.hpp merge_near_points.cpp bg_operators.hpp bg_operators.cpp geos_helpers.hpp geos_helpers.cpp
+trim_paths_tests_SOURCES = trim_paths.hpp trim_paths.cpp trim_paths_tests.cpp boost_unit_test.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp merge_near_points.hpp merge_near_points.cpp bg_operators.hpp bg_operators.cpp geos_helpers.hpp geos_helpers.cpp
+outline_bridges_tests_SOURCES = outline_bridges_tests.cpp outline_bridges.hpp outline_bridges.cpp bg_operators.hpp bg_operators.cpp bg_helpers.hpp bg_helpers.cpp eulerian_paths.hpp eulerian_paths.cpp segmentize.hpp segmentize.cpp boost_unit_test.cpp merge_near_points.hpp merge_near_points.cpp geos_helpers.hpp geos_helpers.cpp
 
 TESTS = $(check_PROGRAMS)
 

bg_helpers.cpp

@@ -1,40 +1,12 @@
 #include "eulerian_paths.hpp"
 #ifdef GEOS_VERSION
-#include <geos/io/WKTReader.h>
-#include <geos/io/WKTReader.inl>
-#include <geos/io/WKTWriter.h>
 #include <geos/operation/buffer/BufferOp.h>
+#include "geos_helpers.hpp"
 #endif // GEOS_VERSION
 
 #include "bg_operators.hpp"
 #include "bg_helpers.hpp"
 
-#ifdef GEOS_VERSION
-template <typename T>
-std::unique_ptr<geos::geom::Geometry> to_geos(const T& mp) {
-  geos::io::WKTReader reader;
-  std::stringstream ss;
-  ss << bg::wkt(mp);
-  return reader.read(ss.str());
-}
-
-multi_polygon_type_fp from_geos(const std::unique_ptr<geos::geom::Geometry>& g) {
-  geos::io::WKTWriter writer;
-  std::string geos_wkt = writer.write(g.get());
-  boost::replace_all(geos_wkt, "EMPTY, ", "");
-  boost::replace_all(geos_wkt, ", EMPTY", "");
-  multi_polygon_type_fp ret;
-  if (strncmp(geos_wkt.c_str(), "MULTIPOLYGON", 12) == 0) {
-    bg::read_wkt(geos_wkt, ret);
-  } else {
-    polygon_type_fp poly;
-    bg::read_wkt(geos_wkt, poly);
-    ret.push_back(poly);
-  }
-  return ret;
-}
-#endif //GEOS_VERSION
-
 namespace bg_helpers {
 
 // The below implementations of buffer are similar to bg::buffer but
@@ -49,7 +21,7 @@ multi_polygon_type_fp buffer(multi_polygon_type_fp const & geometry_in, coordina
   }
 #ifdef GEOS_VERSION
   auto geos_in = to_geos(geometry_in);
-  return from_geos(
+  return from_geos<multi_polygon_type_fp>(
       std::unique_ptr<geos::geom::Geometry>(
           geos::operation::buffer::BufferOp::bufferOp(geos_in.get(), expand_by, points_per_circle/4)));
 #else
@@ -98,7 +70,7 @@ multi_polygon_type_fp buffer(linestring_type_fp const & geometry_in, CoordinateT
   }
 #ifdef GEOS_VERSION
   auto geos_in = to_geos(geometry_in);
-  return from_geos(
+  return from_geos<multi_polygon_type_fp>(
       std::unique_ptr<geos::geom::Geometry>(
           geos::operation::buffer::BufferOp::bufferOp(geos_in.get(), expand_by, points_per_circle/4)));
 #else
@@ -123,7 +95,7 @@ multi_polygon_type_fp buffer(multi_linestring_type_fp const & geometry_in, Coord
   multi_linestring_type_fp mls = eulerian_paths::make_eulerian_paths(geometry_in, true, true);
 #ifdef GEOS_VERSION
   auto geos_in = to_geos(mls);
-  return from_geos(
+  return from_geos<multi_polygon_type_fp>(
       std::unique_ptr<geos::geom::Geometry>(
           geos::operation::buffer::BufferOp::bufferOp(geos_in.get(), expand_by, points_per_circle/4)));
 #else

bg_helpers.hpp

@@ -3,10 +3,6 @@
 
 #include "eulerian_paths.hpp"
 #include <boost/functional/hash/hash.hpp>
-#ifdef GEOS_VERSION
-#include <geos/io/WKTReader.h>
-#include <geos/io/WKTWriter.h>
-#endif // GEOS_VERSION
 
 namespace bg_helpers {
 

bg_operators.cpp

@@ -1,9 +1,17 @@
 #include "geometry.hpp"
 #include "geometry_int.hpp"
 #include "bg_helpers.hpp"
+#ifdef GEOS_VERSION
+#include <geos/operation/union/CascadedUnion.h>
+#include <geos/util/TopologyException.h>
+#include "geos_helpers.hpp"
+#endif // GEOS_VERSION
 
 #include "bg_operators.hpp"
 
+using std::unique_ptr;
+using std::vector;
+
 template <typename polygon_type_t, typename rhs_t>
 bg::model::multi_polygon<polygon_type_t> operator-(
     const bg::model::multi_polygon<polygon_type_t>& lhs,
@@ -160,3 +168,92 @@ bg::model::multi_polygon<polygon_type_t> operator+(const bg::model::multi_polygo
 }
 
 template multi_polygon_type_fp operator+(const multi_polygon_type_fp&, const multi_polygon_type_fp&);
+
+template <typename Addition>
+multi_polygon_type_fp reduce(const std::vector<multi_polygon_type_fp>& mpolys,
+                          const Addition& adder,
+                          const std::vector<box_type_fp>& bboxes) {
+  if (mpolys.size() == 0) {
+    return multi_polygon_type_fp();
+  } else if (mpolys.size() == 1) {
+    return mpolys.front();
+  }
+  size_t current = 0;
+  std::vector<multi_polygon_type_fp> new_mpolys;
+  std::vector<box_type_fp> new_bboxes;
+  if (mpolys.size() % 2 == 1) {
+    new_mpolys.push_back(mpolys[current]);
+    new_bboxes.push_back(bboxes[current]);
+    current++;
+  }
+  // There are at least two and the total number is even.
+  for (; current < mpolys.size(); current += 2) {
+    box_type_fp new_bbox = bboxes[current];
+    bg::expand(new_bbox, bboxes[current+1]);
+    new_bboxes.push_back(new_bbox);
+    if (!bg::intersects(bboxes[current], bboxes[current+1])) {
+      new_mpolys.push_back(mpolys[current]);
+      new_mpolys.back().insert(new_mpolys.back().cend(), mpolys[current+1].cbegin(), mpolys[current+1].cend());
+    } else {
+      new_mpolys.push_back(adder(mpolys[current], mpolys[current+1]));
+    }
+  }
+  return reduce(new_mpolys, adder, new_bboxes);
+}
+
+template <typename Addition>
+multi_polygon_type_fp reduce(const std::vector<multi_polygon_type_fp>& mpolys,
+                             const Addition& adder) {
+  std::vector<box_type_fp> bboxes;
+  bboxes.reserve(mpolys.size());
+  for (const auto& mpoly : mpolys) {
+    bboxes.push_back(bg::return_envelope<box_type_fp>(mpoly));
+  }
+  return reduce(mpolys, adder, bboxes);
+}
+
+multi_polygon_type_fp sum(const std::vector<multi_polygon_type_fp>& mpolys) {
+  if (mpolys.size() == 0) {
+    return {};
+  } else if (mpolys.size() == 1) {
+    return mpolys[0];
+  }
+#ifdef GEOS_VERSION
+  std::vector<geos::geom::Geometry*> geos_mpolys;
+  std::vector<std::unique_ptr<geos::geom::Geometry>> geos_mpolys_tmp;
+  for (const auto& mpoly : mpolys) {
+    if (bg::area(mpoly) == 0) {
+      continue;
+    }
+    geos_mpolys_tmp.push_back(to_geos(mpoly));
+    geos_mpolys.push_back(geos_mpolys_tmp.back().get());
+  }
+  if (geos_mpolys.size() == 0) {
+    return {};
+  }
+  try {
+    std::unique_ptr<geos::geom::Geometry> geos_out(
+        geos::operation::geounion::CascadedUnion::Union(&geos_mpolys));
+    return from_geos<multi_polygon_type_fp>(geos_out);
+  } catch (const geos::util::TopologyException& e) {
+    std::cerr << "\nError: Internal error with libgeos.  Upgrading geos may help." << std::endl;
+    throw;
+  }
+#else // !GEOS_VERSION
+  return reduce(mpolys, operator+<polygon_type_fp, multi_polygon_type_fp>);
+#endif // GEOS_VERSION
+}
+
+multi_polygon_type_fp symdiff(const std::vector<multi_polygon_type_fp>& mpolys) {
+  if (mpolys.size() == 0) {
+    return multi_polygon_type_fp();
+  } else if (mpolys.size() == 1) {
+    return mpolys[0];
+  }
+  std::vector<box_type_fp> bboxes;
+  bboxes.reserve(mpolys.size());
+  for (const auto& mpoly : mpolys) {
+    bboxes.push_back(bg::return_envelope<box_type_fp>(mpoly));
+  }
+  return reduce(mpolys, operator^<polygon_type_fp>);
+}

bg_operators.hpp

@@ -47,6 +47,9 @@ template <typename polygon_type_t, typename rhs_t>
 bg::model::multi_polygon<polygon_type_t> operator+(const bg::model::multi_polygon<polygon_type_t>& lhs,
                                                    const rhs_t& rhs);
 
+multi_polygon_type_fp sum(const std::vector<multi_polygon_type_fp>& mpolys);
+multi_polygon_type_fp symdiff(const std::vector<multi_polygon_type_fp>& mpolys);
+
 // It's not great to insert definitions into the bg namespace but they
 // are useful for sorting and maps.
 

geos_helpers.cpp

@@ -0,0 +1,57 @@
+#ifdef GEOS_VERSION
+
+#include "geometry.hpp"
+#include <geos/io/WKTReader.h>
+#include <geos/io/WKTReader.inl>
+#include <geos/io/WKTWriter.h>
+
+template <typename T>
+std::unique_ptr<geos::geom::Geometry> to_geos(const T& mp) {
+  geos::io::WKTReader reader;
+  std::stringstream ss;
+  ss << bg::wkt(mp);
+  return reader.read(ss.str());
+}
+template std::unique_ptr<geos::geom::Geometry> to_geos(const multi_polygon_type_fp& mp);
+template std::unique_ptr<geos::geom::Geometry> to_geos(const multi_linestring_type_fp& mls);
+template std::unique_ptr<geos::geom::Geometry> to_geos(const linestring_type_fp& mls);
+
+
+template <typename T>
+T from_geos(const std::unique_ptr<geos::geom::Geometry>& g);
+
+template <>
+multi_polygon_type_fp from_geos(const std::unique_ptr<geos::geom::Geometry>& g) {
+  geos::io::WKTWriter writer;
+  std::string geos_wkt = writer.write(g.get());
+  boost::replace_all(geos_wkt, "EMPTY, ", "");
+  boost::replace_all(geos_wkt, ", EMPTY", "");
+  multi_polygon_type_fp ret;
+  if (strncmp(geos_wkt.c_str(), "MULTIPOLYGON", 12) == 0) {
+    bg::read_wkt(geos_wkt, ret);
+  } else {
+    polygon_type_fp poly;
+    bg::read_wkt(geos_wkt, poly);
+    ret.push_back(poly);
+  }
+  return ret;
+}
+
+template <>
+multi_linestring_type_fp from_geos(const std::unique_ptr<geos::geom::Geometry>& g) {
+  geos::io::WKTWriter writer;
+  std::string geos_wkt = writer.write(g.get());
+  boost::replace_all(geos_wkt, "EMPTY, ", "");
+  boost::replace_all(geos_wkt, ", EMPTY", "");
+  multi_linestring_type_fp ret;
+  if (strncmp(geos_wkt.c_str(), "MULTILINESTRING", 12) == 0) {
+    bg::read_wkt(geos_wkt, ret);
+  } else {
+    linestring_type_fp ls;
+    bg::read_wkt(geos_wkt, ls);
+    ret.push_back(ls);
+  }
+  return ret;
+}
+
+#endif //GEOS_VERSION

geos_helpers.hpp

@@ -0,0 +1,17 @@
+#ifndef GEOS_HELPERS_HPP
+#define GEOS_HELPERS_HPP
+
+#ifdef GEOS_VERSION
+
+#include <memory>
+#include <geos/geom/Geometry.h>
+#include "geometry.hpp"
+
+template <typename T>
+std::unique_ptr<geos::geom::Geometry> to_geos(const T& mp);
+template <typename T>
+T from_geos(const std::unique_ptr<geos::geom::Geometry>& g);
+
+#endif // GEOS_VERSION
+
+#endif //GEOS_HELPERS_HPP