Merge pull request #141 from mapbox/multithread-low-zooms

Multithread line simplification and polygon cleaning at low zooms

CHANGELOG.md

@@ -1,3 +1,7 @@
+## 1.6.1
+
+* Use multiple threads for line simplification and polygon cleaning
+
 ## 1.6.0
 
 * Add option of parallelized input when reading from a line-delimited file

geometry.cc

@@ -341,6 +341,163 @@ drawvec close_poly(drawvec &geom) {
 	return out;
 }
 
+static bool inside(draw d, int edge, long long area, long long buffer) {
+	long long clip_buffer = buffer * area / 256;
+
+	switch (edge) {
+	case 0:  // top
+		return d.y > -clip_buffer;
+
+	case 1:  // right
+		return d.x < area + clip_buffer;
+
+	case 2:  // bottom
+		return d.y < area + clip_buffer;
+
+	case 3:  // left
+		return d.x > -clip_buffer;
+	}
+
+	fprintf(stderr, "internal error inside\n");
+	exit(EXIT_FAILURE);
+}
+
+// http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
+static draw get_line_intersection(draw p0, draw p1, draw p2, draw p3) {
+	double s1_x = p1.x - p0.x;
+	double s1_y = p1.y - p0.y;
+	double s2_x = p3.x - p2.x;
+	double s2_y = p3.y - p2.y;
+
+	double t;
+	// s = (-s1_y * (p0.x - p2.x) + s1_x * (p0.y - p2.y)) / (-s2_x * s1_y + s1_x * s2_y);
+	t = (s2_x * (p0.y - p2.y) - s2_y * (p0.x - p2.x)) / (-s2_x * s1_y + s1_x * s2_y);
+
+	return draw(VT_LINETO, p0.x + (t * s1_x), p0.y + (t * s1_y));
+}
+
+static draw intersect(draw a, draw b, int edge, long long area, long long buffer) {
+	long long clip_buffer = buffer * area / 256;
+
+	switch (edge) {
+	case 0:  // top
+		return get_line_intersection(a, b, draw(VT_MOVETO, -clip_buffer, -clip_buffer), draw(VT_MOVETO, area + clip_buffer, -clip_buffer));
+		break;
+
+	case 1:  // right
+		return get_line_intersection(a, b, draw(VT_MOVETO, area + clip_buffer, -clip_buffer), draw(VT_MOVETO, area + clip_buffer, area + clip_buffer));
+		break;
+
+	case 2:  // bottom
+		return get_line_intersection(a, b, draw(VT_MOVETO, area + clip_buffer, area + clip_buffer), draw(VT_MOVETO, -clip_buffer, area + clip_buffer));
+		break;
+
+	case 3:  // left
+		return get_line_intersection(a, b, draw(VT_MOVETO, -clip_buffer, area + clip_buffer), draw(VT_MOVETO, -clip_buffer, -clip_buffer));
+		break;
+	}
+
+	fprintf(stderr, "internal error intersecting\n");
+	exit(EXIT_FAILURE);
+}
+
+// http://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm
+static drawvec clip_poly1(drawvec &geom, int z, int detail, int buffer) {
+	drawvec out = geom;
+
+	long long area = 0xFFFFFFFF;
+	if (z != 0) {
+		area = 1LL << (32 - z);
+	}
+
+	for (int edge = 0; edge < 4; edge++) {
+		if (out.size() > 0) {
+			drawvec in = out;
+			out.resize(0);
+
+			draw S = in[in.size() - 1];
+
+			for (unsigned e = 0; e < in.size(); e++) {
+				draw E = in[e];
+
+				if (inside(E, edge, area, buffer)) {
+					if (!inside(S, edge, area, buffer)) {
+						out.push_back(intersect(S, E, edge, area, buffer));
+					}
+					out.push_back(E);
+				} else if (inside(S, edge, area, buffer)) {
+					out.push_back(intersect(S, E, edge, area, buffer));
+				}
+
+				S = E;
+			}
+		}
+	}
+
+	if (out.size() > 0) {
+		// If the polygon begins and ends outside the edge,
+		// the starting and ending points will be left as the
+		// places where it intersects the edge. Need to add
+		// another point to close the loop.
+
+		if (out[0].x != out[out.size() - 1].x || out[0].y != out[out.size() - 1].y) {
+			out.push_back(out[0]);
+		}
+
+		if (out.size() < 3) {
+			fprintf(stderr, "Polygon degenerated to a line segment\n");
+		}
+
+		out[0].op = VT_MOVETO;
+		for (unsigned i = 1; i < out.size(); i++) {
+			out[i].op = VT_LINETO;
+		}
+	}
+
+	return out;
+}
+
+drawvec simple_clip_poly(drawvec &geom, int z, int detail, int buffer) {
+	if (z == 0) {
+		return geom;
+	}
+
+	drawvec out;
+
+	for (unsigned i = 0; i < geom.size(); i++) {
+		if (geom[i].op == VT_MOVETO) {
+			unsigned j;
+			for (j = i + 1; j < geom.size(); j++) {
+				if (geom[j].op != VT_LINETO) {
+					break;
+				}
+			}
+
+			drawvec tmp;
+			for (unsigned k = i; k < j; k++) {
+				tmp.push_back(geom[k]);
+			}
+			tmp = clip_poly1(tmp, z, detail, buffer);
+			if (tmp.size() > 0) {
+				if (tmp[0].x != tmp[tmp.size() - 1].x || tmp[0].y != tmp[tmp.size() - 1].y) {
+					fprintf(stderr, "Internal error: Polygon ring not closed\n");
+					exit(EXIT_FAILURE);
+				}
+			}
+			for (unsigned k = 0; k < tmp.size(); k++) {
+				out.push_back(tmp[k]);
+			}
+
+			i = j - 1;
+		} else {
+			fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op);
+			exit(EXIT_FAILURE);
+		}
+	}
+
+	return out;
+}
+
 drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double *accum_area) {
 	drawvec out;
 	long long pixel = (1 << (32 - detail - z)) * 2;

geometry.hh

@@ -21,6 +21,7 @@ void to_tile_scale(drawvec &geom, int z, int detail);
 drawvec remove_noop(drawvec geom, int type, int shift);
 drawvec clip_point(drawvec &geom, int z, int detail, long long buffer);
 drawvec clean_or_clip_poly(drawvec &geom, int z, int detail, int buffer, bool clip);
+drawvec simple_clip_poly(drawvec &geom, int z, int detail, int buffer);
 drawvec close_poly(drawvec &geom);
 drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double *accum_area);
 drawvec clip_lines(drawvec &geom, int z, int detail, long long buffer);

