mirror/OpenSpace
1
0
Fork 0
mirror of https://github.com/OpenSpace/OpenSpace.git synced 2026-04-23 20:50:59 -05:00
Code Issues Packages Projects Releases Wiki Activity

Replace GdalDataRegion with more general class PixelRegion

Browse Source
This commit is contained in:
Erik Broberg
2016-06-30 11:33:09 -04:00
parent 161d722623
commit 1b6d17575b
2 changed files with 157 additions and 144 deletions
Download Patch File Download Diff File Expand all files Collapse all files
modules/globebrowsing/tile/tiledataset.cpp
+106 -133
View File
@@ -68,110 +68,7 @@ namespace openspace {
bytesPerPixel = bytesPerDatum * numRasters;
textureFormat = TileDataType::getTextureFormat(numRasters, gdalType);
}
//////////////////////////////////////////////////////////////////////////////////
// GDAL Data Region //
//////////////////////////////////////////////////////////////////////////////////
GdalDataRegion::GdalDataRegion(GDALDataset * dataSet,
const ChunkIndex& chunkIndex, int tileLevelDifference) {
GDALRasterBand* firstBand = dataSet->GetRasterBand(1);
// Assume all raster bands have the same data type
// Level = overviewCount - overview (default, levels may be overridden)
int numOverviews = firstBand->GetOverviewCount();
// Generate a patch from the chunkIndex, extract the bounds which
// are used to calculated where in the GDAL data set to read data.
// pixelStart0 and pixelEnd0 defines the interval in the pixel space
// at overview 0
GeodeticPatch patch = GeodeticPatch(chunkIndex);
glm::uvec2 pixelStart0 = geodeticToPixel(dataSet, patch.getCorner(Quad::NORTH_WEST));
glm::uvec2 pixelEnd0 = geodeticToPixel(dataSet, patch.getCorner(Quad::SOUTH_EAST));
glm::uvec2 numPixels0 = pixelEnd0 - pixelStart0;
// Calculate a suitable overview to choose from the GDAL dataset
int minNumPixels0 = glm::min(numPixels0.x, numPixels0.y);
GDALRasterBand* maxOverview;
if (numOverviews <= 0) {
maxOverview = firstBand;
}
else {
maxOverview = firstBand->GetOverview(numOverviews - 1);
}
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) - tileLevelDifference;
if (numOverviews > 0) {
ov = glm::clamp(ov, 0, numOverviews - 1);
}
// Convert the interval [pixelStart0, pixelEnd0] to pixel space at
// the calculated suitable overview, ov. using a >> b = a / 2^b
int toShift = ov + 1;
// Set member variables
overview = ov;
pixelStart = glm::uvec2(pixelStart0.x >> toShift, pixelStart0.y >> toShift);
pixelEnd = glm::uvec2(pixelEnd0.x >> toShift, pixelEnd0.y >> toShift);
numPixels = pixelEnd - pixelStart;
}
glm::uvec2 GdalDataRegion::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];
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 glm::uvec2(glm::round(P), glm::round(L));
}
@@ -185,6 +82,9 @@ namespace openspace {
const glm::ivec2 TileDataset::tilePixelStartOffset = glm::ivec2(0, 0);
const glm::ivec2 TileDataset::tilePixelSizeDifference = glm::ivec2(0, 0);
const PixelRegion TileDataset::padding = PixelRegion(tilePixelStartOffset, tilePixelSizeDifference);
bool TileDataset::GdalHasBeenInitialized = false;
TileDataset::TileDataset(const std::string& gdalDatasetDesc, int minimumPixelSize,
@@ -231,7 +131,7 @@ namespace openspace {
return log2(minimumPixelSize) - log2(sizeLevel0);
}
const int TileDataset::calculateMaxLevel(int tileLevelDifference) {
int TileDataset::calculateMaxLevel(int tileLevelDifference) {
int numOverviews = _dataset->GetRasterBand(1)->GetOverviewCount();
if (numOverviews <= 0) { // No overviews.
return -tileLevelDifference;
@@ -254,6 +154,78 @@ namespace openspace {
return transform;
}
PixelCoordinate TileDataset::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];
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));
}
PixelRegion TileDataset::gdalPixelRegion(const GeodeticPatch& geodeticPatch) const {
Geodetic2 nwCorner = geodeticPatch.getCorner(Quad::NORTH_WEST);
Geodetic2 swCorner = geodeticPatch.getCorner(Quad::SOUTH_EAST);
PixelCoordinate pixelStart = TileDataset::geodeticToPixel(_dataset, nwCorner);
PixelCoordinate pixelEnd = TileDataset::geodeticToPixel(_dataset, swCorner);
PixelRegion gdalRegion(pixelStart, pixelEnd- pixelStart);
return gdalRegion;
}
int TileDataset::gdalOverview(const PixelCoordinate& 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) - _tileLevelDifference;
ov = glm::clamp(ov, 0, overviews - 1);
return ov;
}
int TileDataset::gdalOverview(const ChunkIndex& chunkIndex) const {
int overviews = _dataset->GetRasterBand(1)->GetOverviewCount();
int ov = _dataset->GetRasterBand(1)->GetOverviewCount() - (chunkIndex.level + _tileLevelDifference + 1);
return glm::clamp(ov, 0, overviews - 1);
}
int TileDataset::getMaximumLevel() const {
return _maxLevel;
}
@@ -263,7 +235,10 @@ namespace openspace {
}
std::shared_ptr<TileIOResult> TileDataset::readTileData(ChunkIndex chunkIndex) {
GdalDataRegion region(_dataset, chunkIndex, _tileLevelDifference);
//GdalDataRegion region(_dataset, chunkIndex, _tileLevelDifference);
PixelRegion region = gdalPixelRegion(chunkIndex);
int overview = gdalOverview(chunkIndex);
region.shrinkPow2(overview + 1);
size_t bytesPerLine = _dataLayout.bytesPerPixel * region.numPixels.x;
size_t totalNumBytes = bytesPerLine * region.numPixels.y;
@@ -273,37 +248,35 @@ namespace openspace {
// Read the data (each rasterband is a separate channel)
for (size_t i = 0; i < _dataLayout.numRasters; i++) {
GDALRasterBand* rasterBand;
int pixelSourceScale = 1;
if (_dataset->GetRasterBand(i + 1)->GetOverviewCount() <= 0){
rasterBand = _dataset->GetRasterBand(i + 1);
pixelSourceScale = pow(2, region.overview + 1);
}
else {
rasterBand = _dataset->GetRasterBand(i + 1)->GetOverview(region.overview);
pixelSourceScale = 1;
}
char* dataDestination = imageData + (i * _dataLayout.bytesPerDatum);
int pixelStartX = region.pixelStart.x * pixelSourceScale + tilePixelStartOffset.x;
int pixelStartY = region.pixelStart.y * pixelSourceScale + tilePixelStartOffset.y;
int pixelWidthX = region.numPixels.x * pixelSourceScale + tilePixelSizeDifference.x;
int pixelWidthY = region.numPixels.y * pixelSourceScale + tilePixelSizeDifference.y;
PixelRegion readRegion(region);
GDALRasterBand* rasterBand;
if (_dataset->GetRasterBand(i + 1)->GetOverviewCount() <= 0){
rasterBand = _dataset->GetRasterBand(i + 1);
int pixelSourceScale = pow(2, overview + 1);
readRegion.scale(pixelSourceScale);
}
else {
rasterBand = _dataset->GetRasterBand(i + 1)->GetOverview(overview);
}
PixelRegion gdalPixelRegion;
gdalPixelRegion.start = glm::ivec2(0);
gdalPixelRegion.numPixels = glm::ivec2(rasterBand->GetXSize(), rasterBand->GetYSize());
readRegion.addPadding(padding);
readRegion.clamp(gdalPixelRegion);
// Clamp to be inside dataset
//pixelStartX = glm::max(pixelStartX, 0);
//pixelStartY = glm::max(pixelStartY, 0);
//pixelWidthX = glm::min(pixelStartX + pixelWidthX, rasterBand->GetXSize()) - pixelStartX;
//pixelWidthY = glm::min(pixelStartY + pixelWidthY, rasterBand->GetYSize()) - pixelStartY;
CPLErr err = rasterBand->RasterIO(
GF_Read,
pixelStartX, // Begin read x
pixelStartY, // Begin read y
pixelWidthX, // width to read x
pixelWidthY, // width to read y
readRegion.start.x, // Begin read x
readRegion.start.y, // Begin read y
readRegion.numPixels.x, // width to read x
readRegion.numPixels.y, // width to read y
dataDestination, // Where to put data
region.numPixels.x, // width to write x in destination
region.numPixels.y, // width to write y in destination
@@ -324,7 +297,7 @@ namespace openspace {
if (_doPreprocessing) {
result->preprocessData = preprocess(imageData, region, _dataLayout);
int success;
auto gdalOverview = _dataset->GetRasterBand(1)->GetOverview(region.overview);
auto gdalOverview = _dataset->GetRasterBand(1)->GetOverview(overview);
double missingDataValue = gdalOverview->GetNoDataValue(&success);
if (!success) {
missingDataValue = 32767; // missing data value
@@ -343,7 +316,7 @@ namespace openspace {
return result;
}
char* TileDataset::getImageDataFlippedY(const GdalDataRegion& region,
char* TileDataset::getImageDataFlippedY(const PixelRegion& region,
const TileDataLayout& dataLayout, const char* imageData)
{
size_t bytesPerLine = dataLayout.bytesPerPixel * region.numPixels.x;
@@ -375,7 +348,7 @@ namespace openspace {
std::shared_ptr<TilePreprocessData> TileDataset::preprocess(const char* imageData,
const GdalDataRegion& region, const TileDataLayout& dataLayout)
const PixelRegion& region, const TileDataLayout& dataLayout)
{
size_t bytesPerLine = dataLayout.bytesPerPixel * region.numPixels.x;
size_t totalNumBytes = bytesPerLine * region.numPixels.y;
modules/globebrowsing/tile/tiledataset.h
+51 -11
View File
@@ -60,23 +60,53 @@ namespace openspace {
TextureFormat textureFormat;
};
struct GdalDataRegion {
GdalDataRegion(GDALDataset* dataSet, const ChunkIndex& chunkIndex, int tileLevelDifference);
typedef glm::ivec2 PixelCoordinate;
glm::uvec2 pixelStart;
glm::uvec2 pixelEnd;
glm::uvec2 numPixels;
struct PixelRegion {
int overview;
PixelRegion() : start(0), numPixels(0) { }
PixelRegion(const PixelRegion& o) : start(o.start), numPixels(o.numPixels) { }
PixelRegion(const PixelCoordinate& pixelStart, const PixelCoordinate& numberOfPixels)
: start(pixelStart), numPixels(numberOfPixels) { }
static glm::uvec2 geodeticToPixel(GDALDataset* dataSet, const Geodetic2& geo);
PixelCoordinate start;
PixelCoordinate numPixels;
void addPadding(const PixelRegion& padding) {
start += padding.start;
numPixels += padding.numPixels;
}
void clamp(const PixelRegion& boundingRegion) {
start = glm::max(start, boundingRegion.start);
numPixels = glm::min(end(), boundingRegion.end()) - start;
}
PixelCoordinate end() const {
return start + numPixels;
}
void scale(const glm::dvec2& s) {
start = PixelCoordinate(glm::round(s * glm::dvec2(start)));
numPixels = PixelCoordinate(glm::round(s * glm::dvec2(numPixels)));
}
void scale(double s) {
scale(glm::dvec2(s));
}
void shrinkPow2(int exponent) {
start.x >>= exponent;
start.y >>= exponent;
numPixels.x >>= exponent;
numPixels.y >>= exponent;
}
};
class TileDataset {
public:
@@ -105,6 +135,11 @@ namespace openspace {
const static glm::ivec2 tilePixelStartOffset;
const static glm::ivec2 tilePixelSizeDifference;
const static PixelRegion padding;
static PixelCoordinate geodeticToPixel(GDALDataset* dataSet, const Geodetic2& geo);
private:
@@ -114,17 +149,22 @@ namespace openspace {
// HELPER FUNCTIONS //
//////////////////////////////////////////////////////////////////////////////////
const int calculateMaxLevel(int tileLevelDifference);
PixelRegion gdalPixelRegion(const GeodeticPatch& geodeticPatch) const;
int gdalOverview(const PixelCoordinate& baseRegionSize) const;
int gdalOverview(const ChunkIndex& chunkIndex) const;
int calculateMaxLevel(int tileLevelDifference);
TileDepthTransform calculateTileDepthTransform();
std::shared_ptr<TilePreprocessData> preprocess(const char* imageData,
const GdalDataRegion& region, const TileDataLayout& dataLayout);
const PixelRegion& region, const TileDataLayout& dataLayout);
static int calculateTileLevelDifference(GDALDataset* dataset, int minimumPixelSize);
static char* getImageDataFlippedY(const GdalDataRegion& region,
static char* getImageDataFlippedY(const PixelRegion& region,
const TileDataLayout& dataLayout, const char* imageData);