From 1ddaa92166372d3e3d3553ec29f159347e0943bf Mon Sep 17 00:00:00 2001 From: Eric Fischer Date: Fri, 24 Aug 2018 15:56:12 -0700 Subject: [PATCH] Work on generalizing to arbitrarily many dimensions --- geometry.cpp | 45 ++++++++++++++++++++++--------- geometry.hpp | 8 +++--- mvt.cpp | 37 +++++++++++++++++--------- mvt.hpp | 4 +-- plugin.cpp | 2 +- read_json.cpp | 10 ++++--- serial.cpp | 7 +++-- tile.cpp | 17 +++++++++--- write_json.cpp | 72 +++++++++++++++++--------------------------------- 9 files changed, 112 insertions(+), 90 deletions(-) diff --git a/geometry.cpp b/geometry.cpp index df61008..f240fd2 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -22,7 +22,7 @@ #include "options.hpp" static int pnpoly(drawvec &vert, size_t start, size_t nvert, long long testx, long long testy); -static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax, bool *changed0, bool *changed1, double *el0, double *el1); +static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax, bool *changed0, bool *changed1, std::vector *el0, std::vector *el1); drawvec decode_geometry(FILE *meta, std::atomic *geompos, int z, unsigned tx, unsigned ty, long long *bbox, unsigned initial_x, unsigned initial_y) { drawvec out; @@ -81,7 +81,13 @@ drawvec decode_geometry(FILE *meta, std::atomic *geompos, int z, unsi d.y = wwy; if (op & VT_NODE_3D) { - deserialize_double_io(meta, &d.elevation, geompos); + unsigned long long n; + deserialize_ulong_long_io(meta, &n, geompos); + + d.elevations.resize(n); + for (size_t i = 0; i < n; i++) { + deserialize_double_io(meta, &d.elevations[i], geompos); + } } if (op & VT_NODE_ATTRIB) { deserialize_string_io(meta, d.attributes, geompos); @@ -666,11 +672,11 @@ drawvec clip_lines(drawvec &geom, int z, long long buffer) { 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 e1 = geom[i - 1].elevation; + std::vector e1 = geom[i - 1].elevations; double x2 = geom[i - 0].x; double y2 = geom[i - 0].y; - double e2 = geom[i - 0].elevation; + std::vector e2 = geom[i - 0].elevations; bool changed0 = false, changed1 = false; int c = clip(&x1, &y1, &x2, &y2, min, min, area, area, &changed0, &changed1, &e1, &e2); @@ -800,11 +806,11 @@ drawvec impose_tile_boundaries(drawvec &geom, long long extent) { if (i > 0 && geom[i].op == VT_LINETO && (geom[i - 1].op == VT_MOVETO || geom[i - 1].op == VT_LINETO)) { double x1 = geom[i - 1].x; double y1 = geom[i - 1].y; - double e1 = geom[i - 1].elevation; + std::vector e1 = geom[i - 1].elevations; double x2 = geom[i - 0].x; double y2 = geom[i - 0].y; - double e2 = geom[i - 0].elevation; + std::vector e2 = geom[i - 0].elevations; bool changed0 = false, changed1 = false; int c = clip(&x1, &y1, &x2, &y2, 0, 0, extent, extent, &changed0, &changed1, &e1, &e2); @@ -1096,7 +1102,7 @@ static int computeOutCode(double x, double y, double xmin, double ymin, double x return code; } -static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax, bool *changed0, bool *changed1, double *e0, double *e1) { +static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, double ymin, double xmax, double ymax, bool *changed0, bool *changed1, std::vector *e0, std::vector *e1) { int outcode0 = computeOutCode(*x0, *y0, xmin, ymin, xmax, ymax); int outcode1 = computeOutCode(*x1, *y1, xmin, ymin, xmax, ymax); int accept = 0; @@ -1111,7 +1117,8 @@ static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, dou } 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, e = *e0; + double x = *x0, y = *y0; + std::vector e = *e0; // At least one endpoint is outside the clip rectangle; pick it. int outcodeOut = outcode0 ? outcode0 : outcode1; @@ -1120,20 +1127,34 @@ static int clip(double *x0, double *y0, double *x1, double *y1, double xmin, dou // 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); - e = *e0 + (*e1 - *e0) * (ymax - *y0) / (*y1 - *y0); y = ymax; + + e.resize(0); + for (size_t i = 0; i < e0->size() && i < e1->size(); i++) { + e.push_back((*e0)[i] + ((*e1)[i] - (*e0)[i]) * (ymax - *y0) / (*y1 - *y0)); + } } else if (outcodeOut & BOTTOM) { // point is below the clip rectangle x = *x0 + (*x1 - *x0) * (ymin - *y0) / (*y1 - *y0); - e = *e0 + (*e1 - *e0) * (ymin - *y0) / (*y1 - *y0); y = ymin; + + e.resize(0); + for (size_t i = 0; i < e0->size() && i < e1->size(); i++) { + e.push_back((*e0)[i] + ((*e1)[i] - (*e0)[i]) * (ymin - *y0) / (*y1 - *y0)); + } } else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle y = *y0 + (*y1 - *y0) * (xmax - *x0) / (*x1 - *x0); - e = *e0 + (*e1 - *e0) * (xmax - *x0) / (*x1 - *x0); x = xmax; + + for (size_t i = 0; i < e0->size() && i < e1->size(); i++) { + e.push_back((*e0)[i] + ((*e1)[i] - (*e0)[i]) * (xmax - *x0) / (*x1 - *x0)); + } } else if (outcodeOut & LEFT) { // point is to the left of clip rectangle y = *y0 + (*y1 - *y0) * (xmin - *x0) / (*x1 - *x0); - e = *e0 + (*e1 - *e0) * (xmin - *x0) / (*x1 - *x0); x = xmin; + + for (size_t i = 0; i < e0->size() && i < e1->size(); i++) { + e.push_back((*e0)[i] + ((*e1)[i] - (*e0)[i]) * (xmin - *x0) / (*x1 - *x0)); + } } // Now we move outside point to intersection point to clip diff --git a/geometry.hpp b/geometry.hpp index 3cb1d68..191e3f1 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -25,7 +25,7 @@ struct draw { long long x : 40; signed char op; long long y : 40; - double elevation; + std::vector elevations; std::string attributes; signed char necessary; @@ -33,15 +33,14 @@ struct draw { : x(nx), op(nop), y(ny), - elevation(NAN), necessary(0) { } - draw(int nop, long long nx, long long ny, double nel) + draw(int nop, long long nx, long long ny, std::vector nel) : x(nx), op(nop), y(ny), - elevation(nel), + elevations(nel), necessary(0) { } @@ -49,7 +48,6 @@ struct draw { : x(0), op(0), y(0), - elevation(NAN), necessary(0) { } diff --git a/mvt.cpp b/mvt.cpp index 7ca6829..c74f983 100644 --- a/mvt.cpp +++ b/mvt.cpp @@ -25,11 +25,11 @@ mvt_geometry::mvt_geometry(int nop, long long nx, long long ny) { this->y = ny; } -mvt_geometry::mvt_geometry(int nop, long long nx, long long ny, double nelevation) { +mvt_geometry::mvt_geometry(int nop, long long nx, long long ny, std::vector nelevations) { this->op = nop; this->x = nx; this->y = ny; - this->elevation = nelevation; + this->elevations = nelevations; } // https://github.com/mapbox/mapnik-vector-tile/blob/master/src/vector_tile_compression.hpp @@ -314,6 +314,7 @@ bool mvt_tile::decode(std::string &message, bool &was_compressed) { protozero::pbf_reader feature_reader(layer_reader.get_message()); mvt_feature feature; std::vector geoms; + size_t dimensions = 2; std::vector elevations; while (feature_reader.next()) { @@ -363,6 +364,12 @@ bool mvt_tile::decode(std::string &message, bool &was_compressed) { break; } + case 7: /* dimensions */ + { + dimensions = feature_reader.get_uint64(); + break; + } + default: feature_reader.skip(); break; @@ -383,11 +390,14 @@ bool mvt_tile::decode(std::string &message, bool &was_compressed) { g += 2; mvt_geometry decoded = mvt_geometry(op, px, py); - if (elevation_index < elevations.size()) { - decoded.elevation = elevations[elevation_index]; - elevation_index++; - } else { - decoded.elevation = NAN; + + for (size_t i = 2; i < dimensions; i++) { + if (elevation_index < elevations.size()) { + decoded.elevations.push_back(elevations[elevation_index]); + elevation_index++; + } else { + decoded.elevations.push_back(NAN); + } } feature.geometry.push_back(decoded); } @@ -505,7 +515,7 @@ std::string mvt_tile::encode() { for (size_t f = 0; f < layers[i].features.size(); f++) { std::string feature_string; protozero::pbf_writer feature_writer(feature_string); - bool has_elevation = false; + size_t dimensions = 2; feature_writer.add_enum(3, layers[i].features[f].type); feature_writer.add_packed_uint32(2, std::begin(layers[i].features[f].tags), std::end(layers[i].features[f].tags)); @@ -552,8 +562,8 @@ std::string mvt_tile::encode() { py = wwy; length++; - if (!std::isnan(geom[g].elevation)) { - has_elevation = true; + if (geom[g].elevations.size() + 2 > dimensions) { + dimensions = geom[g].elevations.size(); } } else if (op == mvt_closepath) { length++; @@ -569,17 +579,20 @@ std::string mvt_tile::encode() { feature_writer.add_packed_uint32(4, std::begin(geometry), std::end(geometry)); - if (has_elevation) { + if (dimensions > 2) { std::vector elevations; for (size_t g = 0; g < geom.size(); g++) { int op = geom[g].op; if (op == mvt_moveto || op == mvt_lineto) { - elevations.push_back(geom[g].elevation); + for (size_t e = 0; e < geom[g].elevations.size(); e++) { + elevations.push_back(geom[g].elevations[e]); + } } } feature_writer.add_packed_double(6, std::begin(elevations), std::end(elevations)); + feature_writer.add_uint64(7, dimensions); } layer_writer.add_message(2, feature_string); diff --git a/mvt.hpp b/mvt.hpp index 99d53cf..3da06b6 100644 --- a/mvt.hpp +++ b/mvt.hpp @@ -29,11 +29,11 @@ struct mvt_geometry { long long x = 0; long long y = 0; int /* mvt_operation */ op = 0; - double elevation = 0; + std::vector elevations; std::vector attributes; mvt_geometry(int op, long long x, long long y); - mvt_geometry(int op, long long x, long long y, double elevation); + mvt_geometry(int op, long long x, long long y, std::vector elevation); bool operator<(mvt_geometry const &s) const { if (y < s.y || (y == s.y && x < s.x)) { diff --git a/plugin.cpp b/plugin.cpp index 192fb92..c72e6dd 100644 --- a/plugin.cpp +++ b/plugin.cpp @@ -77,7 +77,7 @@ static std::vector to_feature(drawvec &geom) { std::vector out; for (size_t i = 0; i < geom.size(); i++) { - out.push_back(mvt_geometry(geom[i].op, geom[i].x, geom[i].y, geom[i].elevation)); + out.push_back(mvt_geometry(geom[i].op, geom[i].x, geom[i].y, geom[i].elevations)); } return out; diff --git a/read_json.cpp b/read_json.cpp index 7fd3acc..2df866f 100644 --- a/read_json.cpp +++ b/read_json.cpp @@ -82,12 +82,13 @@ void parse_geometry(int t, json_object *j, drawvec &out, int op, const char *fna draw d(op, x, y); - size_t maybe_attr = 3; - if (j->length > 2 && j->array[2]->type == JSON_NUMBER) { - d.elevation = j->array[2]->number; - maybe_attr = 4; + for (size_t i = 2; i < j->length; i++) { + if (j->array[i]->type == JSON_NUMBER) { + d.elevations.push_back(j->array[i]->number); + } } +#if 0 if (j->length > maybe_attr) { if (j->array[maybe_attr]->type == JSON_HASH) { char *s = json_stringify(j->array[maybe_attr]); @@ -95,6 +96,7 @@ void parse_geometry(int t, json_object *j, drawvec &out, int op, const char *fna free(s); // stringify } } +#endif out.push_back(d); } else { diff --git a/serial.cpp b/serial.cpp index 57b93ba..772cea8 100644 --- a/serial.cpp +++ b/serial.cpp @@ -206,7 +206,7 @@ static void write_geometry(drawvec const &dv, std::atomic *fpos, FILE for (size_t i = 0; i < dv.size(); i++) { if (dv[i].op == VT_MOVETO || dv[i].op == VT_LINETO) { int op = dv[i].op; - if (!std::isnan(dv[i].elevation)) { + if (dv[i].elevations.size() != 0) { op |= VT_NODE_3D; } if (dv[i].attributes.size() != 0) { @@ -218,7 +218,10 @@ static void write_geometry(drawvec const &dv, std::atomic *fpos, FILE serialize_long_long(out, dv[i].y - wy, fpos, fname); if (op & VT_NODE_3D) { - serialize_double(out, dv[i].elevation, fpos, fname); + serialize_ulong_long(out, dv[i].elevations.size(), fpos, fname); + for (size_t j = 0; j < dv[i].elevations.size(); j++) { + serialize_double(out, dv[i].elevations[j], fpos, fname); + } } if (op & VT_NODE_ATTRIB) { serialize_string(out, dv[i].attributes, fpos, fname); diff --git a/tile.cpp b/tile.cpp index 18f754a..648dd72 100644 --- a/tile.cpp +++ b/tile.cpp @@ -62,7 +62,7 @@ std::vector to_feature(drawvec &geom, mvt_layer &layer) { std::vector out; for (size_t i = 0; i < geom.size(); i++) { - mvt_geometry g(geom[i].op, geom[i].x, geom[i].y, geom[i].elevation); + mvt_geometry g(geom[i].op, geom[i].x, geom[i].y, geom[i].elevations); if (geom[i].attributes.size() != 0) { json_pull *jp = json_begin_string(geom[i].attributes.c_str()); @@ -1891,13 +1891,22 @@ long long write_tile(FILE *geoms, std::atomic *geompos_in, char *meta sf.geometry.size() == 1) { double x = (double) partials[which_partial].geoms[0][0].x * partials[which_partial].clustered; double y = (double) partials[which_partial].geoms[0][0].y * partials[which_partial].clustered; - double e = (double) partials[which_partial].geoms[0][0].elevation * partials[which_partial].clustered; + x += sf.geometry[0].x; y += sf.geometry[0].y; - e += sf.geometry[0].elevation; + partials[which_partial].geoms[0][0].x = x / (partials[which_partial].clustered + 1); partials[which_partial].geoms[0][0].y = y / (partials[which_partial].clustered + 1); - partials[which_partial].geoms[0][0].elevation = e / (partials[which_partial].clustered + 1); + + std::vector e; + for (size_t i = 0; i < partials[which_partial].geoms[0][0].elevations.size() && + i < sf.geometry[0].elevations.size(); i++) { + e.push_back(partials[which_partial].geoms[0][0].elevations[i] * partials[which_partial].clustered); + e[i] += sf.geometry[0].elevations[i]; + e[i] /= (partials[which_partial].clustered + 1); + } + + partials[which_partial].geoms[0][0].elevations = e; } preserve_attributes(arg->attribute_accum, attribute_accum_state, sf, stringpool, pool_off, partials[which_partial]); diff --git a/write_json.cpp b/write_json.cpp index be2e7cf..d6504bc 100644 --- a/write_json.cpp +++ b/write_json.cpp @@ -240,15 +240,15 @@ struct lonlat { double lat; long long x; long long y; - double elevation; + std::vector elevations; - lonlat(int nop, double nlon, double nlat, long long nx, long long ny, double nelevation) + lonlat(int nop, double nlon, double nlat, long long nx, long long ny, std::vector nelevations) : op(nop), lon(nlon), lat(nlat), x(nx), y(ny), - elevation(nelevation) { + elevations(nelevations) { } }; @@ -372,6 +372,19 @@ void stringify_val(std::string &out, mvt_feature const &feature, mvt_layer const } } +void write_coordinates(json_writer &state, lonlat const &p) { + state.json_write_array(); + + state.json_write_float(p.lon); + state.json_write_float(p.lat); + + for (size_t i = 0; i < p.elevations.size(); i++) { + state.json_write_number(p.elevations[i]); + } + + state.json_end_array(); +} + void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y, bool comma, bool name, bool zoom, bool dropped, unsigned long long index, long long sequence, long long extent, bool complain, json_writer &state) { for (size_t f = 0; f < layer.features.size(); f++) { mvt_feature const &feat = layer.features[f]; @@ -479,9 +492,9 @@ void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y double lat, lon; projection->unproject(wx, wy, 32, &lon, &lat); - ops.push_back(lonlat(op, lon, lat, px, py, feat.geometry[g].elevation)); + ops.push_back(lonlat(op, lon, lat, px, py, feat.geometry[g].elevations)); } else { - ops.push_back(lonlat(op, 0, 0, 0, 0, NAN)); + ops.push_back(lonlat(op, 0, 0, 0, 0, std::vector())); } } @@ -491,14 +504,7 @@ void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y state.json_write_string("Point"); state.json_write_string("coordinates"); - - state.json_write_array(); - state.json_write_float(ops[0].lon); - state.json_write_float(ops[0].lat); - if (!std::isnan(ops[0].elevation)) { - state.json_write_float(ops[0].elevation); - } - state.json_end_array(); + write_coordinates(state, ops[0]); } else { state.json_write_string("type"); state.json_write_string("MultiPoint"); @@ -507,13 +513,7 @@ void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y state.json_write_array(); for (size_t i = 0; i < ops.size(); i++) { - state.json_write_array(); - state.json_write_float(ops[i].lon); - state.json_write_float(ops[i].lat); - if (!std::isnan(ops[i].elevation)) { - state.json_write_float(ops[i].elevation); - } - state.json_end_array(); + write_coordinates(state, ops[i]); } state.json_end_array(); @@ -534,13 +534,7 @@ void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y state.json_write_array(); for (size_t i = 0; i < ops.size(); i++) { - state.json_write_array(); - state.json_write_float(ops[i].lon); - state.json_write_float(ops[i].lat); - if (!std::isnan(ops[i].elevation)) { - state.json_write_float(ops[i].elevation); - } - state.json_end_array(); + write_coordinates(state, ops[i]); } state.json_end_array(); @@ -556,37 +550,19 @@ void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y for (size_t i = 0; i < ops.size(); i++) { if (ops[i].op == VT_MOVETO) { if (sstate == 0) { - state.json_write_array(); - state.json_write_float(ops[i].lon); - state.json_write_float(ops[i].lat); - if (!std::isnan(ops[i].elevation)) { - state.json_write_float(ops[i].elevation); - } - state.json_end_array(); + write_coordinates(state, ops[i]); sstate = 1; } else { state.json_end_array(); state.json_write_array(); - state.json_write_array(); - state.json_write_float(ops[i].lon); - state.json_write_float(ops[i].lat); - if (!std::isnan(ops[i].elevation)) { - state.json_write_float(ops[i].elevation); - } - state.json_end_array(); + write_coordinates(state, ops[i]); sstate = 1; } } else { - state.json_write_array(); - state.json_write_float(ops[i].lon); - state.json_write_float(ops[i].lat); - if (!std::isnan(ops[i].elevation)) { - state.json_write_float(ops[i].elevation); - } - state.json_end_array(); + write_coordinates(state, ops[i]); } }