Update ElevationLayer::assembleImage to precompute reprojection grid …

…like ImageLayer does.

src/osgEarth/ElevationLayer

@@ -69,9 +69,6 @@ namespace osgEarth
         void setNoDataPolicy(const ElevationNoDataPolicy& value);
         const ElevationNoDataPolicy& getNoDataPolicy() const;
 
-        //! Override from VisibleLayer
-        //virtual void setVisible(bool value);
-
     public: // methods
 
         /**
@@ -133,10 +130,11 @@ namespace osgEarth
         virtual ~ElevationLayer() { }
 
     private:
-        void assembleHeightField(
-            const TileKey& key,
-            osg::ref_ptr<osg::HeightField>& out_hf,
-            ProgressCallback* progress) const;
+        GeoHeightField assembleHeightField(const TileKey& key, ProgressCallback* progress);
+        //void assembleHeightField(
+        //    const TileKey& key,
+        //    osg::ref_ptr<osg::HeightField>& out_hf,
+        //    ProgressCallback* progress) const;
 
         void normalizeNoDataValues(osg::HeightField* hf) const;
 

src/osgEarth/ElevationLayer.cpp

@@ -22,6 +22,7 @@
 #include <osgEarth/MemCache>
 #include <osgEarth/Metrics>
 #include <osgEarth/NetworkMonitor>
+#include <osgEarth/Math>
 #include <cinttypes>
 
 using namespace osgEarth;
@@ -216,6 +217,145 @@ ElevationLayer::applyProfileOverrides(osg::ref_ptr<const Profile>& inOutProfile)
     }
 }
 
+GeoHeightField
+ElevationLayer::assembleHeightField(const TileKey& key, ProgressCallback* progress)
+{
+    // If we got here, assert that there's a non-null layer profile.
+    if (!getProfile())
+    {
+        setStatus(Status::Error(Status::AssertionFailure, "assembleImage with undefined profile"));
+        return { };
+    }
+
+    GeoHeightField output;
+    unsigned targetLOD = getProfile()->getEquivalentLOD(key.getProfile(), key.getLOD());
+
+    std::vector<TileKey> intersectingKeys;
+    getProfile()->getIntersectingTiles(key, intersectingKeys);
+
+    using KeyedSource = std::pair<TileKey, GeoHeightField>;
+    std::vector<KeyedSource> sources;
+
+    if (intersectingKeys.size() > 0)
+    {
+        bool hasAtLeastOneSourceAtTargetLOD = false;
+
+        for (auto& intersectingKey : intersectingKeys)
+        {
+            TileKey subKey = intersectingKey;
+            GeoHeightField subTile;
+            while (subKey.valid() && !subTile.valid())
+            {
+                subTile = createHeightFieldInKeyProfile(subKey, progress);
+                if (!subTile.valid())
+                    subKey.makeParent();
+
+                if (progress && progress->isCanceled())
+                    return {};
+            }
+
+            if (subTile.valid())
+            {
+                if (subKey.getLOD() == targetLOD)
+                {
+                    hasAtLeastOneSourceAtTargetLOD = true;
+                }
+
+                // got a valid image, so add it to our sources collection:
+                sources.emplace_back(subKey, subTile);
+            }
+        }
+
+        // If we actually got at least one piece of usable data,
+        // move ahead and build a mosaic of all sources.
+        if (hasAtLeastOneSourceAtTargetLOD)
+        {
+            unsigned cols = getTileSize();
+            unsigned rows = getTileSize();
+
+            // sort the sources by LOD (highest first).
+            std::sort(
+                sources.begin(), sources.end(),
+                [](const auto& lhs, const auto& rhs) {
+                    return lhs.first.getLOD() > rhs.first.getLOD();
+                });
+
+            // assume all tiles to mosaic are in the same SRS.
+            auto* key_srs = key.getExtent().getSRS();
+            auto* source_srs = sources[0].second.getSRS();
+
+            // new output:
+            auto mosaic = new osg::HeightField();
+            mosaic->allocate(cols, rows);
+            auto& heights = mosaic->getHeightList();
+            for(int i=0; i<cols*rows; ++i)
+                heights[i] = NO_DATA_VALUE;
+
+            // Working set of points. it's much faster to xform an entire vector all at once.
+            std::vector<osg::Vec3d> points;
+            points.resize(cols * rows);
+
+            double minx, miny, maxx, maxy;
+            key.getExtent().getBounds(minx, miny, maxx, maxy);
+            double dx = (maxx - minx) / (double)(cols-1);
+            double dy = (maxy - miny) / (double)(rows-1);
+
+            // Source bounds and projection grid. We use cols-1 and rows-1
+            // since we want to sample edge to edge, unlike with imagery
+            Bounds sourceBounds;
+            sources[0].second.getSRS()->getBounds(sourceBounds);
+
+            // build a grid of sample points:
+            for (unsigned t = 0; t < rows; ++t)
+            {
+                double y = miny + (dy * (double)t);
+                for (unsigned s = 0; s < cols; ++s)
+                {
+                    double x = minx + (dx * (double)s);
+                    points[t * cols + s] = { x, y, 0.0 };
+                }
+            }
+
+            // transform the sample points to the SRS of our source data tiles.
+            // NOTE: point.z() will hold a vertical offset if the layers' vdatums are different;
+            // we will add it back in later.
+            if (source_srs && key_srs)
+            {
+                key_srs->transform(points, source_srs);
+
+                if (sourceBounds.valid())
+                {
+                    for (auto& point : points)
+                    {
+                        point.x() = clamp(point.x(), sourceBounds.xMin(), sourceBounds.xMax());
+                        point.y() = clamp(point.y(), sourceBounds.yMin(), sourceBounds.yMax());
+                    }
+                }
+            }
+
+            // Mosaic our sources into a single output image.
+            for (int row = 0; row < rows; ++row)
+            {
+                for (int col = 0; col < cols; ++col)
+                {
+                    int i = row * cols + col;
+                    auto& h = mosaic->getHeight(col, row);
+                    for (unsigned k = 0; k < sources.size() && h == NO_DATA_VALUE; ++k)
+                    {
+                        h = sources[k].second.getElevation(points[i].x(), points[i].y(), INTERP_BILINEAR)
+                            - points[i].z(); // apply vdatum offset
+                    }
+                }
+            }
+
+            output = GeoHeightField(mosaic, key.getExtent());
+        }
+    }
+
+    return output;
+}
+
+#if 0
 void
 ElevationLayer::assembleHeightField(const TileKey& key,
                                     osg::ref_ptr<osg::HeightField>& out_hf,
@@ -354,6 +494,7 @@ ElevationLayer::assembleHeightField(const TileKey& key,
         out_hf = 0;
     }
 }
+#endif
 
 GeoHeightField
 ElevationLayer::createHeightField(const TileKey& key)
@@ -494,9 +635,7 @@ ElevationLayer::createHeightFieldInKeyProfile(const TileKey& key, ProgressCallba
             else
             {
                 // If the profiles are different, use a compositing method to assemble the tile.
-                osg::ref_ptr<osg::HeightField> hf;
-                assembleHeightField(key, hf, progress);
-                result = GeoHeightField(hf.get(), key.getExtent());
+                result = assembleHeightField(key, progress);
             }
 
             // Check for cancelation before writing to a cache

src/osgEarth/GeoData

@@ -896,6 +896,14 @@ namespace osgEarth
          */
         const GeoExtent& getExtent() const;
 
+        /**
+         * Shortcut to get the spatial reference system describing
+         * the projection of the image.
+         */
+        const SpatialReference* getSRS() const {
+            return getExtent().getSRS();
+        }
+
         /**
          * The minimum height in the heightfield
          */