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

This commit is contained in:
Kevin Kreiser 2014-10-08 13:47:36 -04:00
parent 00382775a9
commit b7d3a798cd
3 changed files with 204 additions and 1 deletions

View file

@ -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<vertex2d>& 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<double>::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<vertex2d> 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_;

View file

@ -27,6 +27,8 @@
#include <mapnik/geometry.hpp>
#include <mapnik/geometry_container.hpp>
#include <mapnik/wkt/wkt_grammar.hpp>
#include <mapnik/wkt/wkt_generator_grammar.hpp>
#include <boost/spirit/include/karma.hpp>
// stl
#include <string>
@ -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<std::string>;
static const mapnik::wkt::wkt_multi_generator<sink_type, mapnik::geometry_container> generator;
sink_type sink(wkt);
return boost::spirit::karma::generate(sink, generator, paths);
}
}
#endif // MAPNIK_WKT_FACTORY_HPP

View file

@ -0,0 +1,97 @@
#include <boost/detail/lightweight_test.hpp>
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include "utils.hpp"
#include <mapnik/layer.hpp>
#include <mapnik/wkt/wkt_factory.hpp>
#include <mapnik/wkt/wkt_grammar_impl.hpp>
#include <mapnik/wkt/wkt_generator_grammar.hpp>
#include <mapnik/wkt/wkt_generator_grammar_impl.hpp>
#include <mapnik/simplify.hpp>
#include <mapnik/simplify_converter.hpp>
// stl
#include <stdexcept>
// 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<mapnik::geometry_type> 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<std::string> args;
for (int i=1;i<argc;++i)
{
args.push_back(argv[i]);
}
bool quiet = std::find(args.begin(), args.end(), "-q")!=args.end();
BOOST_TEST(set_working_dir(args));
ramer_douglas_peucker( std::string("LineString(0 0,2 2,3 5,4 1,5 0,6 7,7 0)"),
4,
std::string("LineString(0 0,6 7,7 0)"));
ramer_douglas_peucker( std::string("LineString(0 0,2 2,3 5,4 1,5 0,6 7,7 0)"),
2,
std::string("LineString(0 0,3 5,5 0,6 7,7 0)"));
ramer_douglas_peucker( std::string("LineString(10 0,9 -4,7 -7,4 -9,0 -10,-4 -9,-7 -7,-9 -4,-10 0,-9 4,-7 7,-4 9,0 10,4 9,7 7,9 4)"),
4,
std::string("LineString(10 0,0 -10,-10 0,0 10,9 4)"));
ramer_douglas_peucker( std::string("LineString(0 0,1 1,2 2,0 10,0 0)"),
10,
std::string("LineString(0 0,0 0)"));
ramer_douglas_peucker( std::string("LineString(0 0,1 1,2 2,0 10,0 0)"),
8,
std::string("LineString(0 0,0 10,0 0)"));
ramer_douglas_peucker( std::string("LineString(0 0,1 1,2 2,0 10,0 0)"),
1,
std::string("LineString(0 0,2 2,0 10,0 0)"));
ramer_douglas_peucker( std::string("LineString(0 0, 1 -1, 2 2, 0 -10, 0 0, -5 7, 4 6)"),
3,
std::string("LineString(0 0,0 -10,-5 7,4 6)"));
if (!::boost::detail::test_errors())
{
if (quiet) std::clog << "\x1b[1;32m.\x1b[0m";
else std::clog << "C++ simplify conversions: \x1b[1;32m✓ \x1b[0m\n";
::boost::detail::report_errors_remind().called_report_errors_function = true;
}
else
{
return ::boost::report_errors();
}
}