R/rbdr.R

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

#' @title Simulation of Boolean Model of Deterministic Rectangles
#' @export rbdr  bdrcoverageprob  bdrcovar
#' 
#' @description Functions for simulating a Boolean model with grains that are deterministic rectangles.
#' 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.
#' An introduction can be found in (Chiu et al., 2013).
#' Also described in this help file are functions for calculating the coverage probability and covariance.
#' 
#' @param lambda Intensity of the germ process (which is a Poisson point process)
#' @param grain Rectangle object specifying the grain
#' @param win 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.
#' @param xy A raster object that specifies the pixel coordinates of the desired covariance image. \code{xy} works in similar fashion to passing an image or pixel mask through the \code{xy} argument of \code{\link[spatstat.geom]{as.mask}} in \pkg{spatstat}.

#' @return 
#' Depends on the function used (see Functions section).

#' @section WARNING:
#' The returned object of \code{rbdr} is only the foreground of a binary map and thus can have much smaller extent than the simulation window (e.g. when the simulated set is empty).
#' 
#' 
#' @examples 
#' grain <- owin(xrange = c(-5, 5), yrange = c(-5, 5))
#' win <- owin(xrange = c(0, 100), c(0, 100))
#' lambda <- 4.2064E-3
#' xi <- rbdr(lambda, grain, win)
#' 
#' cp_theoretical <- bdrcoverageprob(lambda, grain)
#' xy <- as.mask(dilationAny(win, win), eps = c(1, 1))
#' cvc_theoretical <- bdrcovar(lambda, grain, xy)

#' @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 rbdr Returns an \code{owin} that is a set generated by simulating a Boolean
#'  model with a specified intensity and fixed rectangular grain.
#'  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 spatstat's \code{\link[spatstat.random]{rpoispp}}.
rbdr <- function(lambda, grain, win, seed = NULL){
  grainlib <- solist(grain)
  bufferdist <- 1.1 *
    (diameter.owin(grain)/2 + sqrt(sum(unlist(centroid.owin(grain))^2)))
  
  if (!missing(seed)){set.seed(seed)}
  pp <- rpoispp(lambda = lambda, win = dilation.owin(win, bufferdist), nsim = 1, drop = TRUE) 
  if (pp$n == 0 ){return( setminus.owin(square(r = 1), square(r = 1)) )}
  xibuffer <- placegrainsfromlib(pp, grainlib, w = win)
  xi <- intersect.owin(xibuffer, win)
  return(xi)
}

#' @describeIn rbdr Returns the true coverage probability given the intensity and grain.
bdrcoverageprob  <- function(lambda, grain){
  return (1 - exp(- lambda * area.owin(grain)))
}


#theoretical set covariance of a rectangle
#grain is an owin rectangle object
#v1, v2 is shift vector list for x and y directions equally, or vec is list with x and y components
setcov.vec.rect <- function(rectang, v1 = NULL, v2 = NULL, vec = NULL){
  stopifnot(is.rectangle(rectang))
  if (!is.null(vec)){
    v1 <- vec$x
    v2 <- vec$y
  }
  xwidth <- pmax(0, (rectang$xrange[[2]] - rectang$xrange[[1]]) - abs(v1))
  ywidth <- pmax(0, (rectang$yrange[[2]] - rectang$yrange[[1]]) - abs(v2))
  return(xwidth * ywidth)
}

#' @describeIn rbdr Returns an image of the covariance as calculated from disc radius and intensity.
bdrcovar <- function(lambda, grain, xy){
  expectedsetcovariance <- as.im.funxy(function(x, y) setcov.vec.rect(grain, v1 = x, v2 = y)  ,
  W = Frame(xy), xy = xy)
  p <- bdrcoverageprob(lambda, grain)
  covariance <- eval.im(2 * p - 1 + (1 - p ) ^ 2 * exp(lambda * expectedsetcovariance))
  return(covariance)
}