diff --git a/NEWS.md b/NEWS.md index a1a5e445..c0696fd1 100755 --- a/NEWS.md +++ b/NEWS.md @@ -7,6 +7,7 @@ Fix: regression of the `stars` package makes `rasterize_terrain()` extremely slo New: `catalog_intersects()` support a `SpatExtent` Fix: `lidR` can fully works without `raster` ans `sp` Fix: [#742](https://github.com/r-lidar/lidR/discussions/742) +Fix: local maximum filter `lmf()` with a fixed windows is now 20 times faster. ## lidR v4.1.0 (Release date: 2024-01-31) diff --git a/src/LAS.cpp b/src/LAS.cpp index 4dbd6b8f..cca57515 100644 --- a/src/LAS.cpp +++ b/src/LAS.cpp @@ -358,12 +358,20 @@ NumericVector LAS::compute_range(DataFrame flightlines) void LAS::filter_local_maxima(NumericVector ws, double min_height, bool circular) { + char UKN = 0; + char NLM = 1; + char LMX = 2; + bool abort = false; bool vws = ws.length() > 1; SpatialIndex tree(las); Progress pb(npoints, "Local maximum filter: "); + std::vector state(npoints); + std::fill(state.begin(), state.end(), 0); + for (int i = 0 ; i < npoints ; i++) { if (Z[i] < min_height) state[i] = NLM; } + #pragma omp parallel for num_threads(ncpu) for (unsigned int i = 0 ; i < npoints ; i++) { @@ -373,8 +381,7 @@ void LAS::filter_local_maxima(NumericVector ws, double min_height, bool circular double hws = (vws) ? ws[i]/2 : ws[0]/2; - if (Z[i] < min_height) - continue; + if (state[i] == NLM) continue; // Get the points within a windows centered on the current point std::vector pts; @@ -403,7 +410,13 @@ void LAS::filter_local_maxima(NumericVector ws, double min_height, bool circular if(z > zmax) { is_lm = false; - break; + } + + // If the windows size is fixed, then, if a point is below the central point + // we are sure it is not a LM + if (!vws && z < Z[i]) + { + state[pt.id] = NLM; } // Found one equal. If this one was already tagged LM we can't have two lm @@ -456,6 +469,10 @@ void LAS::filter_local_maxima(NumericVector ws) else Rcpp::stop("C++ unexpected internal error in 'filter_local_maxima': invalid windows."); // # nocov + std::vector state(npoints); + std::fill(state.begin(), state.end(), 0); + for (int i = 0 ; i < npoints ; i++) { if (Z[i] < min_height) state[i] = NLM; } + SpatialIndex tree(las, skip); Progress pb(npoints, "Local maximum filter: ");