R/algorithm-dec.R

In lidR: Airborne LiDAR Data Manipulation and Visualization for Forestry Applications

#' Point Cloud Decimation Algorithm
#'
#' This function is made to be used in \link{decimate_points}. It implements an algorithm that
#' randomly removes points or pulses to reach the desired density over the whole area (see
#' \code{\link[=area]{area}}).
#'
#' @param density numeric. The desired output density.
#'
#' @param use_pulse logical. Decimate by removing random pulses instead of random points (requires running
#' \link{retrieve_pulses} first)
#'
#' @export
#'
#' @family point cloud decimation algorithms
#'
#' @examples
#' LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
#' las = readLAS(LASfile, select = "xyz")
#'
#' # Reach a pulse density of 1 on the overall dataset
#' thinned1 = decimate_points(las, random(1))
#' plot(rasterize_density(las))
#' plot(rasterize_density(thinned1))
#' @name sample_random
random = function(density, use_pulse = FALSE)
{
  assert_is_a_number(density)
  assert_all_are_positive(density)
  assert_is_a_bool(use_pulse)

  density   <- lazyeval::uq(density)
  use_pulse <- lazyeval::uq(use_pulse)

  f = function(las)
  {
    assert_is_valid_context(LIDRCONTEXTDEC, "random")

    if(use_pulse & !"pulseID" %in% names(las))
    {
      warning("No 'pulseID' attribute found. Decimation by points is used.")
      use_pulse <- FALSE
    }

    n <- round(density*area(las))

    if (use_pulse)
      return(.selected_pulses(las@data$pulseID, n))
    else
    {
      if (nrow(las@data) > n)
      {
        idx = sample(1:nrow(las@data), n)
        idx = sort(idx)
        return(idx)
      }
      else
        return(1:nrow(las@data))
    }
  }

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

#' Point Cloud Decimation Algorithm
#'
#' This function is made to be used in \link{decimate_points}. It implements an algorithm that
#' creates a grid with a given resolution and filters the point cloud by randomly selecting some
#' points in each cell. It is designed to produce point clouds that have uniform densities throughout
#' the coverage area. For each cell, the proportion of points or pulses that will be retained is computed
#' using the actual local density and the desired density. If the desired density is greater than the actual
#' density it returns an unchanged set of points (it cannot increase the density). The cell size must be
#' large enough to compute a coherent local density. For example, in a 2 points/m^2 point cloud, 25 square
#' meters would be feasible; however 1 square meter cells would not be feasible because density does
#' not have meaning at this scale.
#'
#' @param density numeric. The desired output density.
#'
#' @param res numeric. The resolution of the grid used to filter the point cloud
#'
#' @param use_pulse logical. Decimate by removing random pulses instead of random points (requires running
#' \link{retrieve_pulses} first)
#'
#' @export
#'
#' @family point cloud decimation algorithms
#'
#' @examples
#' LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
#' las = readLAS(LASfile, select = "xyz")
#'
#' # Select points randomly to reach an homogeneous density of 1
#' thinned <- decimate_points(las, homogenize(1,5))
#' plot(rasterize_density(thinned, 10))
#' @name sample_homogenize
homogenize = function(density, res = 5, use_pulse = FALSE)
{
  assert_is_a_number(density)
  assert_all_are_positive(density)
  assert_is_a_bool(use_pulse)
  assert_is_a_number(res)
  assert_all_are_positive(res)

  density   <- lazyeval::uq(density)
  res       <- lazyeval::uq(res)
  use_pulse <- lazyeval::uq(use_pulse)

  f = function(las)
  {
    assert_is_valid_context(LIDRCONTEXTDEC, "homogenize")

    if (use_pulse & !"pulseID" %in% names(las))
    {
      warning("No 'pulseID' attribute found. Decimation by points is used.")
      use_pulse <- FALSE
    }

    pulseID <- NULL

    n       <- round(density*res^2)
    layout  <- raster_layout(las, res)
    cells   <- get_group(layout, las)

    if (use_pulse)
      return(las@data[, .I[.selected_pulses(pulseID, n)], by = cells]$V1)
    else
      return(las@data[, .I[.selected_pulses(1:.N, n)], by = cells]$V1)
  }

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

#' Point Cloud Decimation Algorithm
#'
#' These functions are made to be used in \link{decimate_points}. They implement algorithms that
#' create a grid with a given resolution and filters the point cloud by selecting the highest/lowest
#' point within each cell.
#'
#' @param res numeric. The resolution of the grid used to filter the point cloud
#'
#' @export
#'
#' @family point cloud decimation algorithms
#'
#' @examples
#' LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
#' las = readLAS(LASfile, select = "xyz")
#'
#' # Select the highest point within each cell of an overlayed grid
#' thinned = decimate_points(las, highest(4))
#' #plot(thinned)
#'
#' # Select the lowest point within each cell of an overlayed grid
#' thinned = decimate_points(las, lowest(4))
#' #plot(thinned)
#' @name sample_maxima
highest = function(res = 1)
{
  assert_is_a_number(res)
  assert_all_are_positive(res)

  res <- lazyeval::uq(res)

  f = function(las)
  {
    assert_is_valid_context(LIDRCONTEXTDEC, "highest")
    layout <- raster_layout(las, res)
    return(C_highest(las, layout))
  }

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

#' @family point cloud decimation algorithms
#' @export
#' @name sample_maxima
lowest = function(res = 1)
{
  assert_is_a_number(res)
  assert_all_are_positive(res)

  res <- lazyeval::uq(res)

  f = function(las)
  {
    assert_is_valid_context(LIDRCONTEXTDEC, "lowest")
    layout <- raster_layout(las, res)
    return(C_lowest(las, layout))
  }

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

#' Point Cloud Decimation Algorithm
#'
#' This functions is made to be used in \link{decimate_points}. It implements an algorithm that
#' creates a 3D grid with a given resolution and filters the point cloud by randomly selecting
#' n points within each voxel
#'
#' @param res numeric. The resolution of the voxel grid used to filter the point cloud
#' @param n integer. The number of points to select
#'
#' @examples
#' LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
#' las <- readLAS(LASfile, select = "xyz")
#' thinned <- decimate_points(las, random_per_voxel(8, 1))
#' #plot(thinned)
#' @family point cloud decimation algorithms
#' @export
#' @name sample_per_voxel
random_per_voxel = function(res = 1, n = 1)
{
  assert_all_are_positive(n)
  assert_all_are_positive(res)
  n <- as.integer(n)
  if (length(res) == 1) res <- c(res, res)

  n <- lazyeval::uq(n)
  res <- lazyeval::uq(res)

  f = function(las)
  {
    by <- group_grid_3d(las$X, las$Y, las$Z, res)
    return(las@data[, .I[.selected_pulses(1:.N, n)], by = by]$V1)
  }

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


.selected_pulses = function(pulseID, n)
{
  p <- unique(pulseID)

  if (n > length(p))
    return(rep(TRUE, length(pulseID)))

  selectedPulses <- sample(p, n)
  selectedPulses <- pulseID %in% selectedPulses
  selectedPulses <- sort(selectedPulses)

  return(selectedPulses)
}