[codex] Optimize capped elevation LOD lookup

src/applications/osgearth_benchmarks/main.cpp

@@ -5,7 +5,10 @@
 
 #include <benchmark/benchmark.h>
 
+#include <osgEarth/ElevationPool>
 #include <osgEarth/GeoData>
+#include <osgEarth/HeightFieldUtils>
+#include <osgEarth/Map>
 #include <osgEarth/SpatialReference>
 #include <osgEarth/StringUtils>
 #include <osgEarth/Cache>
@@ -16,6 +19,62 @@
 using namespace osgEarth;
 namespace fs = std::filesystem;
 
+namespace
+{
+    class ConstantElevationLayer : public ElevationLayer
+    {
+    public:
+        META_LayerNoOptions(osgEarth, ConstantElevationLayer, ElevationLayer, constant_elevation);
+
+        void init() override
+        {
+            ElevationLayer::init();
+            setProfile(Profile::create(Profile::GLOBAL_GEODETIC));
+            setMaxDataLevel(12u);
+            options().tileSize() = 9u;
+        }
+
+    protected:
+        GeoHeightField createHeightFieldImplementation(
+            const TileKey& key,
+            ProgressCallback* progress) const override
+        {
+            osg::ref_ptr<osg::HeightField> hf = HeightFieldUtils::createReferenceHeightField(
+                key.getExtent(),
+                getTileSize(),
+                getTileSize(),
+                0u,
+                false,
+                42.0f);
+
+            return GeoHeightField(hf.get(), key.getExtent());
+        }
+
+        virtual ~ConstantElevationLayer() { }
+    };
+
+    osg::ref_ptr<Map> createElevationBenchmarkMap(unsigned dataExtentCount)
+    {
+        osg::ref_ptr<Map> map = new Map();
+        map->setProfile(Profile::create(Profile::GLOBAL_GEODETIC));
+
+        osg::ref_ptr<ConstantElevationLayer> layer = new ConstantElevationLayer();
+
+        DataExtentList dataExtents;
+        dataExtents.reserve(dataExtentCount);
+        for (unsigned i = 0; i < dataExtentCount; ++i)
+        {
+            dataExtents.emplace_back(map->getProfile()->getExtent(), 0u, 12u);
+        }
+
+        layer->setDataExtents(dataExtents);
+        map->addLayer(layer.get());
+        map->getElevationPool()->setMap(map.get());
+
+        return map;
+    }
+}
+
 static void BM_GeoPointTransform(benchmark::State& state)
 {
     auto wgs84 = osgEarth::SpatialReference::get("wgs84");
@@ -58,6 +117,56 @@ static void BM_GeoExtentIntersects(benchmark::State& state)
 }
 BENCHMARK(BM_GeoExtentIntersects);
 
+static void BM_ElevationPoolSampleMapCoordsFixedResolution(benchmark::State& state)
+{
+    osg::ref_ptr<Map> map = createElevationBenchmarkMap(256u);
+    osg::ref_ptr<ElevationPool> pool = map->getElevationPool();
+
+    std::vector<osg::Vec3d> points;
+    points.reserve(state.range(0));
+    for (int i = 0; i < state.range(0); ++i)
+    {
+        double x = -1.0 + 2.0 * static_cast<double>(i % 64) / 63.0;
+        double y = -1.0 + 2.0 * static_cast<double>((i / 64) % 64) / 63.0;
+        points.emplace_back(x, y, 0.0);
+    }
+
+    const Distance resolution(10000000.0, Units::METERS);
+
+    int warmupCount = pool->sampleMapCoords(
+        points.begin(),
+        points.end(),
+        resolution,
+        nullptr,
+        nullptr);
+
+    if (warmupCount != static_cast<int>(points.size()))
+    {
+        state.SkipWithError("ElevationPool failed to sample all benchmark points");
+        return;
+    }
+
+    for (auto _ : state)
+    {
+        int count = pool->sampleMapCoords(
+            points.begin(),
+            points.end(),
+            resolution,
+            nullptr,
+            nullptr);
+
+        if (count != static_cast<int>(points.size()))
+        {
+            state.SkipWithError("ElevationPool failed to sample all benchmark points");
+            return;
+        }
+
+        benchmark::DoNotOptimize(count);
+        benchmark::ClobberMemory();
+    }
+}
+BENCHMARK(BM_ElevationPoolSampleMapCoordsFixedResolution)->Arg(4096);
+
 
 
 const int NUM_CACHE_IMAGES = 1000;

src/osgEarth/ElevationPool

@@ -240,7 +240,7 @@ namespace osgEarth
             ProgressCallback* progress);        
 
         //! Best LOD this a point for the working set, or -1 if no data in index
-        int getLOD(double x, double y, WorkingSet*);
+        int getLOD(double x, double y, WorkingSet*, int maxLOD);
 
         osg::ref_ptr<ElevationTile> getOrCreateRaster(
             MapData& snapshot,

src/osgEarth/ElevationPool.cpp

@@ -123,7 +123,7 @@ ElevationPool::setMap(const Map* map)
 }
 
 int
-ElevationPool::getLOD(double x, double y, WorkingSet* ws)
+ElevationPool::getLOD(double x, double y, WorkingSet* ws, int maxLOD)
 {
     double point[2] = { x, y };
     int maxiestMaxLevel = -1;
@@ -141,7 +141,9 @@ ElevationPool::getLOD(double x, double y, WorkingSet* ws)
             index->Search(point, point, [&](const unsigned& level)
                 {
                     maxiestMaxLevel = std::max(maxiestMaxLevel, (int)level);
-                    return RTREE_KEEP_SEARCHING;
+                    return maxiestMaxLevel >= maxLOD ?
+                        RTREE_STOP_SEARCHING :
+                        RTREE_KEEP_SEARCHING;
                 });
         }
     }
@@ -332,7 +334,7 @@ ElevationPool::prepareEnvelope(ElevationPool::Envelope& env, const GeoPoint& ref
         resolutionInMapUnits,
         ELEVATION_TILE_SIZE);
 
-    env._lod = std::min(getLOD(refPointMap.x(), refPointMap.y(), ws), (int)maxLOD);
+    env._lod = std::min(getLOD(refPointMap.x(), refPointMap.y(), ws, maxLOD), (int)maxLOD);
 
     // This can happen if the elevation data publishes no data extents
     if (env._lod < 0 && !_mapData.layers.empty())
@@ -662,7 +664,7 @@ ElevationPool::sampleMapCoords(
             resolutionInMapUnits,
             ELEVATION_TILE_SIZE);
 
-        lod = std::min(getLOD(p.x(), p.y(), ws), (int)computedLOD);
+        lod = std::min(getLOD(p.x(), p.y(), ws, computedLOD), (int)computedLOD);
 
         if (lod < 0)
         {
@@ -751,7 +753,7 @@ ElevationPool::getSample(const GeoPoint& p, unsigned maxLOD, WorkingSet* ws, Pro
     maxLOD = std::min(maxLOD, static_cast<unsigned>(std::numeric_limits<int>::max()));
 
     // returns the best LOD for the given point, or -1 if there is no data there
-    int lod = std::min(getLOD(p.x(), p.y(), ws), (int)maxLOD);
+    int lod = std::min(getLOD(p.x(), p.y(), ws, (int)maxLOD), (int)maxLOD);
 
     if (lod >= 0)
     {