R/curl.R

In ambient: A Generator of Multidimensional Noise

#' Generate curl noise
#'
#' One of the use cases for fractal noise is to simulate natural phenomena.
#' perlin/simplex noise are e.g. often used to create flow fields, but this can
#' be problematic as they are not divergence-free (particles will concentrate at
#' sinks/gutters in the field). An approach to avoid this is to take the curl of
#' a field instead. The curl operator is ensured to produce divergence-free
#' output, when supplied with continuous fields such as those generated by
#' simplex and perlin noise. The end result is a field that is incompressible,
#' thus modelling fluid dynamics quite well.
#'
#' @param x,y,z The coordinates to generate the curl for as unquoted expressions
#' @param generator The noise generating function, such as [gen_simplex], or
#' [fracture()]
#' @param ... Further arguments to `generator`
#' @param seed A seed for the generator. For 2D curl the seed is a single
#' integer and for 3D curl it must be a vector of 3 integers. If `NULL` the
#' seeds will be random.
#' @param delta The offset to use for the partial derivative of the `generator`.
#' If `NULL`, it will be set as 1e-4 of the largest range of the dimensions.
#' @param mod A modification function taking the coordinates along with the
#' output of the `generator` call and allow modifications of it prior to
#' calculating the curl. The function will get the coordinates as well as a
#' `value` holding the generator output for each coordinate. If the curl is
#' requested in 2D the value will be a numeric vector and `mod()` should return
#' a numeric vector of the same length. IF the curl is requested in 3D the value
#' is a list of three numeric vectors (x, y, and z) and `mod()` should return a
#' list of three vectors of the same length. Passing NULL will use the generator
#' values unmodified.
#'
#' @export
#'
#' @family derived values
#'
#' @references
#' Bridson, Robert. Hourihan, Jim. Nordenstam, Marcus (2007). *Curl-noise for procedural fluid flow*.
#' ACM Transactions on Graphics 26(3): 46. doi:10.1145/1275808.1276435.
#'
#' @examples
#' grid <- long_grid(seq(0, 1, l = 100), seq(0, 1, l = 100))
#'
#' # Use one of the generators
#' grid$curl <- curl_noise(gen_simplex, x = grid$x, y = grid$y)
#' plot(grid$x, grid$y, type = 'n')
#' segments(grid$x, grid$y, grid$x + grid$curl$x / 100, grid$y + grid$curl$y / 100)
#'
#' # If the curl of fractal noise is needed, pass in `fracture` instead
#' grid$curl <- curl_noise(fracture, x = grid$x, y = grid$y, noise = gen_simplex,
#'                         fractal = fbm, octaves = 4)
#' plot(grid$x, grid$y, type = 'n')
#' segments(grid$x, grid$y, grid$x + grid$curl$x / 500, grid$y + grid$curl$y / 500)
#'
curl_noise <- function(generator, x, y, z = NULL, ..., seed = NULL, delta = NULL, mod = NULL) {
  if (is.null(z) || length(z) == 1) {
    .curl_noise2d(generator, x, y, z = z, ..., seed = seed, delta = delta, mod = mod)
  } else {
    .curl_noise3d(generator, x, y, z, ..., seed = seed, delta = delta, mod = mod)
  }
}

.curl_noise2d <- function(generator, x, y, ..., seed = NULL, delta = NULL, mod = NULL) {
  if (is.null(seed)) {
    seed <- rep(random_seed(), 2)
  } else {
    seed <- rep_len(seed, 2)
  }
  if (is.null(delta)) {
    delta <- max(diff(range(x)), diff(range(y))) * 1e-4
  }

  valx1 <- with_mod2(x = x + delta, y = y, generator, mod, seed = seed[1], ...)
  valx2 <- with_mod2(x = x - delta, y = y, generator, mod, seed = seed[1], ...)
  valy1 <- with_mod2(x = x, y = y + delta, generator, mod, seed = seed[2], ...)
  valy2 <- with_mod2(x = x, y = y - delta, generator, mod, seed = seed[2], ...)

  velocity_x <- -(valy1 - valy2) / (2 * delta)
  velocity_y <- (valx1 - valx2)  / (2 * delta)

  data.frame(x = velocity_x, y = velocity_y)
}
with_mod2 <- function(x, y, gen, mod, ...) {
  value <- gen(x = x, y = y, ...)
  if (is.null(mod)) return(value)
  mod(x = x, y = y, value)
}


.curl_noise3d <- function(generator, x, y, z, ..., seed = NULL, delta = NULL, mod = NULL) {
  seed <- random_seed(3, seed)
  if (is.null(delta)) {
    delta <- max(diff(range(x)), diff(range(y)), diff(range(z))) * 1e-4
  }

  valx1 <- with_mod3(x = x + delta, y = y, z = z, generator, mod, seed = seed, ...)
  valx2 <- with_mod3(x = x - delta, y = y, z = z, generator, mod, seed = seed, ...)
  valy1 <- with_mod3(x = x, y = y + delta, z = z, generator, mod, seed = seed, ...)
  valy2 <- with_mod3(x = x, y = y - delta, z = z, generator, mod, seed = seed, ...)
  valz1 <- with_mod3(x = x, y = y, z = z + delta, generator, mod, seed = seed, ...)
  valz2 <- with_mod3(x = x, y = y, z = z - delta, generator, mod, seed = seed, ...)

  velocity_x <- ((valy1[[3]] - valy2[[3]]) - (valz1[[2]] - valz2[[2]])) / (2 * delta)
  velocity_y <- ((valz1[[1]] - valz2[[1]]) - (valx1[[3]] - valx2[[3]])) / (2 * delta)
  velocity_z <- ((valx1[[2]] - valx2[[2]]) - (valy1[[1]] - valy2[[1]])) / (2 * delta)

  data.frame(x = velocity_x, y = velocity_y, z = velocity_z)
}
with_mod3 <- function(x, y, z, seed, gen, mod, ...) {
  value <- list(
    gen(x = x, y = y, z = z, seed = seed[1], ...),
    gen(x = x, y = y, z = z, seed = seed[2], ...),
    gen(x = x, y = y, z = z, seed = seed[3], ...)
  )
  if (is.null(mod)) return(value)
  mod(x = x, y = y, z = z, value)
}