From fcdc674f95950b6fe7b424955ca50cb5610c1c76 Mon Sep 17 00:00:00 2001 From: Erik Broberg Date: Tue, 5 Jul 2016 13:11:38 -0400 Subject: [PATCH] Implement repeat wrap mode on pixel read regions before performing gdal RasterIO --- modules/globebrowsing/tile/tiledataset.cpp | 264 +++++++++++++-------- modules/globebrowsing/tile/tiledataset.h | 169 ++++++++++++- 2 files changed, 334 insertions(+), 99 deletions(-) diff --git a/modules/globebrowsing/tile/tiledataset.cpp b/modules/globebrowsing/tile/tiledataset.cpp index 1b32832ade..239f2601f3 100644 --- a/modules/globebrowsing/tile/tiledataset.cpp +++ b/modules/globebrowsing/tile/tiledataset.cpp @@ -51,6 +51,10 @@ namespace { namespace openspace { + std::ostream& operator<<(std::ostream& os, const openspace::PixelRegion& pr) { + return os << pr.start.x << ", " << pr.start.y << " with size " << pr.numPixels.x << ", " << pr.numPixels.y; + } + ////////////////////////////////////////////////////////////////////////////////// // Tile Data Layout // ////////////////////////////////////////////////////////////////////////////////// @@ -118,6 +122,9 @@ namespace openspace { delete _dataset; } + + + ////////////////////////////////////////////////////////////////////////////////// // Public interface // ////////////////////////////////////////////////////////////////////////////////// @@ -162,6 +169,9 @@ namespace openspace { } + + + ////////////////////////////////////////////////////////////////////////////////// // Initialization // ////////////////////////////////////////////////////////////////////////////////// @@ -196,10 +206,14 @@ namespace openspace { } + + + ////////////////////////////////////////////////////////////////////////////////// // GDAL helper methods // ////////////////////////////////////////////////////////////////////////////////// + bool TileDataset::gdalHasOverviews() const { return _dataset->GetRasterBand(1)->GetOverviewCount() > 0; } @@ -252,6 +266,9 @@ namespace openspace { } + + + ////////////////////////////////////////////////////////////////////////////////// // ReadTileData helper functions // ////////////////////////////////////////////////////////////////////////////////// @@ -300,11 +317,11 @@ namespace openspace { PixelRegion region = gdalPixelRegion(chunkIndex); // pixel region at overview zero region.downscalePow2(overview + 1); // pixel region at suitable overview - // Create an IORegion based on that overview pixel region + // Create an IORegion based on that overview pixel region IODescription io; io.read.overview = overview; io.read.region = region; - io.write.region = { PixelCoordinate(0, 0), region.numPixels }; + io.write.region = region; // Handle the case where the dataset does not have overviews if (!gdalHasOverviews()) { @@ -315,12 +332,19 @@ namespace openspace { // For correct sampling in height dataset, we need to pad the texture tile io.read.region.pad(padding); io.write.region.pad(padding); + io.write.region.start = PixelCoordinate(0, 0); // write region starts in origin + + // Doing this may cause invalid regions, i.e. having negative pixel coordinates // or being too large etc. For now, just clamp PixelRegion overviewRegion = gdalPixelRegion(gdalRasterBand(overview)); - //io.read.region.clampTo(overviewRegion); - //io.write.region.clampTo(overviewRegion); + + //if (chunkIndex.level < 3) { + // io.read.region.clampTo(overviewRegion); + // io.write.region.clampTo(overviewRegion); + //} + io.write.bytesPerLine = _dataLayout.bytesPerPixel * io.write.region.numPixels.x; io.write.totalNumBytes = io.write.bytesPerLine * io.write.region.numPixels.y; @@ -328,7 +352,7 @@ namespace openspace { return io; } - char* TileDataset::readImageData(const IODescription& io, CPLErr& worstError) const { + char* TileDataset::readImageData(IODescription& io, CPLErr& worstError) const { // allocate memory for the image char* imageData = new char[io.write.totalNumBytes]; @@ -340,24 +364,7 @@ namespace openspace { // for every raster (or data channel, i.e. R in RGB) char* dataDestination = imageData + (i * _dataLayout.bytesPerDatum); - // 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 - dataDestination += (io.write.totalNumBytes - io.write.bytesPerLine); - - CPLErr err = 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 - dataDestination, // 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 + 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); @@ -367,88 +374,157 @@ namespace openspace { return imageData; } - char* TileDataset::readImageData2(const IODescription& io, CPLErr& worstError) const { - std::vector imageDataChannels(_dataLayout.numRasters); - size_t numByterPerChannel = io.write.totalNumBytes / _dataLayout.numRasters; - // Read the data (each rasterband is a separate channel) - for (size_t i = 0; i < _dataLayout.numRasters; i++) { - imageDataChannels[i] = new char[numByterPerChannel]; + CPLErr TileDataset::repeatedRasterIO(GDALRasterBand* rasterBand, const IODescription& fullIO, char* dataDestination) const { + CPLErr worstError = CPLErr::CE_None; - GDALRasterBand* rasterBand = gdalRasterBand(io.read.overview, i + 1); + // 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. - CPLErr err = 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 - imageDataChannels[i], // 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 - 0, // Pixel spacing - 0); // Line spacing - - // CE_None = 0, CE_Debug = 1, CE_Warning = 2, CE_Failure = 3, CE_Fatal = 4 - worstError = std::max(worstError, err); - } - - // Combined image data - char* imageData = new char[io.write.totalNumBytes]; - size_t yx = 0; - size_t c = 0; - for (size_t y = 0; y < io.write.region.numPixels.y; y++) { - for (size_t x = 0; x < io.write.region.numPixels.x; x++) { - for (size_t c = 0; c < _dataLayout.numRasters; c++) { - size_t combinedChannelIndex = (yx * _dataLayout.numRasters + c) * _dataLayout.bytesPerDatum; - size_t separateChannelIndex = yx * _dataLayout.bytesPerDatum; - - ghoul_assert(combinedChannelIndex < io.write.totalNumBytes, "Invalid combined index!"); - ghoul_assert(separateChannelIndex < numByterPerChannel, "invalid single index!"); - - char* channelData = imageDataChannels[c]; - for (size_t b = 0; b < _dataLayout.bytesPerDatum; b++) { - char val = channelData[separateChannelIndex + b]; - imageData[combinedChannelIndex + b] = channelData[separateChannelIndex+b]; - } - } - yx++; - } - } - - for (size_t c = 0; c < _dataLayout.numRasters; c++) { - char * singleChannel = imageDataChannels[c]; - delete[] singleChannel; - } - - // GDAL reads pixel lines top to bottom, we want the opposit - return flipImageYAxis(imageData, io.write); - } + // Make a copy of the full IO desription as we will have to modify it + IODescription io = fullIO; + PixelRegion gdalRegion = gdalPixelRegion(rasterBand); - char* TileDataset::flipImageYAxis(char*& imageData, const IODescription::WriteData& writeData) const { - char* imageDataYflipped = new char[writeData.totalNumBytes]; - for (size_t y = 0; y < writeData.region.numPixels.y; y++) { - size_t yi_flipped = y * writeData.bytesPerLine; - size_t yi = (writeData.region.numPixels.y - 1 - y) * writeData.bytesPerLine; - size_t i = 0; - for (size_t x = 0; x < writeData.region.numPixels.x; x++) { - for (size_t c = 0; c < _dataLayout.numRasters; c++) { - for (size_t b = 0; b < _dataLayout.bytesPerDatum; b++) { - imageDataYflipped[yi_flipped + i] = imageData[yi + i]; - i++; - } + // 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("-"); + LDEBUG("repeated read: " << io.read.region); + LDEBUG("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); + + worstError = std::max(worstError, err); } } } - // Delete the old data and return the new - delete[] imageData; - imageData = nullptr; + - return imageDataYflipped; + CPLErr err = rasterIO(rasterBand, io, dataDestination); + worstError = std::max(worstError, err); + return worstError; } + CPLErr TileDataset::rasterIO(GDALRasterBand* rasterBand, const IODescription& io, char* dataDestination) const { + PixelRegion gdalRegion = gdalPixelRegion(rasterBand); + + LDEBUG("rasterIO read: " << io.read.region); + LDEBUG("rasterIO write: " << io.write.region); + + 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 TileDataset::preprocess(std::shared_ptr result, const PixelRegion& region) const { size_t bytesPerLine = _dataLayout.bytesPerPixel * region.numPixels.x; size_t totalNumBytes = bytesPerLine * region.numPixels.y; diff --git a/modules/globebrowsing/tile/tiledataset.h b/modules/globebrowsing/tile/tiledataset.h index 4fa09b5cc7..347b176af1 100644 --- a/modules/globebrowsing/tile/tiledataset.h +++ b/modules/globebrowsing/tile/tiledataset.h @@ -62,8 +62,14 @@ namespace openspace { typedef glm::ivec2 PixelCoordinate; + + struct PixelRegion { + enum Side { + LEFT = 0, TOP, RIGHT, BOTTOM, + }; + PixelRegion() : start(0), numPixels(0) { } PixelRegion(const PixelRegion& o) : start(o.start), numPixels(o.numPixels) { } PixelRegion(const PixelCoordinate& pixelStart, const PixelCoordinate& numberOfPixels) @@ -83,6 +89,36 @@ namespace openspace { return start + numPixels; } + void resizeToStartAt(const PixelCoordinate& c) { + numPixels -= (c - start); + start = c; + } + + void setLeft(int x) { + numPixels.x += (start.x - x); + start.x = x; + } + + void setTop(int p) { + numPixels.y += (start.y - p); + start.y = p; + } + + void setRight(int x) { + numPixels.x = x - start.x; + } + + void setBottom(int y) { + numPixels.y = y - start.y; + } + + + + void alignLeft(int x) { start.x = x; } + void alignTop(int y) { start.y = y; } + void alignRight(int x) { start.x = x - numPixels.x; } + void alignBottom(int y) { start.y = y - numPixels.y; } + void scale(const glm::dvec2& s) { start = PixelCoordinate(glm::round(s * glm::dvec2(start))); numPixels = PixelCoordinate(glm::round(s * glm::dvec2(numPixels))); @@ -92,6 +128,10 @@ namespace openspace { scale(glm::dvec2(s)); } + int area() { + return numPixels.x * numPixels.y; + } + void downscalePow2(int exponent) { start.x >>= exponent; start.y >>= exponent; @@ -106,12 +146,133 @@ namespace openspace { numPixels.y <<= exponent; } + int edge(Side side) const { + switch (side) { + case LEFT: return start.x; + case TOP: return start.y; + case RIGHT: return start.x + numPixels.x; + case BOTTOM: return start.y + numPixels.y; + } + } + + int edgeDirectionSign(Side side) const { + return side < RIGHT ? -1 : 1; + } + + void align(Side side, int pos) { + switch (side) { + case LEFT: alignLeft(pos); break; + case TOP: alignTop(pos); break; + case RIGHT: alignRight(pos); break; + case BOTTOM: alignBottom(pos); break; + } + } + + void move(Side side, int amount) { + switch (side) { + case LEFT: start.x -= amount; break; + case TOP: start.y -= amount; break; + case RIGHT: start.x += amount; break; + case BOTTOM: start.y += amount; break; + } + } + + void setSide(Side side, int pos) { + switch (side) { + case LEFT: setLeft(pos); break; + case TOP: setTop(pos); break; + case RIGHT: setRight(pos); break; + case BOTTOM: setBottom(pos); break; + } + } + + bool lineIntersect(Side side, int p) { + switch (side) + { + case openspace::PixelRegion::LEFT: + case openspace::PixelRegion::RIGHT: + return start.x <= p && p <= (start.x + numPixels.x); + + case openspace::PixelRegion::TOP: + case openspace::PixelRegion::BOTTOM: + return start.y <= p && p <= (start.y + numPixels.y); + } + } + + bool isInside(const PixelRegion& r) const { + PixelCoordinate e = end(); + PixelCoordinate re = r.end(); + return r.start.x <= start.x && e.x <= re.x + && r.start.y <= start.y && e.y <= re.y; + } + + PixelRegion cut(Side side, int p) { + if (!lineIntersect(side, p)) { + return PixelRegion({ 0, 0 }, {0, 0}); + } + + PixelRegion cutOff(*this); + int cutSize = 0; + switch (side) { + case LEFT: + setLeft(p); + cutOff.setRight(p - cutSize); + break; + case TOP: + setTop(p); + cutOff.setBottom(p - cutSize); + break; + case RIGHT: + setRight(p); + cutOff.setLeft(p + cutSize); + break; + case BOTTOM: + setBottom(p); + cutOff.setTop(p + cutSize); + break; + } + return cutOff; + } + + PixelRegion cutDelta(Side side, int p) { + if (p < 1) { + return PixelRegion({ 0, 0 }, { 0, 0 }); + } + else return cut(side, edge(side) - edgeDirectionSign(side) * p); + } + + bool equals(const PixelRegion& r) const { + return start == r.start && numPixels == r.numPixels; + } + PixelCoordinate start; PixelCoordinate numPixels; }; + + struct IODescription { + IODescription cut(PixelRegion::Side side, int pos) { + PixelRegion readPreCut = read.region; + PixelRegion writePreCut = write.region; + + IODescription whatCameOff = *this; + whatCameOff.read.region = read.region.cut(side, pos); + + PixelCoordinate cutSize = whatCameOff.read.region.numPixels; + int cutWidth = (side % 2 == 0) ? cutSize.x : cutSize.y; + whatCameOff.write.region = write.region.cutDelta(side, cutWidth); + + + if (cutWidth == 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"); + } + + return whatCameOff; + } + struct ReadData { int overview; PixelRegion region; @@ -184,9 +345,9 @@ namespace openspace { PixelCoordinate geodeticToPixel(const Geodetic2& geo) const; IODescription getIODescription(const ChunkIndex& chunkIndex) const; - char* readImageData(const IODescription& io, CPLErr& worstError) const; - char* readImageData2(const IODescription& io, CPLErr& worstError) const; - char* flipImageYAxis(char*& imageData, const IODescription::WriteData& writeData) const; + char* readImageData(IODescription& io, CPLErr& worstError) const; + CPLErr rasterIO(GDALRasterBand* rasterBand, const IODescription& io, char* dst) const; + CPLErr repeatedRasterIO(GDALRasterBand* rasterBand, const IODescription& io, char* dst) const; std::shared_ptr preprocess(std::shared_ptr result, const PixelRegion& region) const; CPLErr postProcessErrorCheck(std::shared_ptr ioResult, const IODescription& io) const; @@ -217,6 +378,4 @@ namespace openspace { - - #endif // __TILE_DATASET_H__ \ No newline at end of file