From 92434d4a0aab8e64d48739e9f55c75546dad33dc Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Thu, 7 Sep 2023 16:48:12 -0700 Subject: [PATCH] Pull clipping and scaling code back out into clip.cpp --- clip.cpp | 743 +++++++++++++++++++++++++++++++++++++++++++++++++- geometry.cpp | 756 +-------------------------------------------------- geometry.hpp | 6 +- serial.cpp | 2 +- tile.cpp | 2 +- 5 files changed, 750 insertions(+), 759 deletions(-) diff --git a/clip.cpp b/clip.cpp index 7d87d67e..0a6d0b43 100644 --- a/clip.cpp +++ b/clip.cpp @@ -8,6 +8,745 @@ #include "compression.hpp" #include "mvt.hpp" +void to_tile_scale(drawvec &geom, int z, int detail) { + if (32 - detail - z < 0) { + for (size_t i = 0; i < geom.size(); i++) { + 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)))); + } + } else { + for (size_t i = 0; i < geom.size(); i++) { + 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 *= (1LL << (32 - detail - z)); + d.y *= (1LL << (32 - detail - z)); + out.push_back(d); + } + return out; +} + +drawvec remove_noop(drawvec geom, int type, int shift) { + // first pass: remove empty linetos + + long long ox = 0, oy = 0; + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + long long nx = std::round((double) geom[i].x / (1LL << shift)); + long long ny = std::round((double) geom[i].y / (1LL << shift)); + + if (geom[i].op == VT_LINETO && nx == ox && ny == oy) { + continue; + } + + if (geom[i].op == VT_CLOSEPATH) { + out.push_back(geom[i]); + } else { /* moveto or lineto */ + out.push_back(geom[i]); + ox = nx; + oy = ny; + } + } + + // second pass: remove unused movetos + + if (type != VT_POINT) { + geom = out; + out.resize(0); + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + if (i + 1 >= geom.size()) { + // followed by end-of-geometry: not needed + continue; + } + + if (geom[i + 1].op == VT_MOVETO) { + // followed by another moveto: not needed + continue; + } + + if (geom[i + 1].op == VT_CLOSEPATH) { + // followed by closepath: not possible + fprintf(stderr, "Shouldn't happen\n"); + i++; // also remove unused closepath + continue; + } + } + + out.push_back(geom[i]); + } + } + + // second pass: remove empty movetos + + if (type == VT_LINE) { + geom = out; + out.resize(0); + + for (size_t i = 0; i < geom.size(); i++) { + if (i > 1 && geom[i].op == VT_MOVETO) { + if (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; + } + } + + out.push_back(geom[i]); + } + } + + return out; +} + +double get_area_scaled(const drawvec &geom, size_t i, size_t j) { + const double max_exact_double = (double) ((1LL << 53) - 1); + + // keep scaling the geometry down until we can calculate its area without overflow + for (long long scale = 2; scale < (1LL << 30); scale *= 2) { + long long bx = geom[i].x; + long long by = geom[i].y; + bool again = false; + + // https://en.wikipedia.org/wiki/Shoelace_formula + double area = 0; + for (size_t k = i; k < j; k++) { + area += (double) ((geom[k].x - bx) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].y - by) / scale); + if (std::fabs(area) >= max_exact_double) { + again = true; + break; + } + area -= (double) ((geom[k].y - by) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].x - bx) / scale); + if (std::fabs(area) >= max_exact_double) { + again = true; + break; + } + } + + if (again) { + continue; + } else { + area /= 2; + return area * scale * scale; + } + } + + fprintf(stderr, "get_area_scaled: can't happen\n"); + exit(EXIT_IMPOSSIBLE); +} + +double get_area(const drawvec &geom, size_t i, size_t j) { + const double max_exact_double = (double) ((1LL << 53) - 1); + + // Coordinates in `geom` are 40-bit integers, so there is no good way + // to multiply them without possible precision loss. Since they probably + // do not use the full precision, shift them nearer to the origin so + // their product is more likely to be exactly representable as a double. + // + // (In practice they are actually 34-bit integers: 32 bits for the + // Mercator world plane, plus another two bits so features can stick + // off either the left or right side. But that is still too many bits + // for the product to fit either in a 64-bit long long or in a + // double where the largest exact integer is 2^53.) + // + // If the intermediate calculation still exceeds 2^53, start trying to + // recalculate the area by scaling down the geometry. This will not + // produce as precise an area, but it will still be close, and the + // sign will be correct, which is more important, since the sign + // determines the winding order of the rings. We can then use that + // sign with this generally more precise area calculation. + + long long bx = geom[i].x; + long long by = geom[i].y; + + // https://en.wikipedia.org/wiki/Shoelace_formula + double area = 0; + bool overflow = false; + for (size_t k = i; k < j; k++) { + area += (double) (geom[k].x - bx) * (double) (geom[i + ((k - i + 1) % (j - i))].y - by); + if (std::fabs(area) >= max_exact_double) { + overflow = true; + } + area -= (double) (geom[k].y - by) * (double) (geom[i + ((k - i + 1) % (j - i))].x - bx); + if (std::fabs(area) >= max_exact_double) { + overflow = true; + } + } + area /= 2; + + if (overflow) { + double scaled_area = get_area_scaled(geom, i, j); + if ((area < 0 && scaled_area > 0) || (area > 0 && scaled_area < 0)) { + area = -area; + } + } + + return area; +} + +double get_mp_area(drawvec &geom) { + double ret = 0; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + ret += get_area(geom, i, j); + i = j - 1; + } + } + + return ret; +} + +static void decode_clipped(mapbox::geometry::multi_polygon &t, drawvec &out, double scale) { + out.clear(); + + for (size_t i = 0; i < t.size(); i++) { + for (size_t j = 0; j < t[i].size(); j++) { + drawvec ring; + + for (size_t k = 0; k < t[i][j].size(); k++) { + ring.push_back(draw((k == 0) ? VT_MOVETO : VT_LINETO, std::round(t[i][j][k].x / scale), std::round(t[i][j][k].y / scale))); + } + + if (ring.size() > 0 && ring[ring.size() - 1] != ring[0]) { + fprintf(stderr, "Had to close ring\n"); + ring.push_back(draw(VT_LINETO, ring[0].x, ring[0].y)); + } + + double area = get_area(ring, 0, ring.size()); + + if ((j == 0 && area < 0) || (j != 0 && area > 0)) { + fprintf(stderr, "Ring area has wrong sign: %f for %zu\n", area, j); + exit(EXIT_IMPOSSIBLE); + } + + for (size_t k = 0; k < ring.size(); k++) { + out.push_back(ring[k]); + } + } + } +} + +drawvec clean_or_clip_poly(drawvec &geom, int z, int buffer, bool clip, bool try_scaling) { + geom = remove_noop(geom, VT_POLYGON, 0); + mapbox::geometry::multi_polygon result; + + double scale = 16.0; + if (!try_scaling) { + scale = 1.0; + } + + bool again = true; + while (again) { + mapbox::geometry::wagyu::wagyu wagyu; + again = false; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + if (j >= i + 4) { + mapbox::geometry::linear_ring lr; + + for (size_t k = i; k < j; k++) { + lr.push_back(mapbox::geometry::point(geom[k].x * scale, geom[k].y * scale)); + } + + if (lr.size() >= 3) { + wagyu.add_ring(lr); + } + } + + i = j - 1; + } + } + + if (clip) { + long long area = 0xFFFFFFFF; + if (z != 0) { + area = 1LL << (32 - z); + } + long long clip_buffer = buffer * area / 256; + + mapbox::geometry::linear_ring lr; + + lr.push_back(mapbox::geometry::point(scale * -clip_buffer, scale * -clip_buffer)); + lr.push_back(mapbox::geometry::point(scale * -clip_buffer, scale * (area + clip_buffer))); + lr.push_back(mapbox::geometry::point(scale * (area + clip_buffer), scale * (area + clip_buffer))); + lr.push_back(mapbox::geometry::point(scale * (area + clip_buffer), scale * -clip_buffer)); + lr.push_back(mapbox::geometry::point(scale * -clip_buffer, scale * -clip_buffer)); + + wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip); + } + + try { + result.clear(); + wagyu.execute(mapbox::geometry::wagyu::clip_type_union, result, mapbox::geometry::wagyu::fill_type_positive, mapbox::geometry::wagyu::fill_type_positive); + } catch (std::runtime_error &e) { + FILE *f = fopen("/tmp/wagyu.log", "w"); + fprintf(f, "%s\n", e.what()); + fprintf(stderr, "%s\n", e.what()); + fprintf(f, "["); + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + if (j >= i + 4) { + mapbox::geometry::linear_ring lr; + + if (i != 0) { + fprintf(f, ","); + } + fprintf(f, "["); + + for (size_t k = i; k < j; k++) { + lr.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); + if (k != i) { + fprintf(f, ","); + } + fprintf(f, "[%lld,%lld]", geom[k].x, geom[k].y); + } + + fprintf(f, "]"); + + if (lr.size() >= 3) { + } + } + + i = j - 1; + } + } + + fprintf(f, "]"); + fprintf(f, "\n\n\n\n\n"); + + fclose(f); + fprintf(stderr, "Internal error: Polygon cleaning failed. Log in /tmp/wagyu.log\n"); + exit(EXIT_IMPOSSIBLE); + } + + if (scale != 1) { + for (auto const &outer : result) { + for (auto const &ring : outer) { + for (auto const &p : ring) { + if (p.x / scale != std::round(p.x / scale) || + p.y / scale != std::round(p.y / scale)) { + scale = 1; + again = true; + break; + } + } + } + } + } + } + + drawvec ret; + decode_clipped(result, ret, scale); + return ret; +} + +drawvec close_poly(drawvec &geom) { + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + if (j - 1 > i) { + if (geom[j - 1].x != geom[i].x || geom[j - 1].y != geom[i].y) { + fprintf(stderr, "Internal error: polygon not closed\n"); + } + } + + for (size_t n = i; n < j - 1; n++) { + out.push_back(geom[n]); + } + out.push_back(draw(VT_CLOSEPATH, 0, 0)); + + i = j - 1; + } + } + + return out; +} + +static bool inside(std::pair d, int edge, long long minx, long long miny, long long maxx, long long maxy) { + switch (edge) { + case 0: // top + return d.second > miny; + + case 1: // right + return d.first < maxx; + + case 2: // bottom + return d.second < maxy; + + case 3: // left + return d.first > minx; + } + + fprintf(stderr, "internal error inside\n"); + exit(EXIT_FAILURE); +} + +static std::pair intersect(std::pair a, std::pair b, int edge, long long minx, long long miny, long long maxx, long long maxy) { + switch (edge) { + case 0: // top + return std::pair((a.first + (double) (b.first - a.first) * (miny - a.second) / (b.second - a.second)), miny); + + case 1: // right + return std::pair(maxx, (a.second + (double) (b.second - a.second) * (maxx - a.first) / (b.first - a.first))); + + case 2: // bottom + return std::pair((a.first + (double) (b.first - a.first) * (maxy - a.second) / (b.second - a.second)), maxy); + + case 3: // left + return std::pair(minx, (a.second + (double) (b.second - a.second) * (minx - a.first) / (b.first - a.first))); + } + + fprintf(stderr, "internal error intersecting\n"); + exit(EXIT_FAILURE); +} + +// http://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm +static std::vector> clip_poly1(std::vector> &geom, + long long minx, long long miny, long long maxx, long long maxy, + long long ax, long long ay, long long bx, long long by, drawvec &edge_nodes, + bool prevent_simplify_shared_nodes) { + std::vector> out = geom; + + for (int edge = 0; edge < 4; edge++) { + if (out.size() > 0) { + std::vector> in = out; + out.resize(0); + + std::pair S = in[in.size() - 1]; + + for (size_t e = 0; e < in.size(); e++) { + std::pair E = in[e]; + + if (!inside(S, edge, minx, miny, maxx, maxy)) { + // was outside the buffer + + if (!inside(E, edge, minx, miny, maxx, maxy)) { + // still outside the buffer + } else if (!inside(E, edge, ax, ay, bx, by)) { + // outside the tile but inside the buffer + out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge + out.push_back(E); + } else { + out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge + if (prevent_simplify_shared_nodes) { + out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary + edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); + } + out.push_back(E); + } + } else if (!inside(S, edge, ax, ay, bx, by)) { + // was inside the buffer but outside the tile edge + + if (!inside(E, edge, minx, miny, maxx, maxy)) { + // now outside the buffer + out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge + } else if (!inside(E, edge, ax, ay, bx, by)) { + // still outside the tile edge but inside the buffer + out.push_back(E); + } else { + // now inside the tile + if (prevent_simplify_shared_nodes) { + out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary + edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); + } + out.push_back(E); + } + } else { + // was inside the tile + + if (!inside(E, edge, minx, miny, maxx, maxy)) { + // now outside the buffer + if (prevent_simplify_shared_nodes) { + out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary + edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); + } + out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge + } else if (!inside(E, edge, ax, ay, bx, by)) { + // now inside the buffer but outside the tile edge + if (prevent_simplify_shared_nodes) { + out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary + edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); + } + out.push_back(E); + } else { + // still inside the tile + out.push_back(E); + } + } + + 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].first != out[out.size() - 1].first || out[0].second != out[out.size() - 1].second) { + out.push_back(out[0]); + } + + if (out.size() < 3) { + // fprintf(stderr, "Polygon degenerated to a line segment\n"); + out.clear(); + return out; + } + } + + return out; +} + +drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy, + long long ax, long long ay, long long bx, long long by, drawvec &edge_nodes, bool prevent_simplify_shared_nodes) { + drawvec out; + if (prevent_simplify_shared_nodes) { + geom = remove_noop(geom, VT_POLYGON, 0); + } + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_MOVETO) { + size_t j; + for (j = i + 1; j < geom.size(); j++) { + if (geom[j].op != VT_LINETO) { + break; + } + } + + std::vector> tmp; + for (size_t k = i; k < j; k++) { + double x = geom[k].x; + double y = geom[k].y; + tmp.emplace_back(x, y); + } + tmp = clip_poly1(tmp, minx, miny, maxx, maxy, ax, ay, bx, by, edge_nodes, prevent_simplify_shared_nodes); + if (tmp.size() > 0) { + if (tmp[0].first != tmp[tmp.size() - 1].first || tmp[0].second != tmp[tmp.size() - 1].second) { + fprintf(stderr, "Internal error: Polygon ring not closed\n"); + exit(EXIT_FAILURE); + } + } + for (size_t k = 0; k < tmp.size(); k++) { + if (k == 0) { + out.push_back(draw(VT_MOVETO, std::round(tmp[k].first), std::round(tmp[k].second))); + } else { + out.push_back(draw(VT_LINETO, std::round(tmp[k].first), std::round(tmp[k].second))); + } + } + + i = j - 1; + } else { + fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op); + exit(EXIT_FAILURE); + } + } + + return out; +} + +drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy, bool prevent_simplify_shared_nodes) { + drawvec dv; + return simple_clip_poly(geom, minx, miny, maxx, maxy, minx, miny, maxx, maxy, dv, prevent_simplify_shared_nodes); +} + +drawvec simple_clip_poly(drawvec &geom, int z, int buffer, drawvec &edge_nodes, bool prevent_simplify_shared_nodes) { + long long area = 1LL << (32 - z); + long long clip_buffer = buffer * area / 256; + + return simple_clip_poly(geom, -clip_buffer, -clip_buffer, area + clip_buffer, area + clip_buffer, + 0, 0, area, area, edge_nodes, prevent_simplify_shared_nodes); +} + +drawvec clip_point(drawvec &geom, int z, long long buffer) { + long long min = 0; + long long area = 1LL << (32 - z); + + min -= buffer * area / 256; + area += buffer * area / 256; + + return clip_point(geom, min, min, area, area); +} + +drawvec clip_point(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].x >= minx && geom[i].y >= miny && geom[i].x <= maxx && geom[i].y <= maxy) { + out.push_back(geom[i]); + } + } + + return out; +} + +drawvec clip_lines(drawvec &geom, int z, long long buffer) { + long long min = 0; + long long area = 1LL << (32 - z); + min -= buffer * area / 256; + area += buffer * area / 256; + + return clip_lines(geom, min, min, area, area); +} + +drawvec clip_lines(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (i > 0 && (geom[i - 1].op == VT_MOVETO || geom[i - 1].op == VT_LINETO) && geom[i].op == VT_LINETO) { + long long x1 = geom[i - 1].x; + long long y1 = geom[i - 1].y; + + long long x2 = geom[i - 0].x; + long long y2 = geom[i - 0].y; + + int c = clip(&x1, &y1, &x2, &y2, minx, miny, maxx, maxy); + + if (c > 1) { // clipped + out.push_back(draw(VT_MOVETO, std::round(x1), std::round(y1))); + out.push_back(draw(VT_LINETO, std::round(x2), std::round(y2))); + out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); + } else if (c == 1) { // unchanged + out.push_back(geom[i]); + } else { // clipped away entirely + out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); + } + } else { + out.push_back(geom[i]); + } + } + + return out; +} + +#define INSIDE 0 +#define LEFT 1 +#define RIGHT 2 +#define BOTTOM 4 +#define TOP 8 + +static int computeOutCode(long long x, long long y, long long xmin, long long ymin, long long xmax, long long ymax) { + int code = INSIDE; + + if (x < xmin) { + code |= LEFT; + } else if (x > xmax) { + code |= RIGHT; + } + + if (y < ymin) { + code |= BOTTOM; + } else if (y > ymax) { + code |= TOP; + } + + return code; +} + +int clip(long long *x0, long long *y0, long long *x1, long long *y1, long long xmin, long long ymin, long long xmax, long long ymax) { + int outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); + int outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); + int accept = 0; + int changed = 0; + + while (1) { + if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop + accept = 1; + break; + } else if (outcode0 & outcode1) { // Bitwise AND is not 0. Trivially reject and get out of loop + break; + } else { + // failed both tests, so calculate the line segment to clip + // from an outside point to an intersection with clip edge + long long x = *x0, y = *y0; + + // At least one endpoint is outside the clip rectangle; pick it. + int outcodeOut = outcode0 ? outcode0 : outcode1; + + // XXX truncating division + + // Now find the intersection point; + // use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0) + if (outcodeOut & TOP) { // point is above the clip rectangle + x = *x0 + (*x1 - *x0) * (ymax - *y0) / (*y1 - *y0); + y = ymax; + } else if (outcodeOut & BOTTOM) { // point is below the clip rectangle + x = *x0 + (*x1 - *x0) * (ymin - *y0) / (*y1 - *y0); + y = ymin; + } else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle + y = *y0 + (*y1 - *y0) * (xmax - *x0) / (*x1 - *x0); + x = xmax; + } else if (outcodeOut & LEFT) { // point is to the left of clip rectangle + y = *y0 + (*y1 - *y0) * (xmin - *x0) / (*x1 - *x0); + x = xmin; + } + + // Now we move outside point to intersection point to clip + // and get ready for next pass. + if (outcodeOut == outcode0) { + *x0 = x; + *y0 = y; + outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); + changed = 1; + } else { + *x1 = x; + *y1 = y; + outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); + changed = 1; + } + } + } + + if (accept == 0) { + return 0; + } else { + return changed + 1; + } +} + std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int ny, int detail, int buffer, std::set const &keep) { mvt_tile tile, outtile; @@ -74,8 +813,8 @@ std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int if (t == VT_LINE) { geom = clip_lines(geom, nz, buffer); } else if (t == VT_POLYGON) { - drawvec dv; - geom = simple_clip_poly(geom, nz, buffer, dv); + drawvec dv; + geom = simple_clip_poly(geom, nz, buffer, dv, false); } else if (t == VT_POINT) { geom = clip_point(geom, nz, buffer); } diff --git a/geometry.cpp b/geometry.cpp index 9994c2bd..362c9d48 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -81,374 +81,6 @@ drawvec decode_geometry(char **meta, int z, unsigned tx, unsigned ty, long long return out; } -// @@@ -void to_tile_scale(drawvec &geom, int z, int detail) { - if (32 - detail - z < 0) { - for (size_t i = 0; i < geom.size(); i++) { - 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)))); - } - } else { - for (size_t i = 0; i < geom.size(); i++) { - 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 *= (1LL << (32 - detail - z)); - d.y *= (1LL << (32 - detail - z)); - out.push_back(d); - } - return out; -} - -drawvec remove_noop(drawvec geom, int type, int shift) { - // first pass: remove empty linetos - - long long ox = 0, oy = 0; - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - long long nx = std::round((double) geom[i].x / (1LL << shift)); - long long ny = std::round((double) geom[i].y / (1LL << shift)); - - if (geom[i].op == VT_LINETO && nx == ox && ny == oy) { - continue; - } - - if (geom[i].op == VT_CLOSEPATH) { - out.push_back(geom[i]); - } else { /* moveto or lineto */ - out.push_back(geom[i]); - ox = nx; - oy = ny; - } - } - - // second pass: remove unused movetos - - if (type != VT_POINT) { - geom = out; - out.resize(0); - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - if (i + 1 >= geom.size()) { - // followed by end-of-geometry: not needed - continue; - } - - if (geom[i + 1].op == VT_MOVETO) { - // followed by another moveto: not needed - continue; - } - - if (geom[i + 1].op == VT_CLOSEPATH) { - // followed by closepath: not possible - fprintf(stderr, "Shouldn't happen\n"); - i++; // also remove unused closepath - continue; - } - } - - out.push_back(geom[i]); - } - } - - // second pass: remove empty movetos - - if (type == VT_LINE) { - geom = out; - out.resize(0); - - for (size_t i = 0; i < geom.size(); i++) { - if (i > 1 && geom[i].op == VT_MOVETO) { - if (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; - } - } - - out.push_back(geom[i]); - } - } - - return out; -} - -double get_area_scaled(const drawvec &geom, size_t i, size_t j) { - const double max_exact_double = (double) ((1LL << 53) - 1); - - // keep scaling the geometry down until we can calculate its area without overflow - for (long long scale = 2; scale < (1LL << 30); scale *= 2) { - long long bx = geom[i].x; - long long by = geom[i].y; - bool again = false; - - // https://en.wikipedia.org/wiki/Shoelace_formula - double area = 0; - for (size_t k = i; k < j; k++) { - area += (double) ((geom[k].x - bx) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].y - by) / scale); - if (std::fabs(area) >= max_exact_double) { - again = true; - break; - } - area -= (double) ((geom[k].y - by) / scale) * (double) ((geom[i + ((k - i + 1) % (j - i))].x - bx) / scale); - if (std::fabs(area) >= max_exact_double) { - again = true; - break; - } - } - - if (again) { - continue; - } else { - area /= 2; - return area * scale * scale; - } - } - - fprintf(stderr, "get_area_scaled: can't happen\n"); - exit(EXIT_IMPOSSIBLE); -} - -double get_area(const drawvec &geom, size_t i, size_t j) { - const double max_exact_double = (double) ((1LL << 53) - 1); - - // Coordinates in `geom` are 40-bit integers, so there is no good way - // to multiply them without possible precision loss. Since they probably - // do not use the full precision, shift them nearer to the origin so - // their product is more likely to be exactly representable as a double. - // - // (In practice they are actually 34-bit integers: 32 bits for the - // Mercator world plane, plus another two bits so features can stick - // off either the left or right side. But that is still too many bits - // for the product to fit either in a 64-bit long long or in a - // double where the largest exact integer is 2^53.) - // - // If the intermediate calculation still exceeds 2^53, start trying to - // recalculate the area by scaling down the geometry. This will not - // produce as precise an area, but it will still be close, and the - // sign will be correct, which is more important, since the sign - // determines the winding order of the rings. We can then use that - // sign with this generally more precise area calculation. - - long long bx = geom[i].x; - long long by = geom[i].y; - - // https://en.wikipedia.org/wiki/Shoelace_formula - double area = 0; - bool overflow = false; - for (size_t k = i; k < j; k++) { - area += (double) (geom[k].x - bx) * (double) (geom[i + ((k - i + 1) % (j - i))].y - by); - if (std::fabs(area) >= max_exact_double) { - overflow = true; - } - area -= (double) (geom[k].y - by) * (double) (geom[i + ((k - i + 1) % (j - i))].x - bx); - if (std::fabs(area) >= max_exact_double) { - overflow = true; - } - } - area /= 2; - - if (overflow) { - double scaled_area = get_area_scaled(geom, i, j); - if ((area < 0 && scaled_area > 0) || (area > 0 && scaled_area < 0)) { - area = -area; - } - } - - return area; -} - -double get_mp_area(drawvec &geom) { - double ret = 0; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - ret += get_area(geom, i, j); - i = j - 1; - } - } - - return ret; -} - -static void decode_clipped(mapbox::geometry::multi_polygon &t, drawvec &out, double scale) { - out.clear(); - - for (size_t i = 0; i < t.size(); i++) { - for (size_t j = 0; j < t[i].size(); j++) { - drawvec ring; - - for (size_t k = 0; k < t[i][j].size(); k++) { - ring.push_back(draw((k == 0) ? VT_MOVETO : VT_LINETO, std::round(t[i][j][k].x / scale), std::round(t[i][j][k].y / scale))); - } - - if (ring.size() > 0 && ring[ring.size() - 1] != ring[0]) { - fprintf(stderr, "Had to close ring\n"); - ring.push_back(draw(VT_LINETO, ring[0].x, ring[0].y)); - } - - double area = get_area(ring, 0, ring.size()); - - if ((j == 0 && area < 0) || (j != 0 && area > 0)) { - fprintf(stderr, "Ring area has wrong sign: %f for %zu\n", area, j); - exit(EXIT_IMPOSSIBLE); - } - - for (size_t k = 0; k < ring.size(); k++) { - out.push_back(ring[k]); - } - } - } -} - -drawvec clean_or_clip_poly(drawvec &geom, int z, int buffer, bool clip, bool try_scaling) { - geom = remove_noop(geom, VT_POLYGON, 0); - mapbox::geometry::multi_polygon result; - - double scale = 16.0; - if (!try_scaling) { - scale = 1.0; - } - - bool again = true; - while (again) { - mapbox::geometry::wagyu::wagyu wagyu; - again = false; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - if (j >= i + 4) { - mapbox::geometry::linear_ring lr; - - for (size_t k = i; k < j; k++) { - lr.push_back(mapbox::geometry::point(geom[k].x * scale, geom[k].y * scale)); - } - - if (lr.size() >= 3) { - wagyu.add_ring(lr); - } - } - - i = j - 1; - } - } - - if (clip) { - long long area = 0xFFFFFFFF; - if (z != 0) { - area = 1LL << (32 - z); - } - long long clip_buffer = buffer * area / 256; - - mapbox::geometry::linear_ring lr; - - lr.push_back(mapbox::geometry::point(scale * -clip_buffer, scale * -clip_buffer)); - lr.push_back(mapbox::geometry::point(scale * -clip_buffer, scale * (area + clip_buffer))); - lr.push_back(mapbox::geometry::point(scale * (area + clip_buffer), scale * (area + clip_buffer))); - lr.push_back(mapbox::geometry::point(scale * (area + clip_buffer), scale * -clip_buffer)); - lr.push_back(mapbox::geometry::point(scale * -clip_buffer, scale * -clip_buffer)); - - wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip); - } - - try { - result.clear(); - wagyu.execute(mapbox::geometry::wagyu::clip_type_union, result, mapbox::geometry::wagyu::fill_type_positive, mapbox::geometry::wagyu::fill_type_positive); - } catch (std::runtime_error &e) { - FILE *f = fopen("/tmp/wagyu.log", "w"); - fprintf(f, "%s\n", e.what()); - fprintf(stderr, "%s\n", e.what()); - fprintf(f, "["); - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - if (j >= i + 4) { - mapbox::geometry::linear_ring lr; - - if (i != 0) { - fprintf(f, ","); - } - fprintf(f, "["); - - for (size_t k = i; k < j; k++) { - lr.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); - if (k != i) { - fprintf(f, ","); - } - fprintf(f, "[%lld,%lld]", geom[k].x, geom[k].y); - } - - fprintf(f, "]"); - - if (lr.size() >= 3) { - } - } - - i = j - 1; - } - } - - fprintf(f, "]"); - fprintf(f, "\n\n\n\n\n"); - - fclose(f); - fprintf(stderr, "Internal error: Polygon cleaning failed. Log in /tmp/wagyu.log\n"); - exit(EXIT_IMPOSSIBLE); - } - - if (scale != 1) { - for (auto const &outer : result) { - for (auto const &ring : outer) { - for (auto const &p : ring) { - if (p.x / scale != std::round(p.x / scale) || - p.y / scale != std::round(p.y / scale)) { - scale = 1; - again = true; - break; - } - } - } - } - } - } - - drawvec ret; - decode_clipped(result, ret, scale); - return ret; -} -// @@@ - /* pnpoly: Copyright (c) 1970-2003, Wm. Randolph Franklin @@ -556,234 +188,6 @@ void check_polygon(drawvec &geom) { } } -// @@@ -drawvec close_poly(drawvec &geom) { - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - if (j - 1 > i) { - if (geom[j - 1].x != geom[i].x || geom[j - 1].y != geom[i].y) { - fprintf(stderr, "Internal error: polygon not closed\n"); - } - } - - for (size_t n = i; n < j - 1; n++) { - out.push_back(geom[n]); - } - out.push_back(draw(VT_CLOSEPATH, 0, 0)); - - i = j - 1; - } - } - - return out; -} - -static bool inside(std::pair d, int edge, long long minx, long long miny, long long maxx, long long maxy) { - switch (edge) { - case 0: // top - return d.second > miny; - - case 1: // right - return d.first < maxx; - - case 2: // bottom - return d.second < maxy; - - case 3: // left - return d.first > minx; - } - - fprintf(stderr, "internal error inside\n"); - exit(EXIT_FAILURE); -} - -static std::pair intersect(std::pair a, std::pair b, int edge, long long minx, long long miny, long long maxx, long long maxy) { - switch (edge) { - case 0: // top - return std::pair((a.first + (double) (b.first - a.first) * (miny - a.second) / (b.second - a.second)), miny); - - case 1: // right - return std::pair(maxx, (a.second + (double) (b.second - a.second) * (maxx - a.first) / (b.first - a.first))); - - case 2: // bottom - return std::pair((a.first + (double) (b.first - a.first) * (maxy - a.second) / (b.second - a.second)), maxy); - - case 3: // left - return std::pair(minx, (a.second + (double) (b.second - a.second) * (minx - a.first) / (b.first - a.first))); - } - - fprintf(stderr, "internal error intersecting\n"); - exit(EXIT_FAILURE); -} - -// http://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm -static std::vector> clip_poly1(std::vector> &geom, - long long minx, long long miny, long long maxx, long long maxy, - long long ax, long long ay, long long bx, long long by, drawvec &edge_nodes) { - std::vector> out = geom; - - for (int edge = 0; edge < 4; edge++) { - if (out.size() > 0) { - std::vector> in = out; - out.resize(0); - - std::pair S = in[in.size() - 1]; - - for (size_t e = 0; e < in.size(); e++) { - std::pair E = in[e]; - - if (!inside(S, edge, minx, miny, maxx, maxy)) { - // was outside the buffer - - if (!inside(E, edge, minx, miny, maxx, maxy)) { - // still outside the buffer - } else if (!inside(E, edge, ax, ay, bx, by)) { - // outside the tile but inside the buffer - out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge - out.push_back(E); - } else { - out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge - if (prevent[P_SIMPLIFY_SHARED_NODES]) { - out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary - edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); - } - out.push_back(E); - } - } else if (!inside(S, edge, ax, ay, bx, by)) { - // was inside the buffer but outside the tile edge - - if (!inside(E, edge, minx, miny, maxx, maxy)) { - // now outside the buffer - out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge - } else if (!inside(E, edge, ax, ay, bx, by)) { - // still outside the tile edge but inside the buffer - out.push_back(E); - } else { - // now inside the tile - if (prevent[P_SIMPLIFY_SHARED_NODES]) { - out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary - edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); - } - out.push_back(E); - } - } else { - // was inside the tile - - if (!inside(E, edge, minx, miny, maxx, maxy)) { - // now outside the buffer - if (prevent[P_SIMPLIFY_SHARED_NODES]) { - out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary - edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); - } - out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy)); // on buffer edge - } else if (!inside(E, edge, ax, ay, bx, by)) { - // now inside the buffer but outside the tile edge - if (prevent[P_SIMPLIFY_SHARED_NODES]) { - out.push_back(intersect(S, E, edge, ax, ay, bx, by)); // on tile boundary - edge_nodes.push_back(draw(VT_MOVETO, std::round(out.back().first), std::round(out.back().second))); - } - out.push_back(E); - } else { - // still inside the tile - out.push_back(E); - } - } - - 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].first != out[out.size() - 1].first || out[0].second != out[out.size() - 1].second) { - out.push_back(out[0]); - } - - if (out.size() < 3) { - // fprintf(stderr, "Polygon degenerated to a line segment\n"); - out.clear(); - return out; - } - } - - return out; -} - -drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy, - long long ax, long long ay, long long bx, long long by, drawvec &edge_nodes) { - drawvec out; - if (prevent[P_SIMPLIFY_SHARED_NODES]) { - geom = remove_noop(geom, VT_POLYGON, 0); - } - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_MOVETO) { - size_t j; - for (j = i + 1; j < geom.size(); j++) { - if (geom[j].op != VT_LINETO) { - break; - } - } - - std::vector> tmp; - for (size_t k = i; k < j; k++) { - double x = geom[k].x; - double y = geom[k].y; - tmp.emplace_back(x, y); - } - tmp = clip_poly1(tmp, minx, miny, maxx, maxy, ax, ay, bx, by, edge_nodes); - if (tmp.size() > 0) { - if (tmp[0].first != tmp[tmp.size() - 1].first || tmp[0].second != tmp[tmp.size() - 1].second) { - fprintf(stderr, "Internal error: Polygon ring not closed\n"); - exit(EXIT_FAILURE); - } - } - for (size_t k = 0; k < tmp.size(); k++) { - if (k == 0) { - out.push_back(draw(VT_MOVETO, std::round(tmp[k].first), std::round(tmp[k].second))); - } else { - out.push_back(draw(VT_LINETO, std::round(tmp[k].first), std::round(tmp[k].second))); - } - } - - i = j - 1; - } else { - fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op); - exit(EXIT_FAILURE); - } - } - - return out; -} - -drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec dv; - return simple_clip_poly(geom, minx, miny, maxx, maxy, minx, miny, maxx, maxy, dv); -} - -drawvec simple_clip_poly(drawvec &geom, int z, int buffer, drawvec &edge_nodes) { - long long area = 1LL << (32 - z); - long long clip_buffer = buffer * area / 256; - - return simple_clip_poly(geom, -clip_buffer, -clip_buffer, area + clip_buffer, area + clip_buffer, - 0, 0, area, area, edge_nodes); -} -// @@@ - drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *still_needs_simplification, bool *reduced_away, double *accum_area, serial_feature *this_feature, serial_feature *tiny_feature) { drawvec out; const double pixel = (1LL << (32 - detail - z)) * (double) tiny_polygon_size; @@ -904,30 +308,6 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *still_needs_sim return out; } -// @@@ -drawvec clip_point(drawvec &geom, int z, long long buffer) { - long long min = 0; - long long area = 1LL << (32 - z); - - min -= buffer * area / 256; - area += buffer * area / 256; - - return clip_point(geom, min, min, area, area); -} - -drawvec clip_point(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].x >= minx && geom[i].y >= miny && geom[i].x <= maxx && geom[i].y <= maxy) { - out.push_back(geom[i]); - } - } - - return out; -} -// @@@ - int quick_check(long long *bbox, int z, long long buffer) { long long min = 0; long long area = 1LL << (32 - z); @@ -966,47 +346,6 @@ bool point_within_tile(long long x, long long y, int z) { return x >= 0 && y >= 0 && x < area && y < area; } -// @@@ -drawvec clip_lines(drawvec &geom, int z, long long buffer) { - long long min = 0; - long long area = 1LL << (32 - z); - min -= buffer * area / 256; - area += buffer * area / 256; - - return clip_lines(geom, min, min, area, area); -} - -drawvec clip_lines(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (i > 0 && (geom[i - 1].op == VT_MOVETO || geom[i - 1].op == VT_LINETO) && geom[i].op == VT_LINETO) { - long long x1 = geom[i - 1].x; - long long y1 = geom[i - 1].y; - - long long x2 = geom[i - 0].x; - long long y2 = geom[i - 0].y; - - int c = clip(&x1, &y1, &x2, &y2, minx, miny, maxx, maxy); - - if (c > 1) { // clipped - out.push_back(draw(VT_MOVETO, std::round(x1), std::round(y1))); - out.push_back(draw(VT_LINETO, std::round(x2), std::round(y2))); - out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); - } else if (c == 1) { // unchanged - out.push_back(geom[i]); - } else { // clipped away entirely - out.push_back(draw(VT_MOVETO, geom[i].x, geom[i].y)); - } - } else { - out.push_back(geom[i]); - } - } - - return out; -} -// @@@ - static long long square_distance_from_line_fp(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) { long long p2x = segB_x - segA_x; long long p2y = segB_y - segA_y; @@ -1515,14 +854,14 @@ std::vector chop_polygon(std::vector &geoms) { if (maxy - miny > maxx - minx) { // printf("clipping y to %lld %lld %lld %lld\n", minx, miny, maxx, midy); - c1 = simple_clip_poly(geoms[i], minx, miny, maxx, midy); + c1 = simple_clip_poly(geoms[i], minx, miny, maxx, midy, prevent[P_SIMPLIFY_EDGE_NODES]); // printf(" and %lld %lld %lld %lld\n", minx, midy, maxx, maxy); - c2 = simple_clip_poly(geoms[i], minx, midy, maxx, maxy); + c2 = simple_clip_poly(geoms[i], minx, midy, maxx, maxy, prevent[P_SIMPLIFY_EDGE_NODES]); } else { // printf("clipping x to %lld %lld %lld %lld\n", minx, miny, midx, maxy); - c1 = simple_clip_poly(geoms[i], minx, miny, midx, maxy); + c1 = simple_clip_poly(geoms[i], minx, miny, midx, maxy, prevent[P_SIMPLIFY_EDGE_NODES]); // printf(" and %lld %lld %lld %lld\n", midx, midy, maxx, maxy); - c2 = simple_clip_poly(geoms[i], midx, miny, maxx, maxy); + c2 = simple_clip_poly(geoms[i], midx, miny, maxx, maxy, prevent[P_SIMPLIFY_EDGE_NODES]); } if (c1.size() >= geoms[i].size()) { @@ -1551,93 +890,6 @@ std::vector chop_polygon(std::vector &geoms) { } #endif -// @@@ -#define INSIDE 0 -#define LEFT 1 -#define RIGHT 2 -#define BOTTOM 4 -#define TOP 8 - -static int computeOutCode(long long x, long long y, long long xmin, long long ymin, long long xmax, long long ymax) { - int code = INSIDE; - - if (x < xmin) { - code |= LEFT; - } else if (x > xmax) { - code |= RIGHT; - } - - if (y < ymin) { - code |= BOTTOM; - } else if (y > ymax) { - code |= TOP; - } - - return code; -} - -int clip(long long *x0, long long *y0, long long *x1, long long *y1, long long xmin, long long ymin, long long xmax, long long ymax) { - int outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); - int outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); - int accept = 0; - int changed = 0; - - while (1) { - if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop - accept = 1; - break; - } else if (outcode0 & outcode1) { // Bitwise AND is not 0. Trivially reject and get out of loop - break; - } else { - // failed both tests, so calculate the line segment to clip - // from an outside point to an intersection with clip edge - long long x = *x0, y = *y0; - - // At least one endpoint is outside the clip rectangle; pick it. - int outcodeOut = outcode0 ? outcode0 : outcode1; - - // XXX truncating division - - // Now find the intersection point; - // use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0) - if (outcodeOut & TOP) { // point is above the clip rectangle - x = *x0 + (*x1 - *x0) * (ymax - *y0) / (*y1 - *y0); - y = ymax; - } else if (outcodeOut & BOTTOM) { // point is below the clip rectangle - x = *x0 + (*x1 - *x0) * (ymin - *y0) / (*y1 - *y0); - y = ymin; - } else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle - y = *y0 + (*y1 - *y0) * (xmax - *x0) / (*x1 - *x0); - x = xmax; - } else if (outcodeOut & LEFT) { // point is to the left of clip rectangle - y = *y0 + (*y1 - *y0) * (xmin - *x0) / (*x1 - *x0); - x = xmin; - } - - // Now we move outside point to intersection point to clip - // and get ready for next pass. - if (outcodeOut == outcode0) { - *x0 = x; - *y0 = y; - outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); - changed = 1; - } else { - *x1 = x; - *y1 = y; - outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); - changed = 1; - } - } - } - - if (accept == 0) { - return 0; - } else { - return changed + 1; - } -} -// @@@ - drawvec stairstep(drawvec &geom, int z, int detail) { drawvec out; double scale = 1 << (32 - detail - z); diff --git a/geometry.hpp b/geometry.hpp index c7bb7c4c..23395b0f 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -66,7 +66,6 @@ drawvec from_tile_scale(drawvec const &geom, int z, int detail); drawvec remove_noop(drawvec geom, int type, int shift); drawvec clip_point(drawvec &geom, int z, long long buffer); drawvec clean_or_clip_poly(drawvec &geom, int z, int buffer, bool clip, bool try_scaling); -drawvec simple_clip_poly(drawvec &geom, int z, int buffer, drawvec &shared_nodes); drawvec close_poly(drawvec &geom); drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *still_needs_simplification, bool *reduced_away, double *accum_area, serial_feature *this_feature, serial_feature *tiny_feature); int clip(long long *x0, long long *y0, long long *x1, long long *y1, long long xmin, long long ymin, long long xmax, long long ymax); @@ -84,9 +83,10 @@ double get_mp_area(drawvec &geom); drawvec polygon_to_anchor(const drawvec &geom); drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigned long long label_point); -drawvec simple_clip_poly(drawvec &geom, long long x1, long long y1, long long x2, long long y2); +drawvec simple_clip_poly(drawvec &geom, int z, int buffer, drawvec &shared_nodes, bool prevent_simplify_shared_nodes); +drawvec simple_clip_poly(drawvec &geom, long long x1, long long y1, long long x2, long long y2, bool prevent_simplify_shared_nodes); drawvec simple_clip_poly(drawvec &geom, long long x1, long long y1, long long x2, long long y2, - long long ax, long long ay, long long bx, long long by, drawvec &shared_nodes); + long long ax, long long ay, long long bx, long long by, drawvec &shared_nodes, bool prevent_simplify_shared_nodes); drawvec clip_lines(drawvec &geom, long long x1, long long y1, long long x2, long long y2); drawvec clip_point(drawvec &geom, long long x1, long long y1, long long x2, long long y2); void visvalingam(drawvec &ls, size_t start, size_t end, double threshold, size_t retain); diff --git a/serial.cpp b/serial.cpp index 76b9adc9..a462c920 100644 --- a/serial.cpp +++ b/serial.cpp @@ -443,7 +443,7 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) { for (auto &c : clipbboxes) { if (sf.t == VT_POLYGON) { - scaled_geometry = simple_clip_poly(scaled_geometry, SHIFT_RIGHT(c.minx), SHIFT_RIGHT(c.miny), SHIFT_RIGHT(c.maxx), SHIFT_RIGHT(c.maxy)); + scaled_geometry = simple_clip_poly(scaled_geometry, SHIFT_RIGHT(c.minx), SHIFT_RIGHT(c.miny), SHIFT_RIGHT(c.maxx), SHIFT_RIGHT(c.maxy), prevent[P_SIMPLIFY_SHARED_NODES]); } else if (sf.t == VT_LINE) { scaled_geometry = clip_lines(scaled_geometry, SHIFT_RIGHT(c.minx), SHIFT_RIGHT(c.miny), SHIFT_RIGHT(c.maxx), SHIFT_RIGHT(c.maxy)); scaled_geometry = remove_noop(scaled_geometry, sf.t, 0); diff --git a/tile.cpp b/tile.cpp index a4e33fae..8ced5476 100644 --- a/tile.cpp +++ b/tile.cpp @@ -1435,7 +1435,7 @@ bool clip_to_tile(serial_feature &sf, int z, long long buffer) { clipped = clip_lines(sf.geometry, z, buffer); } if (sf.t == VT_POLYGON) { - clipped = simple_clip_poly(sf.geometry, z, buffer, sf.edge_nodes); + clipped = simple_clip_poly(sf.geometry, z, buffer, sf.edge_nodes, prevent[P_SIMPLIFY_SHARED_NODES]); } if (sf.t == VT_POINT) { clipped = clip_point(sf.geometry, z, buffer);