Fix vdatum issue with TMS and XYZ where an assembled HF was getting t…

…he vdatum applied more than once; fix missing negative values in int16 elevation dataset

src/osgEarth/ElevationLayer.cpp

@@ -635,6 +635,8 @@ ElevationLayer::createHeightFieldInKeyProfile(const TileKey& key, ProgressCallba
                 return GeoHeightField::INVALID;
             }
 
+            bool applyVerticalDatumTransformation = !key.getExtent().getSRS()->isVertEquivalentTo(profile->getSRS());
+
             if (key.getProfile()->isHorizEquivalentTo(profile.get()))
             {
                 Threading::ScopedReadLock lock(inUseMutex());
@@ -644,6 +646,9 @@ ElevationLayer::createHeightFieldInKeyProfile(const TileKey& key, ProgressCallba
             {
                 // If the profiles are different, use a compositing method to assemble the tile.
                 result = assembleHeightField(key, progress);
+
+                // assembleHeightField already does this, and more efficiently:
+                applyVerticalDatumTransformation = false;
             }
 
             // Check for cancelation before writing to a cache
@@ -671,10 +676,8 @@ ElevationLayer::createHeightFieldInKeyProfile(const TileKey& key, ProgressCallba
 
             // If the result is good, we now have a heightfield but its vertical values
             // are still relative to the source's vertical datum. Convert them.
-            if (hf.valid() && !key.getExtent().getSRS()->isVertEquivalentTo(profile->getSRS()))
+            if (applyVerticalDatumTransformation && hf.valid())
             {
-                OE_PROFILING_ZONE_NAMED("vdatum xform");
-
                 VerticalDatum::transform(
                     profile->getSRS()->getVerticalDatum(),    // from
                     key.getExtent().getSRS()->getVerticalDatum(),  // to

src/osgEarth/ImageToHeightFieldConverter

@@ -19,12 +19,10 @@
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>
 */
-
-#ifndef OSGEARTH_IMAGE_TO_HEIGHTFIELD_CONVERTER
-#define OSGEARTH_IMAGE_TO_HEIGHTFIELD_CONVERTER 1
+#pragma once
 
 #include <osgEarth/Export>
-
+#include <osgEarth/GeoCommon>
 #include <osg/Image>
 #include <osg/Shape>
 
@@ -36,10 +34,7 @@ namespace osgEarth { namespace Util
     class OSGEARTH_EXPORT ImageToHeightFieldConverter
     {
     public:
-        ImageToHeightFieldConverter();
-
-        /** dtor */
-        virtual ~ImageToHeightFieldConverter() { }
+        ImageToHeightFieldConverter() = default;
 
         /**
         * Instruct the converter to detect and replace "no data" values. It will
@@ -71,9 +66,7 @@ namespace osgEarth { namespace Util
         osg::Image* convert16(const osg::HeightField* hf ) const;
         osg::Image* convert32(const osg::HeightField* hf ) const;
 
-        bool  _replace_nodata;
-        float _nodata_value;
+        bool _replace_nodata = false;
+        float _nodata_value = NO_DATA_VALUE;
     };
 } }
-
-#endif //OSGEARTH_IMAGE_TO_HEIGHTFIELD_CONVERTER

src/osgEarth/ImageToHeightFieldConverter.cpp

@@ -37,13 +37,6 @@ namespace
     }
 }
 
-ImageToHeightFieldConverter::ImageToHeightFieldConverter():
-_replace_nodata( false ),
-_nodata_value( 0.0f )
-{
-  //NOP
-}
-
 void
 ImageToHeightFieldConverter::setRemoveNoDataValues( bool which, float f )
 {

src/osgEarth/ImageUtils.cpp

@@ -2204,6 +2204,17 @@ namespace
         }
     };
 
+    //! Select an appropriate reader for the given data type.
+    //! 
+    //! NOTE!!
+    //! Fopr UNSIGNED data types, we are using a SIGNED reader, except in the case
+    //! of GL_UNSIGNED_BYTE. This is because the OSG TIFF reader always declares
+    //! the data to be unsigned even when it's not. That's a bug, but it does not
+    //! hurt us to generate signed data for unsigned input, so this is an acceptable
+    //! workaround.
+    //! 
+    //! The exception is GL_UNSIGNED_BYTE which is commonly normalized into [0..1]
+    //! so we will maintain signed-ness for that.
     template<int GLFormat>
     inline ImageUtils::PixelReader::ReaderFunc
     chooseReader(GLenum dataType)
@@ -2215,13 +2226,13 @@ namespace
         case GL_UNSIGNED_BYTE:
             return &ColorReader<GLFormat, GLubyte>::read;
         case GL_SHORT:
-            return &ColorReader<GLFormat, GLshort>::read;
         case GL_UNSIGNED_SHORT:
-            return &ColorReader<GLFormat, GLushort>::read;
+            return &ColorReader<GLFormat, GLshort>::read;
+            //return &ColorReader<GLFormat, GLushort>::read;
         case GL_INT:
-            return &ColorReader<GLFormat, GLint>::read;
         case GL_UNSIGNED_INT:
-            return &ColorReader<GLFormat, GLuint>::read;
+            return &ColorReader<GLFormat, GLint>::read;
+            //return &ColorReader<GLFormat, GLuint>::read;
         case GL_FLOAT:
             return &ColorReader<GLFormat, GLfloat>::read;
         case GL_UNSIGNED_SHORT_5_5_5_1:
@@ -2235,6 +2246,7 @@ namespace
         }
     }
 
+    //! Selects a reader based on the input pixel format and type.
     inline ImageUtils::PixelReader::ReaderFunc
     getReader( GLenum pixelFormat, GLenum dataType )
     {

src/osgEarth/TMS.cpp

@@ -1416,23 +1416,26 @@ TMSElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC
         {
             const osg::Image* image = geoImage.getImage();
 
+#if 0
+            ImageToHeightFieldConverter converter;
+            osg::ref_ptr<osg::HeightField> hf = converter.convert(image);
+            return GeoHeightField(hf.get(), key.getExtent());
+
+#else
             // Allocate the heightfield.
             osg::HeightField* hf = new osg::HeightField();
             hf->allocate(image->s(), image->t());
 
-            ImageUtils::PixelReader reader(image);
+            ImageUtils::PixelReader read(image);
             osg::Vec4f pixel;
-
-            for (int c = 0; c < image->s(); c++)
-            {
-                for (int r = 0; r < image->t(); r++)
+            read.forEachPixel([&](auto& i)
                 {
-                    reader(pixel, c, r);
-                    hf->setHeight(c, r, decode(pixel));
-                }
-            }
+                    read(pixel, i.s(), i.t());
+                    hf->setHeight(i.s(), i.t(), decode(pixel));
+                });
 
             return GeoHeightField(hf, key.getExtent());
+#endif
         }
 
         return GeoHeightField::INVALID;

src/osgEarth/TileLayer.cpp

@@ -360,7 +360,7 @@ TileLayer::addedToMap(const Map* map)
         !map->getProfile()->getSRS()->isHorizEquivalentTo(getProfile()->getSRS()))
     {
         l2CacheSize = 16u;
-        OE_DEBUG << LC << "Map/Layer profiles differ; requesting L2 cache" << std::endl;
+        OE_INFO << LC << "Map/Layer profiles differ; requesting L2 cache" << std::endl;
     }
 
     // Use the user defined option if it's set.

src/osgEarth/XYZ.cpp

@@ -23,6 +23,7 @@
 #include <osgEarth/ImageToHeightFieldConverter>
 #include "MetaTile"
 #include <osgDB/FileUtils>
+#include <osgDB/WriteFile>
 
 using namespace osgEarth;
 using namespace osgEarth::XYZ;
@@ -311,7 +312,10 @@ XYZElevationLayer::openImplementation()
     if (status.isError())
         return status;
 
-    setProfile(_imageLayer->getProfile());            
+    setProfile(_imageLayer->getProfile());
+
+    DataExtentList de{ DataExtent(getProfile()->getExtent()) };
+    setDataExtents(de);
 
     return Status::NoError;
 }
@@ -416,16 +420,19 @@ XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC
             osg::HeightField* hf = new osg::HeightField();
             hf->allocate(image->s(), image->t());
 
-            ImageUtils::PixelReader reader(image);
+            ImageUtils::PixelReader read(image);
             osg::Vec4f pixel;
+            read.forEachPixel([&](auto& i)
+                {
+                    read(pixel, i.s(), i.t());
+                    hf->setHeight(i.s(), i.t(), decode(pixel));
+                });
 
-            for (int c = 0; c < image->s(); c++)
+            if (key.is(8, 70, 107))
             {
-                for (int r = 0; r < image->t(); r++)
-                {
-                    reader(pixel, c, r);
-                    hf->setHeight(c, r, decode(pixel));
-                }
+                ImageToHeightFieldConverter c;
+                auto out = c.convert(hf, 32);
+                osgDB::writeImageFile(*out, "test.tif");
             }
 
             return GeoHeightField(hf, key.getExtent());