1 #pragma once
2 
3 #include <mapbox/geojsonvt/types.hpp>
4 
5 namespace mapbox {
6 namespace geojsonvt {
7 namespace detail {
8 
9 // square distance from a point to a segment
getSqSegDist(const vt_point & p,const vt_point & a,const vt_point & b)10 inline double getSqSegDist(const vt_point& p, const vt_point& a, const vt_point& b) {
11     double x = a.x;
12     double y = a.y;
13     double dx = b.x - a.x;
14     double dy = b.y - a.y;
15 
16     if ((dx != 0.0) || (dy != 0.0)) {
17 
18         const double t = ((p.x - a.x) * dx + (p.y - a.y) * dy) / (dx * dx + dy * dy);
19 
20         if (t > 1) {
21             x = b.x;
22             y = b.y;
23 
24         } else if (t > 0) {
25             x += dx * t;
26             y += dy * t;
27         }
28     }
29 
30     dx = p.x - x;
31     dy = p.y - y;
32 
33     return dx * dx + dy * dy;
34 }
35 
36 // calculate simplification data using optimized Douglas-Peucker algorithm
simplify(std::vector<vt_point> & points,size_t first,size_t last,double sqTolerance)37 inline void simplify(std::vector<vt_point>& points, size_t first, size_t last, double sqTolerance) {
38     double maxSqDist = sqTolerance;
39     size_t index = 0;
40     const int64_t mid = (last - first) >> 1;
41     int64_t minPosToMid = last - first;
42 
43     for (auto i = first + 1; i < last; i++) {
44         const double sqDist = getSqSegDist(points[i], points[first], points[last]);
45 
46         if (sqDist > maxSqDist) {
47             index = i;
48             maxSqDist = sqDist;
49 
50         } else if (sqDist == maxSqDist) {
51             // a workaround to ensure we choose a pivot close to the middle of the list,
52             // reducing recursion depth, for certain degenerate inputs
53             // https://github.com/mapbox/geojson-vt/issues/104
54             auto posToMid = std::abs(static_cast<int64_t>(i) - mid);
55             if (posToMid < minPosToMid) {
56                 index = i;
57                 minPosToMid = posToMid;
58             }
59         }
60     }
61 
62     if (maxSqDist > sqTolerance) {
63         // save the point importance in squared pixels as a z coordinate
64         points[index].z = maxSqDist;
65         if (index - first > 1)
66             simplify(points, first, index, sqTolerance);
67         if (last - index > 1)
68             simplify(points, index, last, sqTolerance);
69     }
70 }
71 
simplify(std::vector<vt_point> & points,double tolerance)72 inline void simplify(std::vector<vt_point>& points, double tolerance) {
73     const size_t len = points.size();
74 
75     // always retain the endpoints (1 is the max value)
76     points[0].z = 1.0;
77     points[len - 1].z = 1.0;
78 
79     simplify(points, 0, len - 1, tolerance * tolerance);
80 }
81 
82 } // namespace detail
83 } // namespace geojsonvt
84 } // namespace mapbox
85