Fill voids with the average of eight nearest neighbors.
| 601 | |
| 602 | // Fill voids with the average of eight nearest neighbors. |
| 603 | void SMRFilter::knnfill(PointViewPtr view, std::vector<double>& cz) |
| 604 | { |
| 605 | //ABELL - This potentially means moving a lot of data from the raster |
| 606 | // to the temporary view. This can be improved by either |
| 607 | // 1) using some method other than a KD tree to find neighbors |
| 608 | // 2) build a KDtree from the raster data directly, rather than moving it |
| 609 | // to a view. |
| 610 | |
| 611 | // Create a temporary PointView that encodes our raster values so that we |
| 612 | // can construct a 2D KDIndex and perform nearest neighbor searches. |
| 613 | PointViewPtr temp = view->makeNew(); |
| 614 | PointId i(0); |
| 615 | for (int c = 0; c < m_cols; ++c) |
| 616 | { |
| 617 | for (int r = 0; r < m_rows; ++r) |
| 618 | { |
| 619 | size_t cell = c * m_rows + r; |
| 620 | double val = cz[cell]; |
| 621 | if (std::isnan(val)) |
| 622 | continue; |
| 623 | |
| 624 | PointRef p = temp->point(i++); |
| 625 | p.setField(Id::X, m_bounds.minx + (c + 0.5) * m_args->m_cell); |
| 626 | p.setField(Id::Y, m_bounds.miny + (r + 0.5) * m_args->m_cell); |
| 627 | p.setField(Id::Z, val); |
| 628 | } |
| 629 | } |
| 630 | |
| 631 | // https://github.com/PDAL/PDAL/issues/2794#issuecomment-625297062 |
| 632 | if (!temp->size()) |
| 633 | return; |
| 634 | |
| 635 | KD2Index& kdi = temp->build2dIndex(); |
| 636 | |
| 637 | // Where the raster has voids (i.e., NaN), we search for that cell's eight |
| 638 | // nearest neighbors, and fill the void with the average value of the |
| 639 | // neighbors. |
| 640 | for (int c = 0; c < m_cols; ++c) |
| 641 | { |
| 642 | for (int r = 0; r < m_rows; ++r) |
| 643 | { |
| 644 | size_t cell = c * m_rows + r; |
| 645 | if (!std::isnan(cz[cell])) |
| 646 | continue; |
| 647 | |
| 648 | double x = m_bounds.minx + (c + 0.5) * m_args->m_cell; |
| 649 | double y = m_bounds.miny + (r + 0.5) * m_args->m_cell; |
| 650 | const int k = 8; |
| 651 | PointIdList neighbors = kdi.neighbors(x, y, k); |
| 652 | |
| 653 | double M1(0.0); |
| 654 | size_t j(0); |
| 655 | for (auto const& n : neighbors) |
| 656 | { |
| 657 | j++; |
| 658 | double delta = temp->getFieldAs<double>(Id::Z, n) - M1; |
| 659 | M1 += (delta / j); |
| 660 | } |