From b971d52669974ad5982040942924309f5cf98877 Mon Sep 17 00:00:00 2001 From: Vladimir Byko-Ianko Date: Tue, 14 Jun 2016 09:07:22 +0300 Subject: [PATCH] Review fixes. --- .../feature_segments_checker.cpp | 83 ++++++++++++------- 1 file changed, 52 insertions(+), 31 deletions(-) diff --git a/generator/feature_segments_checker/feature_segments_checker.cpp b/generator/feature_segments_checker/feature_segments_checker.cpp index 4632446a72..c5d299dcf7 100644 --- a/generator/feature_segments_checker/feature_segments_checker.cpp +++ b/generator/feature_segments_checker/feature_segments_checker.cpp @@ -15,6 +15,7 @@ #include "platform/platform.hpp" #include "base/logging.hpp" +#include "base/math.hpp" #include "std/cmath.hpp" #include "std/fstream.hpp" @@ -31,6 +32,9 @@ DEFINE_string(mwm_file_path, "", "Path to an mwm file."); namespace { +using TAltitude = generator::SrtmTile::THeight; +using TAltitudeVec = vector; + routing::BicycleModel const & GetBicycleModel() { static routing::BicycleModel const instance; @@ -76,10 +80,32 @@ void WriteCSV(Cont const & cont, string const & fileName) fout << a.first << "," << a.second << endl; } -/// \returns y = k * x + b. It's an expected altitude in meters. -double ExpectedPointAltitude(double k, double b, double distFromStartMeters) +/// \returns y = k * x + b. It's the expected altitude in meters. +double GetY(double k, double b, double x) { return k * x + b; } + +/// \brief Calculates factors |k| and |b| of a line using linear least squares method. +/// \returns false in case of error (e.g. if the line is parallel to the vertical axis) +/// and true otherwise. +bool LinearLeastSquaresFactors(vector const & xs, vector const & ys, double & k, + double & b) { - return k * distFromStartMeters + b; + double constexpr kEpsilon = 1e-6; + size_t const n = xs.size(); + double mx = 0, my = 0, mxy = 0, mx2 = 0; + for (int i = 0; i < n; ++i) + { + mx += xs[i] / n; + my += ys[i] / n; + mxy += xs[i] * ys[i] / n; + mx2 += xs[i] * xs[i] / n; + } + + if (my::AlmostEqualAbs(mx * mx, mx2, kEpsilon)) + return false; + + k = (my * mx - mxy) / (mx * mx - mx2); + b = my - k * mx; + return true; } class Processor @@ -89,7 +115,7 @@ public: set m_uniqueRoadPoints; /// Key is an altitude difference for a feature in meters. If a feature goes up the key is greater then 0. /// Value is a number of features. - map m_altitudeDiffs; + map m_altitudeDiffs; /// Key is a length of a feature in meters. Value is a number of features. map m_featureLength; /// Key is a length of a feature segment in meters. Value is a segment counter. @@ -109,11 +135,11 @@ public: /// Key is number of meters. It shows altitude deviation of intermediate feature points /// from linear model. /// Value is a number of features. - map m_diffFromLinear; + map m_diffFromLinear; /// Key is number of meters. It shows altitude deviation of intermediate feature points /// from line calculated base on least squares method for all feature points. /// Value is a number of features. - map m_leastSquaresDiff; + map m_leastSquaresDiff; /// Number of features for GetBicycleModel().IsRoad(feature) == true. uint32_t m_roadCount; /// Number of features for empty features with GetBicycleModel().IsRoad(feature). @@ -122,8 +148,8 @@ public: uint32_t m_roadPointCount; /// Number of features for GetBicycleModel().IsRoad(feature) != true. uint32_t m_notRoadCount; - generator::SrtmTile::THeight m_minAltitude = generator::SrtmTile::kInvalidHeight; - generator::SrtmTile::THeight m_maxAltitude = generator::SrtmTile::kInvalidHeight; + TAltitude m_minAltitude = generator::SrtmTile::kInvalidHeight; + TAltitude m_maxAltitude = generator::SrtmTile::kInvalidHeight; Processor(generator::SrtmTileManager & manager) : m_srtmManager(manager), m_roadCount(0), m_emptyRoadCount(0), m_roadPointCount(0), m_notRoadCount(0) @@ -155,20 +181,22 @@ public: m_uniqueRoadPoints.insert(RoughPoint(f.GetPoint(i))); // Preparing feature altitude and length. - vector pointAltitudes(numPoints); + TAltitudeVec pointAltitudes(numPoints); vector pointDists(numPoints); double distFromStartMeters = 0; for (uint32_t i = 0; i < numPoints; ++i) { // Feature segment altitude. - pointAltitudes[i] = m_srtmManager.GetHeight(MercatorBounds::ToLatLon(f.GetPoint(i))); + TAltitude altitude = m_srtmManager.GetHeight(MercatorBounds::ToLatLon(f.GetPoint(i))); + pointAltitudes[i] = altitude == generator::SrtmTile::kInvalidHeight ? 0 : altitude; if (i == 0) { pointDists[i] = 0; continue; } // Feature segment length. - double const segmentLengthMeters = MercatorBounds::DistanceOnEarth(f.GetPoint(i - 1), f.GetPoint(i)); + double const segmentLengthMeters = + MercatorBounds::DistanceOnEarth(f.GetPoint(i - 1), f.GetPoint(i)); distFromStartMeters += segmentLengthMeters; pointDists[i] = distFromStartMeters; } @@ -196,8 +224,8 @@ public: m_featureLength[static_cast(realFeatureLengthMeters)]++; // Feature altitude difference. - generator::SrtmTile::THeight const startAltitude = pointAltitudes[0]; - generator::SrtmTile::THeight const endAltitude = pointAltitudes[numPoints - 1]; + TAltitude const startAltitude = pointAltitudes[0]; + TAltitude const endAltitude = pointAltitudes[numPoints - 1]; int16_t const altitudeDiff = endAltitude - startAltitude; m_altitudeDiffs[altitudeDiff]++; @@ -240,40 +268,32 @@ public: return; double const k = (endAltitude - startAltitude) / realFeatureLengthMeters; - for (uint32_t i = 1; i + 1 < numPoints; ++i) + for (TAltitude i = 1; i + 1 < numPoints; ++i) { int32_t const deviation = - static_cast(ExpectedPointAltitude(k, startAltitude, pointDists[i])) - pointAltitudes[i]; + static_cast(GetY(k, startAltitude, pointDists[i])) - pointAltitudes[i]; m_diffFromLinear[deviation]++; } // Linear least squares for feature points. { - long double const distSum = accumulate(pointDists.begin(), pointDists.end(), 0.0); - uint64_t const altSum = accumulate(pointAltitudes.begin(), pointAltitudes.end(), 0); - long double const distSquareSum = accumulate(pointDists.begin(), pointDists.end(), 0.0, - [](double x, double y) { return x + y * y; }); - long double distMultiplyAltSum = 0; - for (size_t i = 0; i < pointDists.size(); ++i) - distMultiplyAltSum += pointDists[i] * pointAltitudes[i]; - - if (distSum == 0.0 || distSquareSum == 0.0) + double k; + double b; + vector const pointAltitudesMeters(pointAltitudes.begin(), pointAltitudes.end()); + if (!LinearLeastSquaresFactors(pointDists, pointAltitudesMeters, k, b)) return; - double const b = (altSum/distSum - distMultiplyAltSum/distSquareSum) / (numPoints/distSum - distSum/distSquareSum); - double const k = (altSum - b * numPoints) / distSum; - for (uint32_t i = 0; i < numPoints; ++i) { - int32_t const deviation = - static_cast(ExpectedPointAltitude(k, b, pointDists[i])) - pointAltitudes[i]; + TAltitude const deviation = + static_cast(GetY(k, b, pointDists[i])) - pointAltitudes[i]; m_leastSquaresDiff[deviation]++; } } } }; -double CalculateEntropy(map const & diffFromLinear) +double CalculateEntropy(map const & diffFromLinear) { uint32_t innerPointCount = 0; for (auto const & f : diffFromLinear) @@ -328,7 +348,8 @@ int main(int argc, char ** argv) PrintCont(processor.m_diffFromLinear, "Altitude deviation of internal feature points from linear model.", " internal feature point(s) deviate from linear model with ", " meter(s)"); - PrintCont(processor.m_leastSquaresDiff, "Altitude deviation of feature points from least squares line.", + PrintCont(processor.m_leastSquaresDiff, + "Altitude deviation of feature points from least squares line.", " internal feature point(s) deviate from linear model with ", " meter(s)"); cout << endl << FLAGS_mwm_file_path << endl;