diff --git a/generator/altitude_generator.cpp b/generator/altitude_generator.cpp index 209c6d828..e30c0a7ce 100644 --- a/generator/altitude_generator.cpp +++ b/generator/altitude_generator.cpp @@ -37,7 +37,7 @@ public: // AltitudeGetter overrides: Altitude GetAltitude(m2::PointD const & p) override { - return m_srtmManager.GetTriangleHeight(mercator::ToLatLon(p)); + return m_srtmManager.GetAltitude(mercator::ToLatLon(p)); } void PrintStatsAndPurge() override diff --git a/generator/feature_segments_checker/feature_segments_checker.cpp b/generator/feature_segments_checker/feature_segments_checker.cpp index 807cddd25..7ac21f90d 100644 --- a/generator/feature_segments_checker/feature_segments_checker.cpp +++ b/generator/feature_segments_checker/feature_segments_checker.cpp @@ -188,7 +188,7 @@ public: for (uint32_t i = 0; i < numPoints; ++i) { // Feature segment altitude. - geometry::Altitude altitude = m_srtmManager.GetTriangleHeight(mercator::ToLatLon(f.GetPoint(i))); + geometry::Altitude altitude = m_srtmManager.GetAltitude(mercator::ToLatLon(f.GetPoint(i))); pointAltitudes[i] = altitude == geometry::kInvalidAltitude ? 0 : altitude; if (i == 0) { diff --git a/generator/generator_tests/srtm_parser_test.cpp b/generator/generator_tests/srtm_parser_test.cpp index ed0c086d5..427e7db59 100644 --- a/generator/generator_tests/srtm_parser_test.cpp +++ b/generator/generator_tests/srtm_parser_test.cpp @@ -87,21 +87,23 @@ UNIT_TEST(SRTM_TileTest) double const len = 1.0 / (sz - 1); TEST_EQUAL(tile.GetHeight({0, 0}), 4, ()); - TEST_EQUAL(tile.GetTriangleHeight({0, 0}), 4, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetTriangleHeight({0, 0}), 4.0, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetBilinearHeight({0, 0}), 4.0, ()); TEST_EQUAL(tile.GetHeight({len, len}), 2, ()); - TEST_EQUAL(tile.GetTriangleHeight({len, len}), 2, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetTriangleHeight({len, len}), 2.0, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetBilinearHeight({len, len}), 2.0, ()); - // Key difference here: GetHeight snaps on the nearest node, while GetTriangleHeight calculates - // from the nearest nodes' triangle. double l = len / 2; Altitude h = tile.GetHeight({l, l}); TEST(h == 4 || h == 2, (h)); - TEST_EQUAL(tile.GetTriangleHeight({l, l}), 3, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetTriangleHeight({l, l}), 3.0, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetBilinearHeight({l, l}), 3.0, ()); l = 3 * len + len / 2; h = tile.GetHeight({l, l}); TEST(h == -4 || h == -2, (h)); - TEST_EQUAL(tile.GetTriangleHeight({l, l}), -3, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetTriangleHeight({l, l}), -3.0, ()); + TEST_ALMOST_EQUAL_ULPS(tile.GetBilinearHeight({l, l}), -3.0, ()); } /* diff --git a/generator/srtm_coverage_checker/srtm_coverage_checker.cpp b/generator/srtm_coverage_checker/srtm_coverage_checker.cpp index 024197fc4..e2d54442d 100644 --- a/generator/srtm_coverage_checker/srtm_coverage_checker.cpp +++ b/generator/srtm_coverage_checker/srtm_coverage_checker.cpp @@ -48,10 +48,10 @@ public: return (routing::IsRoad(types) && !types.Has(m_ferry)); } - geometry::Altitude GetHeight(ms::LatLon const & coord) + geometry::Altitude GetAltitude(ms::LatLon const & coord) { std::lock_guard guard(m_mutex); - return m_manager.GetTriangleHeight(coord); + return m_manager.GetAltitude(coord); } void Purge() @@ -119,7 +119,7 @@ void CheckCoverage(SafeTileManager & manager) for (size_t i = 0; i < ft.GetPointsCount(); ++i) { - auto const height = manager.GetHeight(mercator::ToLatLon(ft.GetPoint(i))); + auto const height = manager.GetAltitude(mercator::ToLatLon(ft.GetPoint(i))); if (height != geometry::kInvalidAltitude) good++; } @@ -150,9 +150,9 @@ void CheckDistance(SafeTileManager & manager) for (size_t i = 1; i < ft.GetPointsCount(); ++i) { auto const ll1 = mercator::ToLatLon(ft.GetPoint(i-1)); - auto const alt1 = manager.GetHeight(ll1); + auto const alt1 = manager.GetAltitude(ll1); auto const ll2 = mercator::ToLatLon(ft.GetPoint(i)); - auto const alt2 = manager.GetHeight(ll2); + auto const alt2 = manager.GetAltitude(ll2); if (alt1 == geometry::kInvalidAltitude || alt2 == geometry::kInvalidAltitude) { @@ -171,7 +171,7 @@ void CheckDistance(SafeTileManager & manager) ms::LatLon const ll(ll2.m_lat * a + ll1.m_lat * (1 - a), ll2.m_lon * a + ll1.m_lon * (1 - a)); // Get diff between approx altitude and real one. - auto const alt = manager.GetHeight(ll); + auto const alt = manager.GetAltitude(ll); if (alt == geometry::kInvalidAltitude) { LOG_SHORT(LWARNING, ("Invalid altitude for the middle point:", ll)); @@ -242,7 +242,10 @@ int main(int argc, char * argv[]) } else { - cout << "H = " << manager.GetHeight({lat, lon}) << "; TrgH = " << manager.GetTriangleHeight({lat, lon}); + auto const & tile = manager.GetTile({lat, lon}); + cout << "H = " << tile.GetHeight({lat, lon}) << + "; Trg = " << tile.GetTriangleHeight({lat, lon}) << + "; Bilinear = " << tile.GetBilinearHeight({lat, lon}); cout << endl; } } diff --git a/generator/srtm_parser.cpp b/generator/srtm_parser.cpp index b989a5ecb..2ee57cdf4 100644 --- a/generator/srtm_parser.cpp +++ b/generator/srtm_parser.cpp @@ -135,7 +135,7 @@ geometry::Altitude SrtmTile::GetHeightRC(size_t row, size_t col) const return ReverseByteOrder(Data()[ix]); } -geometry::Altitude SrtmTile::GetTriangleHeight(ms::LatLon const & coord) const +double SrtmTile::GetTriangleHeight(ms::LatLon const & coord) const { if (!IsValid()) return geometry::kInvalidAltitude; @@ -179,9 +179,32 @@ geometry::Altitude SrtmTile::GetTriangleHeight(ms::LatLon const & coord) const double const a2 = ((p3.y - p1.y) * (ll.m_lon - p3.x) + (p1.x - p3.x) * (ll.m_lat - p3.y)) / det; double const a3 = 1 - a1 - a2; - return static_cast(std::round(a1 * GetHeightRC(p1.y, p1.x) + - a2 * GetHeightRC(p2.y, p2.x) + - a3 * GetHeightRC(p3.y, p3.x))); + return a1 * GetHeightRC(p1.y, p1.x) + a2 * GetHeightRC(p2.y, p2.x) + a3 * GetHeightRC(p3.y, p3.x); +} + +double SrtmTile::GetBilinearHeight(ms::LatLon const & coord) const +{ + if (!IsValid()) + return geometry::kInvalidAltitude; + + auto const ll = GetCoordInSeconds(coord); + + m2::Point const p1(static_cast(ll.m_lon), static_cast(ll.m_lat)); + auto p2 = p1; + if (p2.x < kArcSecondsInDegree) + ++p2.x; + if (p2.y < kArcSecondsInDegree) + ++p2.y; + + // https://en.wikipedia.org/wiki/Bilinear_interpolation + double const denom = (p2.x - p1.x) * (p2.y - p1.y); + if (denom == 0) + return GetHeightRC(p1.y, p1.x); + + return (GetHeightRC(p1.y, p1.x) * (p2.x - ll.m_lon) * (p2.y - ll.m_lat) + + GetHeightRC(p1.y, p2.x) * (ll.m_lon - p1.x) * (p2.y - ll.m_lat) + + GetHeightRC(p2.y, p1.x) * (p2.x - ll.m_lon) * (ll.m_lat - p1.y) + + GetHeightRC(p2.y, p2.x) * (ll.m_lon - p1.x) * (ll.m_lat - p1.y)) / denom; } // static diff --git a/generator/srtm_parser.hpp b/generator/srtm_parser.hpp index f92d063a0..7a86ebc01 100644 --- a/generator/srtm_parser.hpp +++ b/generator/srtm_parser.hpp @@ -30,8 +30,15 @@ public: /// @{ /// Nearest serialized height. geometry::Altitude GetHeight(ms::LatLon const & coord) const; - /// Height from underlying triangle (better than GetHeight). - geometry::Altitude GetTriangleHeight(ms::LatLon const & coord) const; + /// Triangle interpolation. + double GetTriangleHeight(ms::LatLon const & coord) const; + /// Bilinear interpolation. + double GetBilinearHeight(ms::LatLon const & coord) const; + + geometry::Altitude GetAltitude(ms::LatLon const & coord) const + { + return static_cast(std::round(GetBilinearHeight(coord))); + } /// @} using LatLonKey = std::pair; @@ -68,14 +75,9 @@ public: SrtmTile const & GetTile(ms::LatLon const & coord); - geometry::Altitude GetHeight(ms::LatLon const & coord) + geometry::Altitude GetAltitude(ms::LatLon const & coord) { - return GetTile(coord).GetHeight(coord); - } - - geometry::Altitude GetTriangleHeight(ms::LatLon const & coord) - { - return GetTile(coord).GetTriangleHeight(coord); + return GetTile(coord).GetAltitude(coord); } size_t GeTilesNumber() const { return m_tiles.size(); } diff --git a/topography_generator/generator.cpp b/topography_generator/generator.cpp index 95a10e0cc..41cd02cda 100644 --- a/topography_generator/generator.cpp +++ b/topography_generator/generator.cpp @@ -184,7 +184,7 @@ public: , m_bottomLat(bottomLat) {} - /// @todo Should we use the same approach as in SrtmTile::GetTriangleHeight? + /// @todo Should we use the same approach as in SrtmTile::GetTriangleHeight/GetBilinearHeight? /// This function is used in ASTER fiter only. Altitude GetValue(ms::LatLon const & pos) override {