Skip to content

Commit d4efb4d

Browse files
e-n-fe-n-f
committed
Add tiny polygon reduction / dust to overzoom
1 parent 830885f commit d4efb4d

9 files changed

+562
-133
lines changed

‎Makefile

+5
Original file line numberDiff line numberDiff line change
@@ -355,6 +355,11 @@ overzoom-test: tippecanoe-overzoom
355355
./tippecanoe-decode tests/pbf/12-2145-1391-filter2.pbf 12 2145 1391 > tests/pbf/12-2145-1391-filter2.pbf.json.check
356356
cmp tests/pbf/12-2145-1391-filter2.pbf.json.check tests/pbf/12-2145-1391-filter2.pbf.json
357357
rm tests/pbf/12-2145-1391-filter2.pbf.json.check tests/pbf/12-2145-1391-filter2.pbf
358+
# Tiny polygon reduction
359+
./tippecanoe-overzoom --line-simplification=5 --tiny-polygon-size=50 -o tests/pbf/countries-0-0-0.pbf.out tests/pbf/countries-0-0-0.pbf 0/0/0 0/0/0
360+
./tippecanoe-decode tests/pbf/countries-0-0-0.pbf.out 0 0 0 > tests/pbf/countries-0-0-0.pbf.out.json.check
361+
cmp tests/pbf/countries-0-0-0.pbf.out.json.check tests/pbf/countries-0-0-0.pbf.out.json
362+
rm tests/pbf/countries-0-0-0.pbf.out tests/pbf/countries-0-0-0.pbf.out.json.check
358363

359364
join-test: tippecanoe tippecanoe-decode tile-join
360365
./tippecanoe -q -f -z12 -o tests/join-population/tabblock_06001420.mbtiles -YALAND10:'Land area' -L'{"file": "tests/join-population/tabblock_06001420.json", "description": "population"}'

‎clip.cpp

+288-4
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
#include <stack>
12
#include <stdlib.h>
23
#include <mapbox/geometry/point.hpp>
34
#include <mapbox/geometry/multi_polygon.hpp>
@@ -755,10 +756,274 @@ static std::vector<std::pair<double, double>> clip_poly1(std::vector<std::pair<d
755756
return out;
756757
}
757758

