From 6c58ae416235f2b5ee0957049386d94673424946 Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Fri, 27 Oct 2023 22:39:01 -0700 Subject: [PATCH] Add --buffer-polygons-outward --- README.md | 1 + geometry.cpp | 132 +++++++++++++++++++++++++---------------------- geometry.hpp | 2 +- main.cpp | 1 + man/tippecanoe.1 | 2 + options.hpp | 1 + tile.cpp | 20 +++---- 7 files changed, 84 insertions(+), 75 deletions(-) diff --git a/README.md b/README.md index 6a1c8557..ae5a5eb9 100644 --- a/README.md +++ b/README.md @@ -496,6 +496,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate * `-pt` or `--no-tiny-polygon-reduction`: Don't combine the area of very small polygons into small squares that represent their combined area. * `-pT` or `--no-tiny-polygon-reduction-at-maximum-zoom`: Combine the area of very small polygons into small squares that represent their combined area only at zoom levels below the maximum. * `--tiny-polygon-size=`_size_: Use the specified _size_ for tiny polygons instead of the default 2. Anything above 6 or so will lead to visible artifacts with the default tile detail. + * `-aO` or `--buffer-polygons-outward`: Buffer polygons, except those that have been detected to be continous by `--no-simplification-of-shared-nodes`, outward slightly to prevent narrow polygon spindles fom collapsing away. * `-av` or `--visvalingam`: Use Visvalingam's simplification algorithm rather than Douglas-Peucker's. ### Attempts to improve shared polygon boundaries diff --git a/geometry.cpp b/geometry.cpp index 6013efc0..646d2db9 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -197,6 +197,23 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *still_needs_sim double area = get_area(geom, i, j); + long long minx = LLONG_MAX; + long long maxx = LLONG_MIN; + long long miny = LLONG_MAX; + long long maxy = LLONG_MIN; + for (auto const &d : geom) { + minx = std::min(minx, (long long) d.x); + maxx = std::max(maxx, (long long) d.x); + miny = std::min(miny, (long long) d.y); + maxy = std::max(maxy, (long long) d.y); + } + if (area > 0 && area <= pixel * pixel && area < (maxx - minx) * (maxy - miny) / 5) { + // if the polygon doesn't use most of its area, + // don't let it be dust, because the shape is + // probably something weird and interesting. + area = pixel * pixel * 2; + } + // XXX There is an ambiguity here: If the area of a ring is 0 and it is followed by holes, // we don't know whether the area-0 ring was a hole too or whether it was the outer ring // that these subsequent holes are somehow being subtracted from. I hope that if a polygon @@ -485,6 +502,26 @@ drawvec impose_tile_boundaries(drawvec &geom, long long extent) { return out; } +bool is_shared_node(draw d, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos) { + if (nodepos > 0) { + // offset to global + if (z != 0) { + d.x += tx * (1LL << (32 - z)); + d.y += ty * (1LL << (32 - z)); + } + + // to quadkey + struct node n; + n.index = encode_quadkey((unsigned) d.x, (unsigned) d.y); + + if (bsearch(&n, shared_nodes_map, nodepos / sizeof(node), sizeof(node), nodecmp) != NULL) { + return true; + } + } + + return false; +} + drawvec simplify_lines(drawvec &geom, int z, int tx, int ty, int detail, bool mark_tile_bounds, double simplification, size_t retain, drawvec const &shared_nodes, struct node *shared_nodes_map, size_t nodepos) { int res = 1 << (32 - detail - z); long long area = 1LL << (32 - z); @@ -514,21 +551,8 @@ drawvec simplify_lines(drawvec &geom, int z, int tx, int ty, int detail, bool ma geom[i].necessary = true; } - if (nodepos > 0) { - // offset to global - draw d = geom[i]; - if (z != 0) { - d.x += tx * (1LL << (32 - z)); - d.y += ty * (1LL << (32 - z)); - } - - // to quadkey - struct node n; - n.index = encode_quadkey((unsigned) d.x, (unsigned) d.y); - - if (bsearch(&n, shared_nodes_map, nodepos / sizeof(node), sizeof(node), nodecmp) != NULL) { - geom[i].necessary = true; - } + if (is_shared_node(geom[i], z, tx, ty, shared_nodes_map, nodepos)) { + geom[i].necessary = true; } } } @@ -1399,33 +1423,24 @@ drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigne return out; } -static double anglediff(double dx1, double dy1, double dx2, double dy2) { - // https://gist.github.com/e-n-f/417cd85f6d9b00dab887276376e77f30 - // https://twitter.com/enf/status/795798781186830341 +bool is_continous_poly(drawvec const &geom, size_t i, size_t j, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos) { + size_t count = 0; - double cross = dx1 * dy2 - dy1 * dx2; - double sd = sqrt(dx1 * dx1 + dy1 * dy1); - double cd = sqrt(dx2 * dx2 + dy2 * dy2); + // j - 1 to avoid double-counting the duplicate last node + for (; i < j - 1; i++) { + if (is_shared_node(geom[i], z, tx, ty, shared_nodes_map, nodepos)) { + count++; - double dot = dx1 * dx2 + dy1 * dy2; - double ret; - if (dot < 0) { - if (cross < 0) { - ret = -M_PI - asin(cross / (sd * cd)); - } else { - ret = M_PI - asin(cross / (sd * cd)); + if (count >= 4) { + return true; + } } - } else { - ret = asin(cross / (sd * cd)); } - if (ret < 0) { - // make it 0 to 360, not -180 to 180 - ret += 2 * M_PI; - } - return ret; + + return false; } -drawvec buffer_poly(drawvec const &geom, double buffer) { +drawvec buffer_poly(drawvec const &geom, double buffer, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos) { drawvec out = geom; if (buffer != 0) { @@ -1438,38 +1453,33 @@ drawvec buffer_poly(drawvec const &geom, double buffer) { } } - for (size_t k = i; k < j - 1; k++) { - draw p0 = geom[(k + 0 - i) % (j - i - 1) + i]; - draw p1 = geom[(k + 1 - i) % (j - i - 1) + i]; - draw p2 = geom[(k + 2 - i) % (j - i - 1) + i]; + if (!is_continous_poly(geom, i, j, z, tx, ty, shared_nodes_map, nodepos)) { + for (size_t k = i; k < j - 1; k++) { + draw p0 = geom[(k + 0 - i) % (j - i - 1) + i]; + draw p1 = geom[(k + 1 - i) % (j - i - 1) + i]; + draw p2 = geom[(k + 2 - i) % (j - i - 1) + i]; - double adiff = 2 * M_PI - anglediff(p0.x - p1.x, p0.y - p1.y, - p2.x - p1.x, p2.y - p1.y); + double a10 = atan2(p1.y - p0.y, p1.x - p0.x); + double a21 = atan2(p2.y - p1.y, p2.x - p1.x); - double a1 = atan2(p1.y - p0.y, p1.x - p0.x); + double dx = cos(a10 - 90 * M_PI / 180) + cos(a21 - 90 * M_PI / 180); + double dy = sin(a10 - 90 * M_PI / 180) + sin(a21 - 90 * M_PI / 180); - printf("at %zu %zu %zu ", - (k + 0 - i) % (j - i - 1) + i, - (k + 1 - i) % (j - i - 1) + i, - (k + 2 - i) % (j - i - 1) + i); - printf("%lld,%lld to %lld,%lld to %lld,%lld: ", p0.x, p0.y, p1.x, p1.y, p2.x, p2.y); - printf("angle %f then %f\n", - atan2(p1.y - p0.y, p1.x - p0.x) * 180 / M_PI, - atan2(p2.y - p1.y, p2.x - p1.x) * 180 / M_PI); + // the angle halfway between the angles + // perpendicular to a0->a1 (a10) and a1->a2 (a21) + double a2 = atan2(dy, dx); - double a10 = atan2(p1.y - p0.y, p1.x - p0.x); - double a21 = atan2(p2.y - p1.y, p2.x - p1.x); + out[(k + 1 - i) % (j - i - 1) + i].x = std::round(p1.x + buffer * cos(a2)); + out[(k + 1 - i) % (j - i - 1) + i].y = std::round(p1.y + buffer * sin(a2)); + } - double dx = cos(a10 - 90 * M_PI / 180) + cos(a21 - 90 * M_PI / 180); - double dy = sin(a10 - 90 * M_PI / 180) + sin(a21 - 90 * M_PI / 180); - double a2 = atan2(dy, dx); - - out[(k + 1 - i) % (j - i - 1) + i].x = std::round(p1.x + buffer * cos(a2)); - out[(k + 1 - i) % (j - i - 1) + i].y = std::round(p1.y + buffer * sin(a2)); + out[j - 1].x = out[i].x; + out[j - 1].y = out[i].y; + } else { + printf("skipping continuous ring %zu to %zu\n", i, j); } - out[j - 1].x = out[i].x; - out[j - 1].y = out[i].y; + i = j - 1; } } } diff --git a/geometry.hpp b/geometry.hpp index 11ee5d9f..53bc4609 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -104,6 +104,6 @@ std::string overzoom(mvt_tile tile, int oz, int ox, int oy, int nz, int nx, int std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int ny, int detail, int buffer, std::set const &keep, bool do_compress); -drawvec buffer_poly(drawvec const &geom, double buffer); +drawvec buffer_poly(drawvec const &geom, double buffer, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos); #endif diff --git a/main.cpp b/main.cpp index 2c303a99..4da0654b 100644 --- a/main.cpp +++ b/main.cpp @@ -3181,6 +3181,7 @@ int main(int argc, char **argv) { {"no-tiny-polygon-reduction", no_argument, &prevent[P_TINY_POLYGON_REDUCTION], 1}, {"no-tiny-polygon-reduction-at-maximum-zoom", no_argument, &prevent[P_TINY_POLYGON_REDUCTION_AT_MAXZOOM], 1}, {"tiny-polygon-size", required_argument, 0, '~'}, + {"buffer-polygons-outward", no_argument, &additional[A_BUFFER_POLYGONS_OUTWARD], 1}, {"no-simplification-of-shared-nodes", no_argument, &prevent[P_SIMPLIFY_SHARED_NODES], 1}, {"visvalingam", no_argument, &additional[A_VISVALINGAM], 1}, diff --git a/man/tippecanoe.1 b/man/tippecanoe.1 index a1abaaf2..73aa87a8 100644 --- a/man/tippecanoe.1 +++ b/man/tippecanoe.1 @@ -637,6 +637,8 @@ the line or polygon within one tile unit of its proper location. You can probabl .IP \(bu 2 \fB\fC\-\-tiny\-polygon\-size=\fR\fIsize\fP: Use the specified \fIsize\fP for tiny polygons instead of the default 2. Anything above 6 or so will lead to visible artifacts with the default tile detail. .IP \(bu 2 +\fB\fC\-aO\fR or \fB\fC\-\-buffer\-polygons\-outward\fR: Buffer polygons, except those that have been detected to be continous by \fB\fC\-\-no\-simplification\-of\-shared\-nodes\fR, outward slightly to prevent narrow polygon spindles fom collapsing away. +.IP \(bu 2 \fB\fC\-av\fR or \fB\fC\-\-visvalingam\fR: Use Visvalingam's simplification algorithm rather than Douglas\-Peucker's. .RE .SS Attempts to improve shared polygon boundaries diff --git a/options.hpp b/options.hpp index c56c5fcd..bcf1e9f5 100644 --- a/options.hpp +++ b/options.hpp @@ -26,6 +26,7 @@ #define A_HILBERT ((int) 'h') #define A_VISVALINGAM ((int) 'v') #define A_GENERATE_POLYGON_LABEL_POINTS ((int) 'P') +#define A_BUFFER_POLYGONS_OUTWARD ((int) 'O') #define P_SIMPLIFY ((int) 's') #define P_SIMPLIFY_LOW ((int) 'S') diff --git a/tile.cpp b/tile.cpp index f9622fcf..ed59f5c9 100644 --- a/tile.cpp +++ b/tile.cpp @@ -612,7 +612,10 @@ void *partial_feature_worker(void *v) { int out_detail = (*partials)[i].extra_detail; drawvec geom = (*partials)[i].geoms[0]; - drawvec before_scaling = geom; + + if (additional[A_BUFFER_POLYGONS_OUTWARD] && t == VT_POLYGON) { + geom = buffer_poly(geom, (1LL << (32 - z - out_detail)) * sqrt(2) / 2, z, (*partials)[i].tx, (*partials)[i].ty, a->shared_nodes_map, a->nodepos); + } to_tile_scale(geom, z, out_detail); @@ -627,19 +630,10 @@ void *partial_feature_worker(void *v) { check_polygon(geom); } - if (true || geom.size() < 4) { + if (geom.size() < 4) { if (area > 0) { - // Try reviving the polygon by buffering it outward a little - - geom = buffer_poly(before_scaling, (1LL << (32 - z - out_detail)) * 2.0); - to_tile_scale(geom, z, out_detail); - geom = clean_or_clip_poly(geom, 0, 0, false, true); - - if (geom.size() < 4) { - // OK, that didn't work, make a placeholder - // area is in world coordinates, calculated before scaling down - geom = revive_polygon(before, area, z, out_detail); - } + // area is in world coordinates, calculated before scaling down + geom = revive_polygon(before, area, z, out_detail); } else { geom.clear(); }