| 3531 | */ |
| 3532 | |
| 3533 | CPLErr GDALGridContextProcess(GDALGridContext *psContext, double dfXMin, |
| 3534 | double dfXMax, double dfYMin, double dfYMax, |
| 3535 | GUInt32 nXSize, GUInt32 nYSize, |
| 3536 | GDALDataType eType, void *pData, |
| 3537 | GDALProgressFunc pfnProgress, void *pProgressArg) |
| 3538 | { |
| 3539 | CPLAssert(psContext); |
| 3540 | CPLAssert(pData); |
| 3541 | |
| 3542 | if (nXSize == 0 || nYSize == 0) |
| 3543 | { |
| 3544 | CPLError(CE_Failure, CPLE_IllegalArg, |
| 3545 | "Output raster dimensions should have non-zero size."); |
| 3546 | return CE_Failure; |
| 3547 | } |
| 3548 | |
| 3549 | const double dfDeltaX = (dfXMax - dfXMin) / nXSize; |
| 3550 | const double dfDeltaY = (dfYMax - dfYMin) / nYSize; |
| 3551 | |
| 3552 | // For linear, check if we will need to fallback to nearest neighbour |
| 3553 | // by sampling along the edges. If all points on edges are within |
| 3554 | // triangles, then interior points will also be. |
| 3555 | if (psContext->eAlgorithm == GGA_Linear && |
| 3556 | psContext->sExtraParameters.hQuadTree == nullptr) |
| 3557 | { |
| 3558 | bool bNeedNearest = false; |
| 3559 | int nStartLeft = 0; |
| 3560 | int nStartRight = 0; |
| 3561 | const double dfXPointMin = dfXMin + (0 + 0.5) * dfDeltaX; |
| 3562 | const double dfXPointMax = dfXMin + (nXSize - 1 + 0.5) * dfDeltaX; |
| 3563 | for (GUInt32 nYPoint = 0; !bNeedNearest && nYPoint < nYSize; nYPoint++) |
| 3564 | { |
| 3565 | const double dfYPoint = dfYMin + (nYPoint + 0.5) * dfDeltaY; |
| 3566 | |
| 3567 | if (!GDALTriangulationFindFacetDirected( |
| 3568 | psContext->sExtraParameters.psTriangulation, nStartLeft, |
| 3569 | dfXPointMin, dfYPoint, &nStartLeft)) |
| 3570 | { |
| 3571 | bNeedNearest = true; |
| 3572 | } |
| 3573 | if (!GDALTriangulationFindFacetDirected( |
| 3574 | psContext->sExtraParameters.psTriangulation, nStartRight, |
| 3575 | dfXPointMax, dfYPoint, &nStartRight)) |
| 3576 | { |
| 3577 | bNeedNearest = true; |
| 3578 | } |
| 3579 | } |
| 3580 | int nStartTop = 0; |
| 3581 | int nStartBottom = 0; |
| 3582 | const double dfYPointMin = dfYMin + (0 + 0.5) * dfDeltaY; |
| 3583 | const double dfYPointMax = dfYMin + (nYSize - 1 + 0.5) * dfDeltaY; |
| 3584 | for (GUInt32 nXPoint = 1; !bNeedNearest && nXPoint + 1 < nXSize; |
| 3585 | nXPoint++) |
| 3586 | { |
| 3587 | const double dfXPoint = dfXMin + (nXPoint + 0.5) * dfDeltaX; |
| 3588 | |
| 3589 | if (!GDALTriangulationFindFacetDirected( |
| 3590 | psContext->sExtraParameters.psTriangulation, nStartTop, |
no test coverage detected