R/bounding_box.R

In lacunr: Fast 3D Lacunarity for Voxel Data

#' Create voxel array
#'
#' `bounding_box()` takes a table of voxel data and converts it into a
#' 3-dimensional array, with the original voxels arranged in their correct
#' spatial positions inside of a 3-D "box" of empty voxels. This array can be
#' input directly into [`lacunarity()`] to generate lacunarity curves.
#'
#' @details `bounding_box()` relies on the spatial coordinates of the input data
#'   to determine the dimensions of the resulting array. Noisy point cloud data
#'   will often produce "outlier" voxels surrounding points that are far removed
#'   from the bulk of the point cloud. These can drastically alter the output,
#'   creating an array in which the occupied voxels are surrounded by a large
#'   region of empty space. It is highly recommended that users supply a cutoff
#'   value for `threshold` to ideally remove these outliers. More elaborate
#'   filtering tools for trimming the point data before voxelization are
#'   available in the [`lidR`][lidR::lidR] package.
#'
#'
#' @param x A '`lac_voxels`' object created by [`voxelize()`] (preferred), or a
#'   '`lasmetrics3d`' object created by [`lidR::voxel_metrics()`]. Users can
#'   alternatively supply a [`data.table`][`data.table()`] containing X, Y, Z,
#'   and N columns, in which case the argument `edge_length` is required.
#' @param edge_length a numeric vector of length 3, specifying the X, Y, and Z
#'   dimensions of each voxel. This argument should only be necessary when
#'   supplying voxel data generated by a function other than [`voxelize()`] or
#'   [lidR::voxel_metrics()].
#' @param threshold An integer specifying the minimum number of points to use
#'   when determining if a voxel is occupied. The default is `0`.
#'   `bounding_box()` retains only those voxels where `x$N` is greater than
#'   `threshold`.
#'
#' @return A 3-dimensional integer [`array`][`array()`] containing values of `0`
#'   or `1`, representing the occupancy of a given voxel. Occupied voxels (`1`)
#'   are arranged according to their relative positions in 3-dimensional space,
#'   and fully encapsulated within a rectangular volume of unoccupied voxels
#'   (`0`). The XYZ positions of the voxels are retained in the array
#'   [`dimnames`][`dimnames()`].
#'
#' @export
#'
#' @examples
#' # Basic usage ---------------------------------------------------------------
#' # simulate a diagonal line of points with XYZ coordinates
#' pc <- data.frame(X = as.numeric(0:24), 
#'                  Y = as.numeric(0:24), 
#'                  Z = as.numeric(0:24))
#' # convert point data to cubic voxels of length 5
#' vox <- voxelize(pc, edge_length = c(5,5,5))
#' # convert to voxel array
#' box <- bounding_box(vox)
#' 
#' # Using lidR::voxel_metrics -------------------------------------------------
#' \donttest{if (require("lidR")){
#' # reformat point data into rudimentary LAS object
#' las <- suppressMessages(lidR::LAS(pc))
#' # convert to voxels of length 5
#' vox <- lidR::voxel_metrics(las, ~list(N = length(Z)), res = 5)
#' # convert to voxel array
#' box <- bounding_box(vox)
#' }
#' }
bounding_box <- function(x, threshold = 0, edge_length = NULL){
  UseMethod("bounding_box", x)
}

#' @rdname bounding_box
#' @usage NULL
#' @export
bounding_box.default <- function(x, threshold = 0, edge_length = NULL){
  # ------------------------------- Checks -------------------------------------
  
  # check that edge_length is defined
  if (is.null(edge_length)){
    stop("No voxel dimensions found, please enter a value for 'edge_length'")
  }
  
  # check that edge_length is numeric
  if (!is.vector(edge_length, mode = "numeric")){
    stop("'edge_length' argument must be a numeric vector")
  }
  
  # check that edge_length has length 3
  if (length(edge_length) != 3){
    stop("'edge_length' must be a vector of length 3")
  }
  
  # check that the required columns are present
  if (!all(c("X", "Y", "Z", "N") %in% names(x))){
    stop("Required columns not found, input data.table must have columns $X, $Y, $Z, and $N")
  }
  
  # check that threshold argument is valid
  if (!is.vector(threshold, mode = "numeric") | length(threshold) != 1){
    stop("'threshold' argument must be a single numeric value")
  }
  
  # ------------------------- Process data.table -------------------------------
  
  # locally bind variables to allow for clean non-standard evaluation by data.table
  . <- X <- Y <- Z <- N <- Occupied <- NULL
  
  # subset voxel data.table, keeping rows where N is greater than the cutoff,
  # then rounding to hopefully prevent floating-point precision issues, then 
  # removing the N column and adding an occupancy column
  rounded <- x[N > threshold][,(c("X","Y","Z","N","Occupied")) := .(round(X,10), 
                                                                    round(Y,10), 
                                                                    round(Z,10), 
                                                                    NULL,
                                                                    1L)]
  
  # create bounding box and anti-join with the rounded data.table, adding a
  # column for occupancy
  bbox <- data.table(expand.grid(X = round(seq(min(rounded$X),
                                               max(rounded$X),
                                               by = edge_length[[1]]),10),
                                 Y = round(seq(min(rounded$Y),
                                               max(rounded$Y),
                                               by = edge_length[[2]]),10),
                                 Z = round(seq(min(rounded$Z),
                                               max(rounded$Z),
                                               by = edge_length[[3]]),10)))[
                                                 !rounded, on = .(X, Y, Z)
                                               ][, Occupied := 0L]
  
  # concatenate the rounded and bounding box data.tables
  filled <- data.table::rbindlist(list(rounded, bbox))
  
  # reorder the rows to prepare for array conversion
  data.table::setorder(filled, Z, Y, X)
  
  # ----------------------------- Create array ---------------------------------
  
  # convert data.table into 3 dimensional array
  lac_array <- array(data = filled$Occupied, 
                     dim=c(length(unique(filled$X)), 
                           length(unique(filled$Y)), 
                           length(unique(filled$Z))), 
                     dimnames=list(unique(filled$X), 
                                   unique(filled$Y), 
                                   unique(filled$Z)))
  
  
  return(lac_array)
}

#' @rdname bounding_box
#' @usage NULL
#' @export
bounding_box.lac_voxels <- function(x, threshold = 0, edge_length = NULL){
  
  # check that voxel_resolution attribute is correctly formatted
  if (!is.vector(attr(x, "voxel_resolution"), mode = "numeric") | 
      length(attr(x, "voxel_resolution")) != 3){
    stop("'lac_voxels' object has improperly formatted attribute: 'voxel_resolution'")
  }

  # fetch voxel dimensions from the attributes
  res <- attr(x, "voxel_resolution")
  
  # pass to default method
  NextMethod(bounding_box, edge_length = res)
}

#' @rdname bounding_box
#' @usage NULL
#' @export
bounding_box.lasmetrics3d <- function(x, threshold = 0, edge_length = NULL){
  
  # check that res attribute is correctly formatted
  if (!is.vector(attr(x, "res"), mode = "numeric") | 
      length(attr(x, "res")) != 1){
    stop("'lasmetrics3d' object has improperly formatted attribute: 'res'")
  }
  
  # fetch voxel dimensions from the attributes
  res <- rep(attr(x, "res"), 3)
  
  # pass to default method
  NextMethod(bounding_box, edge_length = res)
}