R/weights.R

In spatstat.geom: Geometrical Functionality of the 'spatstat' Family

#
#	weights.S
#
#	Utilities for computing quadrature weights
#
#	$Revision: 4.40 $	$Date: 2020/01/10 06:54:23 $
#
#
# Main functions:
		
#	gridweights()	    Divide the window frame into a regular nx * ny
#			    grid of rectangular tiles. Given an arbitrary
#			    pattern of (data + dummy) points derive the
#			    'counting weights'.
#
#	dirichletWeights()  Compute the areas of the tiles of the
#			    Dirichlet tessellation generated by the 
#			    given pattern of (data+dummy) points,
#			    restricted to the window.
#	
# Auxiliary functions:	
#			
#       countingweights()   compute the counting weights
#                           for a GENERIC tiling scheme and an arbitrary
#			    pattern of (data + dummy) points,
#			    given the tile areas and the information
#			    that point number k belongs to tile number id[k]. 
#
#
#	gridindex()	    Divide the window frame into a regular nx * ny
#			    grid of rectangular tiles. 
#			    Compute tile membership for arbitrary x,y.
#				    
#       grid1index()        1-dimensional analogue of gridindex()
#
#
#-------------------------------------------------------------------
	
countingweights <- function(id, areas, check=TRUE) {
	#
	# id:        cell indices of n points
	#                     (length n, values in 1:k)
	#
	# areas:     areas of k cells 
	#                     (length k)
	#
  id <- as.integer(id)
  fid <- factor(id, levels=seq_along(areas))
  counts <- table(fid)
  w <- areas[id] / counts[id]     # ensures denominator > 0
  w <- as.vector(w)
#	
# that's it; but check for funny business
#
  if(check) {
    zerocount <- (counts == 0)
    zeroarea <- (areas == 0)
    if(any(!zeroarea & zerocount)) {
      lostfrac <- 1 - sum(w)/sum(areas)
      lostpc <- round(100 * lostfrac, 1)
      if(lostpc >= 1) 
        warning(paste("some tiles with positive area",
                      "do not contain any quadrature points:",
                      "relative error =",
                      paste0(lostpc, "%")))
    }
    if(any(!zerocount & zeroarea)) {
	warning("Some tiles with zero area contain quadrature points")
	warning("Some weights are zero")
	attr(w, "zeroes") <- zeroarea[id]
    }
  }
#
  names(w) <- NULL
  return(w)
}

gridindex <- function(x, y, xrange, yrange, nx, ny) {
	#
	# The box with dimensions xrange, yrange is divided
	# into nx * ny cells.
	#
	# For each point (x[i], y[i]) compute the index (ix, iy)
	# of the cell containing the point.
	# 
	ix <- grid1index(x, xrange, nx)
	iy <- grid1index(y, yrange, ny)
	#
	return(list(ix=ix, iy=iy, index=as.integer((iy-1) * nx + ix)))
}

grid1index <- function(x, xrange, nx) {
	i <- ceiling( nx * (x - xrange[1])/diff(xrange))
	i <- pmax.int(1, i)
	i <- pmin.int(i, nx)
	i
}

