From 749044ce9735216ce750b271fb090f4acc39a144 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Wed, 12 Jan 2022 20:21:27 +0100 Subject: [PATCH] fix misnamed turning_circle parameter --- pfaedle.cfg | 8 +- src/pfaedle/config/MotConfigReader.cpp | 4 +- src/shapevl/Collector.cpp | 67 +++++-------- src/shapevl/Collector.h | 12 +-- src/util/geo/Geo.h | 132 ++++++++++++++++++++----- 5 files changed, 142 insertions(+), 81 deletions(-) diff --git a/pfaedle.cfg b/pfaedle.cfg index 443935f..711bc6d 100644 --- a/pfaedle.cfg +++ b/pfaedle.cfg @@ -5,7 +5,7 @@ [tram, bus, coach, subway, rail, gondola, funicular, ferry] routing_transition_penalty_fac: 0.0083 -routing_station_move_penalty_fac: 0.00087 +routing_station_move_penalty_fac: 0.002 # Regular expressions and station comparision is # always case insensitive! @@ -717,8 +717,8 @@ osm_filter_station: highway=bus_stop amenity=bus_station -osm_filter_turning_cycle: - highway=turning_cycle +osm_filter_turning_circle: + highway=turning_circle highway=turning_loop junction=roundabout highway=mini_roundabout @@ -791,7 +791,7 @@ routing_line_station_from_unmatched_time_penalty: 1.1 # routing_no_lines_penalty_fac: 1 # If the station name does not match, add this penalty -routing_station_unmatched_penalty: 0.1 +routing_station_unmatched_penalty: 0.2 # Punishment (in seconds) to add to the distance # function if a vehicle performans a full turn diff --git a/src/pfaedle/config/MotConfigReader.cpp b/src/pfaedle/config/MotConfigReader.cpp index 0db2cf2..e9f3905 100644 --- a/src/pfaedle/config/MotConfigReader.cpp +++ b/src/pfaedle/config/MotConfigReader.cpp @@ -152,9 +152,9 @@ void MotConfigReader::parse(const std::vector& paths, } } - if (p.hasKey(secStr, "osm_filter_turning_cycle")) { + if (p.hasKey(secStr, "osm_filter_turning_circle")) { for (const auto& kvs : - p.getStrArr(sec.first, "osm_filter_turning_cycle", ' ')) { + p.getStrArr(sec.first, "osm_filter_turning_circle", ' ')) { auto fRule = getFRule(kvs); cfg.osmBuildOpts.turnCycleFilter[fRule.kv.first].insert( osm::AttrFlagPair(fRule.kv.second, getFlags(fRule.flags))); diff --git a/src/shapevl/Collector.cpp b/src/shapevl/Collector.cpp index 6cda27b..def02f5 100644 --- a/src/shapevl/Collector.cpp +++ b/src/shapevl/Collector.cpp @@ -62,10 +62,10 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT, double unmatchedSegmentsLength; // total _an. length of unmatched segments std::vector oldDists; - LINE oldL = getWebMercLine(oldS, &oldDists); + LINE oldL = getLine(oldS, &oldDists); std::vector newDists; - LINE newL = getWebMercLine(newS, &newDists); + LINE newL = getLine(newS, &newDists); std::vector> newLenDists; std::vector> oldLenDists; @@ -88,13 +88,8 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT, newLCut.insert(newLCut.end(), newL.begin(), newL.end()); } - // determine the scale factor between the distance in projected - // coordinates and the real-world distance in meters - auto avgY = - (oldSegs.front().front().getY() + oldSegs.back().back().getY()) / 2; - double fac = cos(2 * atan(exp(avgY / 6378137.0)) - 1.5707965); - - double SEGL = 10; + // convert (roughly) to degrees + double SEGL = 15.0 / util::geo::M_PER_DEG; auto old = _dCache.find(oldLCut); if (old != _dCache.end()) { @@ -102,11 +97,11 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT, if (match != old->second.end()) { fd = match->second; } else { - fd = util::geo::accFrechetDistC(oldLCut, newLCut, SEGL / fac) * fac; + fd = util::geo::accFrechetDistCHav(oldLCut, newLCut, SEGL); _dCache[oldLCut][newLCut] = fd; } } else { - fd = util::geo::accFrechetDistC(oldLCut, newLCut, SEGL / fac) * fac; + fd = util::geo::accFrechetDistCHav(oldLCut, newLCut, SEGL); _dCache[oldLCut][newLCut] = fd; } @@ -115,7 +110,7 @@ double Collector::add(const Trip* oldT, const Shape* oldS, const Trip* newT, unmatchedSegmentsLength = dA.second; double totL = 0; - for (auto l : oldSegs) totL += util::geo::len(l) * fac; + for (auto l : oldSegs) totL += util::geo::latLngLen(l); // filter out shapes with a length of under 5 meters - they are most likely // artifacts @@ -171,8 +166,8 @@ std::vector Collector::segmentize( size_t i = 0; for (auto st : t->getStopTimes()) { cuts.push_back(std::pair( - util::geo::latLngToWebMerc(st.getStop()->getLat(), - st.getStop()->getLng()), + {st.getStop()->getLng(), + st.getStop()->getLat()}, st.getShapeDistanceTravelled())); i++; } @@ -204,7 +199,7 @@ std::vector Collector::segmentize( t->getStopTimes()[i - 1].getStop()->getLng(), t->getStopTimes()[i].getStop()->getLat(), t->getStopTimes()[i].getStop()->getLng()); - double len = util::geo::webMercLen(curL); + double len = util::geo::latLngLen(curL); lenDist.push_back({dist, len}); ret.push_back(curL); @@ -215,12 +210,11 @@ std::vector Collector::segmentize( } // _____________________________________________________________________________ -LINE Collector::getWebMercLine(const Shape* s, std::vector* dists) { +LINE Collector::getLine(const Shape* s, std::vector* dists) { LINE ret; for (size_t i = 0; i < s->getPoints().size(); i++) { - ret.push_back(util::geo::latLngToWebMerc(s->getPoints()[i].lat, - s->getPoints()[i].lng)); + ret.push_back({s->getPoints()[i].lng, s->getPoints()[i].lat}); (*dists).push_back(s->getPoints()[i].travelDist); } @@ -233,19 +227,6 @@ const std::set& Collector::getResults() const { return _results; } // _____________________________________________________________________________ double Collector::getAvgDist() const { return _fdSum / _results.size(); } -// _____________________________________________________________________________ -std::vector Collector::getBins(double mind, double maxd, size_t steps) { - double bin = (maxd - mind) / steps; - double curE = mind + bin; - - std::vector ret; - while (curE <= maxd) { - ret.push_back(curE); - curE += bin; - } - return ret; -} - // _____________________________________________________________________________ void Collector::printCsv(std::ostream* os, const std::set& result) const { @@ -404,32 +385,28 @@ std::pair Collector::getDa(const std::vector& a, assert(a.size() == b.size()); std::pair ret{0, 0}; - // euclidean distance on web mercator is in meters on equator, - // and proportional to cos(lat) in both y directions - - double fac = webMercDistFactor(a.front().front()); - - double SEGL = 10; + // convert (roughly) to degrees + double SEGL = 15.0 / util::geo::M_PER_DEG; for (size_t i = 0; i < a.size(); i++) { - double fd = 0; + double fdMeter = 0; auto old = _dACache.find(a[i]); if (old != _dACache.end()) { auto match = old->second.find(b[i]); if (match != old->second.end()) { - fd = match->second; + fdMeter = match->second; } else { - fd = util::geo::frechetDist(a[i], b[i], SEGL / fac) * fac; - _dACache[a[i]][b[i]] = fd; + fdMeter = util::geo::frechetDistHav(a[i], b[i], SEGL); + _dACache[a[i]][b[i]] = fdMeter; } } else { - fd = util::geo::frechetDist(a[i], b[i], SEGL / fac) * fac; - _dACache[a[i]][b[i]] = fd; + fdMeter = util::geo::frechetDistHav(a[i], b[i], SEGL); + _dACache[a[i]][b[i]] = fdMeter; } - if (fd >= 50) { + if (fdMeter >= 50) { ret.first++; - ret.second += util::geo::webMercLen(a[i]) * 100; + ret.second += util::geo::latLngLen(a[i]); } } diff --git a/src/shapevl/Collector.h b/src/shapevl/Collector.h index bf5465a..18474bf 100644 --- a/src/shapevl/Collector.h +++ b/src/shapevl/Collector.h @@ -31,8 +31,10 @@ struct lineCmp { } for (size_t i = 0; i < a.size(); i++) { - if (util::geo::dist(a[i], b[i]) > .1) { - return (a[i].getX() < b[i].getX()) || (a[i].getX() == b[i].getX() && a[i].getY() < b[i].getY());; + if (util::geo::dist(a[i], b[i]) > .000001) { + return (a[i].getX() < b[i].getX()) || + (a[i].getX() == b[i].getX() && a[i].getY() < b[i].getY()); + ; } } @@ -83,7 +85,7 @@ class Collector { // Return the averaged average frechet distance double getAvgDist() const; - static LINE getWebMercLine(const Shape* s, std::vector* dists); + static LINE getLine(const Shape* s, std::vector* dists); double getAcc() const; @@ -112,13 +114,11 @@ class Collector { std::ostream* _reportOut; std::pair getDa(const std::vector& a, - const std::vector& b); + const std::vector& b); static std::vector segmentize( const Trip* t, const LINE& shape, const std::vector& dists, std::vector>& lenDist); - - static std::vector getBins(double mind, double maxd, size_t steps); }; } // namespace eval diff --git a/src/util/geo/Geo.h b/src/util/geo/Geo.h index c11557c..20c5b36 100644 --- a/src/util/geo/Geo.h +++ b/src/util/geo/Geo.h @@ -56,7 +56,6 @@ const static double AVERAGING_STEP = 20; const static double M_PER_DEG = 111319.4; - // _____________________________________________________________________________ template inline Box pad(const Box& box, double padding) { @@ -1778,6 +1777,29 @@ inline RotatedBox getOrientedEnvelopeAvg(MultiLine ml) { return rbox; } +// _____________________________________________________________________________ +template +inline double haversine(T lat1, T lon1, T lat2, T lon2) { + lat1 *= RAD; + lat2 *= RAD; + + const double dLat = lat2 - lat1; + const double dLon = (lon2 - lon1) * RAD; + + const double sDLat = sin(dLat / 2); + const double sDLon = sin(dLon / 2); + + const double a = (sDLat * sDLat) + (sDLon * sDLon) * cos(lat1) * cos(lat2); + return 6378137.0 * 2.0 * asin(sqrt(a)); +} + +// _____________________________________________________________________________ +template +inline double haversine(const Point& a, const Point& b) { + return haversine(a.getY(), a.getX(), b.getY(), b.getX()); +} + + // _____________________________________________________________________________ template inline Line densify(const Line& l, double d) { @@ -1878,6 +1900,81 @@ inline double accFrechetDistC(const Line& a, const Line& b, double d) { return ca[p.size() * q.size() - 1]; } +// _____________________________________________________________________________ +template +inline double frechetDistCHav(size_t i, size_t j, const Line& p, + const Line& q, std::vector& ca) { + // based on Eiter / Mannila + // http://www.kr.tuwien.ac.at/staff/eiter/et-archive/cdtr9464.pdf + + if (ca[i * q.size() + j] > -1) + return ca[i * q.size() + j]; + else if (i == 0 && j == 0) + ca[i * q.size() + j] = haversine(p[0], q[0]); + else if (i > 0 && j == 0) + ca[i * q.size() + j] = + std::max(frechetDistCHav(i - 1, 0, p, q, ca), haversine(p[i], q[0])); + else if (i == 0 && j > 0) + ca[i * q.size() + j] = + std::max(frechetDistCHav(0, j - 1, p, q, ca), haversine(p[0], q[j])); + else if (i > 0 && j > 0) + ca[i * q.size() + j] = + std::max(std::min(std::min(frechetDistCHav(i - 1, j, p, q, ca), + frechetDistCHav(i - 1, j - 1, p, q, ca)), + frechetDistCHav(i, j - 1, p, q, ca)), + haversine(p[i], q[j])); + else + ca[i * q.size() + j] = std::numeric_limits::infinity(); + + return ca[i * q.size() + j]; +} + +// _____________________________________________________________________________ +template +inline double frechetDistHav(const Line& a, const Line& b, double d) { + // based on Eiter / Mannila + // http://www.kr.tuwien.ac.at/staff/eiter/et-archive/cdtr9464.pdf + + auto p = densify(a, d); + auto q = densify(b, d); + + std::vector ca(p.size() * q.size(), -1.0); + double fd = frechetDistCHav(p.size() - 1, q.size() - 1, p, q, ca); + + return fd; +} + +// _____________________________________________________________________________ +template +inline double accFrechetDistCHav(const Line& a, const Line& b, double d) { + auto p = densify(a, d); + auto q = densify(b, d); + + assert(p.size()); + assert(q.size()); + + std::vector ca(p.size() * q.size(), 0); + + for (size_t i = 0; i < p.size(); i++) + ca[i * q.size() + 0] = std::numeric_limits::infinity(); + for (size_t j = 0; j < q.size(); j++) + ca[j] = std::numeric_limits::infinity(); + ca[0] = 0; + + for (size_t i = 1; i < p.size(); i++) { + for (size_t j = 1; j < q.size(); j++) { + float d = + util::geo::haversine(p[i], q[j]) * util::geo::dist(p[i], p[i - 1]); + ca[i * q.size() + j] = + d + std::min(ca[(i - 1) * q.size() + j], + std::min(ca[i * q.size() + (j - 1)], + ca[(i - 1) * q.size() + (j - 1)])); + } + } + + return ca[p.size() * q.size() - 1]; +} + // _____________________________________________________________________________ template inline Point latLngToWebMerc(double lat, double lng) { @@ -1890,7 +1987,7 @@ inline Point latLngToWebMerc(double lat, double lng) { // _____________________________________________________________________________ template -//TODO: rename to lngLat +// TODO: rename to lngLat inline Point latLngToWebMerc(Point lngLat) { return latLngToWebMerc(lngLat.getY(), lngLat.getX()); } @@ -1904,28 +2001,6 @@ inline Point webMercToLatLng(double x, double y) { return Point(lon, lat); } -// _____________________________________________________________________________ -template -inline double haversine(T lat1, T lon1, T lat2, T lon2) { - lat1 *= RAD; - lat2 *= RAD; - - const double dLat = lat2 - lat1; - const double dLon = (lon2 - lon1) * RAD; - - const double sDLat = sin(dLat / 2); - const double sDLon = sin(dLon / 2); - - const double a = (sDLat * sDLat) + (sDLon * sDLon) * cos(lat1) * cos(lat2); - return 6378137.0 * 2.0 * asin(sqrt(a)); -} - -// _____________________________________________________________________________ -template -inline double haversine(const Point& a, const Point& b) { - return haversine(a.getY(), a.getX(), b.getY(), b.getX()); -} - // _____________________________________________________________________________ template inline double webMercMeterDist(const Point& a, const Point& b) { @@ -1959,6 +2034,15 @@ inline double webMercLen(const Line& g) { return ret; } +// _____________________________________________________________________________ +template +inline double latLngLen(const Line& g) { + double ret = 0; + for (size_t i = 1; i < g.size(); i++) ret += haversine(g[i - 1], g[i]); + return ret; +} + + // _____________________________________________________________________________ template inline double webMercDistFactor(const G& a) {