Add various interpolation methods to XYZ elevation and metatile.

src/osgEarth/MetaTile

@@ -138,7 +138,7 @@ namespace osgEarth { namespace Util
 
         //! Read the interpolated value at the geospatial coordinate (x,y)
         //! (Bilinear interpolation is used)
-        inline bool readAtCoord(typename T::pixel_type& output, double x, double y);
+        inline bool readAtCoord(typename T::pixel_type& output, double x, double y, RasterInterpolation = INTERP_BILINEAR);
 
         //! The scale&bias of the tile relative to the key originally passed
         //! to setCenterTileKey
@@ -437,7 +437,7 @@ namespace osgEarth { namespace Util
     }
 
     template<typename T>
-    bool MetaTile<T>::readAtCoord(typename T::pixel_type& output, double x, double y)
+    bool MetaTile<T>::readAtCoord(typename T::pixel_type& output, double x, double y, RasterInterpolation interp)
     {
         // geoextents of center tile:
         double xmin, ymin, xmax, ymax;
@@ -456,6 +456,14 @@ namespace osgEarth { namespace Util
             return false;
         double p00_s = fract(u) * (double)_width, p00_t = fract(v) * (double)_height;
 
+        unsigned s00 = (unsigned)floor(p00_s);
+        unsigned t00 = (unsigned)floor(p00_t);
+
+        if (interp == INTERP_NEAREST)
+        {
+            return p00_tile->_data.read(output, s00, t00);
+        }
+
         // Neighbors:
         typename Tile* p10_tile = p00_tile;
         double p10_s = p00_s + 1, p10_t = p00_t;
@@ -483,8 +491,6 @@ namespace osgEarth { namespace Util
         }
 
         // calculate weights:
-        unsigned s00 = (unsigned)floor(p00_s);
-        unsigned t00 = (unsigned)floor(p00_t);
         double left_weight = 1.0 - (p00_s - (double)s00);
         double bottom_weight = 1.0 - (p00_t - (double)t00);
 
@@ -511,12 +517,24 @@ namespace osgEarth { namespace Util
                 return false;
             }
 
-            // bilinear:
-            typename T::pixel_type p0, p1;
-            p0 = mult(p00, left_weight) + mult(p10, 1.0 - left_weight);
-            p1 = mult(p01, left_weight) + mult(p11, 1.0 - left_weight);
-            output = mult(p0, bottom_weight) + mult(p1, 1.0 - bottom_weight);
-            return true;            
+            if (interp == INTERP_AVERAGE)
+            {
+                output =
+                    mult(p00, left_weight * bottom_weight) +
+                    mult(p10, (1.0 - left_weight) * bottom_weight) +
+                    mult(p01, left_weight * (1.0 - bottom_weight)) +
+                    mult(p11, (1.0 - left_weight) * (1.0 - bottom_weight));
+            }
+
+            else // INTERP_BILINEAR
+            {
+                typename T::pixel_type p0, p1;
+                p0 = mult(p00, left_weight) + mult(p10, 1.0 - left_weight);
+                p1 = mult(p01, left_weight) + mult(p11, 1.0 - left_weight);
+                output = mult(p0, bottom_weight) + mult(p1, 1.0 - bottom_weight);
+            }
+
+            return true;
         }
     }
 

src/osgEarth/XYZ

@@ -24,6 +24,7 @@
 #include <osgEarth/URI>
 #include <osgEarth/Threading>
 #include <osgEarth/Containers>
+#include <osgEarth/GeoCommon>
 
 /**
  * XYZ layers. These are general purpose tiled layers that conform
@@ -86,6 +87,7 @@ namespace osgEarth { namespace XYZ
         OE_OPTION(std::string, format, {});
         OE_OPTION(std::string, elevationEncoding, {});
         OE_OPTION(bool, stitchEdges, false);
+        OE_OPTION(RasterInterpolation, interpolation, INTERP_BILINEAR);
         static Config getMetadata();
         virtual Config getConfig() const;
     private:
@@ -170,7 +172,7 @@ namespace osgEarth
     class OSGEARTH_EXPORT XYZElevationLayer : public ElevationLayer
     {
     public:
-        typedef XYZ::XYZElevationLayerOptions Options;
+        using Options = XYZ::XYZElevationLayerOptions;
 
     public:
         META_Layer(osgEarth, XYZElevationLayer, Options, ElevationLayer, XYZElevation);
@@ -191,21 +193,35 @@ namespace osgEarth
         const std::string& getFormat() const;
 
         //! Encoding encoding type
+        //! @param value Encoding type, one of "" (default), "mapbox", "terrarium"
         void setElevationEncoding(const std::string& value);
         const std::string& getElevationEncoding() const;
 
+        //! Whether to stitch the edges of abutting tiles together to prevent gaps.
+        //! Some elevation data is image based, and the side of one tile does not match
+        //! the value on the correspond side of the abutting tile. This option will
+        //! resample the data to make the edges match, for a small performance hit.
+        //! @param value True to stitch edges (default is false)
+        void setStitchEdges(const bool& value);
+        const bool& getStitchEdges() const;
+
+        //! Interpolation method to use when stitching edges = true.
+        //! Default is INTERP_BILINEAR.
+        void setInterpolation(const RasterInterpolation& value);
+        const RasterInterpolation& getInterpolation() const;
+
     public: // Layer
 
         //! Establishes a connection to the XYZ data
-        virtual Status openImplementation() override;
+        Status openImplementation() override;
 
         //! Creates a heightfield for the given tile key
-        virtual GeoHeightField createHeightFieldImplementation(const TileKey& key, ProgressCallback* progress) const override;
+        GeoHeightField createHeightFieldImplementation(const TileKey& key, ProgressCallback* progress) const override;
 
     protected: // Layer
 
         //! Called by constructors
-        virtual void init();
+        void init() override;
 
     protected:
 

src/osgEarth/XYZ.cpp

@@ -271,6 +271,8 @@ OE_LAYER_PROPERTY_IMPL(XYZElevationLayer, URI, URL, url);
 OE_LAYER_PROPERTY_IMPL(XYZElevationLayer, bool, InvertY, invertY);
 OE_LAYER_PROPERTY_IMPL(XYZElevationLayer, std::string, Format, format);
 OE_LAYER_PROPERTY_IMPL(XYZElevationLayer, std::string, ElevationEncoding, elevationEncoding);
+OE_LAYER_PROPERTY_IMPL(XYZElevationLayer, bool, StitchEdges, stitchEdges);
+OE_LAYER_PROPERTY_IMPL(XYZElevationLayer, RasterInterpolation, Interpolation, interpolation);
 
 void
 XYZElevationLayer::init()
@@ -373,7 +375,7 @@ XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC
             return {};
         }
 
-        auto hf = new osg::HeightField();
+        osg::ref_ptr<osg::HeightField> hf = new osg::HeightField();
         hf->allocate(getTileSize(), getTileSize());
         std::fill(hf->getHeightList().begin(), hf->getHeightList().end(), NO_DATA_VALUE);
 
@@ -383,6 +385,9 @@ XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC
         osg::Vec4f pixel;
         for (int c = 0; c < getTileSize(); c++)
         {
+            if (progress && progress->isCanceled())
+                return {};
+
             double u = (double)c / (double)(getTileSize() - 1);
             double x = xmin + u * (xmax - xmin);
             for (int r = 0; r < getTileSize(); r++)
@@ -396,7 +401,7 @@ XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC
             }
         }
 
-        return GeoHeightField(hf, key.getExtent());
+        return GeoHeightField(hf.release(), key.getExtent());
     }
 
     else