Put all of this back in geometry.cpp for conflict resolution

This commit is contained in:
Erica Fischer
2023-09-07 16:23:30 -07:00
parent e6d05bc317
commit d09bae801c
2 changed files with 576 additions and 560 deletions
-560
View File
@@ -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<long long> min(minx, miny);
mapbox::geometry::point<long long> max(maxx, maxy);
mapbox::geometry::box<long long> 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<long long> ring;
for (size_t k = i; k < j; k++) {
ring.push_back(mapbox::geometry::point<long long>(geom[k].x, geom[k].y));
}
mapbox::geometry::linear_ring<long long> 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<long long> &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<long long> 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<long long> lr;
for (size_t k = i; k < j; k++) {
lr.push_back(mapbox::geometry::point<long long>(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<long long> lr;
lr.push_back(mapbox::geometry::point<long long>(-clip_buffer, -clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(-clip_buffer, area + clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(area + clip_buffer, area + clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(area + clip_buffer, -clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(-clip_buffer, -clip_buffer));
wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip);
}
mapbox::geometry::multi_polygon<long long> 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<long long> lr;
if (i != 0) {
fprintf(f, ",");
}
fprintf(f, "[");
for (size_t k = i; k < j; k++) {
lr.push_back(mapbox::geometry::point<long long>(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<std::string> const &keep) {
mvt_tile tile, outtile;
+576
View File
@@ -12,6 +12,8 @@
#include <sqlite3.h>
#include <mapbox/geometry/point.hpp>
#include <mapbox/geometry/multi_polygon.hpp>
#include <mapbox/geometry/wagyu/wagyu.hpp>
#include <mapbox/geometry/wagyu/quick_clip.hpp>
#include <mapbox/geometry/snap_rounding.hpp>
#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<long long> &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<long long> 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<long long> lr;
for (size_t k = i; k < j; k++) {
lr.push_back(mapbox::geometry::point<long long>(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<long long> lr;
lr.push_back(mapbox::geometry::point<long long>(-clip_buffer, -clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(-clip_buffer, area + clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(area + clip_buffer, area + clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(area + clip_buffer, -clip_buffer));
lr.push_back(mapbox::geometry::point<long long>(-clip_buffer, -clip_buffer));
wagyu.add_ring(lr, mapbox::geometry::wagyu::polygon_type_clip);
}
mapbox::geometry::multi_polygon<long long> 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<long long> lr;
if (i != 0) {
fprintf(f, ",");
}
fprintf(f, "[");
for (size_t k = i; k < j; k++) {
lr.push_back(mapbox::geometry::point<long long>(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<long long> min(minx, miny);
mapbox::geometry::point<long long> max(maxx, maxy);
mapbox::geometry::box<long long> 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<long long> ring;
for (size_t k = i; k < j; k++) {
ring.push_back(mapbox::geometry::point<long long>(geom[k].x, geom[k].y));
}
mapbox::geometry::linear_ring<long long> 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<drawvec> chop_polygon(std::vector<drawvec> &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);