From d09bae801ca2649d105300d083f88292edecc528 Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Thu, 7 Sep 2023 16:23:30 -0700 Subject: [PATCH] Put all of this back in geometry.cpp for conflict resolution --- clip.cpp | 560 ------------------------------------------------- geometry.cpp | 576 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 576 insertions(+), 560 deletions(-) diff --git a/clip.cpp b/clip.cpp index 1abc3129..039ac4ca 100644 --- a/clip.cpp +++ b/clip.cpp @@ -8,566 +8,6 @@ #include "compression.hpp" #include "mvt.hpp" -drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { - drawvec out; - - mapbox::geometry::point min(minx, miny); - mapbox::geometry::point max(maxx, maxy); - mapbox::geometry::box bbox(min, max); - - 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; - } - } - - mapbox::geometry::linear_ring ring; - for (size_t k = i; k < j; k++) { - ring.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); - } - - mapbox::geometry::linear_ring lr = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox); - - if (lr.size() > 0) { - for (size_t k = 0; k < lr.size(); k++) { - if (k == 0) { - out.push_back(draw(VT_MOVETO, lr[k].x, lr[k].y)); - } else { - out.push_back(draw(VT_LINETO, lr[k].x, lr[k].y)); - } - } - - if (lr.size() > 0 && lr[0] != lr[lr.size() - 1]) { - out.push_back(draw(VT_LINETO, lr[0].x, lr[0].y)); - } - } - - i = j - 1; - } else { - fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op); - exit(EXIT_IMPOSSIBLE); - } - } - - return out; -} - -drawvec simple_clip_poly(drawvec &geom, int z, int buffer) { - 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); -} - -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) { - double x1 = geom[i - 1].x; - double y1 = geom[i - 1].y; - - double x2 = geom[i - 0].x; - double 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, x1, y1)); - out.push_back(draw(VT_LINETO, x2, 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(double x, double y, double xmin, double ymin, double xmax, double 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(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double 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 - double x = *x0, y = *y0; - - // At least one endpoint is outside the clip rectangle; pick it. - int outcodeOut = outcode0 ? outcode0 : outcode1; - - // 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; - } -} - -static void decode_clipped(mapbox::geometry::multi_polygon &t, drawvec &out) { - 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, t[i][j][k].x, t[i][j][k].y)); - } - - 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) { - mapbox::geometry::wagyu::wagyu wagyu; - - 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; - } - } - - 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, geom[k].y)); - } - - 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(-clip_buffer, -clip_buffer)); - lr.push_back(mapbox::geometry::point(-clip_buffer, area + clip_buffer)); - lr.push_back(mapbox::geometry::point(area + clip_buffer, area + clip_buffer)); - lr.push_back(mapbox::geometry::point(area + clip_buffer, -clip_buffer)); - lr.push_back(mapbox::geometry::point(-clip_buffer, -clip_buffer)); - - wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip); - } - - mapbox::geometry::multi_polygon result; - try { - 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); - } - - drawvec ret; - decode_clipped(result, ret); - return ret; -} - -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 x = 0, y = 0; - drawvec out; - - for (size_t i = 0; i < geom.size(); i++) { - if (geom[i].op == VT_LINETO && (long long) std::round((double) geom[i].x / (1LL << shift)) == x && (long long) std::round((double) geom[i].y / (1LL << shift)) == y) { - continue; - } - - if (geom[i].op == VT_CLOSEPATH) { - out.push_back(geom[i]); - } else { /* moveto or lineto */ - out.push_back(geom[i]); - x = std::round((double) geom[i].x / (1LL << shift)); - y = std::round((double) geom[i].y / (1LL << shift)); - } - } - - // 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()) { - continue; - } - - if (geom[i + 1].op == VT_MOVETO) { - continue; - } - - if (geom[i + 1].op == VT_CLOSEPATH) { - 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 (geom[i].op == VT_MOVETO) { - if (i > 0 && geom[i - 1].op == VT_LINETO && (long long) std::round((double) geom[i - 1].x / (1LL << shift)) == (long long) std::round((double) geom[i].x / (1LL << shift)) && (long long) std::round((double) geom[i - 1].y / (1LL << shift)) == (long long) 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; -} - -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; -} - 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; diff --git a/geometry.cpp b/geometry.cpp index 8c859e37..a6bf9b5f 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -12,6 +12,8 @@ #include #include #include +#include +#include #include #include "geometry.hpp" #include "projection.hpp" @@ -79,6 +81,342 @@ 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 x = 0, y = 0; + drawvec out; + + for (size_t i = 0; i < geom.size(); i++) { + if (geom[i].op == VT_LINETO && (long long) std::round((double) geom[i].x / (1LL << shift)) == x && (long long) std::round((double) geom[i].y / (1LL << shift)) == y) { + continue; + } + + if (geom[i].op == VT_CLOSEPATH) { + out.push_back(geom[i]); + } else { /* moveto or lineto */ + out.push_back(geom[i]); + x = std::round((double) geom[i].x / (1LL << shift)); + y = std::round((double) geom[i].y / (1LL << shift)); + } + } + + // 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()) { + continue; + } + + if (geom[i + 1].op == VT_MOVETO) { + continue; + } + + if (geom[i + 1].op == VT_CLOSEPATH) { + 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 (geom[i].op == VT_MOVETO) { + if (i > 0 && geom[i - 1].op == VT_LINETO && (long long) std::round((double) geom[i - 1].x / (1LL << shift)) == (long long) std::round((double) geom[i].x / (1LL << shift)) && (long long) std::round((double) geom[i - 1].y / (1LL << shift)) == (long long) 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) { + 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, t[i][j][k].x, t[i][j][k].y)); + } + + 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) { + mapbox::geometry::wagyu::wagyu wagyu; + + 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; + } + } + + 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, geom[k].y)); + } + + 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(-clip_buffer, -clip_buffer)); + lr.push_back(mapbox::geometry::point(-clip_buffer, area + clip_buffer)); + lr.push_back(mapbox::geometry::point(area + clip_buffer, area + clip_buffer)); + lr.push_back(mapbox::geometry::point(area + clip_buffer, -clip_buffer)); + lr.push_back(mapbox::geometry::point(-clip_buffer, -clip_buffer)); + + wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip); + } + + mapbox::geometry::multi_polygon result; + try { + 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); + } + + drawvec ret; + decode_clipped(result, ret); + return ret; +} +// @@@ + /* pnpoly: Copyright (c) 1970-2003, Wm. Randolph Franklin @@ -186,6 +524,94 @@ 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; +} +// @@@ + +// @@@ +drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) { + drawvec out; + + mapbox::geometry::point min(minx, miny); + mapbox::geometry::point max(maxx, maxy); + mapbox::geometry::box bbox(min, max); + + 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; + } + } + + mapbox::geometry::linear_ring ring; + for (size_t k = i; k < j; k++) { + ring.push_back(mapbox::geometry::point(geom[k].x, geom[k].y)); + } + + mapbox::geometry::linear_ring lr = mapbox::geometry::wagyu::quick_clip::quick_lr_clip(ring, bbox); + + if (lr.size() > 0) { + for (size_t k = 0; k < lr.size(); k++) { + if (k == 0) { + out.push_back(draw(VT_MOVETO, lr[k].x, lr[k].y)); + } else { + out.push_back(draw(VT_LINETO, lr[k].x, lr[k].y)); + } + } + + if (lr.size() > 0 && lr[0] != lr[lr.size() - 1]) { + out.push_back(draw(VT_LINETO, lr[0].x, lr[0].y)); + } + } + + i = j - 1; + } else { + fprintf(stderr, "Unexpected operation in polygon %d\n", (int) geom[i].op); + exit(EXIT_IMPOSSIBLE); + } + } + + return out; +} + +drawvec simple_clip_poly(drawvec &geom, int z, int buffer) { + 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); +} +// @@@ + 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; @@ -306,6 +732,30 @@ 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); @@ -339,6 +789,47 @@ 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) { + double x1 = geom[i - 1].x; + double y1 = geom[i - 1].y; + + double x2 = geom[i - 0].x; + double 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, x1, y1)); + out.push_back(draw(VT_LINETO, x2, 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 double square_distance_from_line(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; @@ -774,6 +1265,91 @@ std::vector chop_polygon(std::vector &geoms) { } } +// @@@ +#define INSIDE 0 +#define LEFT 1 +#define RIGHT 2 +#define BOTTOM 4 +#define TOP 8 + +static int computeOutCode(double x, double y, double xmin, double ymin, double xmax, double 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(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double 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 + double x = *x0, y = *y0; + + // At least one endpoint is outside the clip rectangle; pick it. + int outcodeOut = outcode0 ? outcode0 : outcode1; + + // 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);