From b3e02b75723b515e82f9ae17f71017fbd8fb7909 Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Mon, 17 Jul 2023 16:04:09 -0700 Subject: [PATCH] I'm not sure how this ever worked before --- drop.cpp | 70 ++++++++++++++++++++++++++++++++++++++++++-------------- drop.hpp | 28 +++++++++++++++++++---- main.cpp | 15 ------------ unit.cpp | 16 ++++++++++++- 4 files changed, 91 insertions(+), 38 deletions(-) diff --git a/drop.cpp b/drop.cpp index 36269098..37db0a11 100644 --- a/drop.cpp +++ b/drop.cpp @@ -1,36 +1,39 @@ +#include #include "drop.hpp" #include "options.hpp" #include "geometry.hpp" unsigned long long preserve_point_density_threshold = 0; -int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, double gamma) { +int calc_feature_minzoom(struct index *ix, struct drop_state ds[], int maxzoom, double gamma) { int feature_minzoom = 0; if (gamma >= 0 && (ix->t == VT_POINT || (additional[A_LINE_DROP] && ix->t == VT_LINE) || (additional[A_POLYGON_DROP] && ix->t == VT_POLYGON))) { - for (ssize_t i = maxzoom; i >= 0; i--) { - ds[i].seq++; + for (ssize_t i = 0; i <= maxzoom; i++) { + // This zoom level is now lighter on features than it should be. + ds[i].error -= 1.0 / ds[i].interval; + // printf("z%zd: error %f with interval %f\n", i, ds[i].error, ds[i].interval); } - ssize_t chosen = maxzoom + 1; - for (ssize_t i = maxzoom; i >= 0; i--) { - if (ds[i].seq < 0) { - feature_minzoom = i + 1; - // The feature we are pushing out - // appears in zooms i + 1 through maxzoom, - // so track where that was so we can make sure - // not to cluster something else that is *too* - // far away into it. - for (ssize_t j = i + 1; j <= maxzoom; j++) { + ssize_t chosen = maxzoom + 1; + for (ssize_t i = 0; i <= maxzoom; i++) { + if (ds[i].error < 0) { + // this zoom level is too light, so it is time to emit a feature. + feature_minzoom = i; + + // this feature now appears in this zoom level and all higher zoom levels, + // so each of them has this feature as its last feature, and each of them + // is now one feature heavier than before. + for (ssize_t j = i; j <= maxzoom; j++) { ds[j].previndex = ix->ix; + ds[j].error += ds[j].interval / ds[j].interval; + // printf("z%zd: now error %f\n", j, ds[j].error); } - chosen = i + 1; + chosen = i; break; - } else { - ds[i].seq -= ds[i].interval; } } @@ -44,8 +47,12 @@ int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, d if (ix->ix - ds[i].previndex > ((1LL << (32 - i)) / preserve_point_density_threshold) * ((1LL << (32 - i)) / preserve_point_density_threshold)) { feature_minzoom = i; - for (ssize_t j = i; j <= maxzoom; j++) { + // this feature now appears in this zoom level and all higher zoom levels below `chosen`, + // so each of them has this feature as its last feature, and each of them + // is now one feature heavier than before. + for (ssize_t j = i; j < chosen; j++) { ds[j].previndex = ix->ix; + ds[j].error += ds[j].interval / ds[j].interval; } break; @@ -56,3 +63,32 @@ int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, d return feature_minzoom; } + +void prep_drop_states(struct drop_state ds[], int maxzoom, int basezoom, double droprate) { + if (basezoom < 0) { + basezoom = maxzoom; + } + + // Needs to be signed for interval calculation + // printf("prep! max %d, base %d, rate %f\n", maxzoom, basezoom, droprate); + for (ssize_t i = 0; i <= maxzoom; i++) { + ds[i].previndex = 0; + ds[i].interval = 1; // every feature appears in every zoom level at or above the basezoom + + if (i < basezoom) { + // at zoom levels below the basezoom, the fraction of points that are dropped is + // the drop rate to the power of the number of zooms this zoom is below the basezoom + // + // for example: + // basezoom: 1 (droprate ^ 0) + // basezoom - 1: 2.5 (droprate ^ 1) + // basezoom - 2: 6.25 (droprate ^ 2) + // ... + // basezoom - n: (droprate ^ n) + ds[i].interval = std::exp(std::log(droprate) * (basezoom - i)); + // printf("%zd: interval %f\n", i, ds[i].interval); + } + + ds[i].error = 0; + } +} diff --git a/drop.hpp b/drop.hpp index bd13b443..0a1f5bb5 100644 --- a/drop.hpp +++ b/drop.hpp @@ -8,7 +8,7 @@ // Note that the fields are in a specific order so that `segment` and `t` will // packed together with `seq` so that the total structure size will be only 32 bytes // instead of 40. (Could we save a few more, perhaps, by tracking `len` instead of -// `end` and limiting the size of individual features to 32 bits?) +// `end` and limiting the size of individual features to 2^32 bytes?) struct index { // first and last+1 byte of the feature in the geometry temp file @@ -25,7 +25,7 @@ struct index { unsigned short t : 2; // sequence number (sometimes with gaps in numbering) of the feature in the original input file - unsigned long long seq : (64 - 18); // pack with segment and t to stay in 32 bytes + unsigned long long seq : (64 - 16 - 2); // pack with segment and t to stay in 32 bytes index() : t(0), @@ -33,15 +33,33 @@ struct index { } }; +// Each zoom level has a drop_state that is used to account for the fraction of +// point features that are supposed to be dropped in that zoom level. As it goes +// through the spatially-sorted features, it is basically doing a diffusion dither +// to keep the density of features in each vicinity at each zoom level +// approximately correct by including or excluding individual features +// to maintain the balance. + struct drop_state { - double gap; + // the z-index or hilbert index of the last feature that was placed in this zoom level unsigned long long previndex; + + // the preservation rate (1 or more) for features in this zoom level. + // 1 would be to keep all the features; 2 would drop every other feature; + // 4 every fourth feature, and so on. double interval; - double seq; // floating point because interval is + + // the current accumulated error in this zoom level: + // positive if too many features have been dropped; + // negative if not enough features have been dropped. + // + // this is floating-point because the interval is. + double error; }; extern unsigned long long preserve_point_density_threshold; -int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, double gamma); +int calc_feature_minzoom(struct index *ix, struct drop_state ds[], int maxzoom, double gamma); +void prep_drop_states(struct drop_state ds[], int maxzoom, int basezoom, double droprate); #endif diff --git a/main.cpp b/main.cpp index 61330583..ce648851 100644 --- a/main.cpp +++ b/main.cpp @@ -985,21 +985,6 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split } } -void prep_drop_states(struct drop_state *ds, int maxzoom, int basezoom, double droprate) { - // Needs to be signed for interval calculation - for (ssize_t i = 0; i <= maxzoom; i++) { - ds[i].gap = 0; - ds[i].previndex = 0; - ds[i].interval = 0; - - if (i < basezoom) { - ds[i].interval = std::exp(std::log(droprate) * (basezoom - i)); - } - - ds[i].seq = 0; - } -} - static size_t calc_memsize() { size_t mem; diff --git a/unit.cpp b/unit.cpp index d29d93e7..90789ae8 100644 --- a/unit.cpp +++ b/unit.cpp @@ -3,6 +3,8 @@ #include "text.hpp" #include "drop.hpp" +unsigned int additional[256] = {0}; + TEST_CASE("UTF-8 enforcement", "[utf8]") { REQUIRE(check_utf8("") == std::string("")); REQUIRE(check_utf8("hello world") == std::string("")); @@ -24,4 +26,16 @@ TEST_CASE("index structure packing", "[index]") { REQUIRE(sizeof(struct index) == 32); } -unsigned int additional[256] = {0}; +TEST_CASE("prep drop states", "[prep_drop_state]") { + struct drop_state ds[25]; + + prep_drop_states(ds, 24, 16, 2); + REQUIRE(ds[24].interval == 1); + REQUIRE(ds[17].interval == 1); + REQUIRE(ds[16].interval == 1); + REQUIRE(ds[15].interval == 2); + REQUIRE(ds[14].interval == 4); + + // to fix: because of floating point error this is not quite true + // REQUIRE(ds[0].interval == 65536); +}