Skip to content

Commit

Permalink
Merge pull request #9957 from rouault/GDALRasterizeGeometries_doc
Browse files Browse the repository at this point in the history 
GDALRasterizeGeometries(): various robustness fixes
  • Loading branch information
@rouault
rouault committed May 22, 2024
2 parents 698d879 + 089e97d commit 3a3c9d3
Showing 1 changed file with 135 additions and 83 deletions.
218 changes: 135 additions & 83 deletions alg/gdalrasterize.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -800,29 +800,37 @@ static CPLErr GDALRasterizeOptions(CSLConstList papszOptions, int *pbAllTouched,
}

/* -------------------------------------------------------------------- */
/* OPTIM=[AUTO]/RASTER/VECTOR */
/* OPTIM=[AUTO]/RASTER/VECTOR */
/* -------------------------------------------------------------------- */
*peOptim = GRO_Auto;
pszOpt = CSLFetchNameValue(papszOptions, "OPTIM");
if (pszOpt)
{
if (EQUAL(pszOpt, "RASTER"))
{
*peOptim = GRO_Raster;
}
else if (EQUAL(pszOpt, "VECTOR"))
{
*peOptim = GRO_Vector;
}
else if (EQUAL(pszOpt, "AUTO"))
if (peOptim)
{
*peOptim = GRO_Auto;
if (EQUAL(pszOpt, "RASTER"))
{
*peOptim = GRO_Raster;
}
else if (EQUAL(pszOpt, "VECTOR"))
{
*peOptim = GRO_Vector;
}
else if (EQUAL(pszOpt, "AUTO"))
{
*peOptim = GRO_Auto;
}
else
{
CPLError(CE_Failure, CPLE_AppDefined,
"Unrecognized value '%s' for OPTIM.", pszOpt);
return CE_Failure;
}
}
else
{
CPLError(CE_Failure, CPLE_AppDefined,
"Unrecognized value '%s' for OPTIM.", pszOpt);
return CE_Failure;
CPLError(CE_Warning, CPLE_NotSupported,
"Option OPTIM is not supported by this function");
}
}

