Work on generalizing to arbitrarily many dimensions

This commit is contained in:
Eric Fischer 2018-08-24 15:56:12 -07:00
parent 68686e813b
commit 1ddaa92166
9 changed files with 112 additions and 90 deletions

View File

@ -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<double> *el0, std::vector<double> *el1);
drawvec decode_geometry(FILE *meta, std::atomic<long long> *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<long long> *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<double> 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<double> 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<double> 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<double> 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<double> *e0, std::vector<double> *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<double> 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

View File

@ -25,7 +25,7 @@ struct draw {
long long x : 40;
signed char op;
long long y : 40;
double elevation;
std::vector<double> 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<double> 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) {
}

37
mvt.cpp
View File

@ -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<double> 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<uint32_t> geoms;
size_t dimensions = 2;
std::vector<double> 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<double> 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);

View File

@ -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<double> elevations;
std::vector<unsigned long> 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<double> elevation);
bool operator<(mvt_geometry const &s) const {
if (y < s.y || (y == s.y && x < s.x)) {

View File

@ -77,7 +77,7 @@ static std::vector<mvt_geometry> to_feature(drawvec &geom) {
std::vector<mvt_geometry> 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;

View File

@ -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 {

View File

@ -206,7 +206,7 @@ static void write_geometry(drawvec const &dv, std::atomic<long long> *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<long long> *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);

View File

@ -62,7 +62,7 @@ std::vector<mvt_geometry> to_feature(drawvec &geom, mvt_layer &layer) {
std::vector<mvt_geometry> 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<long long> *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<double> 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]);

View File

@ -240,15 +240,15 @@ struct lonlat {
double lat;
long long x;
long long y;
double elevation;
std::vector<double> 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<double> 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<double>()));
}
}
@ -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]);
}
}