From a12bea7129081a7642cfd3529744e64a2aca348d Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Thu, 14 Mar 2024 14:36:51 -0700 Subject: [PATCH] Straighten out hilbert and quadkey encoding --- clip.cpp | 1 + geometry.cpp | 21 ++++++++++++--------- geometry.hpp | 2 -- main.cpp | 11 +++++++---- projection.cpp | 8 ++++++++ projection.hpp | 5 +++++ serial.cpp | 7 ++++--- 7 files changed, 37 insertions(+), 18 deletions(-) diff --git a/clip.cpp b/clip.cpp index 83d9b947..4b5c5fb1 100644 --- a/clip.cpp +++ b/clip.cpp @@ -10,6 +10,7 @@ #include "evaluator.hpp" #include "serial.hpp" #include "attribute.hpp" +#include "projection.hpp" static std::vector> clip_poly1(std::vector> &geom, long long minx, long long miny, long long maxx, long long maxy, diff --git a/geometry.cpp b/geometry.cpp index ac063c5c..f8e4903c 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -497,7 +497,7 @@ drawvec simplify_lines(drawvec &geom, int z, int tx, int ty, int detail, bool ma // to quadkey struct node n; - n.index = encode_quadkey((unsigned) d.x, (unsigned) d.y); + n.index = encode_quadkey(coordinate_to_encodable(d.x), coordinate_to_encodable(d.y)); if (bsearch(&n, shared_nodes_map, nodepos / sizeof(node), sizeof(node), nodecmp) != NULL) { geom[i].necessary = true; @@ -578,8 +578,8 @@ drawvec reorder_lines(const drawvec &geom) { // instead of down and to the right // so that it will coalesce better - unsigned long long l1 = encode_index(geom[0].x, geom[0].y); - unsigned long long l2 = encode_index(geom[geom.size() - 1].x, geom[geom.size() - 1].y); + unsigned long long l1 = encode_index(coordinate_to_encodable(geom[0].x), coordinate_to_encodable(geom[0].y)); + unsigned long long l2 = encode_index(coordinate_to_encodable(geom[geom.size() - 1].x), coordinate_to_encodable(geom[geom.size() - 1].y)); if (l1 > l2) { drawvec out; @@ -1302,6 +1302,9 @@ drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigne unsigned wx, wy; decode_index(label_point, &wx, &wy); + long long wwx = decoded_to_coordinate(wx); + long long wwy = decoded_to_coordinate(wy); + // upper left of tile in world coordinates long long tx1 = 0, ty1 = 0; // lower right of tile in world coordinates; @@ -1340,16 +1343,16 @@ drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigne const long long label_spacing = spiral_dist * (tx2 - tx1); - long long x1 = floor(std::min(bx1 - wx, bx2 - wx) / label_spacing); - long long x2 = ceil(std::max(bx1 - wx, bx2 - wx) / label_spacing); + long long x1 = floor(std::min(bx1 - wwx, bx2 - wwx) / label_spacing); + long long x2 = ceil(std::max(bx1 - wwx, bx2 - wwx) / label_spacing); - long long y1 = floor(std::min(by1 - wy, by2 - wy) / label_spacing - 0.5); - long long y2 = ceil(std::max(by1 - wy, by2 - wy) / label_spacing); + long long y1 = floor(std::min(by1 - wwy, by2 - wwy) / label_spacing - 0.5); + long long y2 = ceil(std::max(by1 - wwy, by2 - wwy) / label_spacing); for (long long lx = x1; lx <= x2; lx++) { for (long long ly = y1; ly <= y2; ly++) { - long long x = lx * label_spacing + wx; - long long y = ly * label_spacing + wy; + long long x = lx * label_spacing + wwx; + long long y = ly * label_spacing + wwy; if (((unsigned long long) lx & 1) == 1) { y += label_spacing / 2; diff --git a/geometry.hpp b/geometry.hpp index cdc455e0..f847afd1 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -20,8 +20,6 @@ #define VT_LINETO 2 #define VT_CLOSEPATH 7 -#define GLOBAL_DETAIL 32 - // The bitfield is to make sizeof(draw) be 16 instead of 24 // at the cost, apparently, of a 0.7% increase in running time // for packing and unpacking. diff --git a/main.cpp b/main.cpp index fdac4777..72c00780 100644 --- a/main.cpp +++ b/main.cpp @@ -2069,7 +2069,7 @@ std::pair read_input(std::vector &sources, char *fname, i #endif struct node n; - n.index = encode_quadkey((unsigned) x, (unsigned) y); + n.index = encode_quadkey(coordinate_to_encodable(x), coordinate_to_encodable(y)); fwrite_check((char *) &n, sizeof(struct node), 1, readers[0].nodefile, &readers[0].nodepos, "vertices"); } @@ -2152,7 +2152,7 @@ std::pair read_input(std::vector &sources, char *fname, i unsigned wx, wy; decode_quadkey(here.index, &wx, &wy); double lon, lat; - tile2lonlat(wx, wy, 32, &lon, &lat); + tile2lonlat(decoded_to_coordinate(wx), decoded_to_coordinate(wy), GLOBAL_DETAIL, &lon, &lat); printf("{\"type\":\"Feature\", \"properties\":{}, \"geometry\":{\"type\":\"Point\", \"coordinates\":[%f,%f]}}\n", lon, lat); #endif } @@ -2505,6 +2505,9 @@ std::pair read_input(std::vector &sources, char *fname, i unsigned xx, yy; decode_index(map[ip].ix, &xx, &yy); + long long gxx = decoded_to_coordinate(xx); + long long gyy = decoded_to_coordinate(yy); + long long nprogress = 100 * ip / indices; if (nprogress != progress) { progress = nprogress; @@ -2520,8 +2523,8 @@ std::pair read_input(std::vector &sources, char *fname, i if (z != 0) { // These are tile numbers, not pixels, // so shift, not round - xxx = xx >> (32 - z); - yyy = yy >> (32 - z); + xxx = gxx >> (GLOBAL_DETAIL - z); + yyy = gyy >> (GLOBAL_DETAIL - z); } double scale = (double) (1LL << (64 - 2 * (z + 8))); diff --git a/projection.cpp b/projection.cpp index f5f71382..eb92b6df 100644 --- a/projection.cpp +++ b/projection.cpp @@ -200,6 +200,14 @@ void decode_quadkey(unsigned long long index, unsigned *wx, unsigned *wy) { } } +unsigned coordinate_to_encodable(long long coord) { + return (unsigned) (coord / (1LL << (GLOBAL_DETAIL - 32))); +} + +long long decoded_to_coordinate(unsigned coord) { + return ((long long) coord) * (1LL << (GLOBAL_DETAIL - 32)); +} + void set_projection_or_exit(const char *optarg) { struct projection *p; for (p = projections; p->name != NULL; p++) { diff --git a/projection.hpp b/projection.hpp index d649ef56..789849ec 100644 --- a/projection.hpp +++ b/projection.hpp @@ -1,6 +1,8 @@ #ifndef PROJECTION_HPP #define PROJECTION_HPP +#define GLOBAL_DETAIL 32 + void lonlat2tile(double lon, double lat, int zoom, long long *x, long long *y); void epsg3857totile(double ix, double iy, int zoom, long long *x, long long *y); void tile2lonlat(long long x, long long y, int zoom, double *lon, double *lat); @@ -26,4 +28,7 @@ void decode_quadkey(unsigned long long index, unsigned *wx, unsigned *wy); unsigned long long encode_hilbert(unsigned int wx, unsigned int wy); void decode_hilbert(unsigned long long index, unsigned *wx, unsigned *wy); +unsigned coordinate_to_encodable(long long coord); +long long decoded_to_coordinate(unsigned coord); + #endif diff --git a/serial.cpp b/serial.cpp index a5c350cc..0bad2db9 100644 --- a/serial.cpp +++ b/serial.cpp @@ -409,7 +409,7 @@ static void add_scaled_node(struct reader *r, serialization_state *sst, draw g) long long y = SHIFT_LEFT(g.y); struct node n; - n.index = encode_quadkey((unsigned) x, (unsigned) y); + n.index = encode_quadkey(coordinate_to_encodable(x), coordinate_to_encodable(y)); fwrite_check((char *) &n, sizeof(struct node), 1, r->nodefile, &r->nodepos, sst->fname); } @@ -704,14 +704,15 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf, std:: midy = SHIFT_LEFT(scaled_geometry[ix].y) & ((1LL << GLOBAL_DETAIL) - 1); } - bbox_index = encode_index(midx, midy); + bbox_index = encode_index(coordinate_to_encodable(midx), coordinate_to_encodable(midy)); if (sf.t == VT_POLYGON && additional[A_GENERATE_POLYGON_LABEL_POINTS]) { drawvec dv = polygon_to_anchor(scaled_geometry); if (dv.size() > 0) { dv[0].x = SHIFT_LEFT(dv[0].x) & ((1LL << GLOBAL_DETAIL) - 1); dv[0].y = SHIFT_LEFT(dv[0].y) & ((1LL << GLOBAL_DETAIL) - 1); - sf.label_point = encode_index(dv[0].x, dv[0].y); + // this could just be serialized as numbers instead of encoding and decoding + sf.label_point = encode_index(coordinate_to_encodable(dv[0].x), coordinate_to_encodable(dv[0].y)); } }