tippecanoe/vt3.cpp

781 lines
20 KiB
C++

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <sqlite3.h>
#include <getopt.h>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <zlib.h>
#include <math.h>
#include <fcntl.h>
#include <dirent.h>
#include <sys/stat.h>
#include <sys/mman.h>
#include <protozero/pbf_reader.hpp>
#include <sys/stat.h>
#include <algorithm>
#include "mvt.hpp"
#include "projection.hpp"
#include "geometry.hpp"
#include "vt3.hpp"
#include "clip.hpp"
bool check_empty_lineto(std::vector<mvt_geometry> &geom) {
for (size_t i = 1; i < geom.size(); i++) {
if (geom[i].op == mvt_lineto && geom[i - 1].x == geom[i].x && geom[i - 1].y == geom[i].y) {
return true;
}
}
return false;
}
std::vector<mvt_geometry> clip_lines(std::vector<mvt_geometry> &geom, long left, long top, long right, long bottom) {
std::vector<mvt_geometry> out;
bool inside = false;
for (size_t i = 0; i < geom.size(); i++) {
if (geom[i].op == mvt_moveto) {
if (geom[i].x >= left && geom[i].x <= right && geom[i].y >= top && geom[i].y <= bottom) {
// First point is inside.
out.push_back(geom[i]);
inside = true;
} else {
inside = false;
}
} else {
if (geom[i].x >= left && geom[i].x <= right && geom[i].y >= top && geom[i].y <= bottom) {
// Next point is inside.
if (inside) {
// Moving from inside to another inside, so just include the new point
out.push_back(geom[i]);
} else {
// Moving from outside to inside, so clip and include the inner portion
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, left, top, right, bottom);
if (c != CLIP_CHANGED_FIRST) {
fprintf(stderr, "Expected first to change\n");
exit(EXIT_FAILURE);
}
mvt_geometry phantom(mvt_moveto, round(x1), round(y1));
phantom.phantom = true;
phantom.id = 0;
out.push_back(phantom);
out.push_back(geom[i]);
if (phantom == geom[i]) {
fprintf(stderr, "clip out to in made empty lineto\n");
}
}
inside = true;
} else {
// Next point is outside.
if (inside) {
// Moving from inside to outside, so clip and include the inner portion
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, left, top, right, bottom);
if (c != CLIP_CHANGED_SECOND) {
fprintf(stderr, "Expected second to change\n");
exit(EXIT_FAILURE);
}
// Inside already drawing, so don't need the start point
// out.push_back(geom[i - 1]);
mvt_geometry phantom(mvt_lineto, round(x2), round(y2));
if (phantom == out[out.size() - 1]) {
if (out[out.size() - 1].id == 0) {
fprintf(stderr, "In to out is clipped away but has no id %lld,%lld\n", phantom.x, phantom.y);
}
} else {
phantom.phantom = true;
phantom.id = 0;
out.push_back(phantom);
}
if (out[out.size() - 2] == out[out.size() - 1]) {
fprintf(stderr, "clip in to out made empty lineto\n");
}
} else {
// Outside to outside, so discardable on both sides
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, left, top, right, bottom);
if (c != CLIP_ELIMINATED) {
mvt_geometry phantom1(mvt_moveto, round(x1), round(y1));
phantom1.phantom = true;
phantom1.id = 0;
out.push_back(phantom1);
mvt_geometry phantom2(mvt_lineto, round(x2), round(y2));
phantom2.phantom = true;
phantom2.id = 0;
out.push_back(phantom2);
if (phantom1 == phantom2) {
fprintf(stderr, "clip in discardable made empty lineto\n");
}
}
}
inside = false;
}
}
}
if (check_empty_lineto(out)) {
fprintf(stderr, "empty lineto in clip\n");
exit(EXIT_FAILURE);
}
return out;
}
static void dump(std::vector<mvt_geometry> geom) {
for (size_t i = 0; i < geom.size(); i++) {
printf("%d %lld,%lld %ld %s\n", geom[i].op, geom[i].x, geom[i].y, geom[i].id, geom[i].phantom ? "true" : "false");
}
}
static std::vector<mvt_geometry> remove_noop(std::vector<mvt_geometry> geom, std::vector<mvt_geometry> orig) {
std::vector<mvt_geometry> out;
for (size_t i = 0; i < geom.size(); i++) {
if (geom[i].op == mvt_moveto && (i + 1 >= geom.size() || geom[i + 1].op == mvt_moveto)) {
if (geom[i].id != 0) {
fprintf(stderr, "Removing a moveto with an id %ld in %zu\n", geom[i].id, geom.size());
dump(geom);
fprintf(stderr, "Orig was\n");
dump(orig);
exit(EXIT_FAILURE);
}
fprintf(stderr, "Discarded something\n");
exit(EXIT_FAILURE);
} else {
out.push_back(geom[i]);
}
}
if (check_empty_lineto(out)) {
fprintf(stderr, "empty lineto in remove noop\n");
exit(EXIT_FAILURE);
}
return out;
}
void clip_lines_to_tile(std::vector<mvt_geometry> geom, mvt_tile &tile, long x, long y, long extent, long buffer, std::vector<mvt_geometry> orig) {
geom = clip_lines(geom, x * extent - buffer, y * extent - buffer, (x + 1) * extent + buffer, (y + 1) * extent + buffer);
geom = remove_noop(geom, orig);
mvt_layer &nl = tile.layers[tile.layers.size() - 1];
mvt_feature &nf = nl.features[nl.features.size() - 1];
nf.geometry = geom;
}
void split_feature(mvt_layer const &layer, mvt_feature const &feature, std::vector<std::vector<mvt_tile>> &subtiles, size_t n) {
static long clipid_pool = 0;
const std::vector<mvt_geometry> &geom = feature.geometry;
long extent = layer.extent;
long nextent = extent / n;
if (nextent * (long) n != extent) {
fprintf(stderr, "Extent %ld doesn't subdivide evenly by %zud\n", extent, n);
exit(EXIT_FAILURE);
}
// Calculate bounding box of feature
long minx = LONG_MAX;
long miny = LONG_MAX;
long maxx = LONG_MIN;
long maxy = LONG_MIN;
for (size_t i = 0; i < geom.size(); i++) {
if (geom[i].op == mvt_moveto || geom[i].op == mvt_lineto) {
if (geom[i].x < minx) {
minx = geom[i].x;
}
if (geom[i].y < miny) {
miny = geom[i].y;
}
if (geom[i].x > maxx) {
maxx = geom[i].x;
}
if (geom[i].y > maxy) {
maxy = geom[i].y;
}
}
}
// Extend bounding box by buffer
long buffer = nextent * 5 / 256; // XXX make configurable
if (buffer <= 0) {
fprintf(stderr, "buffer was eliminated\n");
buffer = 1;
exit(EXIT_FAILURE);
}
minx -= buffer;
miny -= buffer;
maxx += buffer;
maxy += buffer;
long left = floor((double) minx / nextent);
long top = floor((double) miny / nextent);
long right = floor((double) maxx / nextent);
long bottom = floor((double) maxy / nextent);
left = std::max(left, 0L);
top = std::max(top, 0L);
right = std::min(right, (long) (n - 1));
bottom = std::min(bottom, (long) (n - 1));
// Is it bigger than one sub-tile?
// If so, generate an ID for matching
// XXX Is this right for edges on the border?
long nclipid = 0;
if (left != right || top != bottom) {
nclipid = ++clipid_pool;
}
// Set up a corresponding feature within each sub-tile
// (Empty features will have to get torn down later)
for (size_t x = 0; x < n; x++) {
for (size_t y = 0; y < n; y++) {
mvt_feature nf;
nf.tags = feature.tags;
nf.type = feature.type;
nf.id = feature.id;
nf.has_id = feature.has_id;
nf.clipid = nclipid;
mvt_tile &nt = subtiles[x][y];
mvt_layer &nl = nt.layers[nt.layers.size() - 1];
nl.features.push_back(nf);
}
}
if (feature.type == mvt_linestring) {
std::vector<mvt_geometry> ogeom = geom;
std::vector<mvt_geometry> ngeom;
long pointid = 0;
// Part 1: Assign (phantom) point IDs in the middle of any
// segments that cross from one sub-tile to another
for (size_t i = 0; i < ogeom.size(); i++) {
if (i > 0 && ogeom[i].op == mvt_lineto && (floor((double) ogeom[i].x / nextent) != floor((double) ogeom[i - 1].x / nextent))) {
long first = floor((double) ogeom[i - 1].x / nextent) * nextent;
long second = floor((double) ogeom[i].x / nextent) * nextent;
// If <= 2, one endpoint or the other is on the edge, not crossing it
if (llabs(ogeom[i].x - ogeom[i - 1].x) > 2) {
if (second >= first) {
for (long x = first + nextent; x <= second; x += nextent) {
long y = round(ogeom[i - 1].y + (ogeom[i].y - ogeom[i - 1].y) * (double) (x - ogeom[i - 1].x) / (ogeom[i].x - ogeom[i - 1].x));
mvt_geometry p(ogeom[i].op, x, y);
p.id = ++pointid;
p.phantom = true;
ngeom.push_back(p);
}
} else {
for (long x = first; x >= second + nextent; x -= nextent) {
long y = round(ogeom[i - 1].y + (ogeom[i].y - ogeom[i - 1].y) * (double) (x - ogeom[i - 1].x) / (ogeom[i].x - ogeom[i - 1].x));
mvt_geometry p(ogeom[i].op, x, y);
p.id = ++pointid;
p.phantom = true;
ngeom.push_back(p);
}
}
}
}
ngeom.push_back(ogeom[i]);
}
if (check_empty_lineto(ngeom)) {
fprintf(stderr, "empty lineto in adding crossings\n");
exit(EXIT_FAILURE);
}
// Part 1a: same, but for Y axis crossings
ogeom = ngeom;
ngeom.clear();
for (size_t i = 0; i < ogeom.size(); i++) {
if (i > 0 && ogeom[i].op == mvt_lineto && (floor((double) ogeom[i].y / nextent) != floor((double) ogeom[i - 1].y / nextent))) {
long first = floor((double) ogeom[i - 1].y / nextent) * nextent;
long second = floor((double) ogeom[i].y / nextent) * nextent;
// If <= 2, one endpoint or the other is on the edge, not crossing it
if (llabs(ogeom[i].y - ogeom[i - 1].y) > 2) {
if (second >= first) {
for (long y = first + nextent; y <= second; y += nextent) {
long x = round(ogeom[i - 1].x + (ogeom[i].x - ogeom[i - 1].x) * (double) (y - ogeom[i - 1].y) / (ogeom[i].y - ogeom[i - 1].y));
mvt_geometry p(ogeom[i].op, x, y);
p.id = ++pointid;
p.phantom = true;
ngeom.push_back(p);
}
} else {
for (long y = first; y >= second + nextent; y -= nextent) {
long x = round(ogeom[i - 1].x + (ogeom[i].x - ogeom[i - 1].x) * (double) (y - ogeom[i - 1].y) / (ogeom[i].y - ogeom[i - 1].y));
mvt_geometry p(ogeom[i].op, x, y);
p.id = ++pointid;
p.phantom = true;
ngeom.push_back(p);
}
}
}
}
ngeom.push_back(ogeom[i]);
}
if (check_empty_lineto(ngeom)) {
fprintf(stderr, "empty lineto in adding crossings\n");
exit(EXIT_FAILURE);
}
// Part 2: Assign (real) point IDs for any points that are on a sub-tile edge
for (size_t i = 0; i < ngeom.size(); i++) {
if (ngeom[i].x % nextent == 0 || ngeom[i].y % nextent == 0) {
if (ngeom[i].id == 0) {
ngeom[i].id = ++pointid;
ngeom[i].phantom = false;
}
}
}
// Part 3: Clip the geometry to each of the sub-tiles
for (size_t x = 0; x < n; x++) {
for (size_t y = 0; y < n; y++) {
clip_lines_to_tile(ngeom, subtiles[x][y], x, y, nextent, buffer, geom);
}
}
} else {
// deal with other feature types later
mvt_tile &nt = subtiles[0][0];
mvt_layer &nl = nt.layers[nt.layers.size() - 1];
mvt_feature &nf = nl.features[nl.features.size() - 1];
nf.geometry = geom;
}
}
void trim_tile(mvt_tile &tile) {
for (ssize_t i = tile.layers.size() - 1; i >= 0; i--) {
mvt_layer &layer = tile.layers[i];
for (ssize_t j = layer.features.size() - 1; j >= 0; j--) {
mvt_feature &feature = layer.features[j];
if (feature.geometry.size() == 0) {
layer.features.erase(layer.features.begin() + j);
}
}
if (layer.features.size() == 0) {
tile.layers.erase(tile.layers.begin() + i);
}
}
}
mvt_feature *add_to_tile(mvt_feature const &f, mvt_layer const &l, mvt_tile &tile, size_t n) {
size_t k;
for (k = 0; k < tile.layers.size(); k++) {
if (tile.layers[k].name == l.name) {
break;
}
}
if (k == tile.layers.size()) {
mvt_layer nl = mvt_layer();
nl.name = l.name;
nl.extent = l.extent * n;
tile.layers.push_back(nl);
}
mvt_feature nf;
nf.type = f.type;
nf.id = f.id;
nf.has_id = f.has_id;
nf.clipid = f.clipid;
for (size_t kv = 0; kv + 1 < f.tags.size(); kv += 2) {
tile.layers[k].tag(nf, l.keys[f.tags[kv]], l.values[f.tags[kv + 1]]);
}
tile.layers[k].features.push_back(nf);
return &tile.layers[k].features[tile.layers[k].features.size() - 1];
}
struct partial {
long clipid;
const mvt_feature *f;
const mvt_layer *l;
long x;
long y;
bool operator<(const partial &p) const {
return clipid < p.clipid;
}
};
void merge_partials(std::vector<partial> &partials, size_t start, size_t end, mvt_tile &tile, size_t n) {
// Pull out contiguous portions of original geometry
std::vector<std::vector<mvt_geometry>> revised;
for (size_t i = start; i < end; i++) {
std::vector<mvt_geometry> const &geom = partials[i].f->geometry;
for (size_t j = 0; j < geom.size(); j++) {
if (geom[j].op == mvt_moveto) {
size_t k;
std::vector<mvt_geometry> out;
out.push_back(geom[j]);
for (k = j + 1; k < geom.size(); k++) {
if (geom[k].op == mvt_moveto) {
break;
}
out.push_back(geom[k]);
}
if (check_empty_lineto(out)) {
fprintf(stderr, "empty lineto in buffer removal 1\n");
exit(EXIT_FAILURE);
}
// Now need to discard anything that happens in the buffer,
// because it should all be duplicated.
// Do this by location instead of by ID, because geometry is
// only tagged as phantom if it actually touches the tile
// or buffer edge, not if it just happens to be in the buffer
long minx = partials[i].l->extent * partials[i].x;
long maxx = partials[i].l->extent * (partials[i].x + 1);
long miny = partials[i].l->extent * partials[i].y;
long maxy = partials[i].l->extent * (partials[i].y + 1);
std::vector<mvt_geometry> out2;
bool within = false;
for (size_t l = 0; l < out.size(); l++) {
if (out[l].x >= minx && out[l].x <= maxx && out[l].y >= miny && out[l].y <= maxy) {
if (!within) {
out[l].op = mvt_moveto;
}
out2.push_back(out[l]);
within = true;
} else {
within = false;
}
}
if (check_empty_lineto(out2)) {
fprintf(stderr, "empty lineto in buffer removal\n");
exit(EXIT_FAILURE);
}
revised.push_back(out2);
j = k - 1;
}
}
}
mvt_feature *nf = add_to_tile(*partials[start].f, *partials[start].l, tile, n);
#if 0
for (size_t i = 0; i < revised.size(); i++) {
for (size_t j = 0; j < revised[i].size(); j++) {
nf->geometry.push_back(revised[i][j]);
}
}
#else
// Pull out individual arcs from the various geometries
// Multimap because there may be multiple arcs with id 0,
// as the start of multilinestrings. (Also multiple arcs
// with the same id on tile edges, but those duplicates
// need to be eliminated)
std::multimap<long, std::vector<mvt_geometry>> arcs;
for (size_t i = 0; i < revised.size(); i++) {
for (size_t j = 0; j < revised[i].size(); j++) {
std::vector<mvt_geometry> arc;
arc.push_back(revised[i][j]);
long id = revised[i][j].id;
size_t k;
for (k = j + 1; k < revised[i].size(); k++) {
// Or would it be better to keep movetos together here,
// as a way to chain to the next arc?
if (revised[i][k].op == mvt_moveto) {
break;
}
arc.push_back(revised[i][k]);
if (revised[i][k].id != 0) {
revised[i][k].op = mvt_moveto;
break;
}
}
if (check_empty_lineto(arc)) {
fprintf(stderr, "empty lineto in arc making\n");
exit(EXIT_FAILURE);
}
// Don't insert anything that was clipped down to just one moveto
if (arc.size() > 1) {
arcs.insert(std::pair<long, std::vector<mvt_geometry>>(id, arc));
}
j = k - 1;
}
}
#if 0
for (auto a = arcs.begin(); a != arcs.end(); ++a) {
printf("arc %ld:\n", a->first);
dump(a->second);
}
printf("---\n");
#endif
// std::map<long, std::vector<mvt_geometry>> arcs2;
while (arcs.size() > 0) {
std::vector<mvt_geometry> out = arcs.begin()->second;
arcs.erase(arcs.begin());
while (out[out.size() - 1].id != 0) {
long id = out[out.size() - 1].id;
auto found = arcs.find(id);
if (found == arcs.end()) {
break; // why would this happen?
}
std::vector<mvt_geometry> add = found->second;
arcs.erase(found);
while (true) {
auto also = arcs.find(id);
if (also == arcs.end()) {
break;
}
if (also->second != add) {
fprintf(stderr, "Found conflicting arc for %ld\n", id);
dump(add);
printf("vs\n");
dump(also->second);
exit(EXIT_FAILURE);
}
arcs.erase(also);
}
// get rid of artificial joining node if that's what it is
if (out[out.size() - 1].phantom) {
out.pop_back();
}
// skip first because it is the moveto
for (size_t i = 1; i < add.size(); i++) {
out.push_back(add[i]);
}
}
if (check_empty_lineto(out)) {
fprintf(stderr, "empty lineto in reassembly\n");
exit(EXIT_FAILURE);
}
for (size_t i = 0; i < out.size(); i++) {
nf->geometry.push_back(out[i]);
}
// arcs2.insert(std::pair<long, std::vector<mvt_geometry>>(out[0].id, out));
}
#if 0
// Outer loop because removing an arc once it is used invalidates the iterator
bool again = true;
while (again) {
again = false;
for (auto a = arcs.begin(); a != arcs.end(); a++) {
for (size_t j = 0; j < a->second.size(); j++) {
nf->geometry.push_back(a->second[j]);
}
}
}
#endif
#endif
}
mvt_tile reassemble(std::vector<std::vector<mvt_tile>> const &subtiles, size_t n) {
mvt_tile tile;
std::vector<partial> partials;
for (size_t x = 0; x < n; x++) {
for (size_t y = 0; y < n; y++) {
mvt_tile const &t = subtiles[x][y];
for (size_t i = 0; i < t.layers.size(); i++) {
mvt_layer const &l = t.layers[i];
for (size_t j = 0; j < l.features.size(); j++) {
mvt_feature const &f = l.features[j];
if (f.clipid == 0) {
mvt_feature *nf = add_to_tile(f, l, tile, n);
nf->geometry = f.geometry;
} else {
partial p;
p.clipid = f.clipid;
p.f = &f;
p.l = &l;
p.x = x;
p.y = y;
partials.push_back(p);
}
}
}
}
}
// XXX do partials
std::sort(partials.begin(), partials.end());
for (size_t i = 0; i < partials.size(); i++) {
size_t j;
for (j = i; j < partials.size(); j++) {
if (partials[j].clipid != partials[i].clipid) {
break;
}
}
merge_partials(partials, i, j, tile, n);
i = j - 1;
}
return tile;
}
mvt_tile split_and_merge(mvt_tile tile, int tile_zoom) {
// Features will be split into an NxN grid of sub-tiles,
// to be merged back together at the end,
// which should result in the original set of features
// except (perhaps) for their sequence.
size_t n = 1 << tile_zoom;
std::vector<std::vector<mvt_tile>> subtiles;
subtiles.resize(n);
for (size_t i = 0; i < n; i++) {
subtiles[i].resize(n);
}
for (size_t i = 0; i < tile.layers.size(); i++) {
mvt_layer &layer = tile.layers[i];
// Set up a corresponding layer within each sub-tile
// (Empty layers will have to get torn down later)
for (size_t x = 0; x < n; x++) {
for (size_t y = 0; y < n; y++) {
mvt_layer nl;
// Note that for simplicity this is copying *all* the
// keys and values to the sub-layers, not only the ones
// actually needed for the sub-features.
nl.version = layer.version;
nl.extent = layer.extent >> tile_zoom;
nl.name = layer.name;
nl.keys = layer.keys;
nl.values = layer.values;
subtiles[x][y].layers.push_back(nl);
}
}
// Iterate through features, copying them into sub-tiles
for (size_t j = 0; j < layer.features.size(); j++) {
mvt_feature &feature = layer.features[j];
split_feature(layer, feature, subtiles, n);
}
}
// Trim unused features from layers, layers from tiles
for (size_t x = 0; x < n; x++) {
for (size_t y = 0; y < n; y++) {
trim_tile(subtiles[x][y]);
}
}
// Write each tile to PBF
// Decode each tile back from PBF
// Recreate original tile from decoded sub-tiles
tile = reassemble(subtiles, n);
#if 0
tile.layers.clear();
for (size_t x = 0; x < n; x++) {
for (size_t y = 0; y < n; y++) {
for (size_t k = 0; k < subtiles[x][y].layers.size(); k++) {
tile.layers.push_back(subtiles[x][y].layers[k]);
}
}
}
#endif
// Verify that the original data has been recreated
return tile;
}