From dace304182302a8eef550e601d658fdfd3852b27 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Tue, 8 Dec 2015 16:24:17 -0800 Subject: [PATCH 1/5] Offset coordinates for Clipper to keep them positive. Limit very high or low latitudes and longitudes. --- geometry.cc | 22 ++++++++++++---------- projection.c | 17 +++++++++++++++++ 2 files changed, 29 insertions(+), 10 deletions(-) diff --git a/geometry.cc b/geometry.cc index 586bfa46..41f80873 100644 --- a/geometry.cc +++ b/geometry.cc @@ -206,6 +206,8 @@ drawvec shrink_lines(drawvec &geom, int z, int detail, int basezoom, long long * } #endif +#define CLIPPER_OFFSET 0x400000000LL + static void decode_clipped(ClipperLib::PolyNode *t, drawvec &out) { // To make the GeoJSON come out right, we need to do each of the // outer rings followed by its children if any, and then go back @@ -213,19 +215,19 @@ static void decode_clipped(ClipperLib::PolyNode *t, drawvec &out) { ClipperLib::Path p = t->Contour; for (unsigned i = 0; i < p.size(); i++) { - out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X, p[i].Y)); + out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X - CLIPPER_OFFSET, p[i].Y - CLIPPER_OFFSET)); } if (p.size() > 0) { - out.push_back(draw(VT_LINETO, p[0].X, p[0].Y)); + out.push_back(draw(VT_LINETO, p[0].X - CLIPPER_OFFSET, p[0].Y - CLIPPER_OFFSET)); } for (int n = 0; n < t->ChildCount(); n++) { ClipperLib::Path p = t->Childs[n]->Contour; for (unsigned i = 0; i < p.size(); i++) { - out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X, p[i].Y)); + out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X - CLIPPER_OFFSET, p[i].Y - CLIPPER_OFFSET)); } if (p.size() > 0) { - out.push_back(draw(VT_LINETO, p[0].X, p[0].Y)); + out.push_back(draw(VT_LINETO, p[0].X - CLIPPER_OFFSET, p[0].Y - CLIPPER_OFFSET)); } } @@ -252,7 +254,7 @@ drawvec clean_or_clip_poly(drawvec &geom, int z, int detail, int buffer, bool cl drawvec tmp; for (unsigned k = i; k < j; k++) { - path.push_back(ClipperLib::IntPoint(geom[k].x, geom[k].y)); + path.push_back(ClipperLib::IntPoint(geom[k].x + CLIPPER_OFFSET, geom[k].y + CLIPPER_OFFSET)); } if (!clipper.AddPath(path, ClipperLib::ptSubject, true)) { @@ -280,11 +282,11 @@ drawvec clean_or_clip_poly(drawvec &geom, int z, int detail, int buffer, bool cl long long clip_buffer = buffer * area / 256; ClipperLib::Path edge; - edge.push_back(ClipperLib::IntPoint(-clip_buffer, -clip_buffer)); - edge.push_back(ClipperLib::IntPoint(area + clip_buffer, -clip_buffer)); - edge.push_back(ClipperLib::IntPoint(area + clip_buffer, area + clip_buffer)); - edge.push_back(ClipperLib::IntPoint(-clip_buffer, area + clip_buffer)); - edge.push_back(ClipperLib::IntPoint(-clip_buffer, -clip_buffer)); + edge.push_back(ClipperLib::IntPoint(-clip_buffer + CLIPPER_OFFSET, -clip_buffer + CLIPPER_OFFSET)); + edge.push_back(ClipperLib::IntPoint(area + clip_buffer + CLIPPER_OFFSET, -clip_buffer + CLIPPER_OFFSET)); + edge.push_back(ClipperLib::IntPoint(area + clip_buffer + CLIPPER_OFFSET, area + clip_buffer + CLIPPER_OFFSET)); + edge.push_back(ClipperLib::IntPoint(-clip_buffer + CLIPPER_OFFSET, area + clip_buffer + CLIPPER_OFFSET)); + edge.push_back(ClipperLib::IntPoint(-clip_buffer + CLIPPER_OFFSET, -clip_buffer + CLIPPER_OFFSET)); clipper.AddPath(edge, ClipperLib::ptClip, true); } diff --git a/projection.c b/projection.c index 839370ed..ccd4f00d 100644 --- a/projection.c +++ b/projection.c @@ -3,6 +3,23 @@ // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames void latlon2tile(double lat, double lon, int zoom, long long *x, long long *y) { + // Must limit latitude somewhere to prevent overflow. + // 89.9 degrees latitude is 0.621 worlds beyond the edge of the flat earth, + // hopefully far enough out that there are few expectations about the shape. + if (lat < -89.9) { + lat = -89.9; + } + if (lat > 89.9) { + lat = 89.9; + } + + if (lon < -360) { + lon = -360; + } + if (lon > 360) { + lon = 360; + } + double lat_rad = lat * M_PI / 180; unsigned long long n = 1LL << zoom; From 256139b385e775f2a6f41f5705bdfd7005902b93 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Tue, 8 Dec 2015 16:57:04 -0800 Subject: [PATCH 2/5] Clipping is faster with only one duplicate/shifted geometry copy --- tile.cc | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/tile.cc b/tile.cc index 6daa00a8..b21e6b03 100644 --- a/tile.cc +++ b/tile.cc @@ -550,11 +550,17 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi // shifted by 360 degrees, and then make sure both copies get clipped down to size. unsigned n = geom.size(); - for (unsigned i = 0; i < n; i++) { - geom.push_back(draw(geom[i].op, geom[i].x - (1LL << 32), geom[i].y)); + + if (bbox[0] < 0) { + for (unsigned i = 0; i < n; i++) { + geom.push_back(draw(geom[i].op, geom[i].x + (1LL << 32), geom[i].y)); + } } - for (unsigned i = 0; i < n; i++) { - geom.push_back(draw(geom[i].op, geom[i].x + (1LL << 32), geom[i].y)); + + if (bbox[2] > 1LL << 32) { + for (unsigned i = 0; i < n; i++) { + geom.push_back(draw(geom[i].op, geom[i].x - (1LL << 32), geom[i].y)); + } } bbox[0] = 0; From 4bde17f8ff37bb3fff954176ab1839093d87bd69 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Tue, 8 Dec 2015 17:21:59 -0800 Subject: [PATCH 3/5] Remove unnecessary coordinate offsetting. Negative coordinates are OK. --- geometry.cc | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/geometry.cc b/geometry.cc index 41f80873..586bfa46 100644 --- a/geometry.cc +++ b/geometry.cc @@ -206,8 +206,6 @@ drawvec shrink_lines(drawvec &geom, int z, int detail, int basezoom, long long * } #endif -#define CLIPPER_OFFSET 0x400000000LL - static void decode_clipped(ClipperLib::PolyNode *t, drawvec &out) { // To make the GeoJSON come out right, we need to do each of the // outer rings followed by its children if any, and then go back @@ -215,19 +213,19 @@ static void decode_clipped(ClipperLib::PolyNode *t, drawvec &out) { ClipperLib::Path p = t->Contour; for (unsigned i = 0; i < p.size(); i++) { - out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X - CLIPPER_OFFSET, p[i].Y - CLIPPER_OFFSET)); + out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X, p[i].Y)); } if (p.size() > 0) { - out.push_back(draw(VT_LINETO, p[0].X - CLIPPER_OFFSET, p[0].Y - CLIPPER_OFFSET)); + out.push_back(draw(VT_LINETO, p[0].X, p[0].Y)); } for (int n = 0; n < t->ChildCount(); n++) { ClipperLib::Path p = t->Childs[n]->Contour; for (unsigned i = 0; i < p.size(); i++) { - out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X - CLIPPER_OFFSET, p[i].Y - CLIPPER_OFFSET)); + out.push_back(draw((i == 0) ? VT_MOVETO : VT_LINETO, p[i].X, p[i].Y)); } if (p.size() > 0) { - out.push_back(draw(VT_LINETO, p[0].X - CLIPPER_OFFSET, p[0].Y - CLIPPER_OFFSET)); + out.push_back(draw(VT_LINETO, p[0].X, p[0].Y)); } } @@ -254,7 +252,7 @@ drawvec clean_or_clip_poly(drawvec &geom, int z, int detail, int buffer, bool cl drawvec tmp; for (unsigned k = i; k < j; k++) { - path.push_back(ClipperLib::IntPoint(geom[k].x + CLIPPER_OFFSET, geom[k].y + CLIPPER_OFFSET)); + path.push_back(ClipperLib::IntPoint(geom[k].x, geom[k].y)); } if (!clipper.AddPath(path, ClipperLib::ptSubject, true)) { @@ -282,11 +280,11 @@ drawvec clean_or_clip_poly(drawvec &geom, int z, int detail, int buffer, bool cl long long clip_buffer = buffer * area / 256; ClipperLib::Path edge; - edge.push_back(ClipperLib::IntPoint(-clip_buffer + CLIPPER_OFFSET, -clip_buffer + CLIPPER_OFFSET)); - edge.push_back(ClipperLib::IntPoint(area + clip_buffer + CLIPPER_OFFSET, -clip_buffer + CLIPPER_OFFSET)); - edge.push_back(ClipperLib::IntPoint(area + clip_buffer + CLIPPER_OFFSET, area + clip_buffer + CLIPPER_OFFSET)); - edge.push_back(ClipperLib::IntPoint(-clip_buffer + CLIPPER_OFFSET, area + clip_buffer + CLIPPER_OFFSET)); - edge.push_back(ClipperLib::IntPoint(-clip_buffer + CLIPPER_OFFSET, -clip_buffer + CLIPPER_OFFSET)); + edge.push_back(ClipperLib::IntPoint(-clip_buffer, -clip_buffer)); + edge.push_back(ClipperLib::IntPoint(area + clip_buffer, -clip_buffer)); + edge.push_back(ClipperLib::IntPoint(area + clip_buffer, area + clip_buffer)); + edge.push_back(ClipperLib::IntPoint(-clip_buffer, area + clip_buffer)); + edge.push_back(ClipperLib::IntPoint(-clip_buffer, -clip_buffer)); clipper.AddPath(edge, ClipperLib::ptClip, true); } From f04c5e153a0e34ea3a4f72a48df4dfddfe4703f3 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Wed, 9 Dec 2015 15:02:59 -0800 Subject: [PATCH 4/5] Avoid arithmetic overflow in area calculation --- geometry.cc | 62 ++++++++++++++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 25 deletions(-) diff --git a/geometry.cc b/geometry.cc index 586bfa46..36ba691e 100644 --- a/geometry.cc +++ b/geometry.cc @@ -357,41 +357,53 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double double area = 0; for (unsigned k = i; k < j; k++) { - area += geom[k].x * geom[i + ((k - i + 1) % (j - i))].y; - area -= geom[k].y * geom[i + ((k - i + 1) % (j - i))].x; + area += (long double) geom[k].x * (long double) geom[i + ((k - i + 1) % (j - i))].y; + area -= (long double) geom[k].y * (long double) geom[i + ((k - i + 1) % (j - i))].x; } area = area / 2; - if (fabs(area) <= pixel * pixel || (area < 0 && !included_last_outer)) { - // printf("area is only %f vs %lld so using square\n", area, pixel * pixel); + // XXX There is an ambiguity here: If the area of a ring is 0 and it is followed by holes, + // we don't know whether the area-0 ring was a hole too or whether it was the outer ring + // that these subsequent holes are somehow being subtracted from. I hope that if a polygon + // was simplified down to nothing, its holes also became nothing. - *accum_area += area; - if (*accum_area > pixel * pixel) { - // XXX use centroid; + if (area != 0) { + // These are pixel coordinates, so area > 0 for the outer ring. + // If the outer ring of a polygon was reduced to a pixel, its + // inner rings must just have their area de-accumulated rather + // than being drawn since we don't really know where they are. - out.push_back(draw(VT_MOVETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y - pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y + pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y + pixel / 2)); - out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); + if (fabs(area) <= pixel * pixel || (area < 0 && !included_last_outer)) { + // printf("area is only %f vs %lld so using square\n", area, pixel * pixel); - *accum_area -= pixel * pixel; - } + *accum_area += area; + if (area > 0 && *accum_area > pixel * pixel) { + // XXX use centroid; - if (area >= 0) { - included_last_outer = false; - } - } else { - // printf("area is %f so keeping instead of %lld\n", area, pixel * pixel); + out.push_back(draw(VT_MOVETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y - pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x + pixel / 2, geom[i].y + pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y + pixel / 2)); + out.push_back(draw(VT_LINETO, geom[i].x - pixel / 2, geom[i].y - pixel / 2)); - for (unsigned k = i; k <= j && k < geom.size(); k++) { - out.push_back(geom[k]); - } + *accum_area -= pixel * pixel; + } - *reduced = false; + if (area > 0) { + included_last_outer = false; + } + } else { + // printf("area is %f so keeping instead of %lld\n", area, pixel * pixel); - if (area >= 0) { - included_last_outer = true; + for (unsigned k = i; k <= j && k < geom.size(); k++) { + out.push_back(geom[k]); + } + + *reduced = false; + + if (area > 0) { + included_last_outer = true; + } } } From a6be746163f92cfd40fb9943b8d37eb1c591cb90 Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Wed, 9 Dec 2015 15:24:07 -0800 Subject: [PATCH 5/5] Unrelated code formatting correction --- geojson.c | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/geojson.c b/geojson.c index 041492cf..c933657e 100644 --- a/geojson.c +++ b/geojson.c @@ -65,21 +65,21 @@ int CPUS; int TEMP_FILES; void init_cpus() { - CPUS = sysconf(_SC_NPROCESSORS_ONLN); - if (CPUS < 1) { - CPUS = 1; - } + CPUS = sysconf(_SC_NPROCESSORS_ONLN); + if (CPUS < 1) { + CPUS = 1; + } - TEMP_FILES = 64; - struct rlimit rl; - if (getrlimit(RLIMIT_NOFILE, &rl) != 0) { - perror("getrlimit"); - } else { - TEMP_FILES = rl.rlim_cur / 3; - if (TEMP_FILES > CPUS * 4) { - TEMP_FILES = CPUS * 4; - } - } + TEMP_FILES = 64; + struct rlimit rl; + if (getrlimit(RLIMIT_NOFILE, &rl) != 0) { + perror("getrlimit"); + } else { + TEMP_FILES = rl.rlim_cur / 3; + if (TEMP_FILES > CPUS * 4) { + TEMP_FILES = CPUS * 4; + } + } } size_t fwrite_check(const void *ptr, size_t size, size_t nitems, FILE *stream, const char *fname) {