mirror/OpenSpace
1
0
Fork 0
mirror of https://github.com/OpenSpace/OpenSpace.git synced 2026-01-05 11:09:37 -06:00
Code Issues Packages Projects Releases Wiki Activity

Carving out cachable geodetic tiles from GDAL datasets

Browse Source
This commit is contained in:
Erik Broberg
2016-05-02 11:13:41 -04:00
parent 4656d997e4
commit 3a3038e862
12 changed files with 829 additions and 686 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
data/scene/debugglobe/map_service_configs/TERRAIN.wms
View File

@@ -3,4 +3,11 @@
<ServerUrl>http://198.102.45.23/arcgis/rest/services/worldelevation3d/terrain3d?</ServerUrl>
<TiledGroupName>GCS_Elevation</TiledGroupName>
</Service>
<DataWindow>
<UpperLeftX>-180.0</UpperLeftX>
<UpperLeftY>90.0</UpperLeftY>
<LowerRightX>180.0</LowerRightX>
<LowerRightY>-90.0</LowerRightY>
<YOrigin>bottom</YOrigin>
</DataWindow>
</GDAL_WMS>

476
modules/globebrowsing/geodetics/geodetic2.cpp
View File

@@ -23,9 +23,8 @@
****************************************************************************************/
#include <modules/globebrowsing/geodetics/geodetic2.h>
#include <modules/globebrowsing/globes/chunknode.h>
#include <modules/globebrowsing/geodetics/angle.h>
#include <modules/globebrowsing/geodetics/ellipsoid.h>
#include <ghoul/misc/assert.h>
@@ -33,260 +32,285 @@
#include <math.h>
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<Scalar>;
// 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<Scalar>;
// 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<Scalar>;
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<Scalar>;
// 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<Scalar>;
// 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<Scalar>;
// 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

174
modules/globebrowsing/geodetics/geodetic2.h
View File

@@ -31,6 +31,7 @@
#include <ostream>
// 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

9
modules/globebrowsing/globes/chunknode.cpp
View File

@@ -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<ChunkNode>(new ChunkNode(_owner, nwBounds, this));
_children[Quad::NORTH_EAST] = std::unique_ptr<ChunkNode>(new ChunkNode(_owner, neBounds, this));
_children[Quad::SOUTH_WEST] = std::unique_ptr<ChunkNode>(new ChunkNode(_owner, swBounds, this));
_children[Quad::NORTH_EAST] = std::unique_ptr<ChunkNode>(new ChunkNode(_owner, neBounds, this));
_children[Quad::SOUTH_EAST] = std::unique_ptr<ChunkNode>(new ChunkNode(_owner, seBounds, this));
}

41
modules/globebrowsing/globes/chunknode.h
View File

@@ -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<ChunkIndex> 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<ChunkNode> _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<ChunkNode> _children[4];
ChunkLodGlobe& _owner;
GeodeticPatch _patch;
bool _isVisible;
bool _isVisible;
};
} // namespace openspace

219
modules/globebrowsing/other/gdaldataconverter.cpp
View File

