Skip to content

Commit

Permalink
Merge branch 'grafi-tt-fill-polygons-by-tiles'
Browse files Browse the repository at this point in the history
  • Loading branch information
systemed committed Nov 3, 2016
2 parents 413a8b5 + a2d9870 commit a62c455
Show file tree
Hide file tree
Showing 3 changed files with 139 additions and 129 deletions.
76 changes: 64 additions & 12 deletions src/coordinates.cpp
Expand Up @@ -18,9 +18,12 @@ double lat2latp(double lat) { return rad2deg(log(tan(deg2rad(lat+90.0)/2.0))); }
double latp2lat(double latp) { return rad2deg(atan(exp(deg2rad(latp)))*2.0)-90.0; }

// Tile conversions
uint lon2tilex(double lon, uint z) { return scalbn((lon+180.0) * (1/360.0), (int)z); }
uint latp2tiley(double latp, uint z) { return scalbn((180.0-latp) * (1/360.0), (int)z); }
uint lat2tiley(double lat, uint z) { return latp2tiley(lat2latp(lat), z); }
double lon2tilexf(double lon, uint z) { return scalbn((lon+180.0) * (1/360.0), (int)z); }
double latp2tileyf(double latp, uint z) { return scalbn((180.0-latp) * (1/360.0), (int)z); }
double lat2tileyf(double lat, uint z) { return latp2tileyf(lat2latp(lat), z); }
uint lon2tilex(double lon, uint z) { return lon2tilexf(lon, z); }
uint latp2tiley(double latp, uint z) { return latp2tileyf(latp, z); }
uint lat2tiley(double lat, uint z) { return lat2tileyf(lat, z); }
double tilex2lon(uint x, uint z) { return scalbn(x, -(int)z) * 360.0 - 180.0; }
double tiley2latp(uint y, uint z) { return 180.0 - scalbn(y, -(int)z) * 360.0; }
double tiley2lat(uint y, uint z) { return latp2lat(tiley2latp(y, z)); }
Expand All @@ -44,15 +47,64 @@ double meter2degp(double meter, double latp) {
return rad2deg((1/RadiusMeter) * (meter / cos(deg2rad(latp2lat(latp)))));
}

