https://github.com/cran/spatstat
Revision e575bb3736f6d70d2bd2b23ea2e4474cdbae00be authored by Adrian Baddeley on 11 November 2010, 13:09:43 UTC, committed by cran-robot on 11 November 2010, 13:09:43 UTC
1 parent cafe44b
Tip revision: e575bb3736f6d70d2bd2b23ea2e4474cdbae00be authored by Adrian Baddeley on 11 November 2010, 13:09:43 UTC
version 1.21-1
version 1.21-1
Tip revision: e575bb3
weights.R
#
# weights.S
#
# Utilities for computing quadrature weights
#
# $Revision: 4.25 $ $Date: 2010/07/14 01:56:44 $
#
#
# 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.
#
# 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(areas))
counts <- table(fid)
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 <- 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(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 && (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)
}
Computing file changes ...