gridweights <- function(X, ntile=NULL, ..., window=NULL, verbose=FALSE,
                        npix=NULL, areas=NULL) {
	#
	# Compute counting weights based on a regular tessellation of the
	# window frame into ntile[1] * ntile[2] rectangular tiles.
	#
	# Arguments X and (optionally) 'window' are interpreted as a
	# point pattern.
	#
	# The window frame is divided into a regular ntile[1] * ntile[2] grid
	# of rectangular tiles. The counting weights based on this tessellation
	# are computed for the points (x, y) of the pattern.
	#
        # 'npix' determines the dimensions of the pixel raster used to
        # approximate tile areas.
	
	X <- as.ppp(X, window)
	x <- X$x
	y <- X$y
        win <- X$window

        # determine number of tiles
        if(is.null(ntile))
          ntile <- default.ntile(X)
        if(length(ntile) == 1)
          ntile <- rep.int(ntile, 2)
        nx <- ntile[1]
        ny <- ntile[2]

        if(verbose) 
          cat(paste("grid weights for a", nx, "x", ny, "grid of tiles\n"))

        ## determine pixel resolution in case it is required
        if(!is.null(npix)) {
          npix <- ensure2vector(npix)
        } else {
          npix <- pmax(rev(spatstat.options("npixel")),
                       c(nx, ny))
          if(is.mask(win))
            npix <- pmax(npix, rev(dim(win)))
        }
        
        if(is.null(areas)) {
  	  # compute tile areas
          switch(win$type,
                 rectangle = {
                   nxy <- nx * ny
                   tilearea <- area(win)/nxy
                   areas <- rep.int(tilearea, nxy)
                   zeroareas <- rep(FALSE, nxy)
                 },
                 polygonal = {
                   areamat <- polytileareaEngine(win,
                                                 win$xrange, win$yrange,
                                                 nx, ny)
                   ## convert from 'im' to 'gridindex' ordering
                   areas <- as.vector(t(areamat))
                   zeroareas <- (areas == 0)
                   if(verbose)
                     splat("Split polygonal window of area", area(win),
                           "into", nx, "x", ny, "grid of tiles",
                           "of total area", sum(areas))
                 },
                 mask = {
                   win <- as.mask(win, dimyx=rev(npix))
                   if(verbose) 
                     splat("Converting mask dimensions to",
                           npix[1], "x", npix[2], "pixels")
                   ## extract pixel coordinates inside window
                   rxy <- rasterxy.mask(win, drop=TRUE)
                   xx <- rxy$x
                   yy <- rxy$y
                   
                   ## classify all pixels into tiles
                   pixelid <- gridindex(xx, yy, 
                                        win$xrange, win$yrange, nx, ny)$index
                   pixelid <- factor(pixelid, levels=seq_len(nx * ny))
                                
                   ## compute digital areas of tiles
                   tilepixels <- as.vector(table(pixelid))
                   pixelarea <- win$xstep * win$ystep
                   areas <- tilepixels * pixelarea
                   zeroareas <- (tilepixels == 0)
                 }
               )
        } else zeroareas <- (areas == 0)
        
        id <- gridindex(x, y, win$xrange, win$yrange, nx, ny)$index

        if(win$type != "rectangle" && any(uhoh <- zeroareas[id])) {
          # this can happen: the tile has digital area zero
          # but contains a data/dummy point
          slivers <- unique(id[uhoh])
          switch(win$type,
                 mask = {
                   offence <- "digital area zero"
                   epsa <- pixelarea/2
                 },
                 polygonal = {
                   offence <- "very small area"
                   epsa <- min(areas[!zeroareas])/10
                 })
          areas[slivers] <- epsa
          nsliver <- length(slivers)
          extraarea <- nsliver * epsa
          extrafrac <- extraarea/area(win)
          if(verbose || extrafrac > 0.01) {
            splat(nsliver, ngettext(nsliver, "tile", "tiles"),
                  "of", offence,
                  ngettext(nsliver, "was", "were"),
                  "given nominal area", signif(epsa, 3),
                  "increasing the total area by",
                  signif(extraarea, 5), "square units or",
                  paste0(round(100 * extrafrac, 1), "% of total area"))
            if(extrafrac > 0.01)
              warning(paste("Repairing tiles with", offence,
                            "caused a",
                            paste0(round(100 * extrafrac), "%"),
                            "change in total area"),
                      call.=FALSE)
          }
        }
 
	# compute counting weights 
	w <- countingweights(id, areas)

        # attach information about weight construction parameters
        attr(w, "weight.parameters") <-
          list(method="grid", ntile=ntile, npix=npix, areas=areas)
        
	return(w)
}


# dirichlet.weights <- function(...) {
#   .Deprecated("dirichletWeights", package="spatstat")
#   dirichletWeights(...)
# }

dirichletWeights <- function(X, window = NULL, exact=TRUE, ...) {
  #'
  #' Compute weights based on Dirichlet tessellation of the window 
  #' induced by the point pattern X. 
  #' The weights are just the tile areas.
  #'
  #' NOTE:	X should contain both data and dummy points,
  #' if you need these weights for the B-T-B method.
  #'
  #' Arguments X and (optionally) 'window' are interpreted as a
  #' point pattern.
  #'
  #' If the window is a rectangle, we invoke Rolf Turner's "deldir"
  #' package to compute the areas of the tiles of the Dirichlet
  #' tessellation of the window frame induced by the points.
  #' [NOTE: the functionality of deldir to create dummy points
  #' is NOT used. ]
  #'	if exact=TRUE	compute the exact areas, using "deldir"
  #'	if exact=FALSE      compute the digital areas using exactdt()
  #' 
  #' If the window is a mask, we compute the digital area of
  #' each tile of the Dirichlet tessellation by counting pixels.
  #'
  #'
  #' 
  X <- as.ppp(X, window)

  if(!exact && is.polygonal(Window(X)))
    Window(X) <- as.mask(Window(X))
  
  #' compute tile areas
  w <- dirichletAreas(X)

  #' zero areas can occur due to discretisation or weird geometry
  zeroes <- (w == 0)
  if(any(zeroes)) {
    #' compute weights for subset
    nX <- npoints(X)
    Xnew <- X[!zeroes]
    wnew <- dirichletAreas(Xnew)
    w    <- numeric(nX)
    w[!zeroes] <- wnew
    #' map deleted points to nearest retained point
    jj <- nncross(X[zeroes], Xnew, what="which")
    #' map retained points to themselves
    ii <- Xseq <- seq_len(nX)
    ii[zeroes] <- (ii[!zeroes])[jj]
    #' redistribute weights
    nshare <- table(factor(ii, levels=Xseq))
    w <- w[ii]/nshare[ii]
  }
  #' attach information about weight construction parameters
  attr(w, "weight.parameters") <- list(method="dirichlet", exact=exact)
  return(w)
}

default.ntile <- function(X) { 
	# default number of tiles (n x n) for tile weights
        # when data and dummy points are X
  X <- as.ppp(X)
  guess.ngrid <- 10 * floor(sqrt(X$n)/10)
  max(5, guess.ngrid/2)
}