From b7d3a798cd200b172fdcdd11edd9513c28424713 Mon Sep 17 00:00:00 2001 From: Kevin Kreiser Date: Wed, 8 Oct 2014 13:47:36 -0400 Subject: [PATCH] add ramer douglas peucker to the list of simplify converters. also add tests for the simplify converters but so far only testing peucker. also add a method to the wkt factor for converting geometry container to wkt --- include/mapnik/simplify_converter.hpp | 97 +++++++++++++++++++- include/mapnik/wkt/wkt_factory.hpp | 11 +++ tests/cpp_tests/simplify_converters_test.cpp | 97 ++++++++++++++++++++ 3 files changed, 204 insertions(+), 1 deletion(-) create mode 100644 tests/cpp_tests/simplify_converters_test.cpp diff --git a/include/mapnik/simplify_converter.hpp b/include/mapnik/simplify_converter.hpp index 5fa53b815..80a87b5ab 100644 --- a/include/mapnik/simplify_converter.hpp +++ b/include/mapnik/simplify_converter.hpp @@ -171,7 +171,8 @@ private: switch (algorithm_) { case visvalingam_whyatt: - return output_vertex_cached(x, y); + case douglas_peucker: + return output_vertex_cached(x, y); case radial_distance: return output_vertex_distance(x, y); case zhao_saalfeld: @@ -359,6 +360,8 @@ private: return status_ = process; case zhao_saalfeld: return status_ = cache; + case douglas_peucker: + return init_vertices_RDP(); default: throw std::runtime_error("simplification algorithm not yet implemented"); } @@ -433,6 +436,98 @@ private: return status_ = process; } + void RDP(std::vector& vertices, const size_t start, const size_t end) + { + // Squared length of a vector + auto sqlen = [] (const vertex2d& vec) { return vec.x*vec.x + vec.y*vec.y; }; + // Compute square distance of p to a line segment + auto segment_distance = [&sqlen] (const vertex2d& p, const vertex2d& a, const vertex2d& b, const vertex2d& dir, double dir_sq_len) + { + //special case where segment has same start and end point at which point we are just doing a radius check + if(dir_sq_len == 0) + return sqlen(vertex2d(p.x - b.x, p.y - b.y, SEG_END)); + + //project p onto dir by ((p dot dir / dir dot dir) * dir) + double scale = ((p.x - a.x) * dir.x + (p.y - a.y) * dir.y) / dir_sq_len; + double projected_x = dir.x * scale; + double projected_y = dir.y * scale; + double projected_origin_distance = projected_x * projected_x + projected_y * projected_y; + //projected point doesn't lie on the segment + if(projected_origin_distance > dir_sq_len) + { + //projected point lies past the end of the segment + if(scale > 0) + return sqlen(vertex2d(p.x - b.x, p.y - b.y, SEG_END)); + //projected point lies before the beginning of the segment + else + return sqlen(vertex2d(p.x - a.x, p.y - a.y, SEG_END)); + }//projected point lies on the segment + else + return sqlen(vertex2d(p.x - (projected_x + a.x), p.y - (projected_y + a.y), SEG_END)); + }; + + // Compute the directional vector along the segment + vertex2d dir = vertex2d(vertices[end].x - vertices[start].x, vertices[end].y - vertices[start].y, SEG_END); + double dir_sq_len = sqlen(dir); + + // Find the point with the maximum distance from this line segment + double max = std::numeric_limits::min(); + size_t keeper = 0; + for(size_t i = start + 1; i < end; ++i) + { + double d = segment_distance(vertices[i], vertices[start], vertices[end], dir, dir_sq_len); + if (d > max) + { + keeper = i; + max = d; + } + } + + // Split at the vertex that is furthest outside of the tolerance + // NOTE: we work in square distances to avoid sqrt so we sqaure tolerance accordingly + if (max > tolerance_ * tolerance_) + { + // Make sure not to smooth out the biggest outlier (keeper) + if(keeper - start != 1) RDP(vertices, start, keeper); + if(end - keeper != 1) RDP(vertices, keeper, end); + }// Everyone between the start and the end was close enough to the line + else + { + // Mark each of them as discarded + for(size_t i = start + 1; i < end; ++i) + { + vertices[i].cmd = SEG_END; + } + } + } + + status init_vertices_RDP() + { + //slurp out the original vertices + std::vector vertices; + //vertices.reserve(geom_.size()); + vertex2d vtx(vertex2d::no_init); + while ((vtx.cmd = geom_.vertex(&vtx.x, &vtx.y)) != SEG_END) + { + vertices.push_back(vtx); + } + + //run ramer douglas peucker on it + if(vertices.size() > 2) + { + RDP(vertices, 0, vertices.size() - 1); + } + + //slurp the points back out that haven't been marked as discarded + for(vertex2d& vertex : vertices) + { + if(vertex.cmd != SEG_END) + vertices_.emplace_back(vertex); + } + + return status_ = process; + } + Geometry& geom_; double tolerance_; status status_; diff --git a/include/mapnik/wkt/wkt_factory.hpp b/include/mapnik/wkt/wkt_factory.hpp index c546d5953..837574c31 100644 --- a/include/mapnik/wkt/wkt_factory.hpp +++ b/include/mapnik/wkt/wkt_factory.hpp @@ -27,6 +27,8 @@ #include #include #include +#include +#include // stl #include @@ -43,6 +45,15 @@ inline bool from_wkt(std::string const& wkt, mapnik::geometry_container & paths) return qi::phrase_parse(first, last, g, space, paths); } +inline bool to_wkt(const mapnik::geometry_container& paths, std::string& wkt) +{ + using sink_type = std::back_insert_iterator; + static const mapnik::wkt::wkt_multi_generator generator; + sink_type sink(wkt); + return boost::spirit::karma::generate(sink, generator, paths); +} + + } #endif // MAPNIK_WKT_FACTORY_HPP diff --git a/tests/cpp_tests/simplify_converters_test.cpp b/tests/cpp_tests/simplify_converters_test.cpp new file mode 100644 index 000000000..189d22350 --- /dev/null +++ b/tests/cpp_tests/simplify_converters_test.cpp @@ -0,0 +1,97 @@ +#include +#include +#include +#include +#include +#include "utils.hpp" + +#include +#include +#include +#include +#include +#include +#include + +// stl +#include + +// Ramer Douglas Peucker generalization +void ramer_douglas_peucker(const std::string& wkt_in, double tolerance, const std::string& expected) +{ + //grab the geom + mapnik::geometry_container multi_input; + if (!mapnik::from_wkt(wkt_in , multi_input)) + { + throw std::runtime_error("Failed to parse WKT"); + } + //setup the generalization + mapnik::simplify_converter generalizer(multi_input.front()); + generalizer.set_simplify_algorithm(mapnik::simplify_algorithm_from_string("douglas-peucker").get()); + generalizer.set_simplify_tolerance(tolerance); + //suck the vertices back out of it + mapnik::geometry_type* output = new mapnik::geometry_type(multi_input.front().type()); + mapnik::CommandType cmd; + double x, y; + while((cmd = (mapnik::CommandType)generalizer.vertex(&x, &y)) != mapnik::SEG_END) + { + output->push_vertex(x, y, cmd); + } + //construct the answer + mapnik::geometry_container multi_out; + multi_out.push_back(output); + std::string wkt_out; + BOOST_TEST(mapnik::to_wkt(multi_out, wkt_out)); + BOOST_TEST_EQ(wkt_out, expected); +} + +int main(int argc, char** argv) +{ + std::vector args; + for (int i=1;i