From f43d18eb73afbb569e88b103a5c497b5ffe4208b Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Mon, 21 Dec 2015 15:56:21 -0800 Subject: [PATCH 1/4] Bring back the old simple polygon clipping algorithm --- geometry.cc | 157 ++++++++++++++++++++++++++++++++++++++++++++++++++++ geometry.hh | 1 + 2 files changed, 158 insertions(+) diff --git a/geometry.cc b/geometry.cc index 0005a5c..1b71b97 100644 --- a/geometry.cc +++ b/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; diff --git a/geometry.hh b/geometry.hh index 3f057ec..385b13a 100644 --- a/geometry.hh +++ b/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); From 977533e44988b83938b07d3921cf9891c8d61b64 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Tue, 5 Jan 2016 12:29:40 -0800 Subject: [PATCH 2/4] Use multiple threads within a single tile for geometric simplification --- tile.cc | 183 ++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 144 insertions(+), 39 deletions(-) diff --git a/tile.cc b/tile.cc index 600922e..8cf286d 100644 --- a/tile.cc +++ b/tile.cc @@ -465,6 +465,92 @@ void rewrite(drawvec &geom, int z, int nextzoom, int maxzoom, long long *bbox, u } } +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 *partials; + int task; + int tasks; +}; + +void *partial_feature_worker(void *v) { + struct partial_arg *a = (struct partial_arg *) v; + std::vector *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) { 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 partials; std::vector > features; for (i = 0; i < nlayers; i++) { features.push_back(std::vector()); @@ -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 (!reduced) { - if (t == VT_LINE) { - geom = remove_noop(geom, t, 32 - z - line_detail); - } + if (geom.size() > 0) { + partial p; + p.geom = geom; + 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 = CPUS / (1 << z); + 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 - 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); - } + // This is serial because decode_meta() unifies duplicates + for (unsigned i = 0; i < partials.size(); i++) { + drawvec geom = partials[i].geom; + partials[i].geom.clear(); // avoid keeping two copies in memory + long long layer = partials[i].layer; + char *meta = partials[i].meta; + signed char t = partials[i].t; + int segment = partials[i].segment; + long long original_seq = partials[i].original_seq; if (t == VT_POINT || to_feature(geom, NULL)) { struct coalesce c; c.type = t; - if (geom.size() > 0) { - c.index = encode(geom[0].x, geom[0].y); - 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.index = partials[i].index; + c.index2 = partials[i].index2; c.geom = geom; c.metasrc = meta; c.coalesced = false; From c8573634e1d5c29b55fff534e73c13a428db9cad Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Tue, 5 Jan 2016 13:56:36 -0800 Subject: [PATCH 3/4] Track how many threads are active to calculate how many sub-threads to use --- tile.cc | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/tile.cc b/tile.cc index 8cf286d..88312ba 100644 --- a/tile.cc +++ b/tile.cc @@ -551,7 +551,7 @@ void *partial_feature_worker(void *v) { 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) { +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; @@ -810,7 +810,7 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi } } - int tasks = CPUS / (1 << z); + int tasks = ceil((double) CPUS / *running); if (tasks < 1) { tasks = 1; } @@ -1036,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) { @@ -1075,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)); @@ -1108,6 +1109,7 @@ void *run_thread(void *vargs) { } } + arg->running--; return NULL; } @@ -1204,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++) { @@ -1241,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"); From 22293ca6e8acbaeb6083dfa0076e2fd2537d8db2 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Wed, 6 Jan 2016 13:29:59 -0800 Subject: [PATCH 4/4] Bump version number --- CHANGELOG.md | 4 ++++ version.h | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a15b744..9b4a2d2 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/version.h b/version.h index 7ca57b9..0a102b7 100644 --- a/version.h +++ b/version.h @@ -1 +1 @@ -#define VERSION "tippecanoe v1.6.0\n" +#define VERSION "tippecanoe v1.6.1\n"