Reduce exponent again

2 still leaves some tiles very dense;
3 reduces contrast perhaps too far

Makefile

@@ -1,6 +1,6 @@
 PREFIX ?= /usr/local
 
-all: tippecanoe enumerate decode
+all: tippecanoe enumerate decode simplify
 
 install: tippecanoe
 	mkdir -p $(PREFIX)/bin
@@ -25,6 +25,9 @@ enumerate: enumerate.o
 decode: decode.o vector_tile.pb.o projection.o
 	g++ $(PG) $(LIBS) -O3 -g -Wall -o $@ $^ -lm -lz -lprotobuf-lite -lsqlite3
 
+simplify: simplify.o projection.o
+	gcc $(PG) $(LIBS) -O3 -g -Wall -o $@ $^
+
 libjsonpull.a: jsonpull.o
 	ar rc $@ $^
 	ranlib $@

geojson.c

@@ -58,318 +58,6 @@ int mb_geometry[GEOM_TYPES] = {
 	VT_POLYGON,
 };
 
-struct feature_attribute {
-	int type; /* JSON_STRING, JSON_NUMBER, JSON_TRUE, JSON_FALSE, JSON_NULL */
-	char *key;
-	char *value;
-};
-
-struct feature_geometry {
-	int op; /* VT_MOVETO, VT_LINETO, VT_CLOSEPATH */
-	double lat;
-	double lon;
-
-	struct feature_geometry *next;
-};
-
-struct feature {
-	int type; /* VT_POINT, VT_LINE, VT_POLYGON */
-	int nattributes;
-	struct feature_attribute *attributes;
-	struct feature_geometry *geometries;
-};
-
-void geo_free(struct feature *geo) {
-	int i;
-	for (i = 0; i < geo->nattributes; i++) {
-		free(geo->attributes[i].key);
-		free(geo->attributes[i].value);
-	}
-	free(geo->attributes);
-	while (geo->geometries != NULL) {
-		struct feature_geometry *next = geo->geometries->next;
-		free(geo->geometries);
-		geo->geometries = next;
-	}
-	free(geo);
-}
-
-void geo_print(struct feature *geo) {
-	printf("{ \"type\": \"Feature\", \"properties\": {");
-	int i;
-	for (i = 0; i < geo->nattributes; i++) {
-		if (i != 0) {
-			printf(", ");
-		}
-		printf("\"%s\": \"%s\"", geo->attributes[i].key, geo->attributes[i].value);
-	}
-	printf("}, \"geometry\": { \"type\": \"%s\", \"coordinates\": [",
-		geo->type == VT_POLYGON ? "MultiPolygon" :
-		geo->type == VT_LINE ? "MultiLineString" :
-		geo->type == VT_POINT ? "MultiPoint" :
-		"???");
-
-	struct feature_geometry *g;
-	if (geo->type == VT_POINT) {
-		for (g = geo->geometries; g != NULL; g = g->next) {
-			printf(" [ %f, %f ]", g->lon, g->lat);
-			if (g->next != NULL) {
-				printf(",");
-			}
-		}
-	} else if (geo->type == VT_LINE) {
-		printf(" [");
-		for (g = geo->geometries; g != NULL; g = g->next) {
-			printf(" [ %f, %f ]", g->lon, g->lat);
-			if (g->next != NULL && g->next->op == VT_MOVETO) {
-				printf(" ], [");
-			} else if (g->next != NULL) {
-				printf(",");
-			}
-		}
-		printf(" ] ");
-	} else if (geo->type == VT_POLYGON) {
-		printf(" [ [");
-		for (g = geo->geometries; g != NULL; g = g->next) {
-			if (g->op == VT_CLOSEPATH) {
-				if (g->next != NULL) {
-					printf(" ], [");
-				}
-			} else {
-				printf(" [ %f, %f ]", g->lon, g->lat);
-				if (g->next != NULL && g->next->op == VT_MOVETO) {
-					printf(" ], [");
-				} else if (g->next != NULL && g->next->op != VT_CLOSEPATH) {
-					printf(",");
-				}
-			}
-		}
-
-		printf("] ] ");
-	}
-
-	printf("] } }\n");
-}
-
-struct feature_geometry **geo_parse(int t, json_object *j, int op, const char *fname, int line, struct feature_geometry **geom) {
-	if (j == NULL || j->type != JSON_ARRAY) {
-		fprintf(stderr, "%s:%d: expected array for type %d\n", fname, line, t);
-		return geom;
-	}
-
-	int within = geometry_within[t];
-	struct feature_geometry **began = geom;
-	if (within >= 0) {
-		int i;
-		for (i = 0; i < j->length; i++) {
-			if (within == GEOM_POINT) {
-				if (i == 0 || mb_geometry[t] == GEOM_MULTIPOINT) {
-					op = VT_MOVETO;
-				} else {
-					op = VT_LINETO;
-				}
-			}
-
-			geom = geo_parse(within, j->array[i], op, fname, line, geom);
-		}
-	} else {
-		if (j->length >= 2 && j->array[0]->type == JSON_NUMBER && j->array[1]->type == JSON_NUMBER) {
-			unsigned x, y;
-			double lon = j->array[0]->number;
-			double lat = j->array[1]->number;
-			latlon2tile(lat, lon, 32, &x, &y);
-
-			if (j->length > 2) {
-				static int warned = 0;
-
-				if (!warned) {
-					fprintf(stderr, "%s:%d: ignoring dimensions beyond two\n", fname, line);
-					warned = 1;
-				}
-			}
-
-
-			*geom = malloc(sizeof(struct feature_geometry));
-			(*geom)->op = op;
-			(*geom)->lat = lat;
-			(*geom)->lon = lon;
-			(*geom)->next = NULL;
-			geom = &((*geom)->next);
-		} else {
-			fprintf(stderr, "%s:%d: malformed point\n", fname, line);
-		}
-	}
-
-	if (t == GEOM_POLYGON) {
-		if (geom != began) {
-			*geom = malloc(sizeof(struct feature_geometry));
-			(*geom)->op = VT_CLOSEPATH;
-			(*geom)->lat = 0;
-			(*geom)->lon = 0;
-			(*geom)->next = NULL;
-			geom = &((*geom)->next);
-		}
-	}
-
-	return geom;
-}
-
-struct feature *parse_feature(json_object *geometry, json_object *properties, const char *fname, int line) {
-	json_object *geometry_type = json_hash_get(geometry, "type");
-	if (geometry_type == NULL) {
-		static int warned = 0;
-		if (!warned) {
-			fprintf(stderr, "%s:%d: null geometry (additional not reported)\n", fname, line);
-			warned = 1;
-		}
-
-		return NULL;
-	}
-
-	if (geometry_type->type != JSON_STRING) {
-		fprintf(stderr, "%s:%d: geometry without type\n", fname, line);
-		return NULL;
-	}
-
-	if (properties == NULL || (properties->type != JSON_HASH && properties->type != JSON_NULL)) {
-		fprintf(stderr, "%s:%d: feature without properties hash\n", fname, line);
-		return NULL;
-	}
-
-	json_object *coordinates = json_hash_get(geometry, "coordinates");
-	if (coordinates == NULL || coordinates->type != JSON_ARRAY) {
-		fprintf(stderr, "%s:%d: feature without coordinates array\n", fname, line);
-		return NULL;
-	}
-
-	int t;
-	for (t = 0; t < GEOM_TYPES; t++) {
-		if (strcmp(geometry_type->string, geometry_names[t]) == 0) {
-			break;
-		}
-	}
-	if (t >= GEOM_TYPES) {
-		fprintf(stderr, "%s:%d: Can't handle geometry type %s\n", fname, line, geometry_type->string);
-		return NULL;
-	}
-
-	struct feature *g = malloc(sizeof(struct feature));
-	g->type = mb_geometry[t];
-	g->attributes = NULL;
-	g->geometries = NULL;
-	g->nattributes = 0;
-
-	int nprop = 0;
-	if (properties->type == JSON_HASH) {
-		nprop = properties->length;
-	}
-
-	g->attributes = malloc(nprop * sizeof(struct feature_attribute));
-	int m = 0;
-
-	int i;
-	for (i = 0; i < nprop; i++) {
-		if (properties->keys[i]->type == JSON_STRING) {
-			g->attributes[m].key = strdup(properties->keys[i]->string);
-
-			if (properties->values[i] != NULL && properties->values[i]->type == JSON_STRING) {
-				g->attributes[m].type = VT_STRING;
-				g->attributes[m].value = strdup(properties->values[i]->string);
-				m++;
-			} else if (properties->values[i] != NULL && properties->values[i]->type == JSON_NUMBER) {
-				g->attributes[m].type = VT_NUMBER;
-				g->attributes[m].value = strdup(properties->values[i]->string);
-				m++;
-			} else if (properties->values[i] != NULL && (properties->values[i]->type == JSON_TRUE || properties->values[i]->type == JSON_FALSE)) {
-				g->attributes[m].type = VT_BOOLEAN;
-				g->attributes[m].value = strdup(properties->values[i]->string);
-				m++;
-			} else if (properties->values[i] != NULL && (properties->values[i]->type == JSON_NULL)) {
-				;
-			} else {
-				fprintf(stderr, "%s:%d: Unsupported property type for %s\n", fname, line, properties->keys[i]->string);
-				continue;
-			}
-		}
-	}
-
-	g->nattributes = m;
-	geo_parse(t, coordinates, VT_MOVETO, fname, line, &(g->geometries));
-
-	return g;
-}
-
-
-		/* XXX check for any non-features in the outer object */
-
-/*
-			if (bbox != NULL) {
-				if (x < bbox[0]) {
-					bbox[0] = x;
-				}
-				if (y < bbox[1]) {
-					bbox[1] = y;
-				}
-				if (x > bbox[2]) {
-					bbox[2] = x;
-				}
-				if (y > bbox[3]) {
-					bbox[3] = y;
-				}
-			}
-*/
-
-
-
-
-
-void parse_outer(FILE *f, char *fname) {
-	json_pull *jp = json_begin_file(f);
-	long long seq = 0;
-
-	while (1) {
-		json_object *j = json_read(jp);
-		if (j == NULL) {
-			if (jp->error != NULL) {
-				fprintf(stderr, "%s:%d: %s\n", fname, jp->line, jp->error);
-			}
-
-			json_free(jp->root);
-			break;
-		}
-
-		json_object *type = json_hash_get(j, "type");
-		json_object *geometry = json_hash_get(j, "geometry");
-		json_object *properties = json_hash_get(j, "properties");
-
-		if (type != NULL && type->type == JSON_STRING && strcmp(type->string, "Feature") == 0 &&
-		    geometry != NULL && geometry->type == JSON_HASH &&
-		    properties != NULL && properties->type == JSON_HASH) {
-			struct feature *g = parse_feature(geometry, properties, fname, jp->line);
-			json_free(j);
-			geo_print(g);
-			geo_free(g);
-
-			if (seq % 10000 == 0) {
-				fprintf(stderr, "Read %.2f million features\r", seq / 1000000.0);
-			}
-			seq++;
-		}
-	}
-
-	json_end(jp);
-}
-
-int main(int argc, char **argv) {
-	parse_outer(stdin, "standard input");
-}
-
-#define main xmain
-
-
-
-
-
 size_t fwrite_check(const void *ptr, size_t size, size_t nitems, FILE *stream, const char *fname, json_pull *source) {
 	size_t w = fwrite(ptr, size, nitems, stream);
 	if (w != nitems) {

simplify.c

@@ -0,0 +1,148 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <unistd.h>
+#include <sys/mman.h>
+#include "projection.h"
+
+int cmp(const void *v1, const void *v2) {
+	const unsigned long long *p1 = v1;
+	const unsigned long long *p2 = v2;
+
+	if (*p1 < *p2) {
+		return -1;
+	} else if (*p1 > *p2) {
+		return 1;
+	} else {
+		return 0;
+	}
+}
+
+inline double fpow(double val, double e) {
+	if (val == 0) {
+		return 0;
+	} else {
+		return exp(log(val) * e);
+	}
+}
+
+void out(unsigned long long index) {
+	unsigned x, y;
+	decode(index, &x, &y);
+
+	double lat, lon;
+	tile2latlon(x, y, 32, &lat, &lon);
+
+	printf("%.6f,%.6f", lat, lon);
+}
+
+// in other words, a z14 pixel, about 25 feet
+double scale = (double) (1LL << (64 - 2 * (14 + 8)));
+
+void test() {
+	srandomdev();
+	int i;
+	for (i = 0; i < 50000; i++) {
+		unsigned long long r1 = random() << 32;
+		r1 |= random();
+		unsigned long long r2 = r1 + scale;
+
+		unsigned x1, y1, x2, y2;
+		decode(r1, &x1, &y1);
+		decode(r2, &x2, &y2);
+
+		double lat1, lon1, lat2, lon2;
+		tile2latlon(x1, y1, 32, &lat1, &lon1);
+		tile2latlon(x2, y2, 32, &lat2, &lon2);
+
+		double rat = cos(lat1 + M_PI / 180);
+		double latd = lat1 - lat2;
+		double lond = (lon1 - lon2) * rat;
+		double d = sqrt(latd * latd + lond * lond) / .00000274;
+
+		if (d != 0) {
+			printf("%.6f\n", d);
+		}
+	}
+}
+
+int main() {
+	char s[2000];
+	long long seq = 0;
+	long long size = 0;
+
+	char geomname[2000] = "/tmp/geom.XXXXXXXX";
+	int fd = mkstemp(geomname);
+	if (fd < 0) {
+		perror(geomname);
+		exit(EXIT_FAILURE);
+	}
+	FILE *fp = fopen(geomname, "wb");
+	unlink(geomname);
+
+	while (fgets(s, 2000, stdin)) {
+		double lat, lon;
+		if (sscanf(s, "%lf,%lf", &lat, &lon) != 2) {
+			//fprintf(stderr, "Couldn't understand %s", s);
+			continue;
+		}
+
+		if (seq++ % 10000 == 0) {
+			fprintf(stderr, "Read %.2f million points\r", seq / 1000000.0);
+		}
+
+		unsigned x, y;
+		latlon2tile(lat, lon, 32, &x, &y);
+
+		unsigned long long enc = encode(x, y);
+		fwrite(&enc, 1, sizeof(unsigned long long), fp);
+		size++;
+	}
+
+	fclose(fp);
+
+	fprintf(stderr, "sorting %lld points\n", size);
+
+	unsigned long long *geom = mmap(NULL, size * sizeof(unsigned long long), PROT_READ | PROT_WRITE, MAP_PRIVATE, fd, 0);
+	if (geom == MAP_FAILED) {
+		perror("mmap");
+		exit(EXIT_FAILURE);
+	}
+
+	qsort(geom, size, sizeof(unsigned long long), cmp);
+
+	long long i;
+	double exponent = 2.3;
+
+	out(geom[0]);
+	printf("\n");
+
+	for (i = 1; i < size; i++) {
+		long long dist = geom[i] - geom[i - 1];
+
+		if (dist >= scale) {
+			out(geom[i]);
+			printf(" // far enough %lld vs %lf\n", dist, scale);
+		} else if (dist == 0) {
+			// eliminate exact duplicates.
+			// or should we produce them, but fewer?
+		} else {
+			double gap = dist / scale;
+			long long oi = i;
+			long long count = 0;
+
+			while (fpow((geom[i] - geom[oi - 1]) / scale, exponent) <= gap && i < size) {
+				i++;
+				count++;
+			}
+
+			if (i < size) {
+				if ((i > oi) && (seq++ % 2 == 0)) {
+					i--;
+				}
+				out(geom[i]);
+				printf(" // skipped %lld\n", (i - oi));
+			}
+		}
+	}
+}