| 628 | */ |
| 629 | |
| 630 | CPLErr GDALGridMovingAverage(const void *poOptionsIn, GUInt32 nPoints, |
| 631 | const double *padfX, const double *padfY, |
| 632 | const double *padfZ, double dfXPoint, |
| 633 | double dfYPoint, double *pdfValue, |
| 634 | CPL_UNUSED void *hExtraParamsIn) |
| 635 | { |
| 636 | // TODO: For optimization purposes pre-computed parameters should be moved |
| 637 | // out of this routine to the calling function. |
| 638 | |
| 639 | const GDALGridMovingAverageOptions *const poOptions = |
| 640 | static_cast<const GDALGridMovingAverageOptions *>(poOptionsIn); |
| 641 | // Pre-compute search ellipse parameters. |
| 642 | const double dfRadius1Square = poOptions->dfRadius1 * poOptions->dfRadius1; |
| 643 | const double dfRadius2Square = poOptions->dfRadius2 * poOptions->dfRadius2; |
| 644 | const double dfSearchRadius = |
| 645 | std::max(poOptions->dfRadius1, poOptions->dfRadius2); |
| 646 | const double dfR12Square = dfRadius1Square * dfRadius2Square; |
| 647 | |
| 648 | GDALGridExtraParameters *psExtraParams = |
| 649 | static_cast<GDALGridExtraParameters *>(hExtraParamsIn); |
| 650 | const CPLQuadTree *phQuadTree = psExtraParams->hQuadTree; |
| 651 | |
| 652 | // Compute coefficients for coordinate system rotation. |
| 653 | const double dfAngle = TO_RADIANS * poOptions->dfAngle; |
| 654 | const bool bRotated = dfAngle != 0.0; |
| 655 | |
| 656 | const double dfCoeff1 = bRotated ? cos(dfAngle) : 0.0; |
| 657 | const double dfCoeff2 = bRotated ? sin(dfAngle) : 0.0; |
| 658 | |
| 659 | double dfAccumulator = 0.0; |
| 660 | |
| 661 | GUInt32 n = 0; // Used after for. |
| 662 | if (phQuadTree != nullptr) |
| 663 | { |
| 664 | CPLRectObj sAoi; |
| 665 | sAoi.minx = dfXPoint - dfSearchRadius; |
| 666 | sAoi.miny = dfYPoint - dfSearchRadius; |
| 667 | sAoi.maxx = dfXPoint + dfSearchRadius; |
| 668 | sAoi.maxy = dfYPoint + dfSearchRadius; |
| 669 | int nFeatureCount = 0; |
| 670 | GDALGridPoint **papsPoints = reinterpret_cast<GDALGridPoint **>( |
| 671 | CPLQuadTreeSearch(phQuadTree, &sAoi, &nFeatureCount)); |
| 672 | if (nFeatureCount != 0) |
| 673 | { |
| 674 | for (int k = 0; k < nFeatureCount; k++) |
| 675 | { |
| 676 | const int i = papsPoints[k]->i; |
| 677 | double dfRX = padfX[i] - dfXPoint; |
| 678 | double dfRY = padfY[i] - dfYPoint; |
| 679 | |
| 680 | if (bRotated) |
| 681 | { |
| 682 | const double dfRXRotated = |
| 683 | dfRX * dfCoeff1 + dfRY * dfCoeff2; |
| 684 | const double dfRYRotated = |
| 685 | dfRY * dfCoeff1 - dfRX * dfCoeff2; |
| 686 | dfRX = dfRXRotated; |
| 687 | dfRY = dfRYRotated; |
nothing calls this directly
no test coverage detected