R/evaluate-grf.R

In mathiasisaksen/GRFics: Efficient generation of noise functions from approximate Gaussian random fields

#' Interpolate the value of a GRF function in specified locations.
#'
#' @param locations An N x 2 matrix containing the coordinates of the N locations of interest.
#' @param grf.object The GRF object containing the function, generated by \code{\link{generate.grf.object}}.
#' @param function.number The index of the function of interest.
#' @param function.name The name of the function of interest.
#' @param periodic If true, locations outside the rectangular domain of \code{grf.object} are mapped to the corresponding point inside the domain. Otherwise, the locations are mapped to the nearest point on the boundary. It is not recommended to use \code{periodic = FALSE}.
#' @param interpolation.method Specifies the method used for interpolation. Either \code{"nearest"} (which uses the value of the nearest grid cell), \code{"bilinear"} (which performs bilinear interpolation based on the nearest four grid cells), or \code{"bicubic"} (which performs bicubic interpolation based on the nearest 16 neighbors).
#' @param rescale.method Specifies how the interpolated values are rescaled. Either \code{"none"} (no rescaling), \code{"minmax"} (min-max scaling to the interval [0, 1]), or \code{"uniform"} (a scaling to [0, 1] that attempts to make all of the values equally represented).
#' @param return.df If true, the return value is a data frame containing 3 columns: The x- and y-coordinates from \code{locations}, and the interpolated values.
#' If false, only the interpolated values are returned.
#'
#' @return A \code{numeric} of length N containing the interpolated values of the GRF in the specified locations.
#'
#' @note Using \code{interpolation.method = "bicubic"} can lead to interpolated values that are outside the range of the original grid values.
#' If this is combined with \code{rescale.method = "minmax"}, then the resulting values are not guaranteed to be inside the interval [0, 1].
#' This can be fixed by mapping values greater than 1 to 1, and less than 0 to 0.
#'
#' @examples
#' # Comparison of different interpolation and rescaling methods
#' library(GRFics)
#' library(ggplot2)
#'
#' # Create GRF object on a 20 x 20 grid
#' grf.object = generate.grf.object(-1, 1, -1, 1, 20, 20, num.functions = 1,
#'                                  range.parameter = 1,
#'                                  strength.parameter = 1, direction.parameter = pi/4,
#'                                  initial.seed = 123)
#'
#' # Create 500 x 500 grid for interpolation
#' interp.grid = generate.grid.centers(-1, 1, -1, 1, 500, 500)
#' rescale.method = c("minmax", "uniform")
#' comparison.df = data.frame()
#' for (i in 1:2) {
#'   nearest.values = evaluate.grf(interp.grid, grf.object, interpolation.method = "nearest", rescale.method = rescale.method[i])
#'   bilinear.values = evaluate.grf(interp.grid, grf.object, interpolation.method = "bilinear", rescale.method = rescale.method[i])
#'   bicubic.values = evaluate.grf(interp.grid, grf.object, interpolation.method = "bicubic", rescale.method = rescale.method[i])
#'   comparison.df = rbind(comparison.df,
#'                         data.frame(
#'                           x = interp.grid$x, y = interp.grid$y,
#'                           z = c(nearest.values, bilinear.values, bicubic.values),
#'                           interpolation.method = rep(c("1. nearest", "2. bilinear", "3. bicubic"), each = nrow(interp.grid)),
#'                           rescale.method = rescale.method[i]))
#' }
#'
#' # Show comparison plot using facet_grid
#' ggplot()+
#'   theme_bw()+
#'   geom_raster(data = comparison.df, aes(x = x, y = y, fill = z))+
#'   facet_grid(cols = vars(interpolation.method), rows = vars(rescale.method))+
#'   scale_fill_gradientn(colours = c("black", "white"))+
#'   coord_fixed()
#'
#' @export
#' @author Mathias Isaksen  \email{mathiasleanderi@@gmail.com}
evaluate.grf = function(locations, grf.object, function.number = 1, function.name = NULL, periodic = TRUE, interpolation.method = "bicubic", rescale.method = "minmax", return.df = FALSE) {
  if(is.null(grf.object$functions)) {
    stop("No functions have been generated.")
  }
  # Unpack parameters to function scope
  list2env(grf.object$params, environment())
  norm.locs = normalize.locations(locations, grf.object)
  norm.grid = normalize.locations(grf.object$grid, grf.object)
  if (periodic) {
    norm.locs = map.locations.periodic(norm.locs)
  } else {
    norm.locs = map.locations.boundary(norm.locs)
  }

  # Compute separate x- and y-components of grid
  grid.vectors = generate.grid.centers(0, 1, 0, 1, resolution.x, resolution.y, separate = TRUE)
  x.g = grid.vectors$x
  y.g = grid.vectors$y
  h.x = x.g[2] - x.g[1]
  h.y = y.g[2] - y.g[1]

  if (interpolation.method == "bilinear") {
    weight.func = function(x) ifelse(x <= 1 & x >= 0, x, ifelse(x > 1, 1, 0))
  } else if (interpolation.method == "nearest") {
    weight.func = function(x) ifelse(x >= 0.5, 1, 0)
  }

  # Get values of desired function
  function.name = ifelse(is.null(function.name), names(grf.object$functions)[function.number], function.name)
  values = grf.object$functions[[function.name]]

  # Convert onto matrix form for simpler syntax
  value.matrix = matrix(values, nrow = resolution.y, ncol = resolution.x, byrow = TRUE)
  interpolated.values = NA
  # This next part should be moved modularized
  if (interpolation.method == "bicubic") {
    if (periodic) {
      # Bicubic interpolation uses a 4x4 neighborhood. Padding grid in a periodic manner:
      x.g = c(x.g[1] - h.x, x.g, x.g[resolution.x] + h.x)
      y.g = c(y.g[1] - h.y, y.g, y.g[resolution.y] + h.y)
      value.matrix = cbind(value.matrix[, resolution.x], value.matrix, value.matrix[, 1])
      value.matrix = rbind(value.matrix[resolution.y, ], value.matrix, value.matrix[1, ])
    }
    interpolated.values = akima::bicubic(x.g, y.g, t(value.matrix), norm.locs[, 1], norm.locs[, 2])$z

  } else if (interpolation.method %in% c("nearest", "bilinear")) {
    # Finds the indices of the columns that are to the left and right of each location
    left = round(norm.locs[, 1]/h.x)
    right = left + 1

    # Finds the indices of the rows that are below and above each location
    bottom = round(norm.locs[, 2]/h.y)
    top = bottom + 1

    # Periodic mapping
    if (periodic) {
      left = ifelse(left == 0, resolution.x, left)
      right = ifelse(right == (resolution.x+1), 1, right)
      bottom = ifelse(bottom == 0, resolution.y, bottom)
      top = ifelse(top == (resolution.y+1), 1, top)
    } else { # Nearest-boundary mapping
      left = ifelse(left == 0, 1, left)
      right = ifelse(right == (resolution.x+1), resolution.x, right)
      bottom = ifelse(bottom == 0, 1, bottom)
      top = ifelse(top == (resolution.y+1), resolution.y, top)
    }

    # Distance between grid cells, in this case constant
    width.x = x.g[right] - x.g[left]
    width.y = y.g[top] - y.g[bottom]

    #
    weight.x = weight.func((norm.locs[, 1]-x.g[left])/width.x)
    weight.y = weight.func((norm.locs[, 2]-y.g[bottom])/width.y)

    # For each location, compute a weighted mean of the values in the the four closest grid cells
    interpolated.values = (1-weight.x)*(1-weight.y)*value.matrix[cbind(bottom, left)] +
      weight.x*(1-weight.y)*value.matrix[cbind(bottom, right)] +
      (1-weight.x)*weight.y*value.matrix[cbind(top, left)] +
      weight.x*weight.y*value.matrix[cbind(top, right)]
  } else {
    stop("interpolation.method must be nearest, bilinear, or bicubic.")
  }

  if (rescale.method == "uniform") {
    interpolated.values = grf.object$ecdf.functions[[function.name]](interpolated.values)
  } else if (rescale.method == "minmax") {
    interpolated.values = standardize(interpolated.values, lower = min(values), upper = max(values))
  } else if (rescale.method != "none") {
    stop("Invalid value specified for rescale.method.")
  }
  if (!return.df) {
    return(interpolated.values)
  } else {
    return(data.frame(x = locations[, 1], y = locations[, 2], z = interpolated.values))
  }
}