From 5d1c01d4b9dc76df93836e0420c50d17dca69bc0 Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Tue, 20 Feb 2024 22:08:59 -0800 Subject: [PATCH] Detect longitude wraparound only if it reduces self-intersections --- geometry.cpp | 46 ++++++++++++++++++++++++++ geometry.hpp | 3 ++ serial.cpp | 92 +++++++++++++++++++++++++++++++++++++--------------- 3 files changed, 114 insertions(+), 27 deletions(-) diff --git a/geometry.cpp b/geometry.cpp index 31a1c15af..19f293d4f 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -1397,3 +1397,49 @@ drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigne return out; } + +const std::pair SAME_SLOPE = std::make_pair(-INT_MAX, INT_MAX); + +// https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect +// +// beware of +// https://stackoverflow.com/questions/9043805/test-if-two-lines-intersect-javascript-function/16725715#16725715 +// which does not seem to produce correct results. +std::pair get_line_intersection(long long p0_x, long long p0_y, long long p1_x, long long p1_y, + long long p2_x, long long p2_y, long long p3_x, long long p3_y) { + // bounding box reject, x + long long min01x = std::min(p0_x, p1_x); + long long max01x = std::max(p0_x, p1_x); + long long min23x = std::min(p2_x, p3_x); + long long max23x = std::max(p2_x, p3_x); + if (max01x < min23x || max23x < min01x) { + return std::make_pair(-1, -1); + } + + // bounding box reject, y + long long min01y = std::min(p0_y, p1_y); + long long max01y = std::max(p0_y, p1_y); + long long min23y = std::min(p2_y, p3_y); + long long max23y = std::max(p2_y, p3_y); + if (max01y < min23y || max23y < min01y) { + return std::make_pair(-1, -1); + } + + long long d01_x, d01_y, d23_x, d23_y; + d01_x = p1_x - p0_x; + d01_y = p1_y - p0_y; + d23_x = p3_x - p2_x; + d23_y = p3_y - p2_y; + + long long det = (-d23_x * d01_y + d01_x * d23_y); + + if (det != 0) { + double t, s; + t = (d23_x * (p0_y - p2_y) - d23_y * (p0_x - p2_x)) / (double) det; + s = (-d01_y * (p0_x - p2_x) + d01_x * (p0_y - p2_y)) / (double) det; + + return std::make_pair(t, s); + } + + return SAME_SLOPE; +} diff --git a/geometry.hpp b/geometry.hpp index 850f7e798..80caeb171 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -114,4 +114,7 @@ std::string overzoom(const std::string &s, int oz, int ox, int oy, int nz, int n std::unordered_map const &attribute_accum, std::vector const &unidecode_data); +std::pair get_line_intersection(long long p0_x, long long p0_y, long long p1_x, long long p1_y, + long long p2_x, long long p2_y, long long p3_x, long long p3_y); + #endif diff --git a/serial.cpp b/serial.cpp index 83f114496..e82415d59 100644 --- a/serial.cpp +++ b/serial.cpp @@ -298,9 +298,71 @@ serial_feature deserialize_feature(std::string const &geoms, unsigned z, unsigne } static long long scale_geometry(struct serialization_state *sst, long long *bbox, drawvec &geom) { - long long offset = 0; - long long prev = 0; - bool has_prev = false; + if (additional[A_DETECT_WRAPAROUND]) { + drawvec adjusted = geom; + std::vector jumps; + + long long offset = 0; + long long prev = 0; + bool has_prev = false; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_LINETO) { + long long x = geom[i].x; + x += offset; + + if (has_prev) { + // jumps at least 180° in either direction + + if (x - prev > (1LL << 31)) { + offset -= 1LL << 32; + x -= 1LL << 32; + jumps.push_back(i); + } else if (prev - x > (1LL << 31)) { + offset += 1LL << 32; + x += 1LL << 32; + jumps.push_back(i); + } + } + + adjusted[i].x = x; + has_prev = true; + prev = x; + } else { + offset = 0; + prev = geom[i].x; + } + } + + // count whether making these adjustments resulted in + // more or less self-intersections than in the original + + size_t orig_bad = 0; + size_t adjusted_bad = 0; + + for (auto i : jumps) { + for (size_t j = 1; j < geom.size(); j++) { + if (geom[j].op == VT_LINETO && j != i) { + auto before = get_line_intersection(geom[i - 1].x, geom[i - 1].y, geom[i].x, geom[i].y, + geom[j - 1].x, geom[j - 1].y, geom[j].x, geom[j].y); + auto after = get_line_intersection(adjusted[i - 1].x, adjusted[i - 1].y, adjusted[i].x, adjusted[i].y, + adjusted[j - 1].x, adjusted[j - 1].y, adjusted[j].x, adjusted[j].y); + + if (before.first > 0 && before.first < 1) { + orig_bad++; + } + if (after.first > 0 && after.first < 1) { + adjusted_bad++; + } + } + } + } + + if (orig_bad > adjusted_bad) { + geom = std::move(adjusted); + } + } + double scale = 1.0 / (1 << geometry_scale); for (size_t i = 0; i < geom.size(); i++) { @@ -308,30 +370,6 @@ static long long scale_geometry(struct serialization_state *sst, long long *bbox long long x = geom[i].x; long long y = geom[i].y; - if (additional[A_DETECT_WRAPAROUND]) { - if (geom[i].op == VT_LINETO) { - x += offset; - if (has_prev) { - // jumps at least 180° but not exactly 360°, - // which in some data sets is an intentional - // line across the world - if (x - prev > (1LL << 31) && x - prev != (1LL << 32)) { - offset -= 1LL << 32; - x -= 1LL << 32; - } else if (prev - x > (1LL << 31) && prev - x != (1LL << 32)) { - offset += 1LL << 32; - x += 1LL << 32; - } - } - - has_prev = true; - prev = x; - } else { - offset = 0; - prev = x; - } - } - if (x < bbox[0]) { bbox[0] = x; }