Feature Request - Voxel decimation with user selection of retained points

[liamirwin]

Hi JR,

It would be very useful to have the option to decimate a point cloud using voxels but with the added functionality where the user is able to select which point(s) are retained rather than retaining a random point for each voxel like is currently available. This is especially useful when trying to account for differences in point density across multitemporal acquisitions (see below).

Selecting the highest/lowest point, or mean Z point would be great. I am unsure how much more work it would be to code in the flexibility for the user to define which point they are keeping.

Here is my current approach, it is quite slow and relies on LAStools to attribute a voxel ID to each point

Step 1: compute and attribute voxel IDs to each point in the point cloud (using LAStools as I couldn't find a method with lidR)

Link to LAStools output LAZ file: https://drive.google.com/file/d/1re_PDVriFcvcPolURqmDhV1ppHMuabFC/view?usp=sharing

lasvoxel -i G:/Misc/voxel_request_lidR/Megaplot.laz ^
-step_xy 1 ^
-step_z 1 ^
-compute_IDs_and_voxel_table ^
-olaz ^
-odir G:/Misc/voxel_request_lidR ^
-odix 1m_vox

Step 2: Load voxel attributed LAS as data.table and retain desired points (highest n points)

library(lidR)
library(data.table)

#Number of points to retain
n <- 1
# Load LAS with voxel ID attributed to each point
las <- readLAS('G:/Misc/voxel_request_lidR/Megaplot1m_vox.laz')
# Load as data.table
las_df <- las@data
# Rename Voxel ID attribute from LAStools
setnames(las_df, old = 'voxel ID', new = 'voxel_id')
# Slice n highest points per voxel group
las_filt <- las_df[order(-Z), head(.SD, n), by = voxel_id]
# Populate a new LAS object with filtered data
filtered_las <- LAS(data = las_filt,
                    header = LASheader(las_filt),
                    crs = st_crs(las))

[Jean-Romain]

This is to keep 1 point per voxel. For n points per voxel your approach is good, sorting is required. You could mix your approach with mine.

keep_x_per_voxel = function(res = 1, x = "max")
{
  x = match.arg(x, c("min", "max", "mean"))
  x <- lazyeval::uq(x)
  res <- lazyeval::uq(res)

  f = function(las)
  {
    by <- lidR:::group_grid_3d(las$X, las$Y, las$Z, res)
    if (x == "max") return(las@data[, .I[which.max(Z)], by = by]$V1)
    if (x == "min") return(las@data[, .I[which.min(Z)], by = by]$V1)
    if (x == "mean") return(las@data[, .I[which.min(abs(Z - mean(Z)))], by = by]$V1)
  }

  f <- plugin_decimate(f)
  return(f)
}

library(lidR)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile, select = "xyz")

thinned1 <- decimate_points(las, keep_x_per_voxel(2, "max"))

plot(thinned1)