* Inverse distance to a power with nearest neighbor search, with a per-quadrant * search logic. */
| 349 | * search logic. |
| 350 | */ |
| 351 | static CPLErr GDALGridInverseDistanceToAPowerNearestNeighborPerQuadrant( |
| 352 | const void *poOptionsIn, GUInt32 /*nPoints*/, const double *padfX, |
| 353 | const double *padfY, const double *padfZ, double dfXPoint, double dfYPoint, |
| 354 | double *pdfValue, void *hExtraParamsIn) |
| 355 | { |
| 356 | const GDALGridInverseDistanceToAPowerNearestNeighborOptions *const |
| 357 | poOptions = static_cast< |
| 358 | const GDALGridInverseDistanceToAPowerNearestNeighborOptions *>( |
| 359 | poOptionsIn); |
| 360 | const double dfRadius = poOptions->dfRadius; |
| 361 | const double dfSmoothing = poOptions->dfSmoothing; |
| 362 | const double dfSmoothing2 = dfSmoothing * dfSmoothing; |
| 363 | |
| 364 | const GUInt32 nMaxPoints = poOptions->nMaxPoints; |
| 365 | const GUInt32 nMinPointsPerQuadrant = poOptions->nMinPointsPerQuadrant; |
| 366 | const GUInt32 nMaxPointsPerQuadrant = poOptions->nMaxPointsPerQuadrant; |
| 367 | |
| 368 | GDALGridExtraParameters *psExtraParams = |
| 369 | static_cast<GDALGridExtraParameters *>(hExtraParamsIn); |
| 370 | const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree; |
| 371 | CPLAssert(phQuadTree); |
| 372 | |
| 373 | const double dfRPower2 = psExtraParams->dfRadiusPower2PreComp; |
| 374 | const double dfPowerDiv2 = psExtraParams->dfPowerDiv2PreComp; |
| 375 | std::multimap<double, double> oMapDistanceToZValuesPerQuadrant[4]; |
| 376 | |
| 377 | const double dfSearchRadius = dfRadius; |
| 378 | CPLRectObj sAoi; |
| 379 | sAoi.minx = dfXPoint - dfSearchRadius; |
| 380 | sAoi.miny = dfYPoint - dfSearchRadius; |
| 381 | sAoi.maxx = dfXPoint + dfSearchRadius; |
| 382 | sAoi.maxy = dfYPoint + dfSearchRadius; |
| 383 | int nFeatureCount = 0; |
| 384 | GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>( |
| 385 | CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount)); |
| 386 | if (nFeatureCount != 0) |
| 387 | { |
| 388 | for (int k = 0; k < nFeatureCount; k++) |
| 389 | { |
| 390 | const int i = papsPoints[k]->i; |
| 391 | const double dfRX = padfX[i] - dfXPoint; |
| 392 | const double dfRY = padfY[i] - dfYPoint; |
| 393 | |
| 394 | const double dfR2 = dfRX * dfRX + dfRY * dfRY; |
| 395 | // real distance + smoothing |
| 396 | const double dfRsmoothed2 = dfR2 + dfSmoothing2; |
| 397 | if (dfRsmoothed2 < 0.0000000000001) |
| 398 | { |
| 399 | *pdfValue = padfZ[i]; |
| 400 | CPLFree(papsPoints); |
| 401 | return CE_None; |
| 402 | } |
| 403 | // is point within real distance? |
| 404 | if (dfR2 <= dfRPower2) |
| 405 | { |
| 406 | const int iQuadrant = |
| 407 | ((dfRX >= 0) ? 1 : 0) | (((dfRY >= 0) ? 1 : 0) << 1); |
| 408 | oMapDistanceToZValuesPerQuadrant[iQuadrant].insert( |
nothing calls this directly
no test coverage detected