diff --git a/src/osgEarth/Math b/src/osgEarth/Math index 6da88afa99..fc6adaa0c2 100644 --- a/src/osgEarth/Math +++ b/src/osgEarth/Math @@ -16,14 +16,13 @@ * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see */ - -#ifndef OSGEARTH_MATH_H -#define OSGEARTH_MATH_H 1 +#pragma once #include #include #include #include +#include #include #include #include @@ -512,6 +511,12 @@ namespace osgEarth return x + 1; } + template + inline osg::Vec4f mult(const osg::Vec4f& in, T scalar) + { + return osg::Vec4f(in.x()*scalar, in.y()*scalar, in.z()*scalar, in.w()*scalar); + } + // Adapted from Boost - see boost license // https://www.boost.org/users/license.html template inline std::size_t hash_value_unsigned(T val) { @@ -631,7 +636,4 @@ namespace osgEarth const osg::Matrix& m, double& L, double& R, double& B, double& T, double& N, double& F); }; - } -#endif - diff --git a/src/osgEarth/MetaTile b/src/osgEarth/MetaTile index 924c566148..f15108ab0b 100644 --- a/src/osgEarth/MetaTile +++ b/src/osgEarth/MetaTile @@ -27,6 +27,7 @@ #include #include #include +#include #include namespace osgEarth { namespace Util @@ -135,6 +136,10 @@ namespace osgEarth { namespace Util inline const typename T::pixel_type* read(double u, double v); inline const typename T::pixel_type* read(int s, int t); + //! 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, RasterInterpolation = INTERP_BILINEAR); + //! The scale&bias of the tile relative to the key originally passed //! to setCenterTileKey inline const osg::Matrix& getScaleBias() const { @@ -170,6 +175,9 @@ namespace osgEarth { namespace Util Grid _tiles; osg::Matrix _scale_bias; // scale/bias matrix of _centerKey unsigned _width, _height; + + typename Tile* getTile(int tx, int ty); + bool getTileAndPixelCoords(double u, double v, typename Tile*& tile, double& s, double& t); }; template @@ -402,6 +410,135 @@ namespace osgEarth { namespace Util return tile._data.read((unsigned)s, (unsigned)t); } + template + typename MetaTile::Tile* MetaTile::getTile(int tx, int ty) + { + Tile& tile = _tiles(tx, ty); + + // if we already tried to load this tile and failed, bail out + if (tile._failed) + return nullptr; + + // if we still need to load this tile, do so + if (!tile._data.valid() && _createTile != nullptr) + { + TileKey key = _centerKey.createNeighborKey(tx, -ty); + tile._data = _createTile(key, nullptr); + if (!tile._data.valid()) + { + tile._failed = true; + } + } + + if (tile._failed) + return nullptr; + + return &tile; + } + + template + bool MetaTile::readAtCoord(typename T::pixel_type& output, double x, double y, RasterInterpolation interp) + { + // geoextents of center tile: + double xmin, ymin, xmax, ymax; + _centerKey.getExtent().getBounds(xmin, ymin, xmax, ymax); + + // (u, v) ndc relative to center tile: + double u = (x - xmin) / (xmax - xmin); + double v = (y - ymin) / (ymax - ymin); + + int tx = (int)u; + int ty = (int)v; + + // Anchor: + typename Tile* p00_tile = getTile(tx, ty); + if (!p00_tile) + 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; + if (p10_s >= (double)_width) + { + p10_tile = getTile(tx + 1, ty); + p10_s -= (double)_width; + } + + typename Tile* p01_tile = p00_tile; + double p01_s = p00_s, p01_t = p00_t + 1; + if (p01_t >= _height) + { + p01_tile = getTile(tx, ty + 1); + p01_t -= (double)_height; + } + + typename Tile* p11_tile = p00_tile; + double p11_s = p00_s + 1, p11_t = p00_t + 1; + if (p11_s >= (double)_width || p11_t >= (double)_height) + { + p11_tile = getTile(tx + 1, ty + 1); + if (p11_s > (double)_width) p11_s -= (double)_width; + if (p11_t > (double)_height) p11_t -= (double)_height; + } + + // calculate weights: + double left_weight = 1.0 - (p00_s - (double)s00); + double bottom_weight = 1.0 - (p00_t - (double)t00); + + if (!p01_tile || !p10_tile || !p11_tile) + { + return p00_tile->_data.read(output, s00, t00); + } + else + { + typename T::pixel_type p00, p10, p01, p11; + + unsigned s01 = clamp((unsigned)floor(p01_s), 0u, _width-1); + unsigned t01 = clamp((unsigned)floor(p01_t), 0u, _height-1); + unsigned s10 = clamp((unsigned)floor(p10_s), 0u, _width-1); + unsigned t10 = clamp((unsigned)floor(p10_t), 0u, _height-1); + unsigned s11 = clamp((unsigned)floor(p11_s), 0u, _width-1); + unsigned t11 = clamp((unsigned)floor(p11_t), 0u, _height-1); + + if (!p00_tile->_data.read(p00, s00, t00) || + !p01_tile->_data.read(p01, s01, t01) || + !p10_tile->_data.read(p10, s10, t10) || + !p11_tile->_data.read(p11, s11, t11)) + { + return false; + } + + 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; + } + } + + } } #endif // OSGEARTH_METATILE_H diff --git a/src/osgEarth/XYZ b/src/osgEarth/XYZ index c9fab2f744..a89b942148 100644 --- a/src/osgEarth/XYZ +++ b/src/osgEarth/XYZ @@ -23,6 +23,8 @@ #include #include #include +#include +#include /** * XYZ layers. These are general purpose tiled layers that conform @@ -67,8 +69,8 @@ namespace osgEarth { namespace XYZ public: META_LayerOptions(osgEarth, XYZImageLayerOptions, ImageLayer::Options); OE_OPTION(URI, url); - OE_OPTION(bool, invertY); - OE_OPTION(std::string, format); + OE_OPTION(bool, invertY, false); + OE_OPTION(std::string, format, {}); static Config getMetadata(); virtual Config getConfig() const; private: @@ -81,9 +83,11 @@ namespace osgEarth { namespace XYZ public: META_LayerOptions(osgEarth, XYZElevationLayerOptions, ElevationLayer::Options); OE_OPTION(URI, url); - OE_OPTION(bool, invertY); - OE_OPTION(std::string, format); - OE_OPTION(std::string, elevationEncoding); + OE_OPTION(bool, invertY, false); + 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: @@ -168,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); @@ -189,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: @@ -212,6 +230,7 @@ namespace osgEarth private: osg::ref_ptr _imageLayer; + mutable Util::LRUCache _stitchingCache{ true, 64 }; }; } // namespace osgEarth diff --git a/src/osgEarth/XYZ.cpp b/src/osgEarth/XYZ.cpp index 21999a7ddd..d30a06fe16 100644 --- a/src/osgEarth/XYZ.cpp +++ b/src/osgEarth/XYZ.cpp @@ -21,6 +21,7 @@ #include #include #include +#include "MetaTile" #include using namespace osgEarth; @@ -156,6 +157,7 @@ XYZElevationLayerOptions::getConfig() const conf.set("format", _format); conf.set("invert_y", _invertY); conf.set("elevation_encoding", _elevationEncoding); + conf.set("stitch_edges", stitchEdges()); return conf; } @@ -167,6 +169,7 @@ XYZElevationLayerOptions::fromConfig(const Config& conf) conf.get("invert_y", _invertY); conf.get("elevation_encoding", _elevationEncoding); conf.get("interpretation", elevationEncoding()); // compat with QGIS + conf.get("stitch_edges", stitchEdges()); } Config @@ -218,6 +221,7 @@ XYZImageLayer::openImplementation() return parent; DataExtentList dataExtents; + Status status = _driver.open( options().url().get(), options().format().get(), @@ -233,6 +237,11 @@ XYZImageLayer::openImplementation() setProfile(Profile::create("spherical-mercator")); } + if (dataExtents.empty()) + { + dataExtents.push_back(DataExtent(getProfile()->getExtent())); + } + setDataExtents(dataExtents); return Status::NoError; @@ -262,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() @@ -308,39 +319,99 @@ XYZElevationLayer::openImplementation() GeoHeightField XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressCallback* progress) const { - // Make an image, then convert it to a heightfield - GeoImage geoImage = _imageLayer->createImageImplementation(key, progress); - if (geoImage.valid()) + if (!_imageLayer->mayHaveData(key)) + { + return GeoHeightField::INVALID; + } + + std::function decode; + + if (options().elevationEncoding() == "mapbox") + { + decode = [](const osg::Vec4f& p) -> float + { + return -10000.0f + ((p.r() * 255.0f * 256.0f * 256.0f + p.g() * 255.0f * 256.0f + p.b() * 255.0f) * 0.1f); + }; + } + else if (options().elevationEncoding() == "terrarium") + { + decode = [](const osg::Vec4f& p) -> float + { + return (p.r() * 255.0f * 256.0f + p.g() * 255.0f + p.b() * 255.0f / 256.0f) - 32768.0f; + }; + } + else + { + decode = [](const osg::Vec4f& p) -> float + { + return p.r(); + }; + } + + if (options().stitchEdges() == true) { - const osg::Image* image = geoImage.getImage(); + MetaTile metaImage; + metaImage.setCreateTileFunction([this](const TileKey& key, ProgressCallback* progress) + { + Util::LRUCache::Record r; + if (_stitchingCache.get(key, r)) + { + return r.value(); + } + else + { + GeoImage image = _imageLayer->createImage(key, progress); + if (image.valid()) + { + _stitchingCache.insert(key, image); + } + return image; + } + }); + metaImage.setCenterTileKey(key, progress); - if (options().elevationEncoding() == "mapbox") + if (!metaImage.getCenterTile().valid() || !metaImage.getScaleBias().isIdentity()) { - // Allocate the heightfield. - osg::HeightField* hf = new osg::HeightField(); - hf->allocate(image->s(), image->t()); + return {}; + } - ImageUtils::PixelReader reader(image); - osg::Vec4f pixel; + osg::ref_ptr hf = new osg::HeightField(); + hf->allocate(getTileSize(), getTileSize()); + std::fill(hf->getHeightList().begin(), hf->getHeightList().end(), NO_DATA_VALUE); - for (int c = 0; c < image->s(); c++) + double xmin, ymin, xmax, ymax; + key.getExtent().getBounds(xmin, ymin, xmax, ymax); + + 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++) { - for (int r = 0; r < image->t(); r++) + double v = (double)r / (double)(getTileSize() - 1); + double y = ymin + v * (ymax - ymin); + if (metaImage.readAtCoord(pixel, x, y)) { - reader(pixel, c, r); - pixel.r() *= 255.0; - pixel.g() *= 255.0; - pixel.b() *= 255.0; - float h = -10000.0f + ((pixel.r() * 256.0f * 256.0f + pixel.g() * 256.0f + pixel.b()) * 0.1f); - hf->setHeight(c, r, h); + hf->setHeight(c, r, decode(pixel)); } } - - return GeoHeightField(hf, key.getExtent()); } - else if (options().elevationEncoding() == "terrarium") + return GeoHeightField(hf.release(), key.getExtent()); + } + + else + { + // Make an image, then convert it to a heightfield + GeoImage geoImage = _imageLayer->createImageImplementation(key, progress); + if (geoImage.valid()) { + const osg::Image* image = geoImage.getImage(); + // Allocate the heightfield. osg::HeightField* hf = new osg::HeightField(); hf->allocate(image->s(), image->t()); @@ -353,25 +424,13 @@ XYZElevationLayer::createHeightFieldImplementation(const TileKey& key, ProgressC for (int r = 0; r < image->t(); r++) { reader(pixel, c, r); - pixel.r() *= 255.0; - pixel.g() *= 255.0; - pixel.b() *= 255.0; - float h = (pixel.r() * 256.0f + pixel.g() + pixel.b() / 256.0f) - 32768.0f; - hf->setHeight(c, r, h); + hf->setHeight(c, r, decode(pixel)); } } return GeoHeightField(hf, key.getExtent()); } - else - { - ImageToHeightFieldConverter conv; - osg::HeightField* hf = conv.convert(image); - return GeoHeightField(hf, key.getExtent()); - } + return GeoHeightField::INVALID; } - - - return GeoHeightField::INVALID; } diff --git a/tests/mapzen.earth b/tests/mapzen.earth new file mode 100644 index 0000000000..190f2de93e --- /dev/null +++ b/tests/mapzen.earth @@ -0,0 +1,28 @@ + + + + + + + https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png + terrarium + true + + + + + + +