https://github.com/cran/spatstat
Tip revision: 1b6f7ce253524f978d8945645a567ffe709bf90a authored by Adrian Baddeley on 13 October 2016, 01:00:59 UTC
version 1.47-0
version 1.47-0
Tip revision: 1b6f7ce
weights.R
#
# weights.S
#
# Utilities for computing quadrature weights
#
# $Revision: 4.38 $ $Date: 2015/05/16 05:16:46 $
#
#
# 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 times 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)
}