759+
double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) {
760+
long long p2x = segB_x - segA_x;
761+
long long p2y = segB_y - segA_y;
762+
763+
// These calculations must be made in integers instead of floating point
764+
// to make them consistent between x86 and arm floating point implementations.
765+
//
766+
// Coordinates may be up to 34 bits, so their product is up to 68 bits,
767+
// making their sum up to 69 bits. Downshift before multiplying to keep them in range.
768+
double something = ((p2x / 4) * (p2x / 8) + (p2y / 4) * (p2y / 8)) * 32.0;
769+
// likewise
770+
double u = (0 == something) ? 0 : ((point_x - segA_x) / 4 * (p2x / 8) + (point_y - segA_y) / 4 * (p2y / 8)) * 32.0 / (something);
771+
772+
if (u >= 1) {
773+
u = 1;
774+
} else if (u <= 0) {
775+
u = 0;
776+
}
777+
778+
double x = segA_x + u * p2x;
779+
double y = segA_y + u * p2y;
780+
781+
double dx = x - point_x;
782+
double dy = y - point_y;
783+
784+
double out = std::round(sqrt(dx * dx + dy * dy) * 16.0) / 16.0;
785+
return out;
786+
}
787+
788+
// https://github.com/Project-OSRM/osrm-backend/blob/733d1384a40f/Algorithms/DouglasePeucker.cpp
789+
void douglas_peucker(drawvec &geom, int start, int n, double e, size_t kept, size_t retain, bool prevent_simplify_shared_nodes) {
790+
std::stack<int> recursion_stack;
791+
792+
if (!geom[start + 0].necessary || !geom[start + n - 1].necessary) {
793+
fprintf(stderr, "endpoints not marked necessary\n");
794+
exit(EXIT_IMPOSSIBLE);
795+
}
796+
797+
int prev = 0;
798+
for (int here = 1; here < n; here++) {
799+
if (geom[start + here].necessary) {
800+
recursion_stack.push(prev);
801+
recursion_stack.push(here);
802+
prev = here;
803+
804+
if (prevent_simplify_shared_nodes) {
805+
if (retain > 0) {
806+
retain--;
807+
}
808+
}
809+
}
810+
}
811+
// These segments are put on the stack from start to end,
812+
// independent of winding, so note that anything that uses
813+
// "retain" to force it to keep at least N points will
814+
// keep a different set of points when wound one way than
815+
// when wound the other way.
816+
817+
while (!recursion_stack.empty()) {
818+
// pop next element
819+
int second = recursion_stack.top();
820+
recursion_stack.pop();
821+
int first = recursion_stack.top();
822+
recursion_stack.pop();
823+
824+
double max_distance = -1;
825+
int farthest_element_index;
826+
827+
// find index idx of element with max_distance
828+
int i;
829+
if (geom[start + first] < geom[start + second]) {
830+
farthest_element_index = first;
831+
for (i = first + 1; i < second; i++) {
832+
double temp_dist = distance_from_line(geom[start + i].x, geom[start + i].y, geom[start + first].x, geom[start + first].y, geom[start + second].x, geom[start + second].y);
833+
834+
double distance = std::fabs(temp_dist);
835+
836+
if ((distance > e || kept < retain) && (distance > max_distance || (distance == max_distance && geom[start + i] < geom[start + farthest_element_index]))) {
837+
farthest_element_index = i;
838+
max_distance = distance;
839+
}
840+
}
841+
} else {
842+
farthest_element_index = second;
843+
for (i = second - 1; i > first; i--) {
844+
double temp_dist = distance_from_line(geom[start + i].x, geom[start + i].y, geom[start + second].x, geom[start + second].y, geom[start + first].x, geom[start + first].y);
845+
846+
double distance = std::fabs(temp_dist);
847+
848+
if ((distance > e || kept < retain) && (distance > max_distance || (distance == max_distance && geom[start + i] < geom[start + farthest_element_index]))) {
849+
farthest_element_index = i;
850+
max_distance = distance;
851+
}
852+
}
853+
}
854+
855+
if (max_distance >= 0) {
856+
// mark idx as necessary
857+
geom[start + farthest_element_index].necessary = 1;
858+
kept++;
859+
860+
if (geom[start + first] < geom[start + second]) {
861+
if (1 < farthest_element_index - first) {
862+
recursion_stack.push(first);
863+
recursion_stack.push(farthest_element_index);
864+
}
865+
if (1 < second - farthest_element_index) {
866+
recursion_stack.push(farthest_element_index);
867+
recursion_stack.push(second);
868+
}
869+
} else {
870+
if (1 < second - farthest_element_index) {
871+
recursion_stack.push(farthest_element_index);
872+
recursion_stack.push(second);
873+
}
874+
if (1 < farthest_element_index - first) {
875+
recursion_stack.push(first);
876+
recursion_stack.push(farthest_element_index);
877+
}
878+
}
879+
}
880+
}
881+
}
882+
883+
// cut-down version of simplify_lines(), not dealing with shared node preservation
884+
static drawvec simplify_lines_basic(drawvec &geom, int z, int detail, double simplification, size_t retain) {
885+
int res = 1 << (32 - detail - z);
886+
887+
for (size_t i = 0; i < geom.size(); i++) {
888+
if (geom[i].op == VT_MOVETO) {
889+
geom[i].necessary = 1;
890+
} else if (geom[i].op == VT_LINETO) {
891+
geom[i].necessary = 0;
892+
// if this is actually the endpoint, not an intermediate point,
893+
// it will be marked as necessary below
894+
} else {
895+
geom[i].necessary = 1;
896+
}
897+
}
898+
899+
for (size_t i = 0; i < geom.size(); i++) {
900+
if (geom[i].op == VT_MOVETO) {
901+
size_t j;
902+
for (j = i + 1; j < geom.size(); j++) {
903+
if (geom[j].op != VT_LINETO) {
904+
break;
905+
}
906+
}
907+
908+
geom[i].necessary = 1;
909+
geom[j - 1].necessary = 1;
910+
911+
if (j - i > 1) {
912+
douglas_peucker(geom, i, j - i, res * simplification, 2, retain, false);
913+
}
914+
i = j - 1;
915+
}
916+
}
917+
918+
size_t out = 0;
919+
for (size_t i = 0; i < geom.size(); i++) {
920+
if (geom[i].necessary) {
921+
geom[out++] = geom[i];
922+
}
923+
}
924+
geom.resize(out);
925+
return geom;
926+
}
927+
928+
drawvec reduce_tiny_poly(drawvec const &geom, int z, int detail, bool *still_needs_simplification, bool *reduced_away, double *accum_area, double tiny_polygon_size) {
929+
drawvec out;
930+
const double pixel = (1LL << (32 - detail - z)) * (double) tiny_polygon_size;
931+
932+
bool included_last_outer = false;
933+
*still_needs_simplification = false;
934+
*reduced_away = false;
935+
936+
for (size_t i = 0; i < geom.size(); i++) {
937+
if (geom[i].op == VT_MOVETO) {
938+
size_t j;
939+
for (j = i + 1; j < geom.size(); j++) {
940+
if (geom[j].op != VT_LINETO) {
941+
break;
942+
}
943+
}
944+
945+
double area = get_area(geom, i, j);
946+
947+
// XXX There is an ambiguity here: If the area of a ring is 0 and it is followed by holes,
948+
// we don't know whether the area-0 ring was a hole too or whether it was the outer ring
949+
// that these subsequent holes are somehow being subtracted from. I hope that if a polygon
950+
// was simplified down to nothing, its holes also became nothing.
951+
952+
if (area != 0) {
953+
// These are pixel coordinates, so area > 0 for the outer ring.
954+
// If the outer ring of a polygon was reduced to a pixel, its
955+
// inner rings must just have their area de-accumulated rather
956+
// than being drawn since we don't really know where they are.
957+
958+
// i.e., this outer ring is small enough that we are including it
959+
// in a tiny polygon rather than letting it represent itself,
960+
// OR it is an inner ring and we haven't output an outer ring for it to be
961+
// cut out of, so we are just subtracting its area from the tiny polygon
962+
// rather than trying to deal with it geometrically
963+
if ((area > 0 && area <= pixel * pixel) || (area < 0 && !included_last_outer)) {
964+
*accum_area += area;
965+
*reduced_away = true;
966+
967+
if (area > 0 && *accum_area > pixel * pixel) {
968+
// XXX use centroid;
969+
970+
out.emplace_back(VT_MOVETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2);
971+
out.emplace_back(VT_LINETO, geom[i].x - pixel / 2 + pixel, geom[i].y - pixel / 2);
972+
out.emplace_back(VT_LINETO, geom[i].x - pixel / 2 + pixel, geom[i].y - pixel / 2 + pixel);
973+
out.emplace_back(VT_LINETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2 + pixel);
974+
out.emplace_back(VT_LINETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2);
975+
976+
*accum_area -= pixel * pixel;
977+
}
978+
979+
if (area > 0) {
980+
included_last_outer = false;
981+
}
982+
}
983+
// i.e., this ring is large enough that it gets to represent itself
984+
// or it is a tiny hole out of a real polygon, which we are still treating
985+
// as a real geometry because otherwise we can accumulate enough tiny holes
986+
// that we will drop the next several outer rings getting back up to 0.
987+
else {
988+
for (size_t k = i; k < j && k < geom.size(); k++) {
989+
out.push_back(geom[k]);
990+
}
991+
992+
// which means that the overall polygon has a real geometry,
993+
// which means that it gets to be simplified.
994+
*still_needs_simplification = true;
995+
996+
if (area > 0) {
997+
included_last_outer = true;
998+
}
999+
}
1000+
} else {
1001+
// area is 0: doesn't count as either having been reduced away,
1002+
// since it was probably just degenerate from having been clipped,
1003+
// or as needing simplification, since it produces no output.
1004+
}
1005+
1006+
i = j - 1;
1007+
} else {
1008+
fprintf(stderr, "how did we get here with %d in %d?\n", geom[i].op, (int) geom.size());
1009+
1010+
for (size_t n = 0; n < geom.size(); n++) {
1011+
fprintf(stderr, "%d/%lld/%lld ", geom[n].op, geom[n].x, geom[n].y);
1012+
}
1013+
fprintf(stderr, "\n");
1014+
1015+
out.push_back(geom[i]);
1016+
}
1017+
}
1018+
1019+
return out;
1020+
}
1021+
7581022
std::string overzoom(std::vector<input_tile> const &tiles, int nz, int nx, int ny,
7591023
int detail, int buffer, std::set<std::string> const &keep, bool do_compress,
7601024
std::vector<std::pair<unsigned, unsigned>> *next_overzoomed_tiles,
761-
bool demultiply, json_object *filter, bool preserve_input_order, std::unordered_map<std::string, attribute_op> const &attribute_accum, std::vector<std::string> const &unidecode_data) {
1025+
bool demultiply, json_object *filter, bool preserve_input_order, std::unordered_map<std::string, attribute_op> const &attribute_accum, std::vector<std::string> const &unidecode_data, double simplification,
1026+
double tiny_polygon_size) {
7621027
std::vector<source_tile> decoded;
7631028

7641029
for (auto const &t : tiles) {
@@ -784,7 +1049,7 @@ std::string overzoom(std::vector<input_tile> const &tiles, int nz, int nx, int n
7841049
decoded.push_back(out);
7851050
}
7861051

787-
return overzoom(decoded, nz, nx, ny, detail, buffer, keep, do_compress, next_overzoomed_tiles, demultiply, filter, preserve_input_order, attribute_accum, unidecode_data);
1052+
return overzoom(decoded, nz, nx, ny, detail, buffer, keep, do_compress, next_overzoomed_tiles, demultiply, filter, preserve_input_order, attribute_accum, unidecode_data, simplification, tiny_polygon_size);
7881053
}
7891054

7901055
struct tile_feature {
@@ -885,7 +1150,8 @@ static struct preservecmp {
8851150
std::string overzoom(std::vector<source_tile> const &tiles, int nz, int nx, int ny,
8861151
int detail, int buffer, std::set<std::string> const &keep, bool do_compress,
8871152
std::vector<std::pair<unsigned, unsigned>> *next_overzoomed_tiles,
888-
bool demultiply, json_object *filter, bool preserve_input_order, std::unordered_map<std::string, attribute_op> const &attribute_accum, std::vector<std::string> const &unidecode_data) {
1153+
bool demultiply, json_object *filter, bool preserve_input_order, std::unordered_map<std::string, attribute_op> const &attribute_accum, std::vector<std::string> const &unidecode_data, double simplification,
1154+
double tiny_polygon_size) {
8891155
mvt_tile outtile;
8901156
std::shared_ptr<std::string> tile_stringpool = std::make_shared<std::string>();
8911157

@@ -916,6 +1182,7 @@ std::string overzoom(std::vector<source_tile> const &tiles, int nz, int nx, int
9161182
}
9171183

9181184
std::vector<tile_feature> pending_tile_features;
1185+
double accum_area = 0;
9191186

9201187
static const std::string retain_points_multiplier_first = "tippecanoe:retain_points_multiplier_first";
9211188
static const std::string retain_points_multiplier_sequence = "tippecanoe:retain_points_multiplier_sequence";
@@ -1014,6 +1281,23 @@ std::string overzoom(std::vector<source_tile> const &tiles, int nz, int nx, int
10141281
}
10151282
}
10161283

1284+
bool still_need_simplification_after_reduction = false;
1285+
if (t == VT_POLYGON && tiny_polygon_size > 0) {
1286+
bool simplified_away_by_reduction = false;
1287+
1288+
geom = reduce_tiny_poly(geom, nz, detail, &still_need_simplification_after_reduction, &simplified_away_by_reduction, &accum_area, tiny_polygon_size);
1289+
} else {
1290+
still_need_simplification_after_reduction = true;
1291+
}
1292+
1293+
if (simplification > 0 && still_need_simplification_after_reduction) {
1294+
if (t == VT_POLYGON) {
1295+
geom = simplify_lines_basic(geom, nz, detail, simplification, 4);
1296+
} else if (t == VT_LINE) {
1297+
geom = simplify_lines_basic(geom, nz, detail, simplification, 0);
1298+
}
1299+
}
1300+
10171301
// Scale to output tile extent
10181302

10191303
to_tile_scale(geom, nz, det);
@@ -1077,7 +1361,7 @@ std::string overzoom(std::vector<source_tile> const &tiles, int nz, int nx, int
10771361
std::string child = overzoom(sts,
10781362
nz + 1, nx * 2 + x, ny * 2 + y,
10791363
detail, buffer, keep, false, NULL,
1080-
demultiply, filter, preserve_input_order, attribute_accum, unidecode_data);
1364+
demultiply, filter, preserve_input_order, attribute_accum, unidecode_data, simplification, tiny_polygon_size);
10811365
if (child.size() > 0) {
10821366
next_overzoomed_tiles->emplace_back(nx * 2 + x, ny * 2 + y);
10831367
}

0 commit comments

Comments
 (0)