// Add intermediate points so we don't skip tiles on long segments
void insertIntermediateTiles(unordered_set <uint32_t> *tlPtr, int numPoints, LatpLon startLL, LatpLon endLL, uint baseZoom) {
numPoints *= 3; // perhaps overkill, but why not
int32_t dLon = endLL.lon -startLL.lon ;
int32_t dLatp = endLL.latp-startLL.latp;
for (int i=1; i<numPoints-1; i++) {
LatpLon ll = LatpLon { startLL.latp + dLatp/numPoints*i,
startLL.lon + dLon /numPoints*i };
tlPtr->insert(latpLon2index(ll,baseZoom));
template<typename T>
void insertIntermediateTiles(const T &points, uint baseZoom, unordered_set<uint32_t> &tileSet) {
Point p2(0, 0);
for (auto it = points.begin(); it != points.end(); ++it) {
Point p1 = p2;
p2 = *it;

double tileXf2 = lon2tilexf(p2.x(), baseZoom), tileYf2 = latp2tileyf(p2.y(), baseZoom);
uint tileX2 = static_cast<uint>(tileXf2), tileY2 = static_cast<uint>(tileYf2);

// insert vertex
tileSet.insert((tileX2 << 16) + tileY2);
// p1 is not available at the first iteration
if (it == points.begin()) continue;

double tileXf1 = lon2tilexf(p1.x(), baseZoom), tileYf1 = latp2tileyf(p1.y(), baseZoom);
uint tileX1 = static_cast<uint>(tileXf1), tileY1 = static_cast<uint>(tileYf1);
double dx = tileXf2 - tileXf1, dy = tileYf2 - tileYf1;

// insert all X border
if (tileX1 != tileX2) {
double slope = dy / dx;
uint tileXmin = min(tileX1, tileX2);
uint tileXmax = max(tileX1, tileX2);
for (uint tileXcur = tileXmin+1; tileXcur <= tileXmax; tileXcur++) {
uint tileYcur = static_cast<uint>(tileYf1 + (static_cast<double>(tileXcur) - tileXf1) * slope);
tileSet.insert((tileXcur << 16) + tileYcur);
}
}

// insert all Y border
if (tileY1 != tileY2) {
double slope = dx / dy;
uint tileYmin = min(tileY1, tileY2);
uint tileYmax = max(tileY1, tileY2);
for (uint tileYcur = tileYmin+1; tileYcur <= tileYmax; tileYcur++) {
uint tileXcur = static_cast<uint>(tileXf1 + (static_cast<double>(tileYcur) - tileYf1) * slope);
tileSet.insert((tileXcur << 16) + tileYcur);
}
}
}
}

// the range between smallest y and largest y is filled, for each x
void fillCoveredTiles(unordered_set<uint32_t> &tileSet) {
vector<uint32_t> tileList(tileSet.begin(), tileSet.end());
sort(tileList.begin(), tileList.end());

uint prevX = 0, prevY = static_cast<uint>(-2);
for (uint32_t index: tileList) {
uint tileX = index >> 16, tileY = index & 65535;
if (tileX == prevX) {
// this loop has no effect at the first iteration
for (uint fillY = prevY+1; fillY < tileY; fillY++) {
tileSet.insert((tileX << 16) + fillY);
}
}
prevX = tileX, prevY = tileY;
}
}

Expand Down
6 changes: 3 additions & 3 deletions src/osm_store.cpp
Expand Up @@ -180,18 +180,18 @@ struct OSMStore {
if (wayList.outerBegin != wayList.outerEnd) {
// main outer way and inners
Polygon poly;
fillPoints(poly.outer(), ways.at(*wayList.outerBegin++));
if (ways.count(*wayList.outerBegin)==1) { fillPoints(poly.outer(), ways.at(*wayList.outerBegin++)); }
for (auto it = wayList.innerBegin; it != wayList.innerEnd; ++it) {
Ring inner;
fillPoints(inner, ways.at(*it));
if (ways.count(*it)==1) { fillPoints(inner, ways.at(*it)); }
poly.inners().emplace_back(move(inner));
}
mp.emplace_back(move(poly));

// additional outer ways - we don't match them up with inners, that shit is insane
for (auto it = wayList.outerBegin; it != wayList.outerEnd; ++it) {
Polygon outerPoly;
fillPoints(outerPoly.outer(), ways.at(*it));
if (ways.count(*it)==1) { fillPoints(outerPoly.outer(), ways.at(*it)); }
mp.emplace_back(move(outerPoly));
}

Expand Down
186 changes: 72 additions & 114 deletions src/tilemaker.cpp
Expand Up @@ -117,8 +117,6 @@ int main(int argc, char* argv[]) {
map<uint, string> cachedGeometryNames; // | optional names for each one

map< uint, vector<OutputObject> > tileIndex; // objects to be output
map< WayID, vector<OutputObject> > relationOutputObjects; // outputObjects for multipolygons (saved for processing later as ways)
map< WayID, vector<WayID> > wayRelations; // for each way, which relations it's in (therefore we need to keep them)

// ---- Read command-line options

Expand Down Expand Up @@ -362,8 +360,7 @@ int main(int argc, char* argv[]) {
// ---- Read PBF
// note that the order of reading and processing is:
// 1) output nodes -> (remember current position for rewinding to ways) (skip ways) -> (just remember all ways in any relation),
// 2) (for the remembered ways, construct nodeId lists) -> output relations, though the actual output task is delayed until each way's processing
// 3) output ways, with every relation which contains the way
// 2) output ways, and also construct nodeId list for each way in relation -> output relations

cout << "Reading " << inputFile << endl;

Expand All @@ -380,7 +377,6 @@ int main(int argc, char* argv[]) {
uint i,j,k,ct=0;
int64_t nodeId;
bool checkedRelations = false;
bool processedRelations = false;
int wayPosition = -1;
unordered_set<WayID> waysInRelation;

Expand All @@ -390,10 +386,6 @@ int main(int argc, char* argv[]) {
if (infile.eof()) {
if (!checkedRelations) {
checkedRelations = true;
} else if (!processedRelations) {
processedRelations = true;
// NodeId lists for ways were constructed to process relations. Then reset it, because relations processing have ended.
ways.clear();
} else {
break;
}
Expand Down Expand Up @@ -488,35 +480,72 @@ int main(int argc, char* argv[]) {
break;
}

// ---- For the remembered ways, construct nodeId lists
// ---- Read ways

if (!processedRelations && pg.ways_size() > 0) {
if (pg.ways_size() > 0) {
for (j=0; j<pg.ways_size(); j++) {
pbfWay = pg.ways(j);
WayID wayId = pbfWay.id();
if (waysInRelation.count(wayId) > 0) {
// Assemble nodelist
nodeId = 0;
NodeVec nodeVec;
for (k = 0; k < pbfWay.refs_size(); k++) {
nodeId += pbfWay.refs(k);
nodeVec.push_back(static_cast<NodeID>(nodeId));
}
WayID wayId = static_cast<WayID>(pbfWay.id());

// Assemble nodelist
nodeId = 0;
NodeVec nodeVec;
for (k=0; k<pbfWay.refs_size(); k++) {
nodeId += pbfWay.refs(k);
nodeVec.push_back(static_cast<NodeID>(nodeId));
}

osmObject.setWay(&pbfWay, &nodeVec);
// Call Lua to find what layers and tags we want
try { luabind::call_function<int>(luaState, "way_function", &osmObject);
} catch (const luabind::error &er) {
cerr << er.what() << endl << "-- " << lua_tostring(er.state(), -1) << endl;
return -1;
}

if (!osmObject.empty() || waysInRelation.count(wayId)) {
// Store the way's nodes in the global way store
ways.insert_back(wayId, nodeVec);
}

if (!osmObject.empty()) {
// create a list of tiles this way passes through (tileSet)
unordered_set<uint32_t> tileSet;
insertIntermediateTiles(osmStore.nodeListLinestring(nodeVec), baseZoom, tileSet);

// then, for each tile, store the OutputObject for each layer
bool polygonExists = false;
for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
uint32_t index = *it;
for (auto jt = osmObject.outputs.begin(); jt != osmObject.outputs.end(); ++jt) {
if (jt->geomType == POLYGON) {
polygonExists = true;
continue;
}
tileIndex[index].push_back(*jt);
}
}

// for polygon, fill inner tiles
if (polygonExists) {
fillCoveredTiles(tileSet);
for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
uint32_t index = *it;
for (auto jt = osmObject.outputs.begin(); jt != osmObject.outputs.end(); ++jt) {
if (jt->geomType != POLYGON) continue;
tileIndex[index].push_back(*jt);
}
}
}
}
}
continue;
}

if (!processedRelations && !waysInRelation.empty()) {
// forget those ways here to save memory
waysInRelation.clear();
}

// ---- Read relations
// (just multipolygons for now; we should do routes in time)

if (!processedRelations && pg.relations_size() > 0) {
if (pg.relations_size() > 0) {
int typeKey = osmObject.findStringPosition("type");
int mpKey = osmObject.findStringPosition("multipolygon");
int innerKey= osmObject.findStringPosition("inner");
Expand Down Expand Up @@ -551,103 +580,32 @@ int main(int argc, char* argv[]) {
WayID relID = osmObject.osmID;
// Store the relation members in the global relation store
relations.insert_front(relID, outerWayVec, innerWayVec);
// Store this relation in the way->relations map to oblige each way in the relation
// to output it, even if the way is not rendered in its own right.
for (auto it = outerWayVec.cbegin(); it != outerWayVec.cend(); ++it) {
wayRelations[*it].push_back(relID);
}
for (auto it = innerWayVec.cbegin(); it != innerWayVec.cend(); ++it) {
wayRelations[*it].push_back(relID);
}
// Keep output objects
for (auto jt = osmObject.outputs.begin(); jt != osmObject.outputs.end(); ++jt) {
relationOutputObjects[relID].push_back(*jt);
}
}
}
}
continue;
}

if (!processedRelations) {
// Nothing to do
break;
}

// ---- Read ways

if (pg.ways_size() > 0) {
for (j=0; j<pg.ways_size(); j++) {
pbfWay = pg.ways(j);
WayID wayId = static_cast<WayID>(pbfWay.id());

// Assemble nodelist
nodeId = 0;
NodeVec nodeVec;
for (k=0; k<pbfWay.refs_size(); k++) {
nodeId += pbfWay.refs(k);
nodeVec.push_back(static_cast<NodeID>(nodeId));
}

osmObject.setWay(&pbfWay, &nodeVec);
// Call Lua to find what layers and tags we want
try { luabind::call_function<int>(luaState, "way_function", &osmObject);
} catch (const luabind::error &er) {
cerr << er.what() << endl << "-- " << lua_tostring(er.state(), -1) << endl;
return -1;
}

bool inRelation = wayRelations.count(pbfWay.id()) > 0;
if (!osmObject.empty() || inRelation) {
// Store the way's nodes in the global way store
ways.insert_back(wayId, nodeVec);

// create a list of tiles this way passes through (tilelist)
unordered_set <uint32_t> tilelist;
uint lastX, lastY;
for (k=0; k<pbfWay.refs_size(); k++) {
uint tileX = lon2tilex(nodes.at(nodeVec[k]).lon / 10000000.0, baseZoom);
uint tileY = latp2tiley(nodes.at(nodeVec[k]).latp / 10000000.0, baseZoom);
if (k>0) {
// Check we're not skipping any tiles, and insert intermediate nodes if so
// (we should have a simple fill algorithm for polygons, too)
int dx = abs((int)tileX-(int)lastX);
int dy = abs((int)tileY-(int)lastY);
if (dx>1 || dy>1 || (dx==1 && dy==1)) {
insertIntermediateTiles(&tilelist, max(dx,dy), nodes.at(nodeVec[k-1]), nodes.at(nodeVec[k]), baseZoom);
// for each tile the relation may cover, put the output objects.
unordered_set<uint32_t> tileSet;
MultiPolygon mp = osmStore.wayListMultiPolygon(outerWayVec, innerWayVec);
if (mp.size() == 1) {
insertIntermediateTiles(mp[0].outer(), baseZoom, tileSet);
fillCoveredTiles(tileSet);
} else {
for (Polygon poly: mp) {
unordered_set<uint32_t> tileSetTmp;
insertIntermediateTiles(poly.outer(), baseZoom, tileSetTmp);
fillCoveredTiles(tileSetTmp);
tileSet.insert(tileSetTmp.begin(), tileSetTmp.end());
}
}
uint32_t index = tileX * 65536 + tileY;
tilelist.insert( index );
lastX = tileX;
lastY = tileY;
}

// then, for each tile, store the OutputObject for each layer
for (auto it = tilelist.begin(); it != tilelist.end(); ++it) {
uint32_t index = *it;
for (auto jt = osmObject.outputs.begin(); jt != osmObject.outputs.end(); ++jt) {
tileIndex[index].push_back(*jt);
}
}

// if it's in any relations to be output, do the same for each relation
if (inRelation) {
for (auto wt = wayRelations[wayId].begin(); wt != wayRelations[wayId].end(); ++wt) {
WayID relID = *wt;
// relID is now the relation ID
for (auto it = tilelist.begin(); it != tilelist.end(); ++it) {
// index is now the tile index number
uint32_t index = *it;
// add all the OutputObjects for this relation into this tile
for (auto jt = relationOutputObjects[relID].begin(); jt != relationOutputObjects[relID].end(); ++jt) {
tileIndex[index].push_back(*jt);
}
for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
uint32_t index = *it;
for (auto jt = osmObject.outputs.begin(); jt != osmObject.outputs.end(); ++jt) {
tileIndex[index].push_back(*jt);
}
}
}
}
}
continue;
}

// Everything should be ended
Expand Down

0 comments on commit a62c455

Please sign in to comment.