R/pcf_anisotropic_units.R

In antiphon/Kdirectional: Analysis of anisotropy in point patterns using second order statistics

#' Anisotropic pair correlation function, direction vector formulation
#' 
#' Estimate the anisotropic pair correlation function (2d and 3d), as defined in Stoyan 1991, f. 5.2-5.3.
#' This one doesn't use spherical coordinates but direction vectors.
#' 
#' @param x pp, list with $x~coordinates $bbox~bounding box
#' @param directions Matrix of direction vectors
#' @param h widths of epanechnicov kernels, vector of two values, for ranges and angles.
#' @param f If h not given, use h=f/lambda^(1/dim) for range and h=f*pi
#' @param correction "none" or translation. Translation only for rectangle box.
#' @param n_dir Direction vector grid resolution (if units not given). See Details.
#' @param r In case directions not given, use these lengths for generated directions.
#' 
#' @details
#' For example, in 2D the polar vector (r, theta) is now a vector (r*cos(theta), r*sin(theta)).
#' 
#' For 2d, we use a regular grid of n_dir steps on [0,pi], leaving out pi due to antipode-symmetry to 0. 
#' n_dir=7 is the default.
#' 
#' For 3d we use a triangular grid on the expanding sphere. This is generated by
#' subdiving a regular icosahedron n_dir times. n_dir=2 is the default.
#' 
#' @useDynLib Kdirectional
#' @importFrom rgl subdivision3d
#' @export

pcf_directions <- function(x, directions, h, f=0.15, correction="translation", n_dir, r) {
  x <- check_pp(x)
  bbox <- as.matrix(x$bbox)
  dim <- ncol(bbox)
  sidelengths <- bbox_sideLengths(bbox)
  lambda <- nrow(x$x)/prod(sidelengths)
  #
  #
  # Generate a grid of direction vectors
  if(missing(directions)){
    if(missing(r)){ 
      b <- min(sidelengths)*0.2
      r <- seq(0.001, b, length=30)
    }
    if(dim==2){
      if(missing(n_dir)) n_dir <- 7
      dirs <- check_directions(dim=dim, n_dir=n_dir)
      ang <- dirs$theta[[1]]
      directions1 <- cbind(cos(ang), sin(ang))
    }
    else{
      if(missing(n_dir)) n_dir <- 2
      s0 <- icosahedron3d()
      for(i in 1:n_dir) s0 <- subdivision3d(s0)
      vb <- t(s0$vb[1:3,])/s0$vb[4,]
      directions2 <- vb
      # only the top half
      directions1 <- directions2[ directions2[,3]>=0, ]
    }
    # now expand the circles/spheres
    directions <- NULL
    for(ri in r){
      directions <- rbind(directions, ri*directions1)
    }
  }
  #
  #
  # check smoothing parameters
  if(missing(h)) {
    h <- c(f/lambda^(1/dim), f*pi)
  }
  if(length(h) < 2) h <- rep(h, 2)
  # Edge correction
  correction_i <- pmatch(correction, c("none", "translation"))
  #
  # start:
  xc <- as.matrix(x$x)
  res <- c_anisotropic_unit_pcf(xc, directions, h, bbox, correction_i)
  
  # pcf
  g <- res[[1]]/lambda^2
  
  # compile
  nv <- nrow(directions)
  # report also in polar/spherical
  l <- apply(directions, 1, function(x)sqrt(t(x)%*%x))
  phi <- atan2(directions[,2], directions[,1])
  phi[phi<0] <- phi[phi<0] + 2*pi
  r_phi <- cbind(l, phi  )
  if(dim==3) r_phi <- cbind(r_phi, acos(directions[,3]/l)  )
  #
  # ok, done.
  res <- list(est=g, directions=directions, r_phi=r_phi, counts=res[[2]], 
              correction=correction, h=h, dim=dim)
  class(res) <- "pcf_directions"
  res
}