@@ -24,116 +24,159 @@
#include <modules/globebrowsing/other/gdaldataconverter.h>
#include <modules/globebrowsing/other/tileprovider.h>
#include <modules/globebrowsing/geodetics/angle.h>
namespace {
const std::string _loggerCat = "GdalDataConverter";
const std::string _loggerCat = "GdalDataConverter";
}
namespace openspace {
GdalDataConverter::GdalDataConverter()
{
GdalDataConverter::GdalDataConverter()
{
}
}
GdalDataConverter::~GdalDataConverter()
{
GdalDataConverter::~GdalDataConverter()
{
}
}
std::shared_ptr<Texture> GdalDataConverter::convertToOpenGLTexture(
GDALDataset* dataSet,
const GeodeticTileIndex& tileIndex,
int GLType)
{
int nRasters = dataSet->GetRasterCount();
std::shared_ptr<Texture> 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<GDALRasterBand*> 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> texture = std::shared_ptr<Texture>(new Texture(
static_cast<void*>(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<void*>(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> texture = std::shared_ptr<Texture>(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<Scalar>::fromRadians(geo.lat).asDegrees();
Scalar X = Angle<Scalar>::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

46
modules/globebrowsing/other/gdaldataconverter.h
View File

@@ -25,38 +25,46 @@
#ifndef __GDALDATACONVERTER_H__
#define __GDALDATACONVERTER_H__
#include <modules/globebrowsing/other/tileprovider.h>
//#include <modules/globebrowsing/other/tileprovider.h>
#include <ghoul/logging/logmanager.h>
#include <ghoul/opengl/texture.h>
#include <modules/globebrowsing/geodetics/geodetic2.h>
#include "gdal_priv.h"
#include <memory>
namespace openspace {
using namespace ghoul::opengl;
using namespace ghoul::opengl;
class GdalDataConverter
{
public:
GdalDataConverter();
~GdalDataConverter();
// forward declaration
class GeodeticTileIndex;
//class Geodetic2;
std::shared_ptr<Texture> 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<Texture> 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

232
modules/globebrowsing/other/texturetileset.cpp
View File

@@ -39,144 +39,146 @@
#include <glm/glm.hpp>
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<int>(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<int>(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<Texture> TextureTileSet::getTile(const GeodeticPatch& patch) {
return getTile(getTileIndex(patch));
}
std::shared_ptr<Texture> TextureTileSet::getTile(const GeodeticPatch& patch) {
return getTile(getTileIndex(patch));
}
std::shared_ptr<Texture> TextureTileSet::getTile(const GeodeticTileIndex& tileIndex) {
return _testTexture;
}
std::shared_ptr<Texture> 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

223
modules/globebrowsing/other/tileprovider.cpp
View File

@@ -22,7 +22,7 @@
* OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. *
****************************************************************************************/
#include <modules/globebrowsing/geodetics/geodetic2.h>
#include <modules/globebrowsing/other/tileprovider.h>
#include <openspace/engine/downloadmanager.h>
@@ -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<TextureLoadJob*>(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> 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<TextureLoadJob>(new TextureLoadJob(filePath, hashkey));
//_concurrentJobManager.enqueueJob(job);
std::shared_ptr<Texture> 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<DownloadManager::FileFuture> fileFuture = requestTile(tileIndex);
_fileFutureMap.insert_or_assign(hashkey, fileFuture);
}
return nullptr;
}
std::shared_ptr<Texture> TileProvider::loadAndInitTextureDisk(std::string filePath) {
auto textureReader = ghoul::io::TextureReader::ref();
std::shared_ptr<Texture> 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<void*>(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> texture = std::shared_ptr<Texture>(tex);
delete[] imageData;
// Do not free imageData since the texture now has ownership of it
return texture;
}
std::shared_ptr<DownloadManager::FileFuture> 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<DownloadManager::FileFuture> fileFuture = DownloadManager::ref().downloadFile(twmsRequestUrl, localTileFile, overrideFile);
return fileFuture;
}
} // namespace openspace

57
modules/globebrowsing/other/tileprovider.h
View File

@@ -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 <modules/globebrowsing/other/lrucache.h>
#include <modules/globebrowsing/other/concurrentjobmanager.h>
#include "gdal_priv.h"
#include <openspace/engine/downloadmanager.h>
#include <ghoul/logging/logmanager.h>
#include <ghoul/filesystem/filesystem.h> // absPath
#include <ghoul/opengl/texture.h>
#include <ghoul/io/texture/texturereader.h>
#include <openspace/engine/downloadmanager.h>
#include <modules/globebrowsing/geodetics/geodetic2.h>
#include <modules/globebrowsing/other/lrucache.h>
#include <modules/globebrowsing/other/concurrentjobmanager.h>
#include <modules/globebrowsing/other/gdaldataconverter.h>
//////////////////////////////////////////////////////////////////////////////////////
// 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> _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<Texture> getTile(const GeodeticTileIndex& tileIndex);
@@ -111,17 +100,23 @@ namespace openspace {
private:
std::shared_ptr<DownloadManager::FileFuture> requestTile(const GeodeticTileIndex&);
std::shared_ptr<Texture> loadAndInitTextureDisk(std::string filePath);
std::shared_ptr<Texture> getTileInternal(const GeodeticTileIndex& tileIndex, int GLType);
LRUCache<HashKey, std::shared_ptr<Texture>> _tileCache;
std::unordered_map<HashKey, std::shared_ptr<DownloadManager::FileFuture>> _fileFutureMap;
//LRUCache<HashKey, std::shared_ptr<DownloadManager::FileFuture>> _fileFutureCache;
//ConcurrentJobManager<Texture> _concurrentJobManager;
const std::string _filePath;
static bool hasInitializedGDAL;
GDALDataset* _gdalDataSet;
GdalDataConverter _converter;
};
} // namespace openspace
#endif // __TWMS_TILE_PROVIDER_H__
#endif // __TILE_PROVIDER_H__

26
modules/globebrowsing/rendering/patchrenderer.cpp
View File

@@ -77,7 +77,9 @@ namespace openspace {
LatLonPatchRenderer::LatLonPatchRenderer(shared_ptr<Grid> 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<ghoul::opengl::Texture> 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,

5
modules/globebrowsing/rendering/patchrenderer.h
View File

@@ -60,6 +60,7 @@ namespace openspace {
PatchRenderer();
~PatchRenderer();
protected:
unique_ptr<ProgramObject> _programObject;
@@ -102,10 +103,14 @@ namespace openspace {
public:
ClipMapPatchRenderer(shared_ptr<ClipMapGrid> grid);
void update();
void renderPatch(
const Geodetic2& patchSize,
const RenderData& data,
const Ellipsoid& ellipsoid);
private:
shared_ptr<ClipMapGrid> _grid;
};