R/bw.abram.R

In spatstat/spatstat.core: Core Functionality of the 'spatstat' Family

#'
#'   bw.abram.R
#'
#'   Abramson bandwidths
#'
#'   $Revision: 1.8 $ $Date: 2019/04/25 05:34:49 $
#'

bw.abram <- function(X, h0,
                     ..., 
                     at=c("points", "pixels"), 
                     hp=h0, pilot=NULL, trim=5,
                     smoother=density.ppp){
  stopifnot(is.ppp(X))
  at <- match.arg(at)

  if(missing(h0) || is.null(h0)) {
    h0 <- bw.ppl(X)
  } else {
    check.1.real(h0)
    stopifnot(h0 > 0)
  }

  check.1.real(trim)
  stopifnot(trim > 0)
  
  pilot.data <- X
  imwin <- as.im(Window(X), ...)
  
  if(is.im(pilot)){
    if(!compatible.im(imwin,pilot))
      stop("'X' and 'pilot' have incompatible spatial domains", call.=FALSE)
    #' clip the worst small values away
    pilot[pilot<=0] <- min(pilot[pilot>0]) 
  } else if(is.ppp(pilot)){
    if(!compatible.im(imwin,as.im(Window(pilot), ...)))
      stop("'X' and 'pilot' have incompatible spatial domains", call.=FALSE)
    pilot.data <- pilot
  } else if(!is.null(pilot))
    stop("if supplied, 'pilot' must be a pixel image or a point pattern",
         call.=FALSE)
  
  if(!is.im(pilot)) {
    if(is.character(smoother)) {
      smoother <- get(smoother, mode="function")
    } else stopifnot(is.function(smoother))
    pilot <- smoother(pilot.data,sigma=hp,positive=TRUE,...)
  }
  
  pilot <- pilot/integral(pilot) # scale to probability density
  pilotvalues <- safelookup(pilot, pilot.data, warn=FALSE)
  ## geometric mean re-scaler (Silverman, 1986; ch 5).  
  gamma <- exp(mean(log(pilotvalues[pilotvalues > 0])))^(-0.5)

  switch(at,
         points = {
           pilot.X <- safelookup(pilot,X,warn=FALSE)
           bw <- h0 * pmin((pilot.X^(-0.5))/gamma,trim)
         },
         pixels = {
           bw <- eval.im(h0 * pmin((pilot^(-0.5))/gamma, trim))
         })
  
  return(bw)
}