mirror/OpenSpace
1
0
Fork 0
mirror of https://github.com/OpenSpace/OpenSpace.git synced 2026-01-06 11:39:49 -06:00
Code Issues Packages Projects Releases Wiki Activity
Files
c4bdd12e0aa2cfc48a1cafbc522609b5b2c3dc59
OpenSpace/modules/globebrowsing/tile/tiledataset.cpp

810 lines
34 KiB
C++
Raw Blame History

/*****************************************************************************************
* *
* OpenSpace *
* *
* Copyright (c) 2014-2016 *
* *
* Permission is hereby granted, free of charge, to any person obtaining a copy of this *
* software and associated documentation files (the "Software"), to deal in the Software *
* without restriction, including without limitation the rights to use, copy, modify, *
* merge, publish, distribute, sublicense, and/or sell copies of the Software, and to *
* permit persons to whom the Software is furnished to do so, subject to the following *
* conditions: *
* *
* The above copyright notice and this permission notice shall be included in all copies *
* or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, *
* INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A *
* PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT *
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF *
* CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE *
* OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. *
****************************************************************************************/
#include <limits>
#include <ogr_featurestyle.h>
#include <ogr_spatialref.h>
#include <ghoul/logging/logmanager.h>
#include <ghoul/filesystem/filesystem.h> // abspath
#include <ghoul/misc/assert.h>
#include <modules/globebrowsing/tile/tiledataset.h>
#include <modules/globebrowsing/tile/tileprovider/tileprovider.h>
#include <modules/globebrowsing/tile/tileioresult.h>
#include <modules/globebrowsing/geometry/angle.h>
#include <float.h>
#include <sstream>
#include <algorithm>
#include <gdal_priv.h>
#include <openspace/engine/openspaceengine.h>
#include <openspace/engine/configurationmanager.h>
namespace {
const std::string _loggerCat = "TileDataset";
}
namespace openspace {
namespace globebrowsing {
std::ostream& operator<<(std::ostream& os, const PixelRegion& pr) {
return os << pr.start.x << ", " << pr.start.y << " with size " << pr.numPixels.x << ", " << pr.numPixels.y;
}
//////////////////////////////////////////////////////////////////////////////////
// Tile Data Layout //
//////////////////////////////////////////////////////////////////////////////////
TileDataLayout::TileDataLayout() {
}
TileDataLayout::TileDataLayout(GDALDataset* dataSet, GLuint preferredGlType) {
// Assume all raster bands have the same data type
gdalType = preferredGlType != 0 ? TileDataType::getGdalDataType(preferredGlType) : dataSet->GetRasterBand(1)->GetRasterDataType();
glType = TileDataType::getOpenGLDataType(gdalType);
numRasters = dataSet->GetRasterCount();
bytesPerDatum = TileDataType::numberOfBytes(gdalType);
bytesPerPixel = bytesPerDatum * numRasters;
textureFormat = TileDataType::getTextureFormat(numRasters, gdalType);
}
IODescription IODescription::cut(PixelRegion::Side side, int pos) {
PixelRegion readPreCut = read.region;
PixelRegion writePreCut = write.region;
glm::dvec2 ratio;
ratio.x = write.region.numPixels.x / (double) read.region.numPixels.x;
ratio.y = write.region.numPixels.y / (double) read.region.numPixels.y;
double ratioRatio = ratio.x / ratio.y;
//ghoul_assert(glm::abs(ratioRatio - 1.0) < 0.01, "Different read/write aspect ratio!");
IODescription whatCameOff = *this;
whatCameOff.read.region = read.region.globalCut(side, pos);
PixelRange cutSize = whatCameOff.read.region.numPixels;
PixelRange localWriteCutSize = ratio * glm::dvec2(cutSize);
if (cutSize.x == 0 || cutSize.y == 0) {
ghoul_assert(read.region.equals(readPreCut), "Read region should not have been modified");
ghoul_assert(write.region.equals(writePreCut), "Write region should not have been modified");
}
int localWriteCutPos = (side % 2 == 0) ? localWriteCutSize.x : localWriteCutSize.y;
whatCameOff.write.region = write.region.localCut(side, localWriteCutPos);
return whatCameOff;
}
//////////////////////////////////////////////////////////////////////////////////
// Tile Dataset //
//////////////////////////////////////////////////////////////////////////////////
const glm::ivec2 TileDataset::tilePixelStartOffset = glm::ivec2(-2);
const glm::ivec2 TileDataset::tilePixelSizeDifference = glm::ivec2(4);
const PixelRegion TileDataset::padding = PixelRegion(tilePixelStartOffset, tilePixelSizeDifference);
bool TileDataset::GdalHasBeenInitialized = false;
TileDataset::TileDataset(const std::string& gdalDatasetDesc, const Configuration& config)
: _config(config)
, hasBeenInitialized(false)
{
_initData = { "", gdalDatasetDesc, config.minimumTilePixelSize, config.dataType };
ensureInitialized();
_initData.initDirectory = CPLGetCurrentDir();
}
void TileDataset::reset() {
_cached._maxLevel = -1;
if (_dataset != nullptr) {
GDALClose((GDALDatasetH)_dataset);
}
initialize();
}
float TileDataset::noDataValueAsFloat() {
float noDataValue;
if (_dataset && _dataset->GetRasterBand(1))
{
noDataValue = _dataset->GetRasterBand(1)->GetNoDataValue();;
}
else
noDataValue = std::numeric_limits<float>::min();
return noDataValue;
}
void TileDataset::ensureInitialized() {
if (!hasBeenInitialized) {
initialize();
hasBeenInitialized = true;
}
}
void TileDataset::initialize() {
gdalEnsureInitialized();
_dataset = gdalDataset(_initData.gdalDatasetDesc);
//Do any other initialization needed for the TileDataset
_dataLayout = TileDataLayout(_dataset, _initData.dataType);
_depthTransform = calculateTileDepthTransform();
_cached._tileLevelDifference = calculateTileLevelDifference(_initData.minimumPixelSize);
LDEBUG(_initData.gdalDatasetDesc << " - " << _cached._tileLevelDifference);
}
void TileDataset::gdalEnsureInitialized() {
if (!GdalHasBeenInitialized) {
GDALAllRegister();
CPLSetConfigOption("GDAL_DATA", absPath("${MODULE_GLOBEBROWSING}/gdal_data").c_str());
setGdalProxyConfiguration();
GdalHasBeenInitialized = true;
}
}
void TileDataset::setGdalProxyConfiguration() {
ghoul::Dictionary proxySettings;
bool proxyEnabled = OsEng.configurationManager().getValue(ConfigurationManager::KeyHttpProxy, proxySettings);
if (proxyEnabled) {
std::string proxyAddress, proxyPort, proxyUser, proxyPassword, proxyAuth;
bool success = proxySettings.getValue(ConfigurationManager::PartHttpProxyAddress, proxyAddress);
success &= proxySettings.getValue(ConfigurationManager::PartHttpProxyPort, proxyPort);
proxySettings.getValue(ConfigurationManager::PartHttpProxyAuthentication, proxyAuth);
std::string proxyAuthString = "BASIC";
if (proxyAuth == "basic" || proxyAuth == "") {
proxyAuthString = "BASIC";
} else if (proxyAuth == "ntlm") {
proxyAuthString = "NTLM";
} else if (proxyAuth == "digest") {
proxyAuthString = "DIGEST";
} else if (proxyAuth == "any") {
proxyAuthString = "ANY";
} else {
success = false;
}
bool userAndPassword = proxySettings.getValue(ConfigurationManager::PartHttpProxyUser, proxyUser);
userAndPassword &= proxySettings.getValue(ConfigurationManager::PartHttpProxyPassword, proxyPassword);
if (success) {
std::string proxy = proxyAddress + ":" + proxyPort;
CPLSetConfigOption("GDAL_HTTP_PROXY", proxy.c_str());
LDEBUG("Using proxy server " << proxy);
if (userAndPassword) {
std::string proxyUserPwd = proxyUser + ":" + proxyPassword;
CPLSetConfigOption("GDAL_HTTP_PROXYUSERPWD", proxyUserPwd.c_str());
CPLSetConfigOption("GDAL_HTTP_PROXYAUTH", proxyAuthString.c_str());
LDEBUG("Using authentication method: " << proxyAuthString);
}
} else {
LERROR("Invalid proxy settings for GDAL");
}
} else {
LDEBUG("Setting up GDAL without proxy server");
}
}
GDALDataset* TileDataset::gdalDataset(const std::string& gdalDatasetDesc) {
GDALDataset* dataset = (GDALDataset *)GDALOpen(gdalDatasetDesc.c_str(), GA_ReadOnly);
if (!dataset) {
std::string correctedPath = ghoul::filesystem::FileSystem::ref().pathByAppendingComponent(_initData.initDirectory, gdalDatasetDesc);
dataset = (GDALDataset *)GDALOpen(correctedPath.c_str(), GA_ReadOnly);
if (!dataset) {
throw ghoul::RuntimeError("Failed to load dataset:\n" + gdalDatasetDesc);
}
}
// Commenting away the following for now since it is not supported for older
// versions of GDAL. Only used for debug info.
/*
const std::string originalDriverName = dataset->GetDriverName();
if (originalDriverName != "WMS") {
LDEBUG(" " << originalDriverName);
LDEBUG(" " << dataset->GetGCPProjection());
LDEBUG(" " << dataset->GetProjectionRef());
GDALDriver* driver = dataset->GetDriver();
char** metadata = driver->GetMetadata();
for (int i = 0; metadata[i] != nullptr; i++) {
LDEBUG(" " << metadata[i]);
}
const char* in_memory = "";
//GDALDataset* vrtDataset = driver->CreateCopy(in_memory, dataset, false, nullptr, nullptr, nullptr);
}
*/
return dataset;
}
TileDataset::~TileDataset() {
delete _dataset;
}
//////////////////////////////////////////////////////////////////////////////////
// Public interface //
//////////////////////////////////////////////////////////////////////////////////
std::shared_ptr<TileIOResult> TileDataset::readTileData(TileIndex tileIndex) {
ensureInitialized();
IODescription io = getIODescription(tileIndex);
CPLErr worstError = CPLErr::CE_None;
// Build the Tile IO Result from the data we queride
std::shared_ptr<TileIOResult> result = std::make_shared<TileIOResult>();
result->imageData = readImageData(io, worstError);
result->error = worstError;
result->tileIndex = tileIndex;
result->dimensions = glm::uvec3(io.write.region.numPixels, 1);
result->nBytesImageData = io.write.totalNumBytes;
if (_config.doPreProcessing) {
result->preprocessData = preprocess(result, io.write.region);
result->error = std::max(result->error, postProcessErrorCheck(result, io));
}
return result;
}
std::shared_ptr<TileIOResult> TileDataset::defaultTileData() {
ensureInitialized();
PixelRegion pixelRegion = { PixelCoordinate(0, 0), PixelRange(16, 16) };
std::shared_ptr<TileIOResult> result = std::make_shared<TileIOResult>();
result->tileIndex = { 0, 0, 0 };
result->dimensions = glm::uvec3(pixelRegion.numPixels, 1);
result->nBytesImageData = result->dimensions.x * result->dimensions.y * _dataLayout.bytesPerPixel;
result->imageData = new char[result->nBytesImageData];
for (size_t i = 0; i < result->nBytesImageData; ++i) {
result->imageData[i] = 0;
}
result->error = CPLErr::CE_None;
if (_config.doPreProcessing) {
result->preprocessData = preprocess(result, pixelRegion);
//result->error = std::max(result->error, postProcessErrorCheck(result, io));
}
return result;
}
int TileDataset::maxChunkLevel() {
ensureInitialized();
if (_cached._maxLevel < 0) {
int numOverviews = _dataset->GetRasterBand(1)->GetOverviewCount();
_cached._maxLevel = -_cached._tileLevelDifference;
if (numOverviews > 0) {
_cached._maxLevel += numOverviews - 1;
}
}
return _cached._maxLevel;
}
TileDepthTransform TileDataset::getDepthTransform() {
ensureInitialized();
return _depthTransform;
}
const TileDataLayout& TileDataset::getDataLayout() {
ensureInitialized();
return _dataLayout;
}
//////////////////////////////////////////////////////////////////////////////////
// Initialization //
//////////////////////////////////////////////////////////////////////////////////
int TileDataset::calculateTileLevelDifference(int minimumPixelSize) {
GDALRasterBand* firstBand = _dataset->GetRasterBand(1);
GDALRasterBand* maxOverview;
int numOverviews = firstBand->GetOverviewCount();
int sizeLevel0;
if (numOverviews <= 0) { // No overviews. Use first band.
maxOverview = firstBand;
}
else { // Pick the highest overview.
maxOverview = firstBand->GetOverview(numOverviews - 1);
}
sizeLevel0 = maxOverview->GetXSize();
double diff = log2(minimumPixelSize) - log2(sizeLevel0);
return diff;
}
TileDepthTransform TileDataset::calculateTileDepthTransform() {
GDALRasterBand* firstBand = _dataset->GetRasterBand(1);
bool isFloat = (_dataLayout.gdalType == GDT_Float32 || _dataLayout.gdalType == GDT_Float64);
double maximumValue = isFloat ? 1.0 : TileDataType::getMaximumValue(_dataLayout.gdalType);
TileDepthTransform transform;
transform.depthOffset = firstBand->GetOffset();
transform.depthScale = firstBand->GetScale() * maximumValue;
return transform;
}
//////////////////////////////////////////////////////////////////////////////////
// GDAL helper methods //
//////////////////////////////////////////////////////////////////////////////////
bool TileDataset::gdalHasOverviews() const {
return _dataset->GetRasterBand(1)->GetOverviewCount() > 0;
}
int TileDataset::gdalOverview(const PixelRange& regionSizeOverviewZero) const {
GDALRasterBand* firstBand = _dataset->GetRasterBand(1);
int minNumPixels0 = glm::min(regionSizeOverviewZero.x, regionSizeOverviewZero.y);
int overviews = firstBand->GetOverviewCount();
GDALRasterBand* maxOverview = overviews ? firstBand->GetOverview(overviews - 1) : firstBand;
int sizeLevel0 = maxOverview->GetXSize();
// The dataset itself may not have overviews but even if it does not, an overview
// for the data region can be calculated and possibly be used to sample greater
// Regions of the original dataset.
int ov = std::log2(minNumPixels0) - std::log2(sizeLevel0 + 1) - _cached._tileLevelDifference;
ov = glm::clamp(ov, 0, overviews - 1);
return ov;
}
int TileDataset::gdalOverview(const TileIndex& tileIndex) const {
int overviews = _dataset->GetRasterBand(1)->GetOverviewCount();
int ov = overviews - (tileIndex.level + _cached._tileLevelDifference + 1);
return glm::clamp(ov, 0, overviews - 1);
}
int TileDataset::gdalVirtualOverview(const TileIndex& tileIndex) const {
int overviews = _dataset->GetRasterBand(1)->GetOverviewCount();
int ov = overviews - (tileIndex.level + _cached._tileLevelDifference + 1);
return ov;
}
PixelRegion TileDataset::gdalPixelRegion(GDALRasterBand* rasterBand) const {
PixelRegion gdalRegion;
gdalRegion.start.x = 0;
gdalRegion.start.y = 0;
gdalRegion.numPixels.x = rasterBand->GetXSize();
gdalRegion.numPixels.y = rasterBand->GetYSize();
return gdalRegion;
}
PixelRegion TileDataset::gdalPixelRegion(const GeodeticPatch& geodeticPatch) const {
Geodetic2 nwCorner = geodeticPatch.getCorner(Quad::NORTH_WEST);
Geodetic2 swCorner = geodeticPatch.getCorner(Quad::SOUTH_EAST);
PixelCoordinate pixelStart = geodeticToPixel(nwCorner);
PixelCoordinate pixelEnd = geodeticToPixel(swCorner);
PixelRegion gdalRegion(pixelStart, pixelEnd - pixelStart);
return gdalRegion;
}
GDALRasterBand* TileDataset::gdalRasterBand(int overview, int raster) const {
GDALRasterBand* rasterBand = _dataset->GetRasterBand(raster);
int numberOfOverviews = rasterBand->GetOverviewCount();
rasterBand = gdalHasOverviews() ? rasterBand->GetOverview(overview) : rasterBand;
ghoul_assert(rasterBand != nullptr, "Rasterband is null");
return rasterBand;
}
std::array<double, 6> TileDataset::getGeoTransform() const {
std::array<double, 6> padfTransform;
CPLErr err = _dataset->GetGeoTransform(&padfTransform[0]);
if (err == CE_Failure) {
GeodeticPatch globalCoverage(Geodetic2(0,0), Geodetic2(M_PI / 2, M_PI));
padfTransform[1] = Angle<Scalar>::fromRadians(
globalCoverage.size().lon).asDegrees() / _dataset->GetRasterXSize();
padfTransform[5] = -Angle<Scalar>::fromRadians(
globalCoverage.size().lat).asDegrees() / _dataset->GetRasterYSize();
padfTransform[0] = Angle<Scalar>::fromRadians(
globalCoverage.getCorner(Quad::NORTH_WEST).lon).asDegrees();
padfTransform[3] = Angle<Scalar>::fromRadians(
globalCoverage.getCorner(Quad::NORTH_WEST).lat).asDegrees();
padfTransform[2] = 0;
padfTransform[4] = 0;
}
return padfTransform;
}
//////////////////////////////////////////////////////////////////////////////////
// ReadTileData helper functions //
//////////////////////////////////////////////////////////////////////////////////
PixelCoordinate TileDataset::geodeticToPixel(const Geodetic2& geo) const {
std::array<double, 6> padfTransform = getGeoTransform();
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];
ghoul_assert(abs(X - Xp) < 1e-10, "inverse should yield X as before");
ghoul_assert(abs(Y - Yp) < 1e-10, "inverse should yield Y as before");
return PixelCoordinate(glm::round(P), glm::round(L));
}
Geodetic2 TileDataset::pixelToGeodetic(const PixelCoordinate& p) const {
std::array<double, 6> padfTransform = getGeoTransform();
Geodetic2 geodetic;
// Should be using radians and not degrees?
geodetic.lon = padfTransform[0] + p.x * padfTransform[1] + p.y * padfTransform[2];
geodetic.lat = padfTransform[3] + p.x * padfTransform[4] + p.y * padfTransform[5];
return geodetic;
}
IODescription TileDataset::getIODescription(const TileIndex& tileIndex) const {
IODescription io;
io.read.region = gdalPixelRegion(tileIndex);
if (gdalHasOverviews()) {
int overview = gdalOverview(tileIndex);
io.read.overview = overview;
io.read.region.downscalePow2(overview + 1);
io.write.region = io.read.region;
io.read.region.pad(padding);
}
else {
io.read.overview = 0;
io.write.region = io.read.region;
int virtualOverview = gdalVirtualOverview(tileIndex);
io.write.region.downscalePow2(virtualOverview + 1);
PixelRegion scaledPadding = padding;
scaledPadding.upscalePow2(std::max(virtualOverview + 1, 0));
io.read.region.pad(scaledPadding);
}
// For correct sampling in height dataset, we need to pad the texture tile
io.write.region.pad(padding);
PixelRange preRound = io.write.region.numPixels;
io.write.region.roundDownToQuadratic();
io.write.region.roundUpNumPixelToNearestMultipleOf(2);
if (preRound != io.write.region.numPixels) {
LDEBUG(tileIndex << " | " << preRound.x << ", " << preRound.y << " --> " << io.write.region.numPixels.x << ", " << io.write.region.numPixels.y);
}
io.write.region.start = PixelCoordinate(0, 0); // write region starts in origin
io.write.bytesPerLine = _dataLayout.bytesPerPixel * io.write.region.numPixels.x;
io.write.totalNumBytes = io.write.bytesPerLine * io.write.region.numPixels.y;
return io;
}
char* TileDataset::readImageData(IODescription& io, CPLErr& worstError) const {
// allocate memory for the image
char* imageData = new char[io.write.totalNumBytes];
// Read the data (each rasterband is a separate channel)
for (size_t i = 0; i < _dataLayout.numRasters; i++) {
GDALRasterBand* rasterBand = gdalRasterBand(io.read.overview, i + 1);
// The final destination pointer is offsetted by one datum byte size
// for every raster (or data channel, i.e. R in RGB)
char* dataDestination = imageData + (i * _dataLayout.bytesPerDatum);
CPLErr err = repeatedRasterIO(rasterBand, io, dataDestination);
// CE_None = 0, CE_Debug = 1, CE_Warning = 2, CE_Failure = 3, CE_Fatal = 4
worstError = std::max(worstError, err);
}
// GDAL reads pixel lines top to bottom, we want the opposit
return imageData;
}
CPLErr TileDataset::repeatedRasterIO(GDALRasterBand* rasterBand, const IODescription& fullIO, char* dataDestination, int depth) const {
std::string spaces = " ";
std::string indentation = spaces.substr(0, 2 * depth);
CPLErr worstError = CPLErr::CE_None;
// NOTE:
// Ascii graphics illustrates the implementation details of this method, for one
// specific case. Even though the illustrated case is specific, readers can
// hopefully find it useful to get the general idea.
// Make a copy of the full IO desription as we will have to modify it
IODescription io = fullIO;
PixelRegion gdalRegion = gdalPixelRegion(rasterBand);
// Example:
// We have an io description that defines a WRITE and a READ region.
// In this case the READ region extends outside of the defined gdal region,
// meaning we will have to do wrapping
// io.write.region io.read.region
// | |
// V V
// +-------+ +-------+
// | | | |--------+
// | | | | |
// | | | | |
// +-------+ +-------+ |
// | | <-- gdalRegion
// | |
// +--------------+
//LDEBUG(indentation << "-");
//LDEBUG(indentation << "repeated read: " << io.read.region);
//LDEBUG(indentation << "repeated write: " << io.write.region);
bool didCutOff = false;
if (!io.read.region.isInside(gdalRegion)) {
// Loop through each side: left, top, right, bottom
for (int i = 0; i < 4; ++i) {
// Example:
// We are currently considering the left side of the pixel region
PixelRegion::Side side = (PixelRegion::Side) i;
IODescription cutoff = io.cut(side, gdalRegion.edge(side));
// Example:
// We cut off the left part that was outside the gdal region, and we now
// have an additional io description for the cut off region.
// Note that the cut-method used above takes care of the corresponding
// WRITE region for us.
// cutoff.write.region cutoff.read.region
// | io.write.region | io.read.region
// | | | |
// V V V V
// +-+-----+ +-+-----+
// | | | | | |--------+
// | | | | | | |
// | | | | | | |
// +-+-----+ +-+-----+ |
// | | <-- gdalRegion
// | |
// +--------------+
if (cutoff.read.region.area() > 0) {
didCutOff = true;
// Wrap by repeating
PixelRegion::Side oppositeSide = (PixelRegion::Side) ((i + 2) % 4);
cutoff.read.region.align(oppositeSide, gdalRegion.edge(oppositeSide));
// Example:
// The cut off region is wrapped to the opposite side of the region,
// i.e. "repeated". Note that we don't want WRITE region to change,
// we're only wrapping the READ region.
// cutoff.write.region io.read.region cutoff.read.region
// | io.write.region | |
// | | V V
// V V +-----+ +-+
// +-+-----+ | |------| |
// | | | | | | |
// | | | | | | |
// | | | +-----+ +-+
// +-+-----+ | | <-- gdalRegion
// | |
// +--------------+
// Example:
// The cutoff region has been repeated along one of its sides, but
// as we can see in this example, it still has a top part outside the
// defined gdal region. This is handled through recursion.
CPLErr err = repeatedRasterIO(rasterBand, cutoff, dataDestination, depth + 1);
worstError = std::max(worstError, err);
}
}
}
//if (depth == 0) {
//LDEBUG(indentation << "main rasterIO read: " << io.read.region);
//LDEBUG(indentation << "main rasterIO write: " << io.write.region);
//}
else if (worstError > CPLErr::CE_None) {
LDEBUG(indentation << "Error reading padding: " << worstError);
}
CPLErr err = rasterIO(rasterBand, io, dataDestination);
worstError = std::max(worstError, err);
// The return error from a repeated rasterIO is ONLY based on the main region,
// which in the usual case will cover the main area of the patch anyway
return err;
}
CPLErr TileDataset::rasterIO(GDALRasterBand* rasterBand, const IODescription& io, char* dataDestination) const {
PixelRegion gdalRegion = gdalPixelRegion(rasterBand);
ghoul_assert(io.read.region.isInside(gdalRegion), "write region of bounds!");
ghoul_assert(io.write.region.start.x >= 0 && io.write.region.start.y >= 0, "Invalid write region");
PixelCoordinate end = io.write.region.end();
size_t largestIndex = (end.y - 1) * io.write.bytesPerLine + (end.x - 1) * _dataLayout.bytesPerPixel;
ghoul_assert(largestIndex <= io.write.totalNumBytes, "Invalid write region");
char * dataDest = dataDestination;
// OBS! GDAL reads pixels top to bottom, but we want our pixels bottom to top.
// Therefore, we increment the destination pointer to the last line on in the
// buffer, and the we specify in the rasterIO call that we want negative line
// spacing. Doing this compensates the flipped Y axis
dataDest += (io.write.totalNumBytes - io.write.bytesPerLine);
// handle requested write region
dataDest -= io.write.region.start.y * io.write.bytesPerLine; // note -= since flipped y axis
dataDest += io.write.region.start.x * _dataLayout.bytesPerPixel;
return rasterBand->RasterIO(
GF_Read,
io.read.region.start.x, // Begin read x
io.read.region.start.y, // Begin read y
io.read.region.numPixels.x, // width to read x
io.read.region.numPixels.y, // width to read y
dataDest, // Where to put data
io.write.region.numPixels.x, // width to write x in destination
io.write.region.numPixels.y, // width to write y in destination
_dataLayout.gdalType, // Type
_dataLayout.bytesPerPixel, // Pixel spacing
-io.write.bytesPerLine); // Line spacing
}
std::shared_ptr<TilePreprocessData> TileDataset::preprocess(std::shared_ptr<TileIOResult> result, const PixelRegion& region) const {
size_t bytesPerLine = _dataLayout.bytesPerPixel * region.numPixels.x;
size_t totalNumBytes = bytesPerLine * region.numPixels.y;
TilePreprocessData* preprocessData = new TilePreprocessData();
preprocessData->maxValues.resize(_dataLayout.numRasters);
preprocessData->minValues.resize(_dataLayout.numRasters);
preprocessData->hasMissingData.resize(_dataLayout.numRasters);
std::vector<float> noDataValues;
noDataValues.resize(_dataLayout.numRasters);
for (size_t c = 0; c < _dataLayout.numRasters; c++) {
preprocessData->maxValues[c] = -FLT_MAX;
preprocessData->minValues[c] = FLT_MAX;
preprocessData->hasMissingData[c] = false;
noDataValues[c] = _dataset->GetRasterBand(1)->GetNoDataValue();
}
for (size_t y = 0; y < region.numPixels.y; y++) {
size_t yi = (region.numPixels.y - 1 - y) * bytesPerLine;
size_t i = 0;
for (size_t x = 0; x < region.numPixels.x; x++) {
for (size_t c = 0; c < _dataLayout.numRasters; c++) {
float noDataValue = _dataset->GetRasterBand(c + 1)->GetNoDataValue();
float val = TileDataType::interpretFloat(_dataLayout.gdalType, &(result->imageData[yi + i]));
if (val != noDataValue) {
preprocessData->maxValues[c] = std::max(val, preprocessData->maxValues[c]);
preprocessData->minValues[c] = std::min(val, preprocessData->minValues[c]);
}
else {
preprocessData->hasMissingData[c] = true;
}
i += _dataLayout.bytesPerDatum;
}
}
}
for (size_t c = 0; c < _dataLayout.numRasters; c++) {
if (preprocessData->maxValues[c] > 8800.0f) {
//LDEBUG("Bad preprocess data: " << preprocessData->maxValues[c] << " at " << region.tileIndex);
}
}
return std::shared_ptr<TilePreprocessData>(preprocessData);
}
CPLErr TileDataset::postProcessErrorCheck(std::shared_ptr<const TileIOResult> result, const IODescription& io) const{
int success;
double missingDataValue = gdalRasterBand(io.read.overview)->GetNoDataValue(&success);
if (!success) {
missingDataValue = 32767; // missing data value for TERRAIN.wms. Should be specified in xml
}
bool hasMissingData = false;
for (size_t c = 0; c < _dataLayout.numRasters; c++) {
hasMissingData |= result->preprocessData->maxValues[c] == missingDataValue;
}
bool onHighLevel = result->tileIndex.level > 6;
if (hasMissingData && onHighLevel) {
return CE_Fatal;
}
// ugly test for heightmap overlay
if (_dataLayout.textureFormat.ghoulFormat == Texture::Format::RG) {
// check the alpha
if (result->preprocessData->maxValues[1] == 0.0
&& result->preprocessData->minValues[1] == 0.0)
{
//return CE_Warning;
}
}
return CE_None;
}
} // namespace globebrowsing
} // namespace openspace
Reference in New Issue View Git Blame Copy Permalink