diff --git a/data/scene/debugglobe/map_service_configs/TERRAIN.wms b/data/scene/debugglobe/map_service_configs/TERRAIN.wms index d755d1371a..6782144c50 100644 --- a/data/scene/debugglobe/map_service_configs/TERRAIN.wms +++ b/data/scene/debugglobe/map_service_configs/TERRAIN.wms @@ -3,4 +3,11 @@ http://198.102.45.23/arcgis/rest/services/worldelevation3d/terrain3d? GCS_Elevation + + -180.0 + 90.0 + 180.0 + -90.0 + bottom + diff --git a/modules/globebrowsing/geodetics/geodetic2.cpp b/modules/globebrowsing/geodetics/geodetic2.cpp index 20349824a4..8f5420208a 100644 --- a/modules/globebrowsing/geodetics/geodetic2.cpp +++ b/modules/globebrowsing/geodetics/geodetic2.cpp @@ -23,9 +23,8 @@ ****************************************************************************************/ #include - -#include #include +#include #include @@ -33,260 +32,285 @@ #include namespace { - const std::string _loggerCat = "Geodetic2"; + const std::string _loggerCat = "Geodetic2"; } namespace openspace { - ////////////////////////////////////////////////////////////////////////////////////// - // GEODETIC2 // - ////////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////////// + // GEODETIC2 // + ////////////////////////////////////////////////////////////////////////////////////// - Geodetic2::Geodetic2() - : lat(0) - , lon(0) - {} + Geodetic2::Geodetic2() + : lat(0) + , lon(0) + {} - Geodetic2::Geodetic2(Scalar latitude, Scalar longitude) - : lat(latitude) - , lon(longitude) - { - - } + Geodetic2::Geodetic2(Scalar latitude, Scalar longitude) + : lat(latitude) + , lon(longitude) + { + + } - Geodetic2::Geodetic2(const Geodetic2& p) - : Geodetic2(p.lat, p.lon) - { - - } - - Vec2 Geodetic2::toLonLatVec2() const { - return Vec2(lon, lat); - } - - bool Geodetic2::operator==(const Geodetic2& other) const { - return lat == other.lat && lon == other.lon; - } - - Geodetic2 Geodetic2 ::operator+(const Geodetic2& other) const { - return Geodetic2(lat - other.lat, lon - other.lon); - } - - Geodetic2 Geodetic2 ::operator-(const Geodetic2& other) const { - return Geodetic2(lat - other.lat, lon - other.lon); - } - - Geodetic2 Geodetic2::operator*(Scalar scalar) const { - return Geodetic2(lat * scalar, lon * scalar); - } - - Geodetic2 Geodetic2::operator/(Scalar scalar) const { - return Geodetic2(lat / scalar, lon / scalar); - } + Geodetic2::Geodetic2(const Geodetic2& p) + : Geodetic2(p.lat, p.lon) + { + + } - ////////////////////////////////////////////////////////////////////////////////////// - // GEODETICPATCH // - ////////////////////////////////////////////////////////////////////////////////////// + Vec2 Geodetic2::toLonLatVec2() const { + return Vec2(lon, lat); + } - GeodeticPatch::GeodeticPatch( - Scalar centerLat, - Scalar centerLon, - Scalar halfSizeLat, - Scalar halfSizeLon) - : _center(Geodetic2(centerLat, centerLon)) - , _halfSize(Geodetic2(halfSizeLat, halfSizeLon)) - { - - } + bool Geodetic2::operator==(const Geodetic2& other) const { + return lat == other.lat && lon == other.lon; + } - GeodeticPatch::GeodeticPatch( - const Geodetic2& center, - const Geodetic2& halfSize) - : _center(center) - , _halfSize(halfSize) - { - - } + Geodetic2 Geodetic2::operator+(const Geodetic2& other) const { + return Geodetic2(lat + other.lat, lon + other.lon); + } - GeodeticPatch::GeodeticPatch(const GeodeticPatch& patch) - : _center(patch._center) - , _halfSize(patch._halfSize) - { - - } + Geodetic2 Geodetic2::operator-(const Geodetic2& other) const { + return Geodetic2(lat - other.lat, lon - other.lon); + } - void GeodeticPatch::setCenter(const Geodetic2& center) { - _center = center; - } + Geodetic2 Geodetic2::operator*(Scalar scalar) const { + return Geodetic2(lat * scalar, lon * scalar); + } - void GeodeticPatch::setHalfSize(const Geodetic2& halfSize) { - _halfSize = halfSize; - } - - Scalar GeodeticPatch::minimalBoundingRadius(const Ellipsoid& ellipsoid) const { - // TODO: THIS FUNCTION IS CURRENTLY ERROR PRONE SINCE THE PATCH IS NOW COVERING - // A PART OF AN ELLIPSOID AND NOT A SPHERE!MUST CHECK IF THIS FUNCTION IS STILL - // VALID. - const Geodetic2& cornerNearEquator = _center.lat > 0 ? southWestCorner() : northWestCorner(); - return glm::length(ellipsoid.geodetic2ToCartesian(_center) - ellipsoid.geodetic2ToCartesian(cornerNearEquator)); - } - /* - Scalar GeodeticPatch::unitArea() const { - Scalar deltaTheta = 2 * _halfSize.lon; - Scalar phiMin = _center.lat - _halfSize.lat; - Scalar phiMax = _center.lat + _halfSize.lat; - return deltaTheta * (sin(phiMax) - sin(phiMin)); - } - */ - const Geodetic2& GeodeticPatch::center() const { - return _center; - } - - const Geodetic2& GeodeticPatch::halfSize() const { - return _halfSize; - } - - Geodetic2 GeodeticPatch::size() const { - return Geodetic2(2 * _halfSize.lat, 2 * _halfSize.lon); - } - - Geodetic2 GeodeticPatch::northWestCorner() const{ - return Geodetic2(_center.lat + _halfSize.lat, _center.lon - _halfSize.lon); - } - - Geodetic2 GeodeticPatch::northEastCorner() const{ - return Geodetic2(_center.lat + _halfSize.lat, _center.lon + _halfSize.lon); - } - - Geodetic2 GeodeticPatch::southWestCorner() const{ - return Geodetic2(_center.lat - _halfSize.lat, _center.lon - _halfSize.lon); - } - - Geodetic2 GeodeticPatch::southEastCorner() const{ - return Geodetic2(_center.lat - _halfSize.lat, _center.lon + _halfSize.lon); - } + Geodetic2 Geodetic2::operator/(Scalar scalar) const { + return Geodetic2(lat / scalar, lon / scalar); + } - Geodetic2 GeodeticPatch::clamp(const Geodetic2& p) const { - using Ang = Angle; - - // Convert to Angles for normalization - Ang centerLat = Ang::fromRadians(_center.lat); - Ang centerLon = Ang::fromRadians(_center.lon); - Ang pointLat = Ang::fromRadians(p.lat); - Ang pointLon = Ang::fromRadians(p.lon); - - // Normalize w.r.t. the center in order for the clamping to done correctly - // - // Example: - // centerLat = 0 deg, halfSize.lat = 10 deg, pointLat = 330 deg - // --> Just clamping pointLat would be clamp(330, -10, 10) = 10 // WRONG! - // Instead, if we first normalize 330 deg around 0, we get -30 deg - // --> clamp(-30, -10, 10) = -10 // CORRECT! - pointLat.normalizeAround(centerLat); - pointLon.normalizeAround(centerLon); - - // get clamp bounds - Geodetic2 max = northEastCorner(); - Geodetic2 min = southWestCorner(); - - return Geodetic2( - glm::clamp(pointLat.asRadians(), min.lat, max.lat), - glm::clamp(pointLon.asRadians(), min.lon, max.lon) - ); - } - Geodetic2 GeodeticPatch::closestCorner(const Geodetic2& p) const { - using Ang = Angle; - // LatLon vector from patch center to the point - Geodetic2 centerToPoint = p - _center; - - // Normalize the difference angles to be centered around 0. - Ang latDiff = Ang::fromRadians(centerToPoint.lat).normalizeAround(Ang::ZERO); - Ang lonDiff = Ang::fromRadians(centerToPoint.lon).normalizeAround(Ang::ZERO); - - // If latDiff > 0 - // --> point p is north of the patch center - // --> the closest corner to the point must be a northern one - // --> set the corner's latitude coordinate to center.lat + halfSize.lat - // else - // --> set corner's latidude coordinate to center.lat - halfSize.lat - Scalar cornerLat = _center.lat + _halfSize.lat * (latDiff > Ang::ZERO ? 1 : -1); + ////////////////////////////////////////////////////////////////////////////////////// + // TILE INDEX // + ////////////////////////////////////////////////////////////////////////////////////// - // We then assigned the corner's longitude coordinate in a similar fashion - Scalar cornerLon = _center.lon + _halfSize.lon * (lonDiff > Ang::ZERO ? 1 : -1); + HashKey GeodeticTileIndex::hashKey() const{ + return x ^ (y << 16) ^ (level << 21); + } - return Geodetic2(cornerLat, cornerLon); - } - - Geodetic2 GeodeticPatch::closestPoint(const Geodetic2& p) const { - // This method finds the closest point on the patch, to the provided - // point p. As we are deali ng with latitude-longitude patches, distance in this - // context refers to great-circle distance. - // (https://en.wikipedia.org/wiki/Great-circle_distance) - // - // This uses a simple clamping approach to find the closest point on the - // patch. A naive castesian clamp is not sufficient for this purpose, - // as illustrated with an example below. + ////////////////////////////////////////////////////////////////////////////////////// + // GEODETICPATCH // + ////////////////////////////////////////////////////////////////////////////////////// - // Example: (degrees are used for latidude, longitude) - // patchCenter = (0,0), patchHalfSize = (45,45), point = (5, 170) - // Note, the point and the patch are on opposite sides of the sphere - // - // cartesian clamp: - // --> clampedPointLat = clamp(5, -45, 45) = 5 - // --> clampedPointLon = clamp(170, -45, 45) = 45 - // --> result: (5, 45) - // --> closest point is actually (45, 45) - // --> The error is significant - // - // This method simply adds an extra clamp on the latitude in these cases. In the - // above example, that would be the following: - // --> clampedPointLat = clamp(180 - 5, -45, 45) = 45 - // - // Just doing this actually makes points returned from this methods being the - // true closest point, great-circle distance-wise. - + GeodeticPatch::GeodeticPatch( + Scalar centerLat, + Scalar centerLon, + Scalar halfSizeLat, + Scalar halfSizeLon) + : _center(Geodetic2(centerLat, centerLon)) + , _halfSize(Geodetic2(halfSizeLat, halfSizeLon)) + { + + } - using Ang = Angle; + GeodeticPatch::GeodeticPatch( + const Geodetic2& center, + const Geodetic2& halfSize) + : _center(center) + , _halfSize(halfSize) + { + + } - // Convert to Angles for normalization - Ang centerLat = Ang::fromRadians(_center.lat); - Ang centerLon = Ang::fromRadians(_center.lon); - Ang pointLat = Ang::fromRadians(p.lat); - Ang pointLon = Ang::fromRadians(p.lon); + GeodeticPatch::GeodeticPatch(const GeodeticPatch& patch) + : _center(patch._center) + , _halfSize(patch._halfSize) + { + + } - // Normalize point with respect to center. This is done because the point - // will later be clamped. See LatLonPatch::clamp(const LatLon&) for explanation - pointLat.normalizeAround(centerLat); - pointLon.normalizeAround(centerLon); - // Calculate the longitud difference between center and point. We normalize around - // zero because we want the "shortest distance" difference, i.e the difference - // should be in the interval [-180 deg, 180 deg] - Ang centerToPointLon = (centerLon - pointLon).normalizeAround(Ang::ZERO); + GeodeticPatch::GeodeticPatch(const GeodeticTileIndex& tileIndex) { + Scalar deltaLat = (2*M_PI) / ((double)(1 << tileIndex.level)); + Scalar deltaLon = (2*M_PI) / ((double)(1 << tileIndex.level)); + Geodetic2 nwCorner(M_PI / 2 - deltaLat * tileIndex.y, -M_PI + deltaLon * tileIndex.x); + _halfSize = Geodetic2(deltaLat / 2, deltaLon / 2); + _center = Geodetic2(nwCorner.lat - _halfSize.lat, nwCorner.lon + _halfSize.lon); + } - // Calculate the longitudinal distance to the closest patch edge - Ang longitudeDistanceToClosestPatchEdge = centerToPointLon.abs() - Ang::fromRadians(_halfSize.lon); - // get clamp bounds - Geodetic2 max = northEastCorner(); - Geodetic2 min = southWestCorner(); - // If the longitude distance to the closest patch edge is larger than 90 deg - // the latitude will have to be clamped to its closest corner, as explained in - // the example above. - Scalar clampedLat = longitudeDistanceToClosestPatchEdge > Ang::QUARTER ? - clampedLat = glm::clamp((Ang::HALF - pointLat).normalizeAround(centerLat).asRadians(), min.lat, max.lat) : - clampedLat = glm::clamp(pointLat.asRadians(), min.lat, max.lat); - // Longitude is just clamped normally - Scalar clampedLon = glm::clamp(pointLon.asRadians(), min.lon, max.lon); + void GeodeticPatch::setCenter(const Geodetic2& center) { + _center = center; + } - return Geodetic2(clampedLat, clampedLon); - } + void GeodeticPatch::setHalfSize(const Geodetic2& halfSize) { + _halfSize = halfSize; + } + + Scalar GeodeticPatch::minimalBoundingRadius(const Ellipsoid& ellipsoid) const { + // TODO: THIS FUNCTION IS CURRENTLY ERROR PRONE SINCE THE PATCH IS NOW COVERING + // A PART OF AN ELLIPSOID AND NOT A SPHERE!MUST CHECK IF THIS FUNCTION IS STILL + // VALID. + const Geodetic2& cornerNearEquator = _center.lat > 0 ? southWestCorner() : northWestCorner(); + return glm::length(ellipsoid.geodetic2ToCartesian(_center) - ellipsoid.geodetic2ToCartesian(cornerNearEquator)); + } + /* + Scalar GeodeticPatch::unitArea() const { + Scalar deltaTheta = 2 * _halfSize.lon; + Scalar phiMin = _center.lat - _halfSize.lat; + Scalar phiMax = _center.lat + _halfSize.lat; + return deltaTheta * (sin(phiMax) - sin(phiMin)); + } + */ + const Geodetic2& GeodeticPatch::center() const { + return _center; + } + + const Geodetic2& GeodeticPatch::halfSize() const { + return _halfSize; + } + + Geodetic2 GeodeticPatch::size() const { + return Geodetic2(2 * _halfSize.lat, 2 * _halfSize.lon); + } + + Geodetic2 GeodeticPatch::northWestCorner() const{ + return Geodetic2(_center.lat + _halfSize.lat, _center.lon - _halfSize.lon); + } + + Geodetic2 GeodeticPatch::northEastCorner() const{ + return Geodetic2(_center.lat + _halfSize.lat, _center.lon + _halfSize.lon); + } + + Geodetic2 GeodeticPatch::southWestCorner() const{ + return Geodetic2(_center.lat - _halfSize.lat, _center.lon - _halfSize.lon); + } + + Geodetic2 GeodeticPatch::southEastCorner() const{ + return Geodetic2(_center.lat - _halfSize.lat, _center.lon + _halfSize.lon); + } + + + Geodetic2 GeodeticPatch::clamp(const Geodetic2& p) const { + using Ang = Angle; + + // Convert to Angles for normalization + Ang centerLat = Ang::fromRadians(_center.lat); + Ang centerLon = Ang::fromRadians(_center.lon); + Ang pointLat = Ang::fromRadians(p.lat); + Ang pointLon = Ang::fromRadians(p.lon); + + // Normalize w.r.t. the center in order for the clamping to done correctly + // + // Example: + // centerLat = 0 deg, halfSize.lat = 10 deg, pointLat = 330 deg + // --> Just clamping pointLat would be clamp(330, -10, 10) = 10 // WRONG! + // Instead, if we first normalize 330 deg around 0, we get -30 deg + // --> clamp(-30, -10, 10) = -10 // CORRECT! + pointLat.normalizeAround(centerLat); + pointLon.normalizeAround(centerLon); + + // get clamp bounds + Geodetic2 max = northEastCorner(); + Geodetic2 min = southWestCorner(); + + return Geodetic2( + glm::clamp(pointLat.asRadians(), min.lat, max.lat), + glm::clamp(pointLon.asRadians(), min.lon, max.lon) + ); + } + + + Geodetic2 GeodeticPatch::closestCorner(const Geodetic2& p) const { + using Ang = Angle; + + // LatLon vector from patch center to the point + Geodetic2 centerToPoint = p - _center; + + // Normalize the difference angles to be centered around 0. + Ang latDiff = Ang::fromRadians(centerToPoint.lat).normalizeAround(Ang::ZERO); + Ang lonDiff = Ang::fromRadians(centerToPoint.lon).normalizeAround(Ang::ZERO); + + // If latDiff > 0 + // --> point p is north of the patch center + // --> the closest corner to the point must be a northern one + // --> set the corner's latitude coordinate to center.lat + halfSize.lat + // else + // --> set corner's latidude coordinate to center.lat - halfSize.lat + Scalar cornerLat = _center.lat + _halfSize.lat * (latDiff > Ang::ZERO ? 1 : -1); + + // We then assigned the corner's longitude coordinate in a similar fashion + Scalar cornerLon = _center.lon + _halfSize.lon * (lonDiff > Ang::ZERO ? 1 : -1); + + return Geodetic2(cornerLat, cornerLon); + } + + + Geodetic2 GeodeticPatch::closestPoint(const Geodetic2& p) const { + // This method finds the closest point on the patch, to the provided + // point p. As we are deali ng with latitude-longitude patches, distance in this + // context refers to great-circle distance. + // (https://en.wikipedia.org/wiki/Great-circle_distance) + // + // This uses a simple clamping approach to find the closest point on the + // patch. A naive castesian clamp is not sufficient for this purpose, + // as illustrated with an example below. + + // Example: (degrees are used for latidude, longitude) + // patchCenter = (0,0), patchHalfSize = (45,45), point = (5, 170) + // Note, the point and the patch are on opposite sides of the sphere + // + // cartesian clamp: + // --> clampedPointLat = clamp(5, -45, 45) = 5 + // --> clampedPointLon = clamp(170, -45, 45) = 45 + // --> result: (5, 45) + // --> closest point is actually (45, 45) + // --> The error is significant + // + // This method simply adds an extra clamp on the latitude in these cases. In the + // above example, that would be the following: + // --> clampedPointLat = clamp(180 - 5, -45, 45) = 45 + // + // Just doing this actually makes points returned from this methods being the + // true closest point, great-circle distance-wise. + + + using Ang = Angle; + + // Convert to Angles for normalization + Ang centerLat = Ang::fromRadians(_center.lat); + Ang centerLon = Ang::fromRadians(_center.lon); + Ang pointLat = Ang::fromRadians(p.lat); + Ang pointLon = Ang::fromRadians(p.lon); + + // Normalize point with respect to center. This is done because the point + // will later be clamped. See LatLonPatch::clamp(const LatLon&) for explanation + pointLat.normalizeAround(centerLat); + pointLon.normalizeAround(centerLon); + + // Calculate the longitud difference between center and point. We normalize around + // zero because we want the "shortest distance" difference, i.e the difference + // should be in the interval [-180 deg, 180 deg] + Ang centerToPointLon = (centerLon - pointLon).normalizeAround(Ang::ZERO); + + // Calculate the longitudinal distance to the closest patch edge + Ang longitudeDistanceToClosestPatchEdge = centerToPointLon.abs() - Ang::fromRadians(_halfSize.lon); + + // get clamp bounds + Geodetic2 max = northEastCorner(); + Geodetic2 min = southWestCorner(); + + // If the longitude distance to the closest patch edge is larger than 90 deg + // the latitude will have to be clamped to its closest corner, as explained in + // the example above. + Scalar clampedLat = longitudeDistanceToClosestPatchEdge > Ang::QUARTER ? + clampedLat = glm::clamp((Ang::HALF - pointLat).normalizeAround(centerLat).asRadians(), min.lat, max.lat) : + clampedLat = glm::clamp(pointLat.asRadians(), min.lat, max.lat); + + // Longitude is just clamped normally + Scalar clampedLon = glm::clamp(pointLon.asRadians(), min.lon, max.lon); + + return Geodetic2(clampedLat, clampedLon); + } } // namespace openspace diff --git a/modules/globebrowsing/geodetics/geodetic2.h b/modules/globebrowsing/geodetics/geodetic2.h index c1b9b00579..7d8b7c978d 100644 --- a/modules/globebrowsing/geodetics/geodetic2.h +++ b/modules/globebrowsing/geodetics/geodetic2.h @@ -31,6 +31,7 @@ #include + // Using double precision typedef double Scalar; typedef glm::dvec2 Vec2; @@ -38,98 +39,135 @@ typedef glm::dvec3 Vec3; namespace openspace { - // Forward declaration - class Ellipsoid; + // Forward declaration + class Ellipsoid; struct Geodetic2 { - Geodetic2(); - Geodetic2(Scalar latitude, Scalar longitude); - Geodetic2(const Geodetic2& src); - - /* - static Geodetic2 fromCartesian(const Vec3& v); - Vec3 asUnitCartesian() const; - */ + Geodetic2(); + Geodetic2(Scalar latitude, Scalar longitude); + Geodetic2(const Geodetic2& src); - Vec2 toLonLatVec2() const; + + /* + static Geodetic2 fromCartesian(const Vec3& v); + Vec3 asUnitCartesian() const; + */ - inline bool operator==(const Geodetic2& other) const; - inline bool operator!=(const Geodetic2& other) const { return !(*this == (other)); } - inline Geodetic2 operator+(const Geodetic2& other) const; - inline Geodetic2 operator-(const Geodetic2& other) const; - inline Geodetic2 operator*(Scalar scalar) const; - inline Geodetic2 operator/(Scalar scalar) const; - Scalar lat; - Scalar lon; + Vec2 toLonLatVec2() const; + + bool operator==(const Geodetic2& other) const; + bool operator!=(const Geodetic2& other) const { return !(*this == (other)); } + + Geodetic2 operator+(const Geodetic2& other) const; + Geodetic2 operator-(const Geodetic2& other) const; + Geodetic2 operator*(Scalar scalar) const; + Geodetic2 operator/(Scalar scalar) const; + + Scalar lat; + Scalar lon; }; + + struct Geodetic3 { - Geodetic2 geodetic2; - Scalar height; + Geodetic2 geodetic2; + Scalar height; }; + + + + +////////////////////////////////////////////////////////////////////////////////////// +// TILE INDEX // +////////////////////////////////////////////////////////////////////////////////////// + + + +using HashKey = unsigned long; + + +struct GeodeticTileIndex { + + + int x, y, level; + + HashKey hashKey() const; +}; + + + + +////////////////////////////////////////////////////////////////////////////////////// +// GEODETICPATCH // +////////////////////////////////////////////////////////////////////////////////////// class GeodeticPatch { public: - GeodeticPatch( - Scalar centerLat, - Scalar centerLon, - Scalar halfSizeLat, - Scalar halfSizeLon); - GeodeticPatch( - const Geodetic2& center, - const Geodetic2& halfSize); - GeodeticPatch(const GeodeticPatch& patch); + GeodeticPatch( + Scalar centerLat, + Scalar centerLon, + Scalar halfSizeLat, + Scalar halfSizeLon); - void setCenter(const Geodetic2&); - void setHalfSize(const Geodetic2&); + GeodeticPatch( + const Geodetic2& center, + const Geodetic2& halfSize); - /** - Returns the minimal bounding radius that together with the LatLonPatch's - center point represents a sphere in which the patch is completely contained. - - TODO : THIS FUNCTION IS CURRENTLY ERROR PRONE SINCE THE PATCH IS NOW COVERING - A PART OF AN ELLIPSOID AND NOT A SPHERE! MUST CHECK IF THIS FUNCTION IS STILL - VALID. - */ - Scalar minimalBoundingRadius(const Ellipsoid& ellipsoid) const; + GeodeticPatch(const GeodeticPatch& patch); - /** - Returns the area of the patch with unit radius - */ - // Scalar unitArea() const; + GeodeticPatch::GeodeticPatch(const GeodeticTileIndex& tileIndex); - Geodetic2 northWestCorner() const; - Geodetic2 northEastCorner() const; - Geodetic2 southWestCorner() const; - Geodetic2 southEastCorner() const; + void setCenter(const Geodetic2&); + void setHalfSize(const Geodetic2&); - /** - * Clamps a point to the patch region - */ - Geodetic2 clamp(const Geodetic2& p) const; + /** + Returns the minimal bounding radius that together with the LatLonPatch's + center point represents a sphere in which the patch is completely contained. + + TODO : THIS FUNCTION IS CURRENTLY ERROR PRONE SINCE THE PATCH IS NOW COVERING + A PART OF AN ELLIPSOID AND NOT A SPHERE! MUST CHECK IF THIS FUNCTION IS STILL + VALID. + */ + Scalar minimalBoundingRadius(const Ellipsoid& ellipsoid) const; - /** - * Returns the corner of the patch that is closest to the given point p - */ - Geodetic2 closestCorner(const Geodetic2& p) const; + /** + Returns the area of the patch with unit radius + */ + // Scalar unitArea() const; - /** - * Returns a point on the patch that minimizes the great-circle distance to - * the given point p. - */ - Geodetic2 closestPoint(const Geodetic2& p) const; - - const Geodetic2& center() const; - const Geodetic2& halfSize() const; - Geodetic2 size() const; + Geodetic2 northWestCorner() const; + Geodetic2 northEastCorner() const; + Geodetic2 southWestCorner() const; + Geodetic2 southEastCorner() const; + + /** + * Clamps a point to the patch region + */ + Geodetic2 clamp(const Geodetic2& p) const; + + /** + * Returns the corner of the patch that is closest to the given point p + */ + Geodetic2 closestCorner(const Geodetic2& p) const; + + /** + * Returns a point on the patch that minimizes the great-circle distance to + * the given point p. + */ + Geodetic2 closestPoint(const Geodetic2& p) const; + + + const Geodetic2& center() const; + const Geodetic2& halfSize() const; + Geodetic2 size() const; private: - Geodetic2 _center; - Geodetic2 _halfSize; + Geodetic2 _center; + Geodetic2 _halfSize; }; } // namespace openspace diff --git a/modules/globebrowsing/globes/chunknode.cpp b/modules/globebrowsing/globes/chunknode.cpp index 14507683d7..c0493a6253 100644 --- a/modules/globebrowsing/globes/chunknode.cpp +++ b/modules/globebrowsing/globes/chunknode.cpp @@ -125,6 +125,7 @@ void ChunkNode::internalRender(const RenderData& data, ChunkIndex& traverseData) LatLonPatchRenderer& patchRenderer = _owner.getPatchRenderer(); patchRenderer.renderPatch(_patch, data, _owner.ellipsoid(), ti); + //patchRenderer.renderPatch(_patch, data, _owner.ellipsoid()); ChunkNode::renderedPatches++; } } @@ -208,15 +209,17 @@ void ChunkNode::split(int depth) { Geodetic2 qs = Geodetic2(0.5 * hs.lat, 0.5 * hs.lon); // Subdivide bounds + GeodeticPatch nwBounds = GeodeticPatch(Geodetic2(c.lat + qs.lat, c.lon - qs.lon), qs); - GeodeticPatch neBounds = GeodeticPatch(Geodetic2(c.lat - qs.lat, c.lon - qs.lon), qs); - GeodeticPatch swBounds = GeodeticPatch(Geodetic2(c.lat + qs.lat, c.lon + qs.lon), qs); + GeodeticPatch neBounds = GeodeticPatch(Geodetic2(c.lat + qs.lat, c.lon + qs.lon), qs); + GeodeticPatch swBounds = GeodeticPatch(Geodetic2(c.lat - qs.lat, c.lon - qs.lon), qs); GeodeticPatch seBounds = GeodeticPatch(Geodetic2(c.lat - qs.lat, c.lon + qs.lon), qs); // Create new chunk nodes + _children[Quad::NORTH_WEST] = std::unique_ptr(new ChunkNode(_owner, nwBounds, this)); - _children[Quad::NORTH_EAST] = std::unique_ptr(new ChunkNode(_owner, neBounds, this)); _children[Quad::SOUTH_WEST] = std::unique_ptr(new ChunkNode(_owner, swBounds, this)); + _children[Quad::NORTH_EAST] = std::unique_ptr(new ChunkNode(_owner, neBounds, this)); _children[Quad::SOUTH_EAST] = std::unique_ptr(new ChunkNode(_owner, seBounds, this)); } diff --git a/modules/globebrowsing/globes/chunknode.h b/modules/globebrowsing/globes/chunknode.h index 5468fb367f..954ece8ad7 100644 --- a/modules/globebrowsing/globes/chunknode.h +++ b/modules/globebrowsing/globes/chunknode.h @@ -54,6 +54,9 @@ enum Quad { struct ChunkIndex { int x, y, level; + + // TODO : NOT SURE WHAT IS RIGHT BUT I SWITCHED THE TWO MIDDLE INDICES FOR IT TO WORK + std::vector childIndices() const { return { { 2 * x + 0, 2 * y + 0, level + 1 }, @@ -89,28 +92,28 @@ public: private: - void internalRender(const RenderData& data, ChunkIndex&); - bool internalUpdateChunkTree(const RenderData& data, ChunkIndex& traverseData); + void internalRender(const RenderData& data, ChunkIndex&); + bool internalUpdateChunkTree(const RenderData& data, ChunkIndex& traverseData); - /** - Uses horizon culling, frustum culling and distance to camera to determine a - desired level. - In the current implementation of the horizon culling and the distance to the - camera, the closer the ellipsoid is to a - sphere, the better this will make the splitting. Using the minimum radius to - be safe. This means that if the ellipsoid has high difference between radii, - splitting might accur even though it is not needed. - */ - int calculateDesiredLevelAndUpdateIsVisible( - const RenderData& data, - const ChunkIndex& traverseData); - - - ChunkNode* _parent; - std::unique_ptr _children[4]; + /** + Uses horizon culling, frustum culling and distance to camera to determine a + desired level. + In the current implementation of the horizon culling and the distance to the + camera, the closer the ellipsoid is to a + sphere, the better this will make the splitting. Using the minimum radius to + be safe. This means that if the ellipsoid has high difference between radii, + splitting might accur even though it is not needed. + */ + int calculateDesiredLevelAndUpdateIsVisible( + const RenderData& data, + const ChunkIndex& traverseData); + + + ChunkNode* _parent; + std::unique_ptr _children[4]; ChunkLodGlobe& _owner; GeodeticPatch _patch; - bool _isVisible; + bool _isVisible; }; } // namespace openspace diff --git a/modules/globebrowsing/other/gdaldataconverter.cpp b/modules/globebrowsing/other/gdaldataconverter.cpp index 43fa207b0c..e90cd4a030 100644 --- a/modules/globebrowsing/other/gdaldataconverter.cpp +++ b/modules/globebrowsing/other/gdaldataconverter.cpp @@ -24,116 +24,159 @@ #include +#include +#include namespace { - const std::string _loggerCat = "GdalDataConverter"; + const std::string _loggerCat = "GdalDataConverter"; } namespace openspace { - GdalDataConverter::GdalDataConverter() - { + GdalDataConverter::GdalDataConverter() + { - } + } - GdalDataConverter::~GdalDataConverter() - { + GdalDataConverter::~GdalDataConverter() + { - } + } - std::shared_ptr GdalDataConverter::convertToOpenGLTexture( - GDALDataset* dataSet, - const GeodeticTileIndex& tileIndex, - int GLType) - { - int nRasters = dataSet->GetRasterCount(); + std::shared_ptr GdalDataConverter::convertToOpenGLTexture( + GDALDataset* dataSet, + const GeodeticTileIndex& tileIndex, + int GLType) + { + int nRasters = dataSet->GetRasterCount(); - ghoul_assert(nRasters > 0, "Bad dataset. Contains no rasterband."); + ghoul_assert(nRasters > 0, "Bad dataset. Contains no rasterband."); - GDALRasterBand* firstBand = dataSet->GetRasterBand(1); + GDALRasterBand* firstBand = dataSet->GetRasterBand(1); - // Level = overviewCount - overview - int overviewCount = firstBand->GetOverviewCount(); - int overview = overviewCount - tileIndex.level - 1; + // Level = overviewCount - overview + int overviewCount = firstBand->GetOverviewCount(); + int overview = overviewCount - tileIndex.level - 1; - // The output texture will have this size - int xSizelevel0 = firstBand->GetOverview(overviewCount - 1)->GetXSize(); - int ySizelevel0 = firstBand->GetOverview(overviewCount - 1)->GetYSize(); + // The output texture will have this size + int xSizelevel0 = firstBand->GetOverview(overviewCount - 1)->GetXSize(); + int ySizelevel0 = firstBand->GetOverview(overviewCount - 1)->GetYSize(); - // Create all the raster bands (Commonly one for each channel: red, green, blue) - std::vector rasterBands; - for (size_t i = 0; i < nRasters; i++) - { - rasterBands.push_back(dataSet->GetRasterBand(i + 1)->GetOverview(overview)); - } - // The data that the texture should read - GLubyte* imageData = new GLubyte[xSizelevel0 * ySizelevel0 * nRasters]; - // Read the data - for (size_t i = 0; i < nRasters; i++) - { - int xBeginRead = tileIndex.x * pow(2, tileIndex.level) * xSizelevel0; - int yBeginRead = tileIndex.y * pow(2, tileIndex.level) * ySizelevel0; - rasterBands[i]->RasterIO( - GF_Read, - xBeginRead, // Begin read x - yBeginRead, // Begin read y - xSizelevel0, // width to read x - ySizelevel0, // width to read y - imageData + i, // Where to put data - xSizelevel0, // width to read x in destination - ySizelevel0, // width to read y in destination - GDT_Byte, // Type - sizeof(GLubyte) * nRasters, // Pixel spacing - 0); // Line spacing - } - GdalDataConverter::TextureFormat textrureFormat = - getTextureFormatFromRasterCount(nRasters); + // The data that the texture should read + GLubyte* imageData = new GLubyte[xSizelevel0 * ySizelevel0 * nRasters]; - // The texture should take ownership of the data - std::shared_ptr texture = std::shared_ptr(new Texture( - static_cast(imageData), - glm::uvec3(xSizelevel0, ySizelevel0, 1), - textrureFormat.ghoulFormat, - textrureFormat.glFormat, - GL_UNSIGNED_BYTE, - Texture::FilterMode::Linear, - Texture::WrappingMode::Repeat)); + // Read the data (each rasterband is a separate channel) + for (size_t i = 0; i < nRasters; i++) + { + GDALRasterBand* rasterBand = dataSet->GetRasterBand(i + 1)->GetOverview(overview); + + int xBeginRead = tileIndex.x * pow(2, tileIndex.level) * xSizelevel0; + int yBeginRead = tileIndex.y * pow(2, tileIndex.level) * ySizelevel0; + rasterBand->RasterIO( + GF_Read, + xBeginRead, // Begin read x + yBeginRead, // Begin read y + xSizelevel0, // width to read x + ySizelevel0, // width to read y + imageData + i, // Where to put data + xSizelevel0, // width to read x in destination + ySizelevel0, // width to read y in destination + GDT_Byte, // Type + sizeof(GLubyte) * nRasters, // Pixel spacing + 0); // Line spacing + } - // Do not free imageData since the texture now has ownership of it - return texture; - } + GdalDataConverter::TextureFormat textrureFormat = + getTextureFormatFromRasterCount(nRasters); - GdalDataConverter::TextureFormat GdalDataConverter::getTextureFormatFromRasterCount( - int rasterCount) - { - TextureFormat format; - switch (rasterCount) - { - case 1: - format.ghoulFormat = Texture::Format::Red; - format.glFormat = GL_RED; - break; - case 2: - format.ghoulFormat = Texture::Format::RG; - format.glFormat = GL_RG; - break; - case 3: - format.ghoulFormat = Texture::Format::RGB; - format.glFormat = GL_RGB; - break; - case 4: - format.ghoulFormat = Texture::Format::RGBA; - format.glFormat = GL_RGBA; - break; - default: + Texture* tex = new Texture( + static_cast(imageData), + glm::uvec3(xSizelevel0, ySizelevel0, 1), + textrureFormat.ghoulFormat, + textrureFormat.glFormat, + GL_UNSIGNED_BYTE, + Texture::FilterMode::Linear, + Texture::WrappingMode::Repeat); - break; - } - return format; - } + // The texture should take ownership of the data + std::shared_ptr texture = std::shared_ptr(tex); + + + // Do not free imageData since the texture now has ownership of it + return texture; + } + + + + glm::uvec2 GdalDataConverter::geodeticToPixel(GDALDataset* dataSet, const Geodetic2& geo) { + double padfTransform[6]; + CPLErr err = dataSet->GetGeoTransform(padfTransform); + + ghoul_assert(err != CE_Failure, "Failed to get transform"); + + Scalar Y = Angle::fromRadians(geo.lat).asDegrees(); + Scalar X = Angle::fromRadians(geo.lon).asDegrees(); + + // convert from pixel and line to geodetic coordinates + // Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2]; + // Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5]; + + // <=> + double* a = &(padfTransform[0]); + double* b = &(padfTransform[3]); + + // Xp = a[0] + P*a[1] + L*a[2]; + // Yp = b[0] + P*b[1] + L*b[2]; + + // <=> + double divisor = (a[2]*b[1] - a[1]*b[2]); + ghoul_assert(divisor != 0.0, "Division by zero!"); + //ghoul_assert(a[2] != 0.0, "a2 must not be zero!"); + double P = (a[0]*b[2] - a[2]*b[0] + a[2]*Y - b[2]*X) / divisor; + double L = (-a[0]*b[1] + a[1]*b[0] - a[1]*Y + b[1]*X) / divisor; + // ref: https://www.wolframalpha.com/input/?i=X+%3D+a0+%2B+a1P+%2B+a2L,+Y+%3D+b0+%2B+b1P+%2B+b2L,+solve+for+P+and+L + + + double Xp = a[0] + P*a[1] + L*a[2]; + double Yp = b[0] + P*b[1] + L*b[2]; + + return glm::uvec2(P, L); + } + + + + GdalDataConverter::TextureFormat GdalDataConverter::getTextureFormatFromRasterCount( + int rasterCount) + { + TextureFormat format; + + switch (rasterCount) + { + case 1: + format.ghoulFormat = Texture::Format::Red; + format.glFormat = GL_RED; + break; + case 2: + format.ghoulFormat = Texture::Format::RG; + format.glFormat = GL_RG; + break; + case 3: + format.ghoulFormat = Texture::Format::RGB; + format.glFormat = GL_RGB; + break; + case 4: + format.ghoulFormat = Texture::Format::RGBA; + format.glFormat = GL_RGBA; + break; + default: + + break; + } + return format; + } } // namespace openspace diff --git a/modules/globebrowsing/other/gdaldataconverter.h b/modules/globebrowsing/other/gdaldataconverter.h index f27faacac7..4eaaa162e4 100644 --- a/modules/globebrowsing/other/gdaldataconverter.h +++ b/modules/globebrowsing/other/gdaldataconverter.h @@ -25,38 +25,46 @@ #ifndef __GDALDATACONVERTER_H__ #define __GDALDATACONVERTER_H__ -#include +//#include #include #include +#include + #include "gdal_priv.h" #include namespace openspace { - using namespace ghoul::opengl; + using namespace ghoul::opengl; - class GdalDataConverter - { - public: - GdalDataConverter(); - ~GdalDataConverter(); + // forward declaration + class GeodeticTileIndex; + //class Geodetic2; - std::shared_ptr convertToOpenGLTexture( - GDALDataset* dataSet, - const GeodeticTileIndex& tileIndex, - int GLType); + class GdalDataConverter + { + public: + GdalDataConverter(); + ~GdalDataConverter(); - private: - struct TextureFormat - { - Texture::Format ghoulFormat; - GLuint glFormat; - }; - TextureFormat getTextureFormatFromRasterCount(int rasterCount); - }; + std::shared_ptr convertToOpenGLTexture( + GDALDataset* dataSet, + const GeodeticTileIndex& tileIndex, + int GLType); + + struct TextureFormat + { + Texture::Format ghoulFormat; + GLuint glFormat; + }; + + TextureFormat getTextureFormatFromRasterCount(int rasterCount); + + glm::uvec2 geodeticToPixel(GDALDataset* dataSet, const Geodetic2& geo); + }; } // namespace openspace diff --git a/modules/globebrowsing/other/texturetileset.cpp b/modules/globebrowsing/other/texturetileset.cpp index 6080ada64f..7b70094f45 100644 --- a/modules/globebrowsing/other/texturetileset.cpp +++ b/modules/globebrowsing/other/texturetileset.cpp @@ -39,144 +39,146 @@ #include namespace { - const std::string _loggerCat = "TextureTileSet"; + const std::string _loggerCat = "TextureTileSet"; } namespace openspace { - TextureTileSet::TextureTileSet( - Geodetic2 sizeLevel0, - Geodetic2 offsetLevel0, - int depth) - : _sizeLevel0(sizeLevel0) - , _offsetLevel0(offsetLevel0) - , _depth(depth) - { - + TextureTileSet::TextureTileSet( + Geodetic2 sizeLevel0, + Geodetic2 offsetLevel0, + int depth) + : _sizeLevel0(sizeLevel0) + , _offsetLevel0(offsetLevel0) + , _depth(depth) + { + - // Read using GDAL + // Read using GDAL - std::string testFile = absPath("map_service_configs/TERRA_CR_B143_2016-04-12.wms"); + //std::string testFile = absPath("map_service_configs/TERRA_CR_B143_2016-04-12.wms"); + //std::string testFile = absPath("map_service_configs/TERRAIN.wms"); + std::string testFile = absPath("textures/earth_bluemarble.jpg"); - - GDALDataset *poDataset; - GDALAllRegister(); - poDataset = (GDALDataset *)GDALOpen(testFile.c_str(), GA_ReadOnly); - assert(poDataset != NULL); - GdalDataConverter conv; + + GDALDataset *poDataset; + GDALAllRegister(); + poDataset = (GDALDataset *)GDALOpen(testFile.c_str(), GA_ReadOnly); + assert(poDataset != nullptr); + GdalDataConverter conv; - GeodeticTileIndex ti; - ti.x = 0; - ti.y = 0; - ti.level = 0; - _testTexture = conv.convertToOpenGLTexture(poDataset, ti, GL_UNSIGNED_BYTE); + GeodeticTileIndex ti; + ti.x = 0; + ti.y = 0; + ti.level = 0; + _testTexture = conv.convertToOpenGLTexture(poDataset, ti, GL_UNSIGNED_BYTE); - _testTexture->uploadTexture(); - _testTexture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); - - /* - - // Set e texture to test - std::string fileName = "textures/earth_bluemarble.jpg"; - //std::string fileName = "../../../build/tiles/tile5_8_12.png"; - //std::string fileName = "tile5_8_12.png"; - _testTexture = std::move(ghoul::io::TextureReader::ref().loadTexture(absPath(fileName))); - - if (_testTexture) { - LDEBUG("Loaded texture from '" << "textures/earth_bluemarble.jpg" << "'"); - _testTexture->uploadTexture(); + _testTexture->uploadTexture(); + _testTexture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); + + /* + + // Set e texture to test + std::string fileName = "textures/earth_bluemarble.jpg"; + //std::string fileName = "../../../build/tiles/tile5_8_12.png"; + //std::string fileName = "tile5_8_12.png"; + _testTexture = std::move(ghoul::io::TextureReader::ref().loadTexture(absPath(fileName))); + + if (_testTexture) { + LDEBUG("Loaded texture from '" << "textures/earth_bluemarble.jpg" << "'"); + _testTexture->uploadTexture(); - // Textures of planets looks much smoother with AnisotropicMipMap rather than linear - // TODO: AnisotropicMipMap crashes on ATI cards ---abock - //_testTexture->setFilter(ghoul::opengl::Texture::FilterMode::AnisotropicMipMap); - _testTexture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); - } - */ - } + // Textures of planets looks much smoother with AnisotropicMipMap rather than linear + // TODO: AnisotropicMipMap crashes on ATI cards ---abock + //_testTexture->setFilter(ghoul::opengl::Texture::FilterMode::AnisotropicMipMap); + _testTexture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); + } + */ + } - TextureTileSet::~TextureTileSet(){ + TextureTileSet::~TextureTileSet(){ - } + } - GeodeticTileIndex TextureTileSet::getTileIndex(const GeodeticPatch& patch) { - // Calculate the level of the index depanding on the size of the incoming patch. - // The level is as big as possible (as far down as possible) but it can't be - // too big since at maximum four tiles should be used to cover a patch - int level = log2(static_cast(glm::max( - _sizeLevel0.lat / (patch.size().lat), - _sizeLevel0.lon / (patch.size().lon)))); - - // If the depth is not big enough, the level must be clamped. - level = glm::min(level, _depth); - - // Calculate the index in x y where the tile should be positioned - Vec2 tileSize = _sizeLevel0.toLonLatVec2() / pow(2, level); - Vec2 nw = patch.northWestCorner().toLonLatVec2(); - Vec2 offset = _offsetLevel0.toLonLatVec2(); - glm::ivec2 tileIndexXY = (nw - offset) / tileSize; + GeodeticTileIndex TextureTileSet::getTileIndex(const GeodeticPatch& patch) { + // Calculate the level of the index depanding on the size of the incoming patch. + // The level is as big as possible (as far down as possible) but it can't be + // too big since at maximum four tiles should be used to cover a patch + int level = log2(static_cast(glm::max( + _sizeLevel0.lat / (patch.size().lat), + _sizeLevel0.lon / (patch.size().lon)))); + + // If the depth is not big enough, the level must be clamped. + level = glm::min(level, _depth); + + // Calculate the index in x y where the tile should be positioned + Vec2 tileSize = _sizeLevel0.toLonLatVec2() / pow(2, level); + Vec2 nw = patch.northWestCorner().toLonLatVec2(); + Vec2 offset = _offsetLevel0.toLonLatVec2(); + glm::ivec2 tileIndexXY = (nw - offset) / tileSize; - // Flip y since indices increase from top to bottom - tileIndexXY.y *= -1; + // Flip y since indices increase from top to bottom + tileIndexXY.y *= -1; - // Create the tileindex - GeodeticTileIndex tileIndex = { tileIndexXY.x, tileIndexXY.y, level }; - return tileIndex; - } + // Create the tileindex + GeodeticTileIndex tileIndex = { tileIndexXY.x, tileIndexXY.y, level }; + return tileIndex; + } - std::shared_ptr TextureTileSet::getTile(const GeodeticPatch& patch) { - return getTile(getTileIndex(patch)); - } + std::shared_ptr TextureTileSet::getTile(const GeodeticPatch& patch) { + return getTile(getTileIndex(patch)); + } - std::shared_ptr TextureTileSet::getTile(const GeodeticTileIndex& tileIndex) { - return _testTexture; - } + std::shared_ptr TextureTileSet::getTile(const GeodeticTileIndex& tileIndex) { + return _testTexture; + } - GeodeticPatch TextureTileSet::getTilePositionAndScale( - const GeodeticTileIndex& tileIndex) { - Geodetic2 tileSize = Geodetic2( - _sizeLevel0.lat / pow(2, tileIndex.level), - _sizeLevel0.lon / pow(2, tileIndex.level)); - Geodetic2 northWest = Geodetic2( - _offsetLevel0.lat + tileIndex.y * tileSize.lat, - _offsetLevel0.lon + tileIndex.x * tileSize.lon); - - return GeodeticPatch( - Geodetic2(northWest.lat - tileSize.lat / 2, northWest.lon + tileSize.lon / 2), - Geodetic2(tileSize.lat / 2, tileSize.lon / 2)); - } + GeodeticPatch TextureTileSet::getTilePositionAndScale( + const GeodeticTileIndex& tileIndex) { + Geodetic2 tileSize = Geodetic2( + _sizeLevel0.lat / pow(2, tileIndex.level), + _sizeLevel0.lon / pow(2, tileIndex.level)); + Geodetic2 northWest = Geodetic2( + _offsetLevel0.lat + tileIndex.y * tileSize.lat, + _offsetLevel0.lon + tileIndex.x * tileSize.lon); + + return GeodeticPatch( + Geodetic2(northWest.lat - tileSize.lat / 2, northWest.lon + tileSize.lon / 2), + Geodetic2(tileSize.lat / 2, tileSize.lon / 2)); + } - glm::mat3 TextureTileSet::getUvTransformationPatchToTile( - GeodeticPatch patch, - const GeodeticTileIndex& tileIndex) - { - GeodeticPatch tile = getTilePositionAndScale(tileIndex); - return getUvTransformationPatchToTile(patch, tile); - } + glm::mat3 TextureTileSet::getUvTransformationPatchToTile( + GeodeticPatch patch, + const GeodeticTileIndex& tileIndex) + { + GeodeticPatch tile = getTilePositionAndScale(tileIndex); + return getUvTransformationPatchToTile(patch, tile); + } - glm::mat3 TextureTileSet::getUvTransformationPatchToTile( - GeodeticPatch patch, - GeodeticPatch tile) - { - Vec2 posDiff = - patch.southWestCorner().toLonLatVec2() - - tile.southWestCorner().toLonLatVec2(); + glm::mat3 TextureTileSet::getUvTransformationPatchToTile( + GeodeticPatch patch, + GeodeticPatch tile) + { + Vec2 posDiff = + patch.southWestCorner().toLonLatVec2() - + tile.southWestCorner().toLonLatVec2(); - glm::mat3 invTileScale = glm::mat3( - { 1 / (tile.halfSize().lon * 2), 0, 0, - 0, 1 / (tile.halfSize().lat * 2), 0, - 0, 0, 1 }); + glm::mat3 invTileScale = glm::mat3( + { 1 / (tile.size().lon), 0, 0, + 0, 1 / (tile.size().lat), 0, + 0, 0, 1 }); - glm::mat3 globalTranslation = glm::mat3( - { 1, 0, 0, - 0, 1, 0, - posDiff.x, posDiff.y, 1 }); + glm::mat3 globalTranslation = glm::mat3( + { 1, 0, 0, + 0, 1, 0, + posDiff.x, posDiff.y, 1 }); - glm::mat3 patchScale = glm::mat3( - { (patch.halfSize().lon * 2), 0, 0, - 0, (patch.halfSize().lat * 2), 0, - 0, 0, 1 }); + glm::mat3 patchScale = glm::mat3( + { (patch.halfSize().lon * 2), 0, 0, + 0, (patch.halfSize().lat * 2), 0, + 0, 0, 1 }); - return invTileScale * globalTranslation * patchScale; - } + return invTileScale * globalTranslation * patchScale; + } } // namespace openspace diff --git a/modules/globebrowsing/other/tileprovider.cpp b/modules/globebrowsing/other/tileprovider.cpp index 48d84186f3..151e5577aa 100644 --- a/modules/globebrowsing/other/tileprovider.cpp +++ b/modules/globebrowsing/other/tileprovider.cpp @@ -22,7 +22,7 @@ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * ****************************************************************************************/ - +#include #include #include @@ -37,69 +37,43 @@ namespace { - const std::string _loggerCat = "TwmsTileProvider"; + const std::string _loggerCat = "TileProvider"; } + namespace openspace { + bool TileProvider::hasInitializedGDAL = false; - TileProvider::TileProvider(int tileCacheSize) - : _tileCache(tileCacheSize) // setting cache size + TileProvider::TileProvider(const std::string& filePath, int tileCacheSize) + : _filePath(filePath) + , _tileCache(tileCacheSize) // setting cache size { int downloadApplicationVersion = 1; if (!DownloadManager::isInitialized()) { DownloadManager::initialize("../tmp_openspace_downloads/", downloadApplicationVersion); } + if (!hasInitializedGDAL) { + GDALAllRegister(); + hasInitializedGDAL = true; + } + + std::string absFilePath = absPath(filePath); + _gdalDataSet = (GDALDataset *)GDALOpen(absFilePath.c_str(), GA_ReadOnly); + //auto desc = _gdalDataSet->GetDescription(); + + ghoul_assert(_gdalDataSet != nullptr, "Failed to load dataset: " << filePath); + } TileProvider::~TileProvider(){ - + delete _gdalDataSet; } void TileProvider::prerender() { - // Remove filefutures that are inactive - auto it = _fileFutureMap.begin(); - auto end = _fileFutureMap.end(); - for (; it != end;){ - - if (it->second->isAborted || - it->second->secondsRemaining > 20 || - it->second->errorMessage.compare("") != 0 || - it->second->isFinished) - { - it = _fileFutureMap.erase(it); - //LDEBUG("removnig filefuture"); - } - else { - it++; - } - } - - - // Move finished texture tiles to cache - /* - while (_concurrentJobManager.numFinishedJobs() > 0) { - auto finishedJob = _concurrentJobManager.popFinishedJob(); - - TextureLoadJob* finishedTextureLoadJob = reinterpret_cast(finishedJob.get()); - - ghoul_assert(finishedTextureLoadJob != nullptr, "unable to reinterpret cast to TextureLoadJob*"); - - auto texture = finishedTextureLoadJob->product(); - // upload to gpu - texture->uploadTexture(); - texture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); - texture->setWrapping(ghoul::opengl::Texture::WrappingMode::ClampToEdge); - - //LDEBUG("uploaded texture " << finishedTextureLoadJob->hashKey()); - HashKey hashkey = finishedTextureLoadJob->hashKey(); - _tileCache.put(hashkey, texture); - _fileFutureMap.erase(hashkey); - } - */ } @@ -109,77 +83,112 @@ namespace openspace { if (_tileCache.exist(hashkey)) { return _tileCache.get(hashkey); } - else if (_fileFutureMap.find(hashkey) != _fileFutureMap.end()) { - if (_fileFutureMap.find(hashkey)->second->isFinished) { - - std::string fileName = _fileFutureMap.find(hashkey)->second->filePath; - std::string filePath = "tiles/" + fileName; - - std::shared_ptr texture = loadAndInitTextureDisk(filePath); - _tileCache.put(hashkey, texture); - _fileFutureMap.erase(hashkey); + else { + //GeodeticTileIndex ti0 = { 0, 0, 0 }; + //auto texture = _converter.convertToOpenGLTexture(_gdalDataSet, tileIndex, GL_UNSIGNED_BYTE); + auto texture = getTileInternal(tileIndex, GL_UNSIGNED_BYTE); + texture->uploadTexture(); + texture->setWrapping(ghoul::opengl::Texture::WrappingMode::ClampToEdge); + texture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); + _tileCache.put(hashkey, texture); + return texture; + } + } - //LDEBUG("Downloaded " << fileName); - - //auto job = std::shared_ptr(new TextureLoadJob(filePath, hashkey)); - //_concurrentJobManager.enqueueJob(job); + + + + std::shared_ptr TileProvider::getTileInternal(const GeodeticTileIndex& tileIndex, int GLType) { + + + + int nRasters = _gdalDataSet->GetRasterCount(); + + ghoul_assert(nRasters > 0, "Bad dataset. Contains no rasterband."); + + GDALRasterBand* firstBand = _gdalDataSet->GetRasterBand(1); + + // Level = overviewCount - overview + int numOverviews = firstBand->GetOverviewCount(); + + + int xSize0 = firstBand->GetXSize(); + int ySize0 = firstBand->GetYSize(); + + GeodeticPatch patch = GeodeticPatch(tileIndex); + glm::uvec2 pixelStart0 = _converter.geodeticToPixel(_gdalDataSet, patch.northWestCorner()); + glm::uvec2 pixelEnd0 = _converter.geodeticToPixel(_gdalDataSet, patch.southEastCorner()); + glm::uvec2 numberOfPixels0 = pixelEnd0 - pixelStart0; + + + int minNumPixels0 = glm::min(numberOfPixels0.x, numberOfPixels0.y); + + int minNumPixelsRequired = 256; + int ov = log2(minNumPixels0) - log2(minNumPixelsRequired); + + ov = glm::clamp(ov, 0, numOverviews-1); + + glm::uvec2 pixelStart(pixelStart0.x >> (ov+1), pixelStart0.y >> (ov + 1)); + glm::uvec2 numberOfPixels(numberOfPixels0.x >> (ov + 1), numberOfPixels0.y >> (ov + 1)); + + // For testing + /*pixelStart = glm::uvec2(0, 0); + numberOfPixels = glm::uvec2(512, 256); + ov = 15; + */ + // The data that the texture should read + GLubyte* imageData = new GLubyte[numberOfPixels.x * numberOfPixels.y * nRasters]; + + // Read the data (each rasterband is a separate channel) + for (size_t i = 0; i < nRasters; i++) { + GDALRasterBand* rasterBand = _gdalDataSet->GetRasterBand(i + 1)->GetOverview(ov); + + int xSize = rasterBand->GetXSize(); + int ySize = rasterBand->GetYSize(); + + rasterBand->RasterIO( + GF_Read, + pixelStart.x, // Begin read x + pixelStart.y, // Begin read y + numberOfPixels.x, // width to read x + numberOfPixels.y, // width to read y + imageData + i, // Where to put data + numberOfPixels.x, // width to write x in destination + numberOfPixels.y, // width to write y in destination + GDT_Byte, // Type + sizeof(GLubyte) * nRasters, // Pixel spacing + 0); // Line spacing + } + + GLubyte* imageDataYflipped = new GLubyte[numberOfPixels.x * numberOfPixels.y * nRasters]; + for (size_t y = 0; y < numberOfPixels.y; y++) { + for (size_t x = 0; x < numberOfPixels.x; x++) { + imageDataYflipped[x + y * numberOfPixels.x] = imageData[x + (numberOfPixels.y - y) * numberOfPixels.x]; } } - else if(_fileFutureMap.size() < 100){ - std::shared_ptr fileFuture = requestTile(tileIndex); - _fileFutureMap.insert_or_assign(hashkey, fileFuture); - } - return nullptr; - } - std::shared_ptr TileProvider::loadAndInitTextureDisk(std::string filePath) { - auto textureReader = ghoul::io::TextureReader::ref(); - std::shared_ptr texture = std::move(textureReader.loadTexture(absPath(filePath))); + GdalDataConverter::TextureFormat textrureFormat = _converter.getTextureFormatFromRasterCount(nRasters); - // upload to gpu - texture->uploadTexture(); - texture->setFilter(ghoul::opengl::Texture::FilterMode::Linear); - texture->setWrapping(ghoul::opengl::Texture::WrappingMode::ClampToEdge); + + Texture* tex = new Texture( + static_cast(imageDataYflipped), + glm::uvec3(numberOfPixels.x, numberOfPixels.y, 1), + textrureFormat.ghoulFormat, + textrureFormat.glFormat, + GL_UNSIGNED_BYTE, + Texture::FilterMode::Linear, + Texture::WrappingMode::Repeat); + + // The texture should take ownership of the data + std::shared_ptr texture = std::shared_ptr(tex); + + delete[] imageData; + + // Do not free imageData since the texture now has ownership of it return texture; - } - - - - std::shared_ptr TileProvider::requestTile(const GeodeticTileIndex& tileIndex) { - // download tile - std::stringstream ss; - //std::string baseUrl = "https://map1c.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi?TIME=2016-04-17&layer=MODIS_Terra_CorrectedReflectance_TrueColor&tilematrixset=EPSG4326_250m&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image%2Fjpeg"; - //ss << baseUrl; - //ss << "&TileMatrix=" << tileIndex.level; - //ss << "&TileCol=" << tileIndex.x; - //ss << "&TileRow=" << tileIndex.y; - // https://map1c.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi?TIME=2016-04-17&layer=MODIS_Terra_CorrectedReflectance_TrueColor&tilematrixset=EPSG4326_250m&Service=WMTS&Request=GetTile&Version=1.0.0&Format=image%2Fjpeg&TileMatrix=0&TileCol=0&TileRow=0 - - std::string baseUrl = "http://mesonet.agron.iastate.edu/cache/tile.py/1.0.0/nexrad-n0q-900913"; - ss << baseUrl; - ss << "/" << tileIndex.level; - ss << "/" << tileIndex.x; - ss << "/" << tileIndex.y; - ss << ".png?1461277159335"; - // http://mesonet.agron.iastate.edu/cache/tile.py/1.0.0/nexrad-n0q-900913/5/8/12.png? - std::string twmsRequestUrl = ss.str(); - - ss = std::stringstream(); - ss << tileIndex.level; - ss << "_" << tileIndex.x; - ss << "_" << tileIndex.y; - std::string filePath = "tiles/tile" + ss.str() + ".png"; - - using ghoul::filesystem::File; - File localTileFile(filePath); - bool overrideFile = true; - - std::shared_ptr fileFuture = DownloadManager::ref().downloadFile(twmsRequestUrl, localTileFile, overrideFile); - return fileFuture; } - } // namespace openspace diff --git a/modules/globebrowsing/other/tileprovider.h b/modules/globebrowsing/other/tileprovider.h index de843344b6..c1bcc791d7 100644 --- a/modules/globebrowsing/other/tileprovider.h +++ b/modules/globebrowsing/other/tileprovider.h @@ -22,39 +22,30 @@ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * ****************************************************************************************/ -#ifndef __TWMS_TILE_PROVIDER_H__ -#define __TWMS_TILE_PROVIDER_H__ +#ifndef __TILE_PROVIDER_H__ +#define __TILE_PROVIDER_H__ -#include -#include +#include "gdal_priv.h" +#include #include #include // absPath #include #include -#include +#include +#include +#include +#include -////////////////////////////////////////////////////////////////////////////////////// -// TILE INDEX // -////////////////////////////////////////////////////////////////////////////////////// - -namespace openspace { - using HashKey = unsigned long; - - struct GeodeticTileIndex { - int x, y, level; - - HashKey hashKey() const { - return x ^ (y << 16) ^ (level << 21); - } - }; -} + + +/* namespace openspace { using namespace ghoul::opengl; @@ -87,13 +78,11 @@ namespace openspace { HashKey _hashkey; std::shared_ptr _texture; }; - - } - +*/ ////////////////////////////////////////////////////////////////////////////////////////// -// TWMS TILE PROVIDER // +// TILE PROVIDER // ////////////////////////////////////////////////////////////////////////////////////////// @@ -102,7 +91,7 @@ namespace openspace { class TileProvider { public: - TileProvider(int tileCacheSize); + TileProvider(const std::string& fileName, int tileCacheSize); ~TileProvider(); std::shared_ptr getTile(const GeodeticTileIndex& tileIndex); @@ -111,17 +100,23 @@ namespace openspace { private: - std::shared_ptr requestTile(const GeodeticTileIndex&); - std::shared_ptr loadAndInitTextureDisk(std::string filePath); + + std::shared_ptr getTileInternal(const GeodeticTileIndex& tileIndex, int GLType); + LRUCache> _tileCache; - std::unordered_map> _fileFutureMap; - //LRUCache> _fileFutureCache; - //ConcurrentJobManager _concurrentJobManager; + const std::string _filePath; + + static bool hasInitializedGDAL; + GDALDataset* _gdalDataSet; + GdalDataConverter _converter; + + + }; } // namespace openspace -#endif // __TWMS_TILE_PROVIDER_H__ \ No newline at end of file +#endif // __TILE_PROVIDER_H__ \ No newline at end of file diff --git a/modules/globebrowsing/rendering/patchrenderer.cpp b/modules/globebrowsing/rendering/patchrenderer.cpp index 17827fe080..074fad5d14 100644 --- a/modules/globebrowsing/rendering/patchrenderer.cpp +++ b/modules/globebrowsing/rendering/patchrenderer.cpp @@ -77,7 +77,9 @@ namespace openspace { LatLonPatchRenderer::LatLonPatchRenderer(shared_ptr grid) : PatchRenderer() , _grid(grid) - , tileProvider(5000) + , tileProvider("map_service_configs/TERRAIN.wms" , 5000) + //, tileProvider("map_service_configs/frmt_wms_virtualearth.xml", 5000) + //, tileProvider("textures/earth_bluemarble.jpg", 5000) { _programObject = OsEng.renderEngine().buildRenderProgram( "LatLonSphereMappingProgram", @@ -99,7 +101,6 @@ namespace openspace { // Get the textures that should be used for rendering GeodeticTileIndex ti = _tileSet.getTileIndex(patch); - renderPatch(patch, data, ellipsoid, ti); } @@ -112,8 +113,11 @@ namespace openspace { using namespace glm; + // PATCH DID NOT MATCH THE TILEINDEX SO I CREATED A NEW PATCH FROM THE INDEX + GeodeticPatch newPatch = patch;//GeodeticPatch(tileIndex); + //GeodeticTileIndex tmpTileIndex = _tileSet.getTileIndex(patch); // TODO : Model transform should be fetched as a matrix directly. @@ -129,18 +133,14 @@ namespace openspace { // Get the textures that should be used for rendering std::shared_ptr tile00; bool usingTile = true; - GeodeticTileIndex ti; - ti.level =tileIndex.level; - ti.x = tileIndex.y; - ti.y = tileIndex.x; - tile00 = tileProvider.getTile(ti); + tile00 = tileProvider.getTile(tileIndex); if (tile00 == nullptr) { tile00 = _tileSet.getTile(tileIndex); usingTile = false; } - glm::mat3 uvTransform = usingTile ? glm::mat3(1) : _tileSet.getUvTransformationPatchToTile(patch, tileIndex); + glm::mat3 uvTransform = usingTile ? glm::mat3(1) : _tileSet.getUvTransformationPatchToTile(newPatch, tileIndex); // Bind and use the texture ghoul::opengl::TextureUnit texUnit; @@ -149,11 +149,11 @@ namespace openspace { _programObject->setUniform("textureSampler", texUnit); _programObject->setUniform("uvTransformPatchToTile", uvTransform); - Geodetic2 swCorner = patch.southWestCorner(); + Geodetic2 swCorner = newPatch.southWestCorner(); _programObject->setUniform("segmentsPerPatch", _grid->xSegments()); _programObject->setUniform("modelViewProjectionTransform", modelViewProjectionTransform); _programObject->setUniform("minLatLon", vec2(swCorner.toLonLatVec2())); - _programObject->setUniform("lonLatScalingFactor", vec2(patch.size().toLonLatVec2())); + _programObject->setUniform("lonLatScalingFactor", vec2(newPatch.size().toLonLatVec2())); _programObject->setUniform("radiiSquared", vec3(ellipsoid.radiiSquared())); glEnable(GL_DEPTH_TEST); @@ -186,6 +186,12 @@ namespace openspace { _programObject->setIgnoreSubroutineUniformLocationError(IgnoreError::Yes); } + void ClipMapPatchRenderer::update() { + + } + + + void ClipMapPatchRenderer::renderPatch( const Geodetic2& patchSize, const RenderData& data, diff --git a/modules/globebrowsing/rendering/patchrenderer.h b/modules/globebrowsing/rendering/patchrenderer.h index e9f9a121f8..8a433ead05 100644 --- a/modules/globebrowsing/rendering/patchrenderer.h +++ b/modules/globebrowsing/rendering/patchrenderer.h @@ -60,6 +60,7 @@ namespace openspace { PatchRenderer(); ~PatchRenderer(); + protected: unique_ptr _programObject; @@ -102,10 +103,14 @@ namespace openspace { public: ClipMapPatchRenderer(shared_ptr grid); + void update(); + void renderPatch( const Geodetic2& patchSize, const RenderData& data, const Ellipsoid& ellipsoid); + + private: shared_ptr _grid; };