Non self-intersecting simplify (#239)

CMakeLists.txt

@@ -77,7 +77,24 @@ ADD_CUSTOM_COMMAND(OUTPUT osmformat.pb.cc osmformat.pb.h
 		   ARGS --cpp_out ${CMAKE_BINARY_DIR} -I ${CMAKE_SOURCE_DIR}/include ${CMAKE_SOURCE_DIR}/include/osmformat.proto)
 
 file(GLOB tilemaker_src_files
-    "src/*.cpp"
+	src/attribute_store.cpp
+	src/geom.cpp
+	src/mbtiles.cpp
+	src/osm_mem_tiles.cpp
+	src/output_object.cpp
+	src/read_pbf.cpp
+	src/shared_data.cpp
+	src/tile_data.cpp
+	src/tile_worker.cpp
+	src/coordinates.cpp
+	src/helpers.cpp
+	src/osm_lua_processing.cpp
+	src/osm_store.cpp
+	src/pbf_blocks.cpp
+	src/read_shp.cpp
+	src/shp_mem_tiles.cpp
+	src/tilemaker.cpp
+	src/write_geometry.cpp
   )
 add_executable(tilemaker vector_tile.pb.cc osmformat.pb.cc ${tilemaker_src_files})
 target_link_libraries(tilemaker ${PROTOBUF_LIBRARY} ${LIBSHP_LIBRARIES} ${SQLITE3_LIBRARIES} ${LUAJIT_LIBRARY} ${LUA_LIBRARIES} ${ZLIB_LIBRARY} ${THREAD_LIB} ${CMAKE_DL_LIBS}
@@ -86,5 +103,5 @@ target_link_libraries(tilemaker ${PROTOBUF_LIBRARY} ${LIBSHP_LIBRARIES} ${SQLITE
 if(MSVC)
     target_link_libraries(tilemaker unofficial::sqlite3::sqlite3)
 endif()
-	
+
 install(TARGETS tilemaker RUNTIME DESTINATION bin)

Makefile

@@ -69,7 +69,7 @@ INC := -I/usr/local/include -isystem ./include -I./src $(LUA_CFLAGS)
 
 all: tilemaker
 
-tilemaker: include/osmformat.pb.o include/vector_tile.pb.o src/mbtiles.o src/pbf_blocks.o src/coordinates.o src/osm_store.o src/helpers.o src/output_object.o src/read_shp.o src/read_pbf.o src/osm_lua_processing.o src/write_geometry.o src/shared_data.o src/tile_worker.o src/tile_data.o src/osm_mem_tiles.o src/shp_mem_tiles.o src/attribute_store.o src/tilemaker.o
+tilemaker: include/osmformat.pb.o include/vector_tile.pb.o src/mbtiles.o src/pbf_blocks.o src/coordinates.o src/osm_store.o src/helpers.o src/output_object.o src/read_shp.o src/read_pbf.o src/osm_lua_processing.o src/write_geometry.o src/shared_data.o src/tile_worker.o src/tile_data.o src/osm_mem_tiles.o src/shp_mem_tiles.o src/attribute_store.o src/tilemaker.o src/geom.o
 	$(CXX) $(CXXFLAGS) -o tilemaker $^ $(INC) $(LIB) $(LDFLAGS)
 
 %.o: %.cpp

include/coordinates.h

@@ -3,7 +3,7 @@
 #define _COORDINATES_H
 
 #include <iostream>
-#include "geomtypes.h"
+#include "geom.h"
 #include <utility>
 #include <unordered_set>
 
@@ -86,85 +86,8 @@ double degp2meter(double degp, double latp);
 
 double meter2degp(double meter, double latp);
 
-template<typename T>
-void insertIntermediateTiles(const T &points, uint baseZoom, std::unordered_set<TileCoordinates> &tileSet) {
-	Point p2(0, 0);
-	for (auto it = points.begin(); it != points.end(); ++it) {
-		// Line is from p1 to p2
-		Point p1 = p2;
-		p2 = *it;
-
-		// Calculate p2 tile, and mark it
-		double tileXf2 = lon2tilexf(p2.x(), baseZoom), tileYf2 = latp2tileyf(p2.y(), baseZoom);
-		TileCoordinate tileX2 = static_cast<TileCoordinate>(tileXf2), tileY2 = static_cast<TileCoordinate>(tileYf2);
-		tileSet.insert(TileCoordinates(tileX2, tileY2));
-		if (it == points.begin()) continue;	// first point, so no line
-
-		// Calculate p1 tile
-		double tileXf1 = lon2tilexf(p1.x(), baseZoom), tileYf1 = latp2tileyf(p1.y(), baseZoom);
-		TileCoordinate tileX1 = static_cast<TileCoordinate>(tileXf1), tileY1 = static_cast<TileCoordinate>(tileYf1);
-		tileSet.insert(TileCoordinates(tileX1,tileY1));
-
-		// Supercover line algorithm from http://eugen.dedu.free.fr/projects/bresenham/
-		int i;                       // loop counter
-		int ystep, xstep;            // the step on y and x axis
-		int error;                   // the error accumulated during the increment
-		int errorprev;               // *vision the previous value of the error variable
-		int y = tileY1, x = tileX1;  // the line points
-		int ddy, ddx;                // compulsory variables: the double values of dy and dx
-		int dx = tileX2 - tileX1;
-		int dy = tileY2 - tileY1;
-
-		if (dy < 0) { ystep = -1; dy = -dy; } else { ystep = 1; }
-		if (dx < 0) { xstep = -1; dx = -dx; } else { xstep = 1; }
-
-		ddy = 2 * dy;  // work with double values for full precision
-		ddx = 2 * dx;
-		if (ddx >= ddy) {  // first octant (0 <= slope <= 1)
-			// compulsory initialization (even for errorprev, needed when dx==dy)
-			errorprev = error = dx;  // start in the middle of the square
-			for (i=0 ; i < dx ; i++) {  // do not use the first point (already done)
-				x += xstep;
-				error += ddy;
-				if (error > ddx){  // increment y if AFTER the middle ( > )
-					y += ystep;
-					error -= ddx;
-					// three cases (octant == right->right-top for directions below):
-					if (error + errorprev < ddx)  // bottom square also
-						tileSet.insert(TileCoordinates(x, y-ystep));
-					else if (error + errorprev > ddx)  // left square also
-						tileSet.insert(TileCoordinates(x-xstep, y));
-					else {  // corner: bottom and left squares also
-						tileSet.insert(TileCoordinates(x, y-ystep));
-						tileSet.insert(TileCoordinates(x-xstep, y));
-					}
-				}
-				tileSet.insert(TileCoordinates(x, y));
-				errorprev = error;
-			}
-		} else {  // the same as above
-			errorprev = error = dy;
-			for (i=0 ; i < dy ; i++){
-				y += ystep;
-				error += ddx;
-				if (error > ddy){
-					x += xstep;
-					error -= ddy;
-					if (error + errorprev < ddy)
-						tileSet.insert(TileCoordinates(x-xstep, y));
-					else if (error + errorprev > ddy)
-						tileSet.insert(TileCoordinates(x, y-ystep));
-					else{
-						tileSet.insert(TileCoordinates(x-xstep, y));
-						tileSet.insert(TileCoordinates(x, y-ystep));
-					}
-				}
-				tileSet.insert(TileCoordinates(x, y));
-				errorprev = error;
-			}
-		}
-	}
-}
+void insertIntermediateTiles(Linestring const &points, uint baseZoom, std::unordered_set<TileCoordinates> &tileSet);
+void insertIntermediateTiles(Ring const &points, uint baseZoom, std::unordered_set<TileCoordinates> &tileSet);
 
 // the range between smallest y and largest y is filled, for each x
 void fillCoveredTiles(std::unordered_set<TileCoordinates> &tileSet);
@@ -184,9 +107,15 @@ public:
 	TileBbox(TileCoordinates i, uint z);
 
 	std::pair<int,int> scaleLatpLon(double latp, double lon) const;
+	std::pair<double, double> floorLatpLon(double latp, double lon) const;
+
 	Box getTileBox() const;
 	Box getExtendBox() const;
 };
 
+// Round coordinates to integer coordinates of bbox
+// TODO: This should be self-intersection aware!!
+MultiPolygon round_coordinates(TileBbox const &bbox, MultiPolygon const &mp);
+
 #endif //_COORDINATES_H
 

include/geom.h

@@ -10,6 +10,7 @@ using uint = unsigned int;
 #include <limits>
 
 // boost::geometry
+#define BOOST_GEOMETRY_INCLUDE_SELF_TURNS
 #include <boost/geometry.hpp>
 #include <boost/geometry/algorithms/intersection.hpp>
 #include <boost/geometry/geometries/geometries.hpp>
@@ -145,5 +146,35 @@ typedef uint64_t WayID;
 typedef std::vector<NodeID> NodeVec;
 typedef std::vector<WayID> WayVec;
 
+// Perform self-intersection aware simplification of geometry types
+Linestring simplify(Linestring const &ls, double max_distance);
+Polygon simplify(Polygon const &p, double max_distance);
+MultiPolygon simplify(MultiPolygon const &mp, double max_distance);
+
+// Combine overlapping elements by performing a union
+template<typename C, typename T>
+void simplify_combine(C &result, T &&new_element)
+{
+    result.push_back(new_element);
+
+   	for(std::size_t i = 0; i < result.size() - 1; ) {
+        if(!boost::geometry::intersects(result[i], result.back())) {
+            ++i;
+            continue;
+        }
+
+        std::vector<T> union_result;
+        boost::geometry::union_(result[i], result.back(), union_result);
+
+        if(union_result.size() != 1) {
+			++i;
+			continue;
+		}
+
+       	result.back() = std::move(union_result[0]);
+       	result.erase(result.begin() + i);
+    } 
+}
+
 #endif //_GEOM_TYPES_H
 

include/helpers.h

@@ -3,7 +3,7 @@
 #define _HELPERS_H
 
 #include <zlib.h>
-#include "geomtypes.h"
+#include "geom.h"
 
 // General helper routines
 

include/osm_lua_processing.h

@@ -6,7 +6,7 @@
 #include <string>
 #include <sstream>
 #include <map>
-#include "geomtypes.h"
+#include "geom.h"
 #include "osm_store.h"
 #include "shared_data.h"
 #include "output_object.h"
@@ -115,19 +115,22 @@ public:
 	// ----	Requests from Lua to write this way/node to a vector tile's Layer
 
     template<class GeometryT>
-    void CorrectGeometry(GeometryT &geom)
+    bool CorrectGeometry(GeometryT &geom)
     {
         geom::correct(geom); // fix wrong orientation
 #if BOOST_VERSION >= 105800
         geom::validity_failure_type failure;
         if (isRelation && !geom::is_valid(geom,failure)) {
             if (verbose) std::cout << "Relation " << originalOsmID << " has " << boost_validity_error(failure) << std::endl;
-            if (failure==10) return; // too few points
         } else if (isWay && !geom::is_valid(geom,failure)) {
             if (verbose) std::cout << "Way " << originalOsmID << " has " << boost_validity_error(failure) << std::endl;
-            if (failure==10) return; // too few points
         }
+        if (failure==boost::geometry::failure_self_intersections || failure == boost::geometry::failure_few_points) 
+			return false;
+		if (failure==boost::geometry::failure_spikes)
+			geom::remove_spikes(geom);
 #endif
+		return true;
     }
 
 	// Add layer

include/osm_store.h

@@ -5,7 +5,7 @@
 #include <utility>
 #include <vector>
 #include <iostream>
-#include "geomtypes.h"
+#include "geom.h"
 #include "coordinates.h"
 namespace geom = boost::geometry;
 

include/output_object.h

@@ -6,7 +6,7 @@
 #include <string>
 #include <map>
 #include <memory>
-#include "geomtypes.h"
+#include "geom.h"
 #include "coordinates.h"
 #include "attribute_store.h"
 #include "osm_store.h"

include/read_shp.h

@@ -6,7 +6,7 @@
 #include <string>
 #include <vector>
 #include <map>
-#include "geomtypes.h"
+#include "geom.h"
 #include "output_object.h"
 #include "osm_lua_processing.h"
 #include "kaguya.hpp"

src/coordinates.cpp

@@ -86,11 +86,16 @@ TileBbox::TileBbox(TileCoordinates i, uint z) {
 }
 
 pair<int,int> TileBbox::scaleLatpLon(double latp, double lon) const {
-	int x = round( (lon -   minLon) / xscale );
-	int y = round( (maxLatp - latp) / yscale );
+	int x = floor( (lon -   minLon) / xscale );
+	int y = floor( (maxLatp - latp) / yscale );
 	return pair<int,int>(x,y);
 }
 
+pair<double, double> TileBbox::floorLatpLon(double latp, double lon) const {
+	auto p = scaleLatpLon(latp, lon);
+	return std::make_pair( -(p.second * yscale - maxLatp), p.first * xscale + minLon);
+}
+
 Box TileBbox::getTileBox() const {
 	double xmargin = (maxLon -minLon )/8192.0;
 	double ymargin = (maxLatp-minLatp)/8192.0;
@@ -103,3 +108,114 @@ Box TileBbox::getExtendBox() const {
 		geom::make<Point>( maxLon+(maxLon-minLon)*(8191.0/8192.0), maxLatp+(maxLatp-minLatp)*2.0));
 }
 
+MultiPolygon round_coordinates(TileBbox const &bbox, MultiPolygon const &mp) 
+{
+	MultiPolygon combined_mp;
+	for(auto new_p: mp) {
+		for(auto &i: new_p.outer()) {
+			auto round_i = bbox.floorLatpLon(i.y(), i.x());
+			i = Point(round_i.second, round_i.first);
+		}
+
+		for(auto &r: new_p.inners()) {
+			for(auto &i: r) {
+				auto round_i = bbox.floorLatpLon(i.y(), i.x());
+				i = Point(round_i.second, round_i.first);
+			}
+		}
+
+		boost::geometry::remove_spikes(new_p);
+		simplify_combine(combined_mp, std::move(new_p));
+	}
+	return combined_mp;
+}
+
+template<typename T>
+void impl_insertIntermediateTiles(T const &points, uint baseZoom, std::unordered_set<TileCoordinates> &tileSet) {
+	Point p2(0, 0);
+	for (auto it = points.begin(); it != points.end(); ++it) {
+		// Line is from p1 to p2
+		Point p1 = p2;
+		p2 = *it;
+
+		// Calculate p2 tile, and mark it
+		double tileXf2 = lon2tilexf(p2.x(), baseZoom), tileYf2 = latp2tileyf(p2.y(), baseZoom);
+		TileCoordinate tileX2 = static_cast<TileCoordinate>(tileXf2), tileY2 = static_cast<TileCoordinate>(tileYf2);
+		tileSet.insert(TileCoordinates(tileX2, tileY2));
+		if (it == points.begin()) continue;	// first point, so no line
+
+		// Calculate p1 tile
+		double tileXf1 = lon2tilexf(p1.x(), baseZoom), tileYf1 = latp2tileyf(p1.y(), baseZoom);
+		TileCoordinate tileX1 = static_cast<TileCoordinate>(tileXf1), tileY1 = static_cast<TileCoordinate>(tileYf1);
+		tileSet.insert(TileCoordinates(tileX1,tileY1));
+
+		// Supercover line algorithm from http://eugen.dedu.free.fr/projects/bresenham/
+		int i;                       // loop counter
+		int ystep, xstep;            // the step on y and x axis
+		int error;                   // the error accumulated during the increment
+		int errorprev;               // *vision the previous value of the error variable
+		int y = tileY1, x = tileX1;  // the line points
+		int ddy, ddx;                // compulsory variables: the double values of dy and dx
+		int dx = tileX2 - tileX1;
+		int dy = tileY2 - tileY1;
+
+		if (dy < 0) { ystep = -1; dy = -dy; } else { ystep = 1; }
+		if (dx < 0) { xstep = -1; dx = -dx; } else { xstep = 1; }
+
+		ddy = 2 * dy;  // work with double values for full precision
+		ddx = 2 * dx;
+		if (ddx >= ddy) {  // first octant (0 <= slope <= 1)
+			// compulsory initialization (even for errorprev, needed when dx==dy)
+			errorprev = error = dx;  // start in the middle of the square
+			for (i=0 ; i < dx ; i++) {  // do not use the first point (already done)
+				x += xstep;
+				error += ddy;
+				if (error > ddx){  // increment y if AFTER the middle ( > )
+					y += ystep;
+					error -= ddx;
+					// three cases (octant == right->right-top for directions below):
+					if (error + errorprev < ddx)  // bottom square also
+						tileSet.insert(TileCoordinates(x, y-ystep));
+					else if (error + errorprev > ddx)  // left square also
+						tileSet.insert(TileCoordinates(x-xstep, y));
+					else {  // corner: bottom and left squares also
+						tileSet.insert(TileCoordinates(x, y-ystep));
+						tileSet.insert(TileCoordinates(x-xstep, y));
+					}
+				}
+				tileSet.insert(TileCoordinates(x, y));
+				errorprev = error;
+			}
+		} else {  // the same as above
+			errorprev = error = dy;
+			for (i=0 ; i < dy ; i++){
+				y += ystep;
+				error += ddx;
+				if (error > ddy){
+					x += xstep;
+					error -= ddy;
+					if (error + errorprev < ddy)
+						tileSet.insert(TileCoordinates(x-xstep, y));
+					else if (error + errorprev > ddy)
+						tileSet.insert(TileCoordinates(x, y-ystep));
+					else{
+						tileSet.insert(TileCoordinates(x-xstep, y));
+						tileSet.insert(TileCoordinates(x, y-ystep));
+					}
+				}
+				tileSet.insert(TileCoordinates(x, y));
+				errorprev = error;
+			}
+		}
+	}
+}
+
+void insertIntermediateTiles(Linestring const &points, uint baseZoom, std::unordered_set<TileCoordinates> &tileSet)
+{
+	impl_insertIntermediateTiles(points, baseZoom, tileSet);
+}
+
+void insertIntermediateTiles(Ring const &points, uint baseZoom, std::unordered_set<TileCoordinates> &tileSet)
+{
+	impl_insertIntermediateTiles(points, baseZoom, tileSet);
+}

src/geom.cpp

@@ -0,0 +1,148 @@
+#include "geom.h"
+
+#include <boost/geometry/geometries/segment.hpp>
+#include <boost/geometry/index/rtree.hpp>
+
+typedef boost::geometry::model::segment<Point> simplify_segment;
+typedef boost::geometry::index::rtree<simplify_segment, boost::geometry::index::quadratic<16>> simplify_rtree;
+
+struct simplify_rtree_counter
+{
+	using value_type = simplify_segment;
+	std::size_t count = 0;
+	void push_back(value_type const &) { ++count; }
+	std::size_t size() const { return count; }
+};
+
+template<typename GeometryType>
+void simplify(GeometryType const &input, GeometryType &output, double max_distance, simplify_rtree const &outer_rtree = simplify_rtree())
+{        
+	simplify_rtree rtree;
+
+    std::deque<std::size_t> nodes(input.size());
+    for(std::size_t i = 0; i < input.size(); ++i) {
+		rtree.insert({ input[i], input[i + 1] });    
+        nodes[i] = i;
+	}
+        
+    std::priority_queue<std::size_t, std::vector<size_t>> pq;
+    for(std::size_t i = 0; i < input.size() - 2; ++i) 
+        pq.push(i);      
+        
+    while(!pq.empty()) {
+        auto entry = pq.top();
+        pq.pop();
+        
+        auto start = nodes[entry];
+        auto middle = nodes[entry + 1];
+        auto end = nodes[entry + 2];                   
+                
+        simplify_segment line(input[start], input[end]);
+        double distance = 0.0;
+        for(auto i = start + 1; i < end; ++i) 
+            distance = std::max(distance, boost::geometry::distance(line, input[i]));          
+    
+        if(boost::geometry::distance(input[start], input[end]) < 2 * max_distance || distance < max_distance) {
+            simplify_rtree_counter result;
+            boost::geometry::index::query(rtree, boost::geometry::index::intersects(line), std::back_inserter(result));
+            boost::geometry::index::query(outer_rtree, boost::geometry::index::intersects(line), std::back_inserter(result));
+
+			std::vector<simplify_segment> nearest;
+			constexpr std::size_t nearest_query_size = 5;
+			boost::geometry::index::query(rtree, boost::geometry::index::nearest(line, nearest_query_size), std::back_inserter(nearest));
+			boost::geometry::index::query(outer_rtree, boost::geometry::index::nearest(line, nearest_query_size), std::back_inserter(nearest));
+
+			double min_distance = std::numeric_limits<double>::max();
+			for(auto const &i: nearest) 
+				min_distance = std::min(min_distance, boost::geometry::distance(line, i));
+
+            std::size_t query_expected = ((start == 0 || end == input.size() - 1) ? 2 : 4);
+            if(result.size() == query_expected && min_distance > max_distance) {
+                nodes.erase(nodes.begin() + entry + 1);
+                rtree.remove(simplify_segment(input[start], input[middle]));
+                rtree.remove(simplify_segment(input[middle], input[end]));
+                rtree.insert(line);
+        
+                if(entry + 2 < nodes.size()) {
+                    pq.push(start);             
+                }
+            }
+        }
+    }
+    
+    for(auto i: nodes)
+        boost::geometry::append(output, input[i]);
+}
+
+Polygon simplify(Polygon const &p, double max_distance) 
+{
+	simplify_rtree outer_rtree;
+	for(std::size_t j = 0; j < p.outer().size() - 1; ++j) 
+		outer_rtree.insert({ p.outer()[j], p.outer()[j + 1] });    
+
+	std::vector<Ring> combined_inners;
+	for(size_t i = 0; i < p.inners().size(); ++i) {
+		Ring new_inner = p.inners()[i];
+		if(boost::geometry::area(new_inner) < 0) {
+			std::reverse(new_inner.begin(), new_inner.end());
+			simplify_combine(combined_inners, std::move(new_inner));
+		}
+	}
+
+	std::vector<Ring> new_inners;
+	for(size_t i = 0; i < combined_inners.size(); ++i) {
+		Ring new_inner;
+		simplify(combined_inners[i], new_inner, max_distance, outer_rtree);
+
+		if(boost::geometry::area(new_inner) > max_distance * max_distance) {
+			simplify_combine(new_inners, std::move(new_inner));
+		}
+	}
+
+	simplify_rtree inners_rtree;
+	for(auto const &inner: new_inners) {
+		for(std::size_t z = 0; z < inner.size() - 1; ++z) 
+			inners_rtree.insert({ inner[z], inner[z + 1] });    
+	} 
+
+	Polygon result;
+	simplify(p.outer(), result.outer(), max_distance, inners_rtree);
+	if(boost::geometry::area(result.outer()) < max_distance * max_distance) {
+		return Polygon();
+	}
+
+	for(auto&& r: new_inners) {
+		std::reverse(r.begin(), r.end());
+		result.inners().push_back(r);
+	} 
+
+	return result;
+}
+
+Linestring simplify(Linestring const &ls, double max_distance) 
+{
+	Linestring result;
+	simplify(ls, result, max_distance);
+	return result;
+}
+
+MultiPolygon simplify(MultiPolygon const &mp, double max_distance) 
+{
+	MultiPolygon combined_mp;
+	for(auto const &p: mp) {
+    	if(!p.outer().empty()) {
+			simplify_combine(combined_mp, Polygon(p));
+		}
+	}
+
+	MultiPolygon result_mp;
+	for(auto const &p: combined_mp) {
+		Polygon new_p = simplify(p, max_distance);
+    	if(!new_p.outer().empty()) {
+			simplify_combine(result_mp, std::move(new_p));
+		}
+	}
+
+	return result_mp;
+}
+

src/osm_lua_processing.cpp

@@ -314,7 +314,7 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) {
 			LatpLon pt = indexStore->nodes_at(osmID);
 			Point p = Point(pt.lon, pt.latp);
 
-            CorrectGeometry(p);
+            if(!CorrectGeometry(p)) return;
 
         	OutputObjectRef oo(new OutputObjectOsmStorePoint(geomType, false,
 	    	    			layers.layerMap[layerName],
@@ -344,7 +344,7 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) {
 				mp.push_back(p);
 			}
 
-            CorrectGeometry(mp);
+            if(!CorrectGeometry(mp)) return;
 
             OutputObjectRef oo(new OutputObjectOsmStoreMultiPolygon(geomType, false,
                             layers.layerMap[layerName],
@@ -355,7 +355,7 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) {
 			// linestring
 			Linestring ls = linestringCached();
 
-            CorrectGeometry(ls);
+            if(!CorrectGeometry(ls)) return;
 
     	    OutputObjectRef oo(new OutputObjectOsmStoreLinestring(geomType, false,
 		    			layers.layerMap[layerName],

src/tilemaker.cpp

@@ -30,7 +30,7 @@
 #include <sys/resource.h>
 #endif
 
-#include "geomtypes.h"
+#include "geom.h"
 
 // Tilemaker code
 #include "helpers.h"

src/write_geometry.cpp

@@ -1,6 +1,10 @@
 #include "write_geometry.h"
 #include <iostream>
 #include "helpers.h"
+
+#include <boost/geometry/geometries/segment.hpp>
+#include <boost/geometry/index/rtree.hpp>
+
 using namespace std;
 namespace geom = boost::geometry;
 extern bool verbose;
@@ -26,20 +30,25 @@ void WriteGeometryVisitor::operator()(const Point &p) const {
 void WriteGeometryVisitor::operator()(const MultiPolygon &mp) const {
 	MultiPolygon current;
 	if (simplifyLevel>0) {
-		// Note that Boost simplify can sometimes produce invalid polygons, resulting in broken coastline etc.
-		// In that case, we just revert to the unsimplified version (at the cost of a larger geometry)
-		// See comments in https://github.com/boostorg/geometry/pull/460
-		// When/if dissolve is merged into Boost.Geometry, we can use that to fix the self-intersections
-		bool v = geom::is_valid(mp);
-		geom::simplify(mp, current, simplifyLevel);
-		if (v && !geom::is_valid(current)) { current=mp; }
+		current = simplify(round_coordinates(*bboxPtr, mp), simplifyLevel);
 	} else {
 		current = mp;
 	}
 
 #if BOOST_VERSION >= 105800
 	geom::validity_failure_type failure;
-	if (verbose && !geom::is_valid(current, failure)) { cout << "Output multipolygon has " << boost_validity_error(failure) << endl; }
+	if (verbose && !geom::is_valid(current, failure)) { 
+		cout << "output multipolygon has " << boost_validity_error(failure) << endl; 
+
+		if (!geom::is_valid(mp, failure)) 
+			cout << "input multipolygon has " << boost_validity_error(failure) << endl; 
+		else
+			cout << "input multipolygon valid" << endl;
+	}
+#else	
+	if (verbose && !geom::is_valid(current)) { 
+		cout << "Output multipolygon is invalid " << endl; 
+	}
 #endif
 
 	pair<int,int> lastPos(0,0);
@@ -71,7 +80,9 @@ void WriteGeometryVisitor::operator()(const MultiPolygon &mp) const {
 void WriteGeometryVisitor::operator()(const MultiLinestring &mls) const {
 	MultiLinestring current;
 	if (simplifyLevel>0) {
-		geom::simplify(mls, current, simplifyLevel);
+		for(auto const &ls: mls) {
+			current.push_back(simplify(ls, simplifyLevel));
+		}
 	} else {
 		current = mls;
 	}
@@ -92,7 +103,7 @@ void WriteGeometryVisitor::operator()(const MultiLinestring &mls) const {
 void WriteGeometryVisitor::operator()(const Linestring &ls) const { 
 	Linestring current;
 	if (simplifyLevel>0) {
-		geom::simplify(ls, current, simplifyLevel);
+		current = simplify(ls, simplifyLevel);
 	} else {
 		current = ls;
 	}