https://github.com/cran/spatstat
Tip revision: 94229838c61febbc52fa9a23dca6216f3321f76a authored by Adrian Baddeley on 01 September 2004, 16:21:11 UTC
version 1.5-3
version 1.5-3
Tip revision: 9422983
dummy.S
#
# dummy.S
#
# Utilities for generating patterns of dummy points
#
# $Revision: 4.11 $ $Date: 2004/06/22 02:42:33 $
#
# corners() corners of window
# gridcenters() points of a rectangular grid
# stratrand() random points in each tile of a rectangular grid
# spokes() Rolf's 'spokes' arrangement
#
# concatxy() concatenate any lists of x, y coordinates
#
# default.dummy() Default action to create a dummy pattern
#
corners <- function(window) {
window <- as.owin(window)
x <- window$xrange[c(1,2,1,2)]
y <- window$yrange[c(1,1,2,2)]
return(list(x=x, y=y))
}
gridcenters <-
gridcentres <- function(window, nx, ny) {
window <- as.owin(window)
xr <- window$xrange
yr <- window$yrange
x <- seq(xr[1], xr[2], length = 2 * nx + 1)[2 * (1:nx)]
y <- seq(yr[1], yr[2], length = 2 * ny + 1)[2 * (1:ny)]
x <- rep(x, ny)
y <- rep(y, rep(nx, ny))
return(list(x=x, y=y))
}
stratrand <- function(window,nx,ny, k=1) {
# divide window into an nx * ny grid of tiles
# and place k points at random in each tile
window <- as.owin(window)
wide <- diff(window$xrange)/(nx - 1)
high <- diff(window$yrange)/(ny - 1)
cent <- gridcentres(window, nx, ny)
cx <- rep(cent$x, k)
cy <- rep(cent$y, k)
n <- nx * ny * k
x <- cx + runif(n, min = -wide/2, max = wide/2)
y <- cy + runif(n, min = -high/2, max = high/2)
return(list(x=x,y=y))
}
tilecentroids <- function (W, nx, ny)
{
W <- as.owin(W)
if(W$type == "rectangle")
return(gridcentres(W, nx, ny))
else {
# approximate
W <- as.mask(W)
xx <- as.vector(raster.x(W)[W$m])
yy <- as.vector(raster.y(W)[W$m])
pid <- gridindex(xx,yy,W$xrange,W$yrange,nx,nx)$index
x <- tapply(xx,pid,mean)
y <- tapply(yy,pid,mean)
return(list(x=x,y=y))
}
}
tilemiddles <- function (W, nx, ny)
{
if(W$type == "rectangle")
return(gridcentres(W, nx, ny))
# pixel approximation to window
# This matches the pixel approximation used to compute tile areas
# and ensures that dummy points are generated only inside those tiles
# that have nonzero digital area
M <- as.mask(W)
xx <- as.vector(raster.x(M))[M$m]
yy <- as.vector(raster.y(M))[M$m]
# label the pixels by tile index
pid <- gridindex(xx,yy,W$xrange,W$yrange,nx,ny)$index
pidvalues <- levels(as.factor(pid))
# Latter is for comparison with the output of tapply( ..., pid, ...)
######## 1st try: centroids of tiles
# (usually inside tile, e.g. if it's convex, but not always inside)
cx <- tapply(xx, pid, mean)
cy <- tapply(yy, pid, mean)
# centroids which are inside window
cok <- inside.owin(cx, cy, W)
x <- cx[cok]
y <- cy[cok]
# tiles for which this strategy did not work
nbg <- pidvalues[!cok]
######### 2nd try: middle point in list of pixels in each tile
# (always inside tile, by construction)
todo <- pid %in% nbg
xx <- xx[todo]
yy <- yy[todo]
pid <- pid[todo]
middle <- function(v) { n <- length(v); mid <- ceiling(n/2); v[mid]}
mx <- tapply(xx,pid,middle)
my <- tapply(yy,pid,middle)
x <- c( x, mx)
y <- c( y, my)
return(list(x=x,y=y))
}
spokes <- function(x, y, nrad = 3, nper = 3, fctr = 1.5, Mdefault=1) {
#
# Rolf Turner's "spokes" arrangement
#
# Places dummy points on radii of circles
# emanating from each data point x[i], y[i]
#
# nrad: number of radii from each data point
# nper: number of dummy points per radius
# fctr: length of largest radius = fctr * M
# where M is mean nearest neighbour distance in data
#
pat <- inherits(x,"ppp")
if(pat) w <- x$w
if(checkfields(x,c("x","y"))) {
y <- x$y
x <- x$x
}
M <- if(length(x) > 1) mean(nndist(x,y)) else Mdefault
lrad <- fctr * M / nper
theta <- 2 * pi * (1:nrad)/nrad
cs <- cos(theta)
sn <- sin(theta)
xt <- lrad * as.vector((1:nper) %o% cs)
yt <- lrad * as.vector((1:nper) %o% sn)
xd <- as.vector(outer(x, xt, "+"))
yd <- as.vector(outer(y, yt, "+"))
tmp <- list(x = xd, y = yd)
if(pat) return(as.ppp(tmp,W=w)[,w]) else return(tmp)
}
# concatenate any number of list(x,y) into a list(x,y)
concatxy <- function(...) {
x <- unlist(lapply(list(...), function(w) {w$x}))
y <- unlist(lapply(list(...), function(w) {w$y}))
if(length(x) != length(y))
stop("Internal error: lengths of x and y unequal")
return(list(x=x,y=y))
}
#------------------------------------------------------------
default.dummy <- function(X, nd=NULL, random=FALSE, ntile=NULL, ..., verbose=FALSE) {
# default action to create dummy points.
# regular grid of nd[1] * nd[2] points
# plus corner points of window frame,
# all clipped to window.
#
X <- as.ppp(X)
win <- X$window
# default dimensions
if(is.null(nd))
nd <- if(!is.null(ntile)) ntile else default.ngrid(X)
if(length(nd) == 1)
nd <- rep(nd, 2)
if(verbose)
cat(paste("dummy point grid", nd[1], "x", nd[2], "\n"))
# make dummy points
dummy <- if(random)
stratrand(win, nd[1], nd[2], 1)
else
tilemiddles(win, nd[1], nd[2])
# add corner points
dummy <- concatxy(dummy, corners(win))
# restrict to window
ok <- inside.owin(dummy$x, dummy$y, win)
if(sum(ok) == 0)
stop("None of the dummy points lies inside the window")
x <- dummy$x[ok]
y <- dummy$y[ok]
X <- ppp(x, y, window=win)
# pass parameters for computing weights
if(is.null(ntile))
ntile <- nd
attr(X, "dummy.parameters") <- list(nd=nd, verbose=verbose)
attr(X, "weight.parameters") <-
append(list(...), list(ntile=ntile, verbose=verbose))
return(X)
}
default.ngrid <- function(X) {
# default dimensions of rectangular grid of dummy points
# for data pattern X
X <- as.ppp(X)
max(30, 10 * ceiling(2 * sqrt(X$n)/10))
}