Expand Down Expand Up @@ -890,6 +898,13 @@ static CPLErr GDALRasterizeGeometriesInternal(
* used. Default size will be estimated based on the GDAL cache buffer size
* using formula: cache_size_bytes/scanline_size_bytes, so the chunk will
* not exceed the cache. Not used in OPTIM=RASTER mode.</li>
* <li>"OPTIM": May be set to "AUTO", "RASTER", "VECTOR". Force the algorithm
* used (results are identical). The raster mode is used in most cases and
* optimise read/write operations. The vector mode is useful with a decent
* amount of input features and optimize the CPU use. That mode has to be used
* with tiled images to be efficient. The auto mode (the default) will chose
* the algorithm based on input and output properties.
* </li>
* </ul>
* @param pfnProgress the progress function to report completion.
* @param pProgressArg callback data for progress function.
Expand Down Expand Up @@ -917,7 +932,7 @@ static CPLErr GDALRasterizeGeometriesInternal(
* // or create it as follows
* // GDALDriverH hMemDriver = GDALGetDriverByName("MEM");
* // GDALDatasetH hMemDset = GDALCreate(hMemDriver, "", nBufXSize,
*nBufYSize, nBandCount, eType, NULL);
* nBufYSize, nBandCount, eType, NULL);
*
* double adfGeoTransform[6];
* // Assign GeoTransform parameters,Omitted here.
Expand All @@ -930,9 +945,10 @@ static CPLErr GDALRasterizeGeometriesInternal(
*
* int bandList[3] = { 1, 2, 3};
* std::vector<double> geomBurnValue(nGeomCount*nBandCount,255.0);
* CPLErr err = GDALRasterizeGeometries(hMemDset, nBandCount, bandList,
* nGeomCount, pahGeoms, pfnTransformer,
*pTransformArg, geomBurnValue.data(), papszOptions, pfnProgress, pProgressArg);
* CPLErr err = GDALRasterizeGeometries(
* hMemDset, nBandCount, bandList, nGeomCount, pahGeoms, pfnTransformer,
* pTransformArg, geomBurnValue.data(), papszOptions,
* pfnProgress, pProgressArg);
* if( err != CE_None )
* {
* // Do something ...
Expand Down Expand Up @@ -1113,14 +1129,24 @@ static CPLErr GDALRasterizeGeometriesInternal(
const GDALDataType eType =
GDALGetNonComplexDataType(poBand->GetRasterDataType());

const int nScanlineBytes = nBandCount * poDS->GetRasterXSize() *
GDALGetDataTypeSizeBytes(eType);
const uint64_t nScanlineBytes = static_cast<uint64_t>(nBandCount) *
poDS->GetRasterXSize() *
GDALGetDataTypeSizeBytes(eType);

int nYChunkSize = 0;
const char *pszYChunkSize =
CSLFetchNameValue(papszOptions, "CHUNKYSIZE");
if (pszYChunkSize == nullptr ||
((nYChunkSize = atoi(pszYChunkSize))) == 0)
#if SIZEOF_VOIDP < 8
// Only on 32-bit systems and in pathological cases
if (nScanlineBytes > std::numeric_limits<size_t>::max())
{
CPLError(CE_Failure, CPLE_OutOfMemory, "Too big raster");
if (bNeedToFreeTransformer)
GDALDestroyTransformer(pTransformArg);
return CE_Failure;
}
#endif

int nYChunkSize =
atoi(CSLFetchNameValueDef(papszOptions, "CHUNKYSIZE", "0"));
if (nYChunkSize <= 0)
{
const GIntBig nYChunkSize64 = GDALGetCacheMax64() / nScanlineBytes;
const int knIntMax = std::numeric_limits<int>::max();
Expand All @@ -1138,8 +1164,8 @@ static CPLErr GDALRasterizeGeometriesInternal(
(poDS->GetRasterYSize() + nYChunkSize - 1) / nYChunkSize,
nYChunkSize);

pabyChunkBuf = static_cast<unsigned char *>(
VSI_MALLOC2_VERBOSE(nYChunkSize, nScanlineBytes));
pabyChunkBuf = static_cast<unsigned char *>(VSI_MALLOC2_VERBOSE(
nYChunkSize, static_cast<size_t>(nScanlineBytes)));
if (pabyChunkBuf == nullptr)
{
if (bNeedToFreeTransformer)
Expand Down Expand Up @@ -1176,10 +1202,13 @@ static CPLErr GDALRasterizeGeometriesInternal(
OGRGeometry::FromHandle(pahGeometries[iShape]),
eBurnValueType,
padfGeomBurnValues
? padfGeomBurnValues + iShape * nBandCount
? padfGeomBurnValues +
static_cast<size_t>(iShape) * nBandCount
: nullptr,
panGeomBurnValues
? panGeomBurnValues +
static_cast<size_t>(iShape) * nBandCount
: nullptr,
panGeomBurnValues ? panGeomBurnValues + iShape * nBandCount
: nullptr,
eBurnValueSource, eMergeAlg, pfnTransformer, pTransformArg);
}

Expand Down Expand Up @@ -1228,26 +1257,41 @@ static CPLErr GDALRasterizeGeometriesInternal(
std::min(static_cast<GIntBig>(knIntMax / nPixelSize / nYBlockSize /
nXBlockSize),
nbMaxBlocks64));
const int nbBlocsX = std::max(
const int nbBlocksX = std::max(
1,
std::min(static_cast<int>(sqrt(static_cast<double>(nbMaxBlocks))),
nXBlocks));
const int nbBlocsY =
std::max(1, std::min(nbMaxBlocks / nbBlocsX, nYBlocks));
const int nbBlocksY =
std::max(1, std::min(nbMaxBlocks / nbBlocksX, nYBlocks));

const uint64_t nChunkSize = static_cast<uint64_t>(nXBlockSize) *
nbBlocksX * nYBlockSize * nbBlocksY;

const int nScanblocks = nXBlockSize * nbBlocsX * nYBlockSize * nbBlocsY;
#if SIZEOF_VOIDP < 8
// Only on 32-bit systems and in pathological cases
if (nChunkSize > std::numeric_limits<size_t>::max())
{
CPLError(CE_Failure, CPLE_OutOfMemory, "Too big raster");
if (bNeedToFreeTransformer)
GDALDestroyTransformer(pTransformArg);
return CE_Failure;
}
#endif

pabyChunkBuf = static_cast<unsigned char *>(
VSI_MALLOC2_VERBOSE(nPixelSize, nScanblocks));
VSI_MALLOC2_VERBOSE(nPixelSize, static_cast<size_t>(nChunkSize)));
if (pabyChunkBuf == nullptr)
{
if (bNeedToFreeTransformer)
GDALDestroyTransformer(pTransformArg);
return CE_Failure;
}

int *panSuccessTransform =
static_cast<int *>(CPLCalloc(sizeof(int), 2));
OGREnvelope sRasterEnvelope;
sRasterEnvelope.MinX = 0;
sRasterEnvelope.MinY = 0;
sRasterEnvelope.MaxX = poDS->GetRasterXSize();
sRasterEnvelope.MaxY = poDS->GetRasterYSize();

/* --------------------------------------------------------------------
*/
Expand All @@ -1262,56 +1306,66 @@ static CPLErr GDALRasterizeGeometriesInternal(
OGRGeometry::FromHandle(pahGeometries[iShape]);
if (poGeometry == nullptr || poGeometry->IsEmpty())
continue;
/* --------------------------------------------------------------------
*/
/* get the envelope of the geometry and transform it to pixels
* coo */
/* --------------------------------------------------------------------
*/
OGREnvelope psGeomEnvelope;
poGeometry->getEnvelope(&psGeomEnvelope);
/* ------------------------------------------------------------ */
/* get the envelope of the geometry and transform it to */
/* pixels coordinates. */
/* ------------------------------------------------------------ */
OGREnvelope sGeomEnvelope;
poGeometry->getEnvelope(&sGeomEnvelope);
if (pfnTransformer != nullptr)
{
int anSuccessTransform[2] = {0};
double apCorners[4];
apCorners[0] = psGeomEnvelope.MinX;
apCorners[1] = psGeomEnvelope.MaxX;
apCorners[2] = psGeomEnvelope.MinY;
apCorners[3] = psGeomEnvelope.MaxY;
// TODO: need to add all appropriate error checking
pfnTransformer(pTransformArg, FALSE, 2, &(apCorners[0]),
&(apCorners[2]), nullptr, panSuccessTransform);
psGeomEnvelope.MinX = std::min(apCorners[0], apCorners[1]);
psGeomEnvelope.MaxX = std::max(apCorners[0], apCorners[1]);
psGeomEnvelope.MinY = std::min(apCorners[2], apCorners[3]);
psGeomEnvelope.MaxY = std::max(apCorners[2], apCorners[3]);
apCorners[0] = sGeomEnvelope.MinX;
apCorners[1] = sGeomEnvelope.MaxX;
apCorners[2] = sGeomEnvelope.MinY;
apCorners[3] = sGeomEnvelope.MaxY;

if (!pfnTransformer(pTransformArg, FALSE, 2, &(apCorners[0]),
&(apCorners[2]), nullptr,
anSuccessTransform) ||
!anSuccessTransform[0] || !anSuccessTransform[1])
{
continue;
}
sGeomEnvelope.MinX = std::min(apCorners[0], apCorners[1]);
sGeomEnvelope.MaxX = std::max(apCorners[0], apCorners[1]);
sGeomEnvelope.MinY = std::min(apCorners[2], apCorners[3]);
sGeomEnvelope.MaxY = std::max(apCorners[2], apCorners[3]);
}

int minBlockX = std::max(0, int(psGeomEnvelope.MinX) / nXBlockSize);
int minBlockY = std::max(0, int(psGeomEnvelope.MinY) / nYBlockSize);
int maxBlockX = std::min(
nXBlocks - 1, int(psGeomEnvelope.MaxX + 1) / nXBlockSize);
int maxBlockY = std::min(
nYBlocks - 1, int(psGeomEnvelope.MaxY + 1) / nYBlockSize);

/* --------------------------------------------------------------------
*/
/* loop over the blocks concerned by the geometry */
/* (by packs of nbBlocsX x nbBlocsY) */
/* --------------------------------------------------------------------
*/

for (int xB = minBlockX; xB <= maxBlockX; xB += nbBlocsX)
if (!sGeomEnvelope.Intersects(sRasterEnvelope))
continue;
sGeomEnvelope.Intersect(sRasterEnvelope);
CPLAssert(sGeomEnvelope.MinX >= 0 &&
sGeomEnvelope.MinX <= poDS->GetRasterXSize());
CPLAssert(sGeomEnvelope.MinY >= 0 &&
sGeomEnvelope.MinY <= poDS->GetRasterYSize());
CPLAssert(sGeomEnvelope.MaxX >= 0 &&
sGeomEnvelope.MaxX <= poDS->GetRasterXSize());
CPLAssert(sGeomEnvelope.MaxY >= 0 &&
sGeomEnvelope.MaxY <= poDS->GetRasterYSize());
const int minBlockX = int(sGeomEnvelope.MinX) / nXBlockSize;
const int minBlockY = int(sGeomEnvelope.MinY) / nYBlockSize;
const int maxBlockX = int(sGeomEnvelope.MaxX + 1) / nXBlockSize;
const int maxBlockY = int(sGeomEnvelope.MaxY + 1) / nYBlockSize;

/* ------------------------------------------------------------ */
/* loop over the blocks concerned by the geometry */
/* (by packs of nbBlocksX x nbBlocksY) */
/* ------------------------------------------------------------ */

for (int xB = minBlockX; xB <= maxBlockX; xB += nbBlocksX)
{
for (int yB = minBlockY; yB <= maxBlockY; yB += nbBlocsY)
for (int yB = minBlockY; yB <= maxBlockY; yB += nbBlocksY)
{

/* --------------------------------------------------------------------
*/
/* ensure to stay in the image */
/* --------------------------------------------------------------------
*/
int remSBX = std::min(maxBlockX - xB + 1, nbBlocsX);
int remSBY = std::min(maxBlockY - yB + 1, nbBlocsY);
int remSBX = std::min(maxBlockX - xB + 1, nbBlocksX);
int remSBY = std::min(maxBlockY - yB + 1, nbBlocksY);
int nThisXChunkSize = nXBlockSize * remSBX;
int nThisYChunkSize = nYBlockSize * remSBY;
if (xB * nXBlockSize + nThisXChunkSize >
Expand Down Expand Up @@ -1343,10 +1397,12 @@ static CPLErr GDALRasterizeGeometriesInternal(
OGRGeometry::FromHandle(pahGeometries[iShape]),
eBurnValueType,
padfGeomBurnValues
? padfGeomBurnValues + iShape * nBandCount
? padfGeomBurnValues +
static_cast<size_t>(iShape) * nBandCount
: nullptr,
panGeomBurnValues
? panGeomBurnValues + iShape * nBandCount
? panGeomBurnValues +
static_cast<size_t>(iShape) * nBandCount
: nullptr,
eBurnValueSource, eMergeAlg, pfnTransformer,
pTransformArg);
Expand All @@ -1369,8 +1425,6 @@ static CPLErr GDALRasterizeGeometriesInternal(
}
}

CPLFree(panSuccessTransform);

if (!pfnProgress(1., "", pProgressArg))
{
CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
Expand Down Expand Up @@ -1483,9 +1537,8 @@ CPLErr GDALRasterizeLayers(GDALDatasetH hDS, int nBandCount, int *panBandList,
int bAllTouched = FALSE;
GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
GDALRasterizeOptim eOptim = GRO_Auto;
if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
&eMergeAlg, &eOptim) == CE_Failure)
&eMergeAlg, nullptr) == CE_Failure)
{
return CE_Failure;
}
Expand Down Expand Up @@ -1890,9 +1943,8 @@ CPLErr GDALRasterizeLayersBuf(
int bAllTouched = FALSE;
GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
GDALRasterMergeAlg eMergeAlg = GRMA_Replace;
GDALRasterizeOptim eOptim = GRO_Auto;
if (GDALRasterizeOptions(papszOptions, &bAllTouched, &eBurnValueSource,
&eMergeAlg, &eOptim) == CE_Failure)
&eMergeAlg, nullptr) == CE_Failure)
{
return CE_Failure;
}
Expand Down

0 comments on commit 3a3c9d3

Please sign in to comment.