Round coordinates instead of truncating them (#60)

* Round coordinates instead of truncating them

* Update all the tests for coordinate rounding changes

* Curses, integer division still truncates

* Fix tests

* Also round instead of shifting when looking for no-op linetos

* Should I worry that the same change for moveto doesn't change any tests?

* Also round instead of shifting when scaling down to maxzoom resolution

* Replace another explicit shift, for origin point

* Round instead of shift when writing clipped geometries to the next zoom

* Fix low-zoom gridding and smaller-than-a-pixel checks

* Don't guess an excessively large maxzoom when there is only one feature

* Add a test for guessing the maxzoom of a single point

* Explicitly sort by index if no other order distinguishes features

* Another affected test

CHANGELOG.md

@@ -1,3 +1,7 @@
+## 2.20.0
+
+* Round coordinates instead of truncating them, for better precision when overzooming
+
 ## 2.19.0
 
 * Don't guess an excessively large maxzoom when there is only one feature

Makefile

@@ -238,7 +238,7 @@ pbf-test:
 enumerate-test:
 	./tippecanoe -q -z5 -f -o tests/ne_110m_admin_0_countries/out/enum.mbtiles tests/ne_110m_admin_0_countries/in.json.gz
 	./tippecanoe-enumerate tests/ne_110m_admin_0_countries/out/enum.mbtiles > tests/ne_110m_admin_0_countries/out/enum.check
-	cmp tests/ne_110m_admin_0_countries/out/enum tests/ne_110m_admin_0_countries/out/enum.check
+	cmp tests/ne_110m_admin_0_countries/out/enum.check tests/ne_110m_admin_0_countries/out/enum
 	rm tests/ne_110m_admin_0_countries/out/enum.mbtiles tests/ne_110m_admin_0_countries/out/enum.check
 
 join-test: tile-join

geometry.cpp

@@ -87,17 +87,17 @@ drawvec decode_geometry(FILE *meta, std::atomic<long long> *geompos, int z, unsi
 
 void to_tile_scale(drawvec &geom, int z, int detail) {
 	for (size_t i = 0; i < geom.size(); i++) {
-		geom[i].x >>= (32 - detail - z);
-		geom[i].y >>= (32 - detail - z);
+		geom[i].x = std::round((double) geom[i].x / (1LL << (32 - detail - z)));
+		geom[i].y = std::round((double) geom[i].y / (1LL << (32 - detail - z)));
 	}
 }
 
 drawvec from_tile_scale(drawvec const &geom, int z, int detail) {
 	drawvec out;
 	for (size_t i = 0; i < geom.size(); i++) {
 		draw d = geom[i];
-		d.x <<= (32 - detail - z);
-		d.y <<= (32 - detail - z);
+		d.x *= (1LL << (32 - detail - z));
+		d.y *= (1LL << (32 - detail - z));
 		out.push_back(d);
 	}
 	return out;
@@ -110,16 +110,16 @@ drawvec remove_noop(drawvec geom, int type, int shift) {
 	drawvec out;
 
 	for (size_t i = 0; i < geom.size(); i++) {
-		if (geom[i].op == VT_LINETO && (geom[i].x >> shift) == x && (geom[i].y >> shift) == y) {
+		if (geom[i].op == VT_LINETO && std::round((double) geom[i].x / (1LL << shift)) == x && std::round((double) geom[i].y / (1LL << shift)) == y) {
 			continue;
 		}
 
 		if (geom[i].op == VT_CLOSEPATH) {
 			out.push_back(geom[i]);
 		} else { /* moveto or lineto */
 			out.push_back(geom[i]);
-			x = geom[i].x >> shift;
-			y = geom[i].y >> shift;
+			x = std::round((double) geom[i].x / (1LL << shift));
+			y = std::round((double) geom[i].y / (1LL << shift));
 		}
 	}
 
@@ -158,7 +158,7 @@ drawvec remove_noop(drawvec geom, int type, int shift) {
 
 		for (size_t i = 0; i < geom.size(); i++) {
 			if (geom[i].op == VT_MOVETO) {
-				if (i > 0 && geom[i - 1].op == VT_LINETO && (geom[i - 1].x >> shift) == (geom[i].x >> shift) && (geom[i - 1].y >> shift) == (geom[i].y >> shift)) {
+				if (i > 0 && geom[i - 1].op == VT_LINETO && std::round((double) geom[i - 1].x / (1LL << shift)) == std::round((double) geom[i].x / (1LL << shift)) && std::round((double) geom[i - 1].y / (1LL << shift)) == std::round((double) geom[i].y / (1LL << shift))) {
 					continue;
 				}
 			}
@@ -1323,8 +1323,8 @@ drawvec stairstep(drawvec &geom, int z, int detail) {
 	double scale = 1 << (32 - detail - z);
 
 	for (size_t i = 0; i < geom.size(); i++) {
-		geom[i].x = std::floor(geom[i].x / scale);
-		geom[i].y = std::floor(geom[i].y / scale);
+		geom[i].x = std::round(geom[i].x / scale);
+		geom[i].y = std::round(geom[i].y / scale);
 	}
 
 	for (size_t i = 0; i < geom.size(); i++) {

main.cpp

@@ -2276,6 +2276,8 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 			for (z = 0; z <= MAX_ZOOM; z++) {
 				unsigned xxx = 0, yyy = 0;
 				if (z != 0) {
+					// These are tile numbers, not pixels,
+					// so shift, not round
 					xxx = xx >> (32 - z);
 					yyy = yy >> (32 - z);
 				}

projection.cpp

@@ -2,6 +2,7 @@
 #include <string.h>
 #include <stdlib.h>
 #include <math.h>
+#include <cmath>
 #include <atomic>
 #include "projection.hpp"
 #include "errors.hpp"
@@ -87,8 +88,8 @@ void epsg3857totile(double ix, double iy, int zoom, long long *x, long long *y)
 	*y = ((1LL << 32) - 1) - (iy * (1LL << 31) / 6378137.0 / M_PI + (1LL << 31));
 
 	if (zoom != 0) {
-		*x >>= (32 - zoom);
-		*y >>= (32 - zoom);
+		*x = std::round((double) *x / (1LL << (32 - zoom)));
+		*y = std::round((double) *y / (1LL << (32 - zoom)));
 	}
 }
 

serial.cpp

@@ -24,7 +24,7 @@
 
 // Offset coordinates to keep them positive
 #define COORD_OFFSET (4LL << 32)
-#define SHIFT_RIGHT(a) ((((a) + COORD_OFFSET) >> geometry_scale) - (COORD_OFFSET >> geometry_scale))
+#define SHIFT_RIGHT(a) ((long long) std::round((double)(a) / (1LL << geometry_scale)))
 #define SHIFT_LEFT(a) ((((a) + (COORD_OFFSET >> geometry_scale)) << geometry_scale) - COORD_OFFSET)
 
 size_t fwrite_check(const void *ptr, size_t size, size_t nitems, FILE *stream, const char *fname) {
@@ -383,8 +383,8 @@ static long long scale_geometry(struct serialization_state *sst, long long *bbox
 					*(sst->initial_x) = 1LL << 31;
 					*(sst->initial_y) = 1LL << 31;
 				} else {
-					*(sst->initial_x) = (((x + COORD_OFFSET) >> geometry_scale) << geometry_scale) - COORD_OFFSET;
-					*(sst->initial_y) = (((y + COORD_OFFSET) >> geometry_scale) << geometry_scale) - COORD_OFFSET;
+					*(sst->initial_x) = SHIFT_LEFT(SHIFT_RIGHT(x));
+					*(sst->initial_y) = SHIFT_LEFT(SHIFT_RIGHT(y));
 				}
 
 				*(sst->initialized) = 1;