Added SrtmTileManager::GetBilinearHeight.

Signed-off-by: Viktor Govako <viktor.govako@gmail.com>

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

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)
       {

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, ());
 }
 
 /*

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;
       }
     }

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<geometry::Altitude>(std::round(a1 * GetHeightRC(p1.y, p1.x) +
+  return a1 * GetHeightRC(p1.y, p1.x) + a2 * GetHeightRC(p2.y, p2.x) + a3 * GetHeightRC(p3.y, p3.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<int> const p1(static_cast<int>(ll.m_lon), static_cast<int>(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

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).
+  /// Triangle interpolation.
-  geometry::Altitude GetTriangleHeight(ms::LatLon const & coord) const;
+  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<geometry::Altitude>(std::round(GetBilinearHeight(coord)));
+  }
   /// @}
 
   using LatLonKey = std::pair<int32_t, int32_t>;
@@ -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);
+    return GetTile(coord).GetAltitude(coord);
-  }
-
-  geometry::Altitude GetTriangleHeight(ms::LatLon const & coord)
-  {
-    return GetTile(coord).GetTriangleHeight(coord);
   }
 
   size_t GeTilesNumber() const { return m_tiles.size(); }

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
   {