Skip to content

Commit

Permalink
Calculate distance metrics for line features
Browse files Browse the repository at this point in the history 
This patch enables line metrics tracking for LineString/MultiLineString features.
See mapbox/geojson-vt#93 for details.

Based on patch from @lbud (Lauren Budorick).
  • Loading branch information
@pozdnyakov
pozdnyakov committed Aug 3, 2018
1 parent 2304c52 commit 5cc008e
Show file tree
Hide file tree
Showing 11 changed files with 281 additions and 72 deletions.
29 changes: 16 additions & 13 deletions include/mapbox/geojsonvt.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,9 @@ struct TileOptions {

// tile buffer on each side
uint16_t buffer = 64;

// enable line metrics tracking for LineString/MultiLineString features
bool lineMetrics = false;
};

struct Options : TileOptions {
Expand Down Expand Up @@ -72,15 +75,15 @@ inline const Tile geoJSONToTile(const geojson& geojson_,
auto tolerance = (options.tolerance / options.extent) / z2;
auto features = detail::convert(features_, tolerance);
if (wrap) {
features = detail::wrap(features, double(options.buffer) / options.extent);
features = detail::wrap(features, double(options.buffer) / options.extent, options.lineMetrics);
}
if (clip) {
if (clip || options.lineMetrics) {
const double p = options.buffer / options.extent;

const auto left = detail::clip<0>(features, (x - p) / z2, (x + 1 + p) / z2, -1, 2);
features = detail::clip<1>(left, (y - p) / z2, (y + 1 + p) / z2, -1, 2);
const auto left = detail::clip<0>(features, (x - p) / z2, (x + 1 + p) / z2, -1, 2, options.lineMetrics);
features = detail::clip<1>(left, (y - p) / z2, (y + 1 + p) / z2, -1, 2, options.lineMetrics);
}
return detail::InternalTile({ features, z, x, y, options.extent, tolerance }).tile;
return detail::InternalTile({ features, z, x, y, options.extent, tolerance, options.lineMetrics }).tile;
}

class GeoJSONVT {
Expand All @@ -94,7 +97,7 @@ class GeoJSONVT {
const uint32_t z2 = 1u << options.maxZoom;

auto converted = detail::convert(features_, (options.tolerance / options.extent) / z2);
auto features = detail::wrap(converted, double(options.buffer) / options.extent);
auto features = detail::wrap(converted, double(options.buffer) / options.extent, options.lineMetrics);

splitTile(features, 0, 0, 0);
}
Expand Down Expand Up @@ -186,7 +189,7 @@ class GeoJSONVT {

it = tiles
.emplace(id,
detail::InternalTile{ features, z, x, y, options.extent, tolerance })
detail::InternalTile{ features, z, x, y, options.extent, tolerance, options.lineMetrics })
.first;
stats[z] = (stats.count(z) ? stats[z] + 1 : 1);
total++;
Expand Down Expand Up @@ -230,19 +233,19 @@ class GeoJSONVT {
const auto& min = tile.bbox.min;
const auto& max = tile.bbox.max;

const auto left = detail::clip<0>(features, (x - p) / z2, (x + 0.5 + p) / z2, min.x, max.x);
const auto left = detail::clip<0>(features, (x - p) / z2, (x + 0.5 + p) / z2, min.x, max.x, options.lineMetrics);

splitTile(detail::clip<1>(left, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y), z + 1,
splitTile(detail::clip<1>(left, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y, options.lineMetrics), z + 1,
x * 2, y * 2, cz, cx, cy);
splitTile(detail::clip<1>(left, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y), z + 1,
splitTile(detail::clip<1>(left, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y, options.lineMetrics), z + 1,
x * 2, y * 2 + 1, cz, cx, cy);

const auto right =
detail::clip<0>(features, (x + 0.5 - p) / z2, (x + 1 + p) / z2, min.x, max.x);
detail::clip<0>(features, (x + 0.5 - p) / z2, (x + 1 + p) / z2, min.x, max.x, options.lineMetrics);

splitTile(detail::clip<1>(right, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y), z + 1,
splitTile(detail::clip<1>(right, (y - p) / z2, (y + 0.5 + p) / z2, min.y, max.y, options.lineMetrics), z + 1,
x * 2 + 1, y * 2, cz, cx, cy);
splitTile(detail::clip<1>(right, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y), z + 1,
splitTile(detail::clip<1>(right, (y + 0.5 - p) / z2, (y + 1 + p) / z2, min.y, max.y, options.lineMetrics), z + 1,
x * 2 + 1, y * 2 + 1, cz, cx, cy);

// if we sliced further down, no need to keep source geometry
Expand Down
128 changes: 93 additions & 35 deletions include/mapbox/geojsonvt/clip.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,12 @@ namespace detail {
template <uint8_t I>
class clipper {
public:
clipper(double k1_, double k2_, bool lineMetrics_ = false)
: k1(k1_), k2(k2_), lineMetrics(lineMetrics_) {}

const double k1;
const double k2;
const bool lineMetrics;

vt_geometry operator()(const vt_point& point) const {
return point;
Expand Down Expand Up @@ -81,76 +85,110 @@ class clipper {
}

private:
vt_line_string newSlice(vt_multi_line_string& parts, vt_line_string& slice, double dist) const {
if (!slice.empty()) {
slice.dist = dist;
parts.push_back(std::move(slice));
vt_line_string newSlice(const vt_line_string& line) const {
vt_line_string slice;
slice.dist = line.dist;
if (lineMetrics) {
slice.segStart = line.segStart;
slice.segEnd = line.segEnd;
}
return {};
return slice;
}

void clipLine(const vt_line_string& line, vt_multi_line_string& slices) const {

const double dist = line.dist;
const size_t len = line.size();
double lineLen = line.segStart;
double segLen = 0.0;
double t = 0.0;

if (len < 2)
return;

vt_line_string slice;
vt_line_string slice = newSlice(line);

for (size_t i = 0; i < (len - 1); ++i) {
const auto& a = line[i];
const auto& b = line[i + 1];
const double ak = get<I>(a);
const double bk = get<I>(b);

if (lineMetrics) segLen = std::hypot((b.x - a.x), (b.y - a.y));

if (ak < k1) {
if (bk > k2) { // ---|-----|-->
slice.push_back(intersect<I>(a, b, k1));
slice.push_back(intersect<I>(a, b, k2));
slice = newSlice(slices, slice, dist);
t = calc_progress<I>(a, b, k1);
slice.push_back(intersect<I>(a, b, k1, t));
if (lineMetrics) slice.segStart = lineLen + segLen * t;

t = calc_progress<I>(a, b, k2);
slice.push_back(intersect<I>(a, b, k2, t));
if (lineMetrics) slice.segEnd = lineLen + segLen * t;
slices.push_back(std::move(slice));

slice = newSlice(line);

} else if (bk >= k1) { // ---|--> |
slice.push_back(intersect<I>(a, b, k1));
t = calc_progress<I>(a, b, k1);
slice.push_back(intersect<I>(a, b, k1, t));
if (lineMetrics) slice.segStart = lineLen + segLen * t;

if (i == len - 2)
slice.push_back(b); // last point
}
} else if (ak >= k2) {
if (bk < k1) { // <--|-----|---
slice.push_back(intersect<I>(a, b, k2));
slice.push_back(intersect<I>(a, b, k1));
slice = newSlice(slices, slice, dist);
t = calc_progress<I>(a, b, k2);
slice.push_back(intersect<I>(a, b, k2, t));
if (lineMetrics) slice.segStart = lineLen + segLen * t;

t = calc_progress<I>(a, b, k1);
slice.push_back(intersect<I>(a, b, k1, t));
if (lineMetrics) slice.segEnd = lineLen + segLen * t;

slices.push_back(std::move(slice));

slice = newSlice(line);
} else if (bk < k2) { // | <--|---
slice.push_back(intersect<I>(a, b, k2));
t = calc_progress<I>(a, b, k2);
slice.push_back(intersect<I>(a, b, k2, t));
if (lineMetrics) slice.segStart = lineLen + segLen * t;

if (i == len - 2)
slice.push_back(b); // last point
}
} else {
slice.push_back(a);

if (bk < k1) { // <--|--- |
slice.push_back(intersect<I>(a, b, k1));
slice = newSlice(slices, slice, dist);
t = calc_progress<I>(a, b, k1);
slice.push_back(intersect<I>(a, b, k1, t));
if (lineMetrics) slice.segEnd = lineLen + segLen * t;
slices.push_back(std::move(slice));
slice = newSlice(line);

} else if (bk > k2) { // | ---|-->
slice.push_back(intersect<I>(a, b, k2));
slice = newSlice(slices, slice, dist);
t = calc_progress<I>(a, b, k2);
slice.push_back(intersect<I>(a, b, k2, t));
if (lineMetrics) slice.segEnd = lineLen + segLen * t;
slices.push_back(std::move(slice));
slice = newSlice(line);

} else if (i == len - 2) { // | --> |
slice.push_back(b);
}
}

if (lineMetrics) lineLen += segLen;
}

// add the final slice
newSlice(slices, slice, dist);
if (!slice.empty()) { // add the final slice
slice.segEnd = lineLen;
slices.push_back(std::move(slice));
}
}

vt_linear_ring clipRing(const vt_linear_ring& ring) const {
const size_t len = ring.size();

vt_linear_ring slice;
slice.area = ring.area;

Expand All @@ -165,27 +203,31 @@ class clipper {

if (ak < k1) {
if (bk >= k1) {
slice.push_back(intersect<I>(a, b, k1)); // ---|--> |
if (bk > k2) // ---|-----|-->
slice.push_back(intersect<I>(a, b, k2));
// ---|--> |
slice.push_back(intersect<I>(a, b, k1, calc_progress<I>(a, b, k1)));
if (bk > k2)
// ---|-----|-->
slice.push_back(intersect<I>(a, b, k2, calc_progress<I>(a, b, k2)));
else if (i == len - 2)
slice.push_back(b); // last point
}
} else if (ak >= k2) {
if (bk < k2) { // | <--|---
slice.push_back(intersect<I>(a, b, k2));
slice.push_back(intersect<I>(a, b, k2, calc_progress<I>(a, b, k2)));
if (bk < k1) // <--|-----|---
slice.push_back(intersect<I>(a, b, k1));
slice.push_back(intersect<I>(a, b, k1, calc_progress<I>(a, b, k1)));
else if (i == len - 2)
slice.push_back(b); // last point
}
} else {
slice.push_back(a);
if (bk < k1) // <--|--- |
slice.push_back(intersect<I>(a, b, k1));
else if (bk > k2) // | ---|-->
slice.push_back(intersect<I>(a, b, k2));
// | --> |
slice.push_back(a);
if (bk < k1)
// <--|--- |
slice.push_back(intersect<I>(a, b, k1, calc_progress<I>(a, b, k1)));
else if (bk > k2)
// | ---|-->
slice.push_back(intersect<I>(a, b, k2, calc_progress<I>(a, b, k2)));
}
}

Expand Down Expand Up @@ -214,7 +256,8 @@ inline vt_features clip(const vt_features& features,
const double k1,
const double k2,
const double minAll,
const double maxAll) {
const double maxAll,
const bool lineMetrics) {

if (minAll >= k1 && maxAll < k2) // trivial accept
return features;
Expand All @@ -239,7 +282,22 @@ inline vt_features clip(const vt_features& features,
continue;

} else {
clipped.emplace_back(vt_geometry::visit(geom, clipper<I>{ k1, k2 }), props, id);
const auto& clippedGeom = vt_geometry::visit(geom, clipper<I>{ k1, k2, lineMetrics });

clippedGeom.match(
[&](const auto&) {
clipped.emplace_back(clippedGeom, props, id);
},
[&](const vt_multi_line_string& result) {
if (lineMetrics) {
for (const auto& segment : result) {
clipped.emplace_back(segment, props, id);
}
} else {
clipped.emplace_back(clippedGeom, props, id);
}
}
);
}
}

Expand Down
7 changes: 4 additions & 3 deletions include/mapbox/geojsonvt/convert.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,14 @@ struct project {
for (size_t i = 0; i < len - 1; ++i) {
const auto& a = result[i];
const auto& b = result[i + 1];
// use Manhattan distance instead of Euclidian to avoid expensive square root
// computation
result.dist += std::abs(b.x - a.x) + std::abs(b.y - a.y);
result.dist += std::hypot((b.x - a.x), (b.y - a.y));
}

simplify(result, tolerance);

result.segStart = 0;
result.segEnd = result.dist;

return result;
}

Expand Down
18 changes: 14 additions & 4 deletions include/mapbox/geojsonvt/tile.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ class InternalTile {
const double z2;
const double tolerance;
const double sq_tolerance;
const bool lineMetrics;

vt_features source_features;
mapbox::geometry::box<double> bbox = { { 2, 1 }, { -1, 0 } };
Expand All @@ -36,14 +37,16 @@ class InternalTile {
const uint32_t x_,
const uint32_t y_,
const uint16_t extent_,
const double tolerance_)
const double tolerance_,
const bool lineMetrics_)
: extent(extent_),
z(z_),
x(x_),
y(y_),
z2(std::pow(2, z)),
tolerance(tolerance_),
sq_tolerance(tolerance_ * tolerance_) {
sq_tolerance(tolerance_ * tolerance_),
lineMetrics(lineMetrics_) {

for (const auto& feature : source) {
const auto& geom = feature.geometry;
Expand Down Expand Up @@ -74,8 +77,15 @@ class InternalTile {
const property_map& props,
const optional<identifier>& id) {
const auto new_line = transform(line);
if (!new_line.empty())
tile.features.push_back({ std::move(new_line), props, id });
if (!new_line.empty()) {
if (lineMetrics) {
property_map newProps = props;
newProps["mapbox_clip_start"] = line.segStart / line.dist;
newProps["mapbox_clip_end"] = line.segEnd / line.dist;
tile.features.push_back({ std::move(new_line), newProps, id });
} else
tile.features.push_back({ std::move(new_line), props, id });
}
}

void addFeature(const vt_polygon& polygon,
Expand Down
Loading

0 comments on commit 5cc008e

Please sign in to comment.