R/rbdd.R

In lacunaritycovariance: Gliding Box Lacunarity and Other Metrics for 2D Random Closed Sets

#' @title Simulation of Boolean Model of Deterministic Discs
#' @export rbdd  bddcoverageprob bddcovar.iso bddcovar bddlambda bdddiscr
#' @importFrom stats fft
#'
#' @description Functions for simulating a Boolean model with grains that are
#'   discs of fixed constant radius (the abbreviation 'bdd' is short for Boolean
#'   model with Deterministic Discs). A Boolean model is a two stage model,
#'   first the locations (called germs) of grains are randomly distributed
#'   according to a Poisson point process, then a random grain is placed on each
#'   germ independently. Introductions to Boolean models are available in many
#'   stochastic geometry books (Chiu et al., 2013). Also described here are
#'   functions for calculating the coverage probability, germ intensity, and
#'   covariance from model parameters for a Boolean model with deterministic discs.
#'
#' @param lambda Intensity of the germ process (which is a Poisson point
#'   process)
#' @param discr Radius of the discs
#' @param coverp Coverage probability of the Boolean model
#' @param window The window to simulate in (an \code{owin} object)
#' @param seed Optional input (default in NULL). Is an integer passed to
#'   \code{\link[base]{set.seed}}. Used to reproduce patterns exactly.

#' @return 
#' See Functions section.

#' @section WARNING:
#' The returned object of \code{rbdd} is an \code{owin} specifying the realisation of the Boolean model within the simulation window. The simulation window is not included, thus the object returned by \code{rbdd} can have much smaller extent than the simulation window (e.g. when the simulated set is empty).
#' 
#' 
#' @examples 
#' # Simulate Boolean model with discs of radius 10.
#' # The coverage probability is very close to 0.5.
#' discr <- 10
#' w <- owin(xrange = c(0, 100), c(0, 100))
#' lambda <- 2.2064E-3 
#' xi <- rbdd(lambda, discr, w)
#' 
#' # Compute properties of Boolean model from parameters
#' cp <- bddcoverageprob(lambda, discr)
#' cvc <- bddcovar(c(-10, 10), c(-10, 10), c(0.2, 0.2), lambda, discr)

#' @references 
#' Chiu, S.N., Stoyan, D., Kendall, W.S. and Mecke, J. (2013) \emph{Stochastic Geometry and Its Applications}, 3rd ed. Chichester, United Kingdom: John Wiley & Sons.
#' @keywords spatial datagen



#' @describeIn rbdd Returns an \code{owin} that is a set generated by simulating a Boolean
#'  model with specified intensity and disc radius.
#'  The window information is not contained in this object.
#'  If the simulated set is empty then an empty \code{owin} object is returned.
#' The point process of germs is generated using \pkg{spatstat}'s \code{\link[spatstat.random]{rpoispp}}.
rbdd <- function(lambda, discr, window, seed = NULL){
  grainlib <- solist(disc(radius = discr))
  bufferdist <- 1.1 * discr

  if (!missing(seed)){set.seed(seed)}
  pp <- rpoispp(lambda = lambda, win = dilation(window, bufferdist), nsim = 1, drop = TRUE) #lambda from B\"{o}m (2002) - chosen to make coverage probability very close to 0.5
  if (pp$n == 0 ){return(complement.owin(window))}
  xibuffer <- placegrainsfromlib(pp, grainlib)
  xi <- intersect.owin(xibuffer, window)
  return(xi)
}

#' @describeIn rbdd Returns the true coverage probability given the intensity and disc radius.
bddcoverageprob  <- function(lambda, discr){
  return (1 - exp(-pi * discr ^ 2 * lambda))
}

#' @describeIn rbdd Returns the germ intensity given coverage probability and disc radius.
bddlambda <- function(coverp, discr){
  return(log(1 - coverp)/ (-pi * discr ^2))
}

#' @describeIn rbdd Returns the disc radius given coverage probability and germ intensity.
bdddiscr <- function(coverp, lambda){
  return(sqrt(log(1 - coverp)/ (-pi * lambda)))
}

#theoretical set covariance of a disc
# @param r is the radius to calculate set covariance (can be a vector)
# @param discr is the radius of disc
setcovdisc <- function(r, discr){
  setcovariance <- r*0
  rsubset <- r[r < 2 * discr]
  setcovariance[r < 2 * discr] <- 2 * discr ^ 2 * acos(rsubset / (2 * discr)) - (rsubset / 2) * sqrt(4 * discr ^ 2 - rsubset ^ 2)
  return(setcovariance)
}

#' @describeIn rbdd Returns the true covariance of points separated by a distance \code{r} given the intensity, \code{lambda} and disc radius \code{discr} of the model.
#' @param r is the radius to calculate covariance
# @param lambda is the intensity of the germ process (Poisson point process)
# @param discr is the radius of the discs.
bddcovar.iso <- function(r, lambda, discr){
  expectedsetcovariance <- setcovdisc(r, discr)
  p <- 1 - exp(-pi * discr ^ 2 * lambda)
  covariance <- 2 * p - 1 + (1 - p ) ^ 2 * exp(lambda * expectedsetcovariance)
  return(covariance)
}
# covariance as a function of vectors given in X, Y columns.
bddcovar.vec <- function(X, Y, lambda, discr){
  rlist <- sqrt(X ^ 2 + Y ^ 2)
  covar <- vector(length(rlist), mode = "numeric")
  for (i in 1:length(rlist)){
    covar[i] <- bddcovar.iso(rlist[i], lambda = lambda, discr = discr)
  }
  return(covar)
}

#' @describeIn rbdd Returns an image of the covariance as calculated from disc radius and intensity.
#' @param xrange range of x values for \code{bddcovar}
#' @param yrange range of y values for \code{bddcovar}
#' @param eps list of length 2 of the steps between samples points in x and y respectively for \code{bddcovar}.
#' If eps is of length 1 then the steps between sample points in the x and y directions will both be equal to eps.
bddcovar <- function(xrange, yrange, eps, lambda, discr){
  if (length(eps) == 1){ eps <- c(eps, eps) }
  xpts <- seq(from = xrange[1], to = xrange[2], by = eps[1])
  ypts <- seq(from = yrange[1], to = yrange[2], by = eps[2])
  mat <- outer(xpts, ypts, FUN = "bddcovar.vec", lambda = lambda, discr = discr) #rows correspond to xstep - just a quirk of outer!
  mat <- t(mat) #now columns correspond to x vals.
  return(im(mat, xcol = xpts, yrow = ypts))
}