bool ok = repair_multi_polygon(mp);

if (!ok && verbose)

std::cerr << ("multipolygon repair failed for object " + std::to_string(objectID) + "\n");

return ok;

:

includeID(includeID),

z6OffsetDivisor(indexZoom >= CLUSTER_ZOOM ? (1 << (indexZoom - CLUSTER_ZOOM)) : 1),

objectsMutex(threadNum * 4),

objects(CLUSTER_ZOOM_AREA),

lowZoomObjects(CLUSTER_ZOOM_AREA),

objectsWithIds(CLUSTER_ZOOM_AREA),

lowZoomObjectsWithIds(CLUSTER_ZOOM_AREA),

indexZoom(indexZoom),

pointStores(threadNum),

linestringStores(threadNum),

multilinestringStores(threadNum),

multipolygonStores(threadNum),

multiPolygonClipCache(ClipCache<MultiPolygon>(threadNum, indexZoom)),

multiLinestringClipCache(ClipCache<MultiLinestring>(threadNum, indexZoom))

{

// TileDataSource can only index up to zoom 14. The caller is responsible for

// ensuring it does not use a higher zoom.

if (indexZoom > 14)

throw std::out_of_range("TileDataSource: indexZoom cannot be higher than 14, but was " + std::to_string(indexZoom));

shardBits = 0;

numShards = 1;

while(numShards < threadNum) {

shardBits++;

numShards *= 2;

}

for (int i = 0; i < threadNum; i++) {

availablePointStoreLeases.push_back(std::make_pair(i, &pointStores[i]));

availableLinestringStoreLeases.push_back(std::make_pair(i, &linestringStores[i]));

availableMultiLinestringStoreLeases.push_back(std::make_pair(i, &multilinestringStores[i]));

availableMultiPolygonStoreLeases.push_back(std::make_pair(i, &multipolygonStores[i]));

}

uint64_t finalized = 0;

for (const auto& vec : pendingSmallIndexObjects) {

for (const auto& tuple : vec) {

finalized++;

addObjectToSmallIndexUnsafe(std::get<0>(tuple), std::get<1>(tuple), std::get<2>(tuple));

}

}

std::cout << "indexed " << finalized << " contended objects" << std::endl;

finalizeObjects<OutputObjectXY>(name(), threadNum, indexZoom, objects.begin(), objects.end(), lowZoomObjects);

finalizeObjects<OutputObjectXYID>(name(), threadNum, indexZoom, objectsWithIds.begin(), objectsWithIds.end(), lowZoomObjectsWithIds);

// Pick the z6 index

const size_t z6x = index.x / z6OffsetDivisor;

const size_t z6y = index.y / z6OffsetDivisor;

if (z6x >= 64 || z6y >= 64) {

if (verbose) std::cerr << "ignoring OutputObject with invalid z" << indexZoom << " coordinates " << index.x << ", " << index.y << " (id: " << id << ")" << std::endl;

return;

}

const size_t z6index = z6x * CLUSTER_ZOOM_WIDTH + z6y;

auto& mutex = objectsMutex[z6index % objectsMutex.size()];

if (mutex.try_lock()) {

addObjectToSmallIndexUnsafe(index, oo, id);

mutex.unlock();

} else {

// add to tlsPendingSmallIndexObjects

if (tlsPendingSmallIndexObjects == nullptr) {

std::lock_guard<std::mutex> lock(objectsMutex[0]);

pendingSmallIndexObjects.push_back(std::vector<std::tuple<TileCoordinates, OutputObject, uint64_t>>());

tlsPendingSmallIndexObjects = &pendingSmallIndexObjects.back();

}

tlsPendingSmallIndexObjects->push_back(std::make_tuple(index, oo, id));

}

// Pick the z6 index

const size_t z6x = index.x / z6OffsetDivisor;

const size_t z6y = index.y / z6OffsetDivisor;

const size_t z6index = z6x * CLUSTER_ZOOM_WIDTH + z6y;

if (id == 0 || !includeID)

objects[z6index].push_back({

oo,

(Z6Offset)(index.x - (z6x * z6OffsetDivisor)),

(Z6Offset)(index.y - (z6y * z6OffsetDivisor))

});

else

objectsWithIds[z6index].push_back({

oo,

(Z6Offset)(index.x - (z6x * z6OffsetDivisor)),

(Z6Offset)(index.y - (z6y * z6OffsetDivisor)),

id

});

// Scan through all shards. Convert to base zoom, then convert to the requested zoom.

collectTilesWithObjectsAtZoomTemplate<OutputObjectXY>(indexZoom, objects.begin(), objects.size(), zooms);

collectTilesWithObjectsAtZoomTemplate<OutputObjectXYID>(indexZoom, objectsWithIds.begin(), objectsWithIds.size(), zooms);

size_t maxZoom = zooms.size() - 1;

int scale = pow(2, indexZoom - maxZoom);

TileCoordinate minx = box.min_corner().x() / scale;

TileCoordinate maxx = box.max_corner().x() / scale;

TileCoordinate miny = box.min_corner().y() / scale;

TileCoordinate maxy = box.max_corner().y() / scale;

for (int x=minx; x<=maxx; x++) {

for (int y=miny; y<=maxy; y++) {

size_t zx = x, zy = y;

for (int zoom = maxZoom; zoom >= 0; zoom--) {

zooms[zoom]->set(zx, zy);

zx /= 2;

zy /= 2;

}

}

}

for(auto const &result: boxRtree)

addCoveredTilesToOutput(indexZoom, zooms, result.first);

for(auto const &result: boxRtreeWithIds)

addCoveredTilesToOutput(indexZoom, zooms, result.first);

uint zoom,

TileCoordinates dstIndex,

std::vector<OutputObjectID>& output

) {

if (zoom < CLUSTER_ZOOM) {

collectLowZoomObjectsForTile<OutputObjectXY>(indexZoom, lowZoomObjects, zoom, dstIndex, output);

collectLowZoomObjectsForTile<OutputObjectXYID>(indexZoom, lowZoomObjectsWithIds, zoom, dstIndex, output);

return;

}

size_t iStart = 0;

size_t iEnd = objects.size();

if (zoom >= CLUSTER_ZOOM) {

TileCoordinate z6x = dstIndex.x / (1 << (zoom - CLUSTER_ZOOM));

TileCoordinate z6y = dstIndex.y / (1 << (zoom - CLUSTER_ZOOM));

if (z6x >= 64 || z6y >= 64) {

if (verbose) std::cerr << "collectObjectsForTile: invalid tile z" << zoom << "/" << dstIndex.x << "/" << dstIndex.y << std::endl;

return;

}

iStart = z6x * CLUSTER_ZOOM_WIDTH + z6y;

iEnd = iStart + 1;

}

collectObjectsForTileTemplate<OutputObjectXY>(indexZoom, objects.begin(), iStart, iEnd, zoom, dstIndex, output);

collectObjectsForTileTemplate<OutputObjectXYID>(indexZoom, objectsWithIds.begin(), iStart, iEnd, zoom, dstIndex, output);

uint zoom,

TileCoordinates dstIndex,

std::vector<OutputObjectID>& output

) {

unsigned int clampedZoom = zoom;

while (clampedZoom > indexZoom) {

clampedZoom--;

dstIndex.x /= 2;

dstIndex.y /= 2;

}

int scale = pow(2, indexZoom - clampedZoom);

TileCoordinates srcIndex1( dstIndex.x *scale , dstIndex.y *scale );

TileCoordinates srcIndex2((dstIndex.x+1)*scale-1, (dstIndex.y+1)*scale-1);

Box box = Box(geom::make<Point>(srcIndex1.x, srcIndex1.y),

geom::make<Point>(srcIndex2.x, srcIndex2.y));

for(auto const& result: boxRtree | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(box))) {

if (result.second.minZoom <= zoom)

output.push_back({result.second, 0});

}

for(auto const& result: boxRtreeWithIds | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(box))) {

if (result.second.oo.minZoom <= zoom)

output.push_back({result.second.oo, result.second.id});

}

NodeID const objectID, const TileBbox &bbox) {

switch(geomType) {

case POINT_: {

throw std::runtime_error("unexpected geomType in buildWayGeometry");

}

case LINESTRING_: {

auto const &ls = retrieveLinestring(objectID);

MultiLinestring out;

if(ls.empty())

return out;

Box extBox = bbox.getExtendBox();

const double minX = extBox.min_corner().x(), maxX = extBox.max_corner().x();

const double minY = extBox.min_corner().y(), maxY = extBox.max_corner().y();

auto pointInsideExtBox = [minX, maxX, minY, maxY](const Point& p) {

return p.x() >= minX && p.x() <= maxX && p.y() >= minY && p.y() <= maxY;

};

bool needsIntersection = !pointInsideExtBox(ls[0]);

Linestring current_ls;

geom::append(current_ls, ls[0]);

for(size_t i = 1; i < ls.size(); ++i) {

boost::geometry::model::segment<Point> segment(ls[i-1], ls[i]);

if(!geom::intersects(segment, bbox.clippingBox)) {

if(current_ls.size() > 1)

out.push_back(std::move(current_ls));

current_ls.clear();

}

geom::append(current_ls, ls[i]);

if (!needsIntersection)

needsIntersection = !pointInsideExtBox(ls[i]);

}

if(current_ls.size() > 1)

out.push_back(std::move(current_ls));

if (!needsIntersection)

return out;

MultiLinestring result;

geom::intersection(out, extBox, result);

return result;

}

case MULTILINESTRING_: {

// Look for a previously clipped version at z-1, z-2, ...

std::shared_ptr<MultiLinestring> cachedClip = multiLinestringClipCache.get(bbox.zoom, bbox.index.x, bbox.index.y, objectID);

MultiLinestring uncached;

if (cachedClip == nullptr) {

const auto& input = retrieveMultiLinestring(objectID);

boost::geometry::assign(uncached, input);

}

const auto &mls = cachedClip == nullptr ? uncached : *cachedClip;

// investigate whether filtering the constituent linestrings improves performance

MultiLinestring result;

geom::intersection(mls, bbox.getExtendBox(), result);

multiLinestringClipCache.add(bbox, objectID, result);

return result;

}

case POLYGON_: {

// Look for a previously clipped version at z-1, z-2, ...

std::shared_ptr<MultiPolygon> cachedClip = multiPolygonClipCache.get(bbox.zoom, bbox.index.x, bbox.index.y, objectID);

MultiPolygon uncached;

if (cachedClip == nullptr) {

// The cached multipolygon uses a non-standard allocator, so copy it

populateMultiPolygon(uncached, objectID);

}

const auto &input = cachedClip == nullptr ? uncached : *cachedClip;

Box box = bbox.clippingBox;

if (bbox.endZoom) {

for(auto const &p: input) {

for(auto const &inner: p.inners()) {

for(std::size_t i = 0; i < inner.size() - 1; ++i)

{

Point p1 = inner[i];

Point p2 = inner[i + 1];

if(geom::within(p1, bbox.clippingBox) != geom::within(p2, bbox.clippingBox)) {

box.min_corner() = Point(

std::min(box.min_corner().x(), std::min(p1.x(), p2.x())),

std::min(box.min_corner().y(), std::min(p1.y(), p2.y())));

box.max_corner() = Point(

std::max(box.max_corner().x(), std::max(p1.x(), p2.x())),

std::max(box.max_corner().y(), std::max(p1.y(), p2.y())));

}

}

}

for(std::size_t i = 0; i < p.outer().size() - 1; ++i) {

Point p1 = p.outer()[i];

Point p2 = p.outer()[i + 1];

if(geom::within(p1, bbox.clippingBox) != geom::within(p2, bbox.clippingBox)) {

box.min_corner() = Point(

std::min(box.min_corner().x(), std::min(p1.x(), p2.x())),

std::min(box.min_corner().y(), std::min(p1.y(), p2.y())));

box.max_corner() = Point(

std::max(box.max_corner().x(), std::max(p1.x(), p2.x())),

std::max(box.max_corner().y(), std::max(p1.y(), p2.y())));

}

}

}

Box extBox = bbox.getExtendBox();

box.min_corner() = Point(

std::max(box.min_corner().x(), extBox.min_corner().x()),

std::max(box.min_corner().y(), extBox.min_corner().y()));

box.max_corner() = Point(

std::min(box.max_corner().x(), extBox.max_corner().x()),

std::min(box.max_corner().y(), extBox.max_corner().y()));

}

MultiPolygon mp;

if (cachedClip == nullptr)

mp = std::move(uncached);

else

geom::assign(mp, input);

fast_clip(mp, box);

geom::correct(mp);

geom::validity_failure_type failure = geom::validity_failure_type::no_failure;

bool valid = geom::is_valid(mp,failure);

if (!valid) {

if (verbose) {

// Build the whole line first and emit it with a single stream write:

// tilemaker runs multi-threaded and chained operator<< calls are not

// atomic, so per-token writes interleave into unreadable output.

std::ostringstream msg;

msg << "invalid multipolygon for object " << objectID

<< " at z" << bbox.zoom << " " << bbox.index.x << "/" << bbox.index.y

<< ": " << validityFailureName(failure) << "\n";

std::cerr << msg.str();

}

if (failure==geom::failure_spikes) {

geom::remove_spikes(mp);

failure = geom::validity_failure_type::no_failure;

valid = geom::is_valid(mp,failure);

}

if (!valid && (failure==geom::failure_self_intersections || failure==geom::failure_intersecting_interiors)) {

// fast_clip can introduce self-intersections; redo the clip with the

// slower but robust Boost intersection against the original geometry.

MultiPolygon output;

if (cachedClip == nullptr) {

MultiPolygon original;

populateMultiPolygon(original, objectID);

geom::intersection(original, box, output);

} else {

geom::intersection(input, box, output);

}

geom::correct(output);

// The intersection result can itself still be invalid for very complex

// multipolygons (e.g. large reservoirs), which previously produced

// dropped or holey tiles. Repair it before returning.

repairMultiPolygon(output, objectID);

multiPolygonClipCache.add(bbox, objectID, output);

return output;

} else if (!valid) {

// occasionally also wrong_topological_dimension, disconnected_interior:

// defects geom::correct cannot mend. Repair mp in place; on failure it

// is left unchanged so behaviour never regresses.

repairMultiPolygon(mp, objectID);

}

}

multiPolygonClipCache.add(bbox, objectID, mp);

return mp;

}

default:

throw std::runtime_error("Invalid output geometry");

}

auto p = retrievePoint(objectID);

LatpLon out;

out.latp = p.y();

out.lon = p.x();

return out;

size_t points = 0, linestrings = 0, polygons = 0;

for (const auto& store : pointStores)

points += store.size();

for (const auto& store : linestringStores)

linestrings += store.size();

for (const auto& store : multilinestringStores)

linestrings += store.size();

for (const auto& store : multipolygonStores)

polygons += store.size();

std::cout << "Generated points: " << (points - 1) << ", lines: " << (linestrings - 2) << ", polygons: " << (polygons - 1) << std::endl;

const std::vector<class TileDataSource *>& sources,

std::vector<std::shared_ptr<TileCoordinatesSet>>& zooms

) {

if (zooms.size() > 15)

throw std::out_of_range("populateTilesAtZoom: expected at most z14 zooms (15), but found " + std::to_string(zooms.size()) + " vectors");

for(size_t i=0; i<sources.size(); i++) {

sources[i]->collectTilesWithObjectsAtZoom(zooms);

sources[i]->collectTilesWithLargeObjectsAtZoom(zooms);

}

const std::vector<bool>& sortOrders,

unsigned int zoom,

TileCoordinates coordinates

) {

std::vector<OutputObjectID> data;

collectObjectsForTile(zoom, coordinates, data);

collectLargeObjectsForTile(zoom, coordinates, data);

sortOutputObjectIDs(sortOrders, data);

data.erase(unique(data.begin(), data.end()), data.end());

return data;

const Linestring& geom,

const std::vector<OutputObject>& outputs,

const uint64_t id

) {

unordered_set<TileCoordinates> tileSet;

try {

insertIntermediateTiles(geom, indexZoom, tileSet);

bool polygonExists = false;

TileCoordinate minTileX = std::numeric_limits<TileCoordinate>::max(), maxTileX = 0, minTileY = std::numeric_limits<TileCoordinate>::max(), maxTileY = 0;

for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {

TileCoordinates index = *it;

minTileX = std::min(index.x, minTileX);

minTileY = std::min(index.y, minTileY);

maxTileX = std::max(index.x, maxTileX);

maxTileY = std::max(index.y, maxTileY);

for (const auto& output : outputs) {

if (output.geomType == POLYGON_) {

polygonExists = true;

continue;

}

addObjectToSmallIndex(index, output, id); // not a polygon

}

}

// for polygon, fill inner tiles

if (polygonExists) {

bool tilesetFilled = false;

uint size = (maxTileX - minTileX + 1) * (maxTileY - minTileY + 1);

for (const auto& output : outputs) {

if (output.geomType != POLYGON_) continue;

if (size>= 16) {

// Larger objects - add to rtree

Box box = Box(geom::make<Point>(minTileX, minTileY),

geom::make<Point>(maxTileX, maxTileY));

addObjectToLargeIndex(box, output, id);

} else {

// Smaller objects - add to each individual tile index

if (!tilesetFilled) { fillCoveredTiles(tileSet); tilesetFilled = true; }

for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {

TileCoordinates index = *it;

addObjectToSmallIndex(index, output, id);

}

}

}

}

} catch(std::out_of_range &err) {

cerr << "Error calculating intermediate tiles: " << err.what() << endl;

}

const MultiLinestring& geom,

const std::vector<OutputObject>& outputs,

const uint64_t id

) {

for (const auto& ls : geom) {

unordered_set<TileCoordinates> tileSet;

insertIntermediateTiles(ls, indexZoom, tileSet);

for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {

TileCoordinates index = *it;

for (const auto& output : outputs) {

addObjectToSmallIndex(index, output, id);

}

}

}

const MultiPolygon& geom,

std::vector<OutputObject>& outputs,

const uint64_t id

) {

unordered_set<TileCoordinates> tileSet;

bool singleOuter = geom.size()==1;

for (const auto& poly : geom) {

unordered_set<TileCoordinates> tileSetTmp;

insertIntermediateTiles(poly.outer(), indexZoom, tileSetTmp);

fillCoveredTiles(tileSetTmp);

if (singleOuter) {

tileSet = std::move(tileSetTmp);

} else {

tileSet.insert(tileSetTmp.begin(), tileSetTmp.end());

}

}

TileCoordinate minTileX = std::numeric_limits<TileCoordinate>::max(), maxTileX = 0, minTileY = std::numeric_limits<TileCoordinate>::max(), maxTileY = 0;

for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {

TileCoordinates index = *it;

minTileX = std::min(index.x, minTileX);

minTileY = std::min(index.y, minTileY);

maxTileX = std::max(index.x, maxTileX);

maxTileY = std::max(index.y, maxTileY);

}

const size_t tileSetSize = tileSet.size();

for (auto& output : outputs) {

if (tileSetSize >= 16) {

// Larger objects - add to rtree

// note that the bbox is currently the envelope of the entire multipolygon,

// which is suboptimal in shapes like (_) ...... (_) where the outers are significantly disjoint

Box box = Box(geom::make<Point>(minTileX, minTileY),

geom::make<Point>(maxTileX, maxTileY));

addObjectToLargeIndex(box, output, id);

} else {

// Smaller objects - add to each individual tile index

for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {

TileCoordinates index = *it;

addObjectToSmallIndex(index, output, id);

}

}

}

const auto& store = pointStore.get(this);

NodeID offset = store.second->size();

store.second->emplace_back(input);

NodeID rv = (store.first << (TILE_DATA_ID_SIZE - shardBits)) + offset;

return rv;

const auto& store = linestringStore.get(this);

linestring_t dst(src.begin(), src.end());

NodeID offset = store.second->size();

store.second->emplace_back(std::move(dst));

NodeID rv = (store.first << (TILE_DATA_ID_SIZE - shardBits)) + offset;

return rv;

const auto& store = multipolygonStore.get(this);

multi_polygon_t dst;

dst.resize(src.size());

for(std::size_t i = 0; i < src.size(); ++i) {

dst[i].outer().resize(src[i].outer().size());

boost::geometry::assign(dst[i].outer(), src[i].outer());

dst[i].inners().resize(src[i].inners().size());

for(std::size_t j = 0; j < src[i].inners().size(); ++j) {

dst[i].inners()[j].resize(src[i].inners()[j].size());

boost::geometry::assign(dst[i].inners()[j], src[i].inners()[j]);

}

}

NodeID offset = store.second->size();

store.second->emplace_back(std::move(dst));

NodeID rv = (store.first << (TILE_DATA_ID_SIZE - shardBits)) + offset;

return rv;

const auto& store = multilinestringStore.get(this);

multi_linestring_t dst;

dst.resize(src.size());

for (std::size_t i=0; i<src.size(); ++i) {

boost::geometry::assign(dst[i], src[i]);

}

NodeID offset = store.second->size();

store.second->emplace_back(std::move(dst));

NodeID rv = (store.first << (TILE_DATA_ID_SIZE - shardBits)) + offset;

return rv;

const auto &input = retrieveMultiPolygon(objectID);

dst.resize(input.size());

for(std::size_t i = 0; i < input.size(); ++i) {

dst[i].outer().resize(input[i].outer().size());

boost::geometry::assign(dst[i].outer(), input[i].outer());

dst[i].inners().resize(input[i].inners().size());

for(std::size_t j = 0; j < input[i].inners().size(); ++j) {

dst[i].inners()[j].resize(input[i].inners()[j].size());

boost::geometry::assign(dst[i].inners()[j], input[i].inners()[j]);

}

}