tile.cc

@@ -465,7 +465,93 @@ void rewrite(drawvec &geom, int z, int nextzoom, int maxzoom, long long *bbox, u
 	}
 }
 
-long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsigned tx, unsigned ty, int detail, int min_detail, int basezoom, struct pool **file_keys, char **layernames, sqlite3 *outdb, double droprate, int buffer, const char *fname, FILE **geomfile, int minzoom, int maxzoom, double todo, char *geomstart, volatile long long *along, double gamma, int nlayers, char *prevent, char *additional, int child_shards, long long *meta_off, long long *pool_off, unsigned *initial_x, unsigned *initial_y) {
+struct partial {
+	drawvec geom;
+	long long layer;
+	char *meta;
+	signed char t;
+	int segment;
+	long long original_seq;
+	bool reduced;
+	unsigned long long index;
+	unsigned long long index2;
+	int z;
+	int line_detail;
+	char *prevent;
+	char *additional;
+};
+
+struct partial_arg {
+	std::vector<struct partial> *partials;
+	int task;
+	int tasks;
+};
+
+void *partial_feature_worker(void *v) {
+	struct partial_arg *a = (struct partial_arg *) v;
+	std::vector<struct partial> *partials = a->partials;
+
+	for (unsigned i = a->task; i < (*partials).size(); i += a->tasks) {
+		drawvec geom = (*partials)[i].geom;
+		(*partials)[i].geom.clear(); // avoid keeping two copies in memory
+		signed char t = (*partials)[i].t;
+		int z = (*partials)[i].z;
+		int line_detail = (*partials)[i].line_detail;
+		char *prevent = (*partials)[i].prevent;
+		char *additional = (*partials)[i].additional;
+
+		if ((t == VT_LINE || t == VT_POLYGON) && !prevent['s' & 0xFF]) {
+			if (1 /* !reduced */) { // XXX why did this not simplify if reduced?
+				if (t == VT_LINE) {
+					geom = remove_noop(geom, t, 32 - z - line_detail);
+				}
+
+				geom = simplify_lines(geom, z, line_detail);
+			}
+		}
+
+#if 0
+		if (t == VT_LINE && z != basezoom) {
+			geom = shrink_lines(geom, z, line_detail, basezoom, &along);
+		}
+#endif
+
+		if (t == VT_LINE && additional['r' & 0xFF]) {
+			geom = reorder_lines(geom);
+		}
+
+		to_tile_scale(geom, z, line_detail);
+
+		if (t == VT_POLYGON) {
+			// Scaling may have made the polygon degenerate.
+			// Give Clipper a chance to try to fix it.
+			geom = clean_or_clip_poly(geom, 0, 0, 0, false);
+			geom = close_poly(geom);
+		}
+
+		// Worth skipping this if not coalescing anyway?
+		if (geom.size() > 0) {
+			(*partials)[i].index = encode(geom[0].x, geom[0].y);
+			(*partials)[i].index2 = encode(geom[geom.size() - 1].x, geom[geom.size() - 1].y);
+
+			// Anything numbered below the start of the line
+			// can't possibly be the next feature.
+			// We want lowest-but-not-under.
+			if ((*partials)[i].index2 < (*partials)[i].index) {
+				(*partials)[i].index2 = ~0LL;
+			}
+		} else {
+			(*partials)[i].index = 0;
+			(*partials)[i].index2 = 0;
+		}
+
+		(*partials)[i].geom = geom;
+	}
+
+	return NULL;
+}
+
+long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsigned tx, unsigned ty, int detail, int min_detail, int basezoom, struct pool **file_keys, char **layernames, sqlite3 *outdb, double droprate, int buffer, const char *fname, FILE **geomfile, int minzoom, int maxzoom, double todo, char *geomstart, volatile long long *along, double gamma, int nlayers, char *prevent, char *additional, int child_shards, long long *meta_off, long long *pool_off, unsigned *initial_x, unsigned *initial_y, volatile int *running) {
 	int line_detail;
 	double fraction = 1;
 
@@ -523,6 +609,7 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi
 		long long original_features = 0;
 		long long unclipped_features = 0;
 
+		std::vector<struct partial> partials;
 		std::vector<std::vector<coalesce> > features;
 		for (i = 0; i < nlayers; i++) {
 			features.push_back(std::vector<coalesce>());
@@ -621,7 +708,7 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi
 					geom = clip_lines(geom, z, line_detail, buffer);
 				}
 				if (t == VT_POLYGON) {
-					geom = clean_or_clip_poly(geom, z, line_detail, buffer, true);
+					geom = simple_clip_poly(geom, z, line_detail, buffer);
 				}
 				if (t == VT_POINT) {
 					geom = clip_point(geom, z, line_detail, buffer);
@@ -706,53 +793,71 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi
 				geom = reduce_tiny_poly(geom, z, line_detail, &reduced, &accum_area);
 			}
 
-			if ((t == VT_LINE || t == VT_POLYGON) && !prevent['s' & 0xFF]) {
+			if (geom.size() > 0) {
-				if (!reduced) {
+				partial p;
-					if (t == VT_LINE) {
+				p.geom = geom;
-						geom = remove_noop(geom, t, 32 - z - line_detail);
+				p.layer = layer;
-					}
+				p.meta = meta;
+				p.t = t;
+				p.segment = segment;
+				p.original_seq = original_seq;
+				p.reduced = reduced;
+				p.z = z;
+				p.line_detail = line_detail;
+				p.prevent = prevent;
+				p.additional = additional;
+				partials.push_back(p);
+			}
+		}
 
-					geom = simplify_lines(geom, z, line_detail);
+		int tasks = ceil((double) CPUS / *running);
+		if (tasks < 1) {
+			tasks = 1;
+		}
+
+		pthread_t pthreads[tasks];
+		partial_arg args[tasks];
+		for (int i = 0; i < tasks; i++) {
+			args[i].task = i;
+			args[i].tasks = tasks;
+			args[i].partials = &partials;
+
+			if (tasks > 1) {
+				if (pthread_create(&pthreads[i], NULL, partial_feature_worker, &args[i]) != 0) {
+					perror("pthread_create");
+					exit(EXIT_FAILURE);
+				}
+			} else {
+				partial_feature_worker(&args[i]);
+			}
+		}
+
+		if (tasks > 1) {
+			for (int i = 0; i < tasks; i++) {
+				void *retval;
+
+				if (pthread_join(pthreads[i], &retval) != 0) {
+					perror("pthread_join");
 				}
 			}
+		}
 
-#if 0
+		// This is serial because decode_meta() unifies duplicates
-			if (t == VT_LINE && z != basezoom) {
+		for (unsigned i = 0; i < partials.size(); i++) {
-				geom = shrink_lines(geom, z, line_detail, basezoom, &along);
+			drawvec geom = partials[i].geom;
-			}
+			partials[i].geom.clear(); // avoid keeping two copies in memory
-#endif
+			long long layer = partials[i].layer;
-
+			char *meta = partials[i].meta;
-			if (t == VT_LINE && additional['r' & 0xFF]) {
+			signed char t = partials[i].t;
-				geom = reorder_lines(geom);
+			int segment = partials[i].segment;
-			}
+			long long original_seq = partials[i].original_seq;
-
-			to_tile_scale(geom, z, line_detail);
-
-			if (t == VT_POLYGON) {
-				// Scaling may have made the polygon degenerate.
-				// Give Clipper a chance to try to fix it.
-				geom = clean_or_clip_poly(geom, 0, 0, 0, false);
-				geom = close_poly(geom);
-			}
 
 			if (t == VT_POINT || to_feature(geom, NULL)) {
 				struct coalesce c;
 
 				c.type = t;
-				if (geom.size() > 0) {
+				c.index = partials[i].index;
-					c.index = encode(geom[0].x, geom[0].y);
+				c.index2 = partials[i].index2;
-					c.index2 = encode(geom[geom.size() - 1].x, geom[geom.size() - 1].y);
-
-					// Anything numbered below the start of the line
-					// can't possibly be the next feature.
-					// We want lowest-but-not-under.
-					if (c.index2 < c.index) {
-						c.index2 = ~0LL;
-					}
-				} else {
-					c.index = 0;
-					c.index2 = 0;
-				}
 				c.geom = geom;
 				c.metasrc = meta;
 				c.coalesced = false;
@@ -931,6 +1036,7 @@ struct write_tile_args {
 	long long *pool_off;
 	unsigned *initial_x;
 	unsigned *initial_y;
+	volatile int *running;
 };
 
 void *run_thread(void *vargs) {
@@ -970,7 +1076,7 @@ void *run_thread(void *vargs) {
 
 			// fprintf(stderr, "%d/%u/%u\n", z, x, y);
 
-			long long len = write_tile(&geom, arg->metabase, arg->stringpool, z, x, y, z == arg->maxzoom ? arg->full_detail : arg->low_detail, arg->min_detail, arg->basezoom, arg->file_keys, arg->layernames, arg->outdb, arg->droprate, arg->buffer, arg->fname, arg->geomfile, arg->minzoom, arg->maxzoom, arg->todo, geomstart, arg->along, arg->gamma, arg->nlayers, arg->prevent, arg->additional, arg->child_shards, arg->meta_off, arg->pool_off, arg->initial_x, arg->initial_y);
+			long long len = write_tile(&geom, arg->metabase, arg->stringpool, z, x, y, z == arg->maxzoom ? arg->full_detail : arg->low_detail, arg->min_detail, arg->basezoom, arg->file_keys, arg->layernames, arg->outdb, arg->droprate, arg->buffer, arg->fname, arg->geomfile, arg->minzoom, arg->maxzoom, arg->todo, geomstart, arg->along, arg->gamma, arg->nlayers, arg->prevent, arg->additional, arg->child_shards, arg->meta_off, arg->pool_off, arg->initial_x, arg->initial_y, arg->running);
 
 			if (len < 0) {
 				int *err = (int *) malloc(sizeof(int));
@@ -1003,6 +1109,7 @@ void *run_thread(void *vargs) {
 		}
 	}
 
+	arg->running--;
 	return NULL;
 }
 
@@ -1099,6 +1206,7 @@ int traverse_zooms(int *geomfd, off_t *geom_size, char *metabase, char *stringpo
 
 		pthread_t pthreads[threads];
 		write_tile_args args[threads];
+		int running = threads;
 
 		int thread;
 		for (thread = 0; thread < threads; thread++) {
@@ -1136,6 +1244,7 @@ int traverse_zooms(int *geomfd, off_t *geom_size, char *metabase, char *stringpo
 			args[thread].initial_y = initial_y;
 
 			args[thread].tasks = dispatches[thread].tasks;
+			args[thread].running = &running;
 
 			if (pthread_create(&pthreads[thread], NULL, run_thread, &args[thread]) != 0) {
 				perror("pthread_create");

version.h

@@ -1 +1 @@
-#define VERSION "tippecanoe v1.6.0\n"
+#define VERSION "tippecanoe v1.6.1\n"