https://github.com/cran/spatstat
Tip revision: cbcb03635eaa80aad79cfed5de45794dab5eb73e authored by Adrian Baddeley on 22 August 2006, 14:23:25 UTC
version 1.9-5
version 1.9-5
Tip revision: cbcb036
weights.S
#
# weights.S
#
# Utilities for computing quadrature weights
#
# $Revision: 4.18 $ $Date: 2006/05/31 08:22:08 $
#
#
# 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'.
#
# dirichlet.weights() 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.
#
# discretise() 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 <- factor(id, levels=seq(areas))
counts <- table(id)
w <- areas[id] / counts[id] # ensures denominator > 0
w <- as.vector(w)
#
# that's it; but check for funny business
#
zerocount <- (counts == 0)
zeroarea <- (areas == 0)
if(any(!zeroarea & zerocount))
warning("some tiles with positive area do not contain any points")
if(any(!zerocount & zeroarea)) {
warning("Some tiles with zero area contain points")
warning("Some weights are zero")
attr(w, "zeroes") <- zeroarea[id]
}
#
names(w) <- NULL
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 <- discretise(x, xrange, nx)
iy <- discretise(y, yrange, ny)
#
return(list(ix=ix, iy=iy, index=(iy-1) * nx + ix))
}
discretise <- function(x, xrange, nx) {
i <- ceiling( nx * (x - xrange[1])/diff(xrange))
i <- pmax(1, i)
i <- pmin(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(ntile, 2)
nx <- ntile[1]
ny <- ntile[2]
if(verbose)
cat(paste("grid weights for a", nx, "x", ny, "grid of tiles\n"))
if(is.null(areas)) {
# compute tile areas
if(win$type == "rectangle") {
tilearea <- area.owin(win)/(nx * ny)
areas <- rep(tilearea, nx * ny)
} else {
# convert window to mask
win <- as.mask(win, dimyx=rev(npix))
if(verbose) {
if(!is.null(npix))
np <- npix
else {
np <- spatstat.options("npixel")
if(length(np) == 1) np <- rep(np, 2)
}
cat(paste("Approximating window by mask (",
np[1], " x ", np[2], " pixels)\n", sep=""))
}
# extract pixel coordinates inside window
xx <- as.vector(raster.x(win)[win$m])
yy <- as.vector(raster.y(win)[win$m])
# classify all pixels into tiles
pixelid <- gridindex(xx, yy,
win$xrange, win$yrange, nx, ny)$index
pixelid <- factor(pixelid, levels=seq(nx * ny))
# compute digital areas of tiles
tilepixels <- as.vector(table(pixelid))
pixelarea <- win$xstep * win$ystep
areas <- tilepixels * pixelarea
}
}
# classify each point according to its tile
if(win$type == "mask") {
# first move each data point to nearest pixel
Xapprox <- nearest.raster.point(X$x, X$y, win, indices=FALSE)
x <- Xapprox$x
y <- Xapprox$y
}
id <- gridindex(x, y, win$xrange, win$yrange, nx, ny)$index
# 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(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)
x <- X$x
y <- X$y
win <- X$window
if (exact) {
# check that deldir is available
old.op <- options(warn=-1)
on.exit(options(old.op))
delthere <- require(deldir,quietly=TRUE)
options(old.op)
if(!delthere) {
warning("\'deldir\' package not found; using discrete approximation\n")
exact <- FALSE
}
}
if(exact && (win$type == "rectangle")) {
rw <- c(win$xrange, win$yrange)
# invoke deldir() with NO DUMMY POINTS
tessellation <- deldir(x, y, dpl=NULL, rw=rw)
# extract tile areas
w <- tessellation$summary[, 'dir.area']
} else {
# Compute digital areas of Dirichlet tiles.
win <- as.mask(win)
X$window <- win
#
# Nearest data point to each pixel:
tileid <- exactdt(X)$i
#
if(win$type == "mask")
# Restrict to window (result is a vector - OK)
tileid <- tileid[win$m]
# Count pixels in each tile
id <- factor(tileid, levels=seq(X$n))
counts <- table(id)
# turn off the christmas lights
class(counts) <- NULL
names(counts) <- NULL
dimnames(counts) <- NULL
# Convert to digital area
pixelarea <- win$xstep * win$ystep
w <- pixelarea * counts
# Check for zero pixel counts
zeroes <- (counts == 0)
if(any(zeroes)) {
warning("some Dirichlet tiles have zero digital area")
attr(w, "zeroes") <- zeroes
}
}
# 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)
}