https://github.com/cran/spatstat
Tip revision: 4712eff0e57ddd871eee94f8a5fd24be181fdf44 authored by Adrian Baddeley on 17 January 2013, 07:33:07 UTC
version 1.31-0
version 1.31-0
Tip revision: 4712eff
Kmeasure.R
#
# Kmeasure.R
#
# $Revision: 1.41 $ $Date: 2011/04/19 00:54:13 $
#
# Kmeasure() compute an estimate of the second order moment measure
#
# Kest.fft() use Kmeasure() to form an estimate of the K-function
#
# second.moment.calc() underlying algorithm
#
Kmeasure <- function(X, sigma, edge=TRUE, ..., varcov=NULL) {
stopifnot(is.ppp(X))
sigma.given <- !missing(sigma) && !is.null(sigma)
varcov.given <- !is.null(varcov)
ngiven <- sigma.given + varcov.given
if(ngiven == 2)
stop(paste("Give only one of the arguments",
sQuote("sigma"), "and", sQuote("varcov")))
if(ngiven == 0)
stop(paste("Please specify smoothing bandwidth", sQuote("sigma"),
"or", sQuote("varcov")))
if(varcov.given) {
stopifnot(is.matrix(varcov) && nrow(varcov) == 2 && ncol(varcov)==2 )
sigma <- NULL
} else {
stopifnot(is.numeric(sigma))
stopifnot(length(sigma) %in% c(1,2))
stopifnot(all(sigma > 0))
if(length(sigma) == 2) {
varcov <- diag(sigma^2)
sigma <- NULL
}
}
second.moment.calc(X, sigma=sigma, edge, "Kmeasure", varcov=varcov)
}
second.moment.calc <- function(x, sigma=NULL, edge=TRUE,
what="Kmeasure", debug=FALSE, ...,
varcov=NULL, expand=FALSE) {
if(is.null(sigma) && is.null(varcov))
stop("must specify sigma or varcov")
choices <- c("kernel", "smooth", "Kmeasure", "Bartlett", "edge", "all", "smoothedge")
if(!(what %in% choices))
stop(paste("Unknown choice: what = ", sQuote(what),
"; available options are:",
paste(sQuote(choices), collapse=", ")))
sig <- if(!is.null(sigma)) sigma else max(c(diag(varcov), sqrt(det(varcov))))
win <- rescue.rectangle(as.owin(x))
rec <- as.rectangle(win)
across <- min(diff(rec$xrange), diff(rec$yrange))
# determine whether to expand window
if(!expand || (6 * sig < across))
result <- second.moment.engine(x, sigma=sigma, edge=edge,
what=what, debug=debug, ..., varcov=varcov)
else {
# need to expand window
bigger <- grow.rectangle(rec, (7 * sig - across)/2)
if(is.ppp(x)) {
# pixellate first (to preserve pixel resolution)
X <- pixellate(x, ..., padzero=TRUE)
np <- npoints(x)
} else if(is.im(x)) {
X <- x
np <- NULL
} else stop("x should be an image or a point pattern")
# now expand
X <- rebound.im(X, bigger)
X <- na.handle.im(X, 0)
out <- second.moment.engine(X, sigma=sigma, edge=edge,
what=what, debug=debug, ...,
obswin=win, varcov=varcov, npts=np)
# now clip it
fbox <- shift(rec, origin="midpoint")
result <- switch(what,
kernel = out[fbox],
smooth = out[win],
Kmeasure = out[fbox],
Bartlett = out[fbox],
edge = out[win],
smoothedge = list(smooth=out$smooth[win],
edge =out$edge[win]),
all =
list(kernel=out$kernel[fbox],
smooth=out$smooth[win],
Kmeasure=out$Kmeasure[fbox],
Bartlett=out$Bartlett[fbox],
edge=out$edge[win]))
}
return(result)
}
second.moment.engine <- function(x, sigma=NULL, edge=TRUE,
what="Kmeasure", debug=FALSE, ...,
obswin = as.owin(x), varcov=NULL,
npts=NULL)
{
stopifnot(what %in% c("kernel", "smooth", "Kmeasure",
"Bartlett", "edge", "all", "smoothedge"))
if(is.ppp(x))
# convert list of points to mass distribution
X <- pixellate(x, ..., padzero=TRUE)
else if(is.im(x))
X <- x
else stop("internal error: unrecognised format for x")
# ensure obswin has same bounding frame as X
if(!missing(obswin))
obswin <- rebound.owin(obswin, as.rectangle(X))
#
if(is.null(npts) && is.ppp(x))
npts <- npoints(x)
# go to work
Y <- X$v
xw <- X$xrange
yw <- X$yrange
# pad with zeroes
nr <- nrow(Y)
nc <- ncol(Y)
Ypad <- matrix(0, ncol=2*nc, nrow=2*nr)
Ypad[1:nr, 1:nc] <- Y
lengthYpad <- 4 * nc * nr
# corresponding coordinates
xw.pad <- xw[1] + 2 * c(0, diff(xw))
yw.pad <- yw[1] + 2 * c(0, diff(yw))
xcol.pad <- X$xcol[1] + X$xstep * (0:(2*nc-1))
yrow.pad <- X$yrow[1] + X$ystep * (0:(2*nr-1))
# set up Gauss kernel
xcol.G <- X$xstep * c(0:(nc-1),-(nc:1))
yrow.G <- X$ystep * c(0:(nr-1),-(nr:1))
xx <- matrix(xcol.G[col(Ypad)], ncol=2*nc, nrow=2*nr)
yy <- matrix(yrow.G[row(Ypad)], ncol=2*nc, nrow=2*nr)
if(!is.null(sigma)) {
# if(max(abs(diff(xw)),abs(diff(yw))) < 6 * sigma)
# warning("sigma is too large for this window")
Kern <- exp(-(xx^2 + yy^2)/(2 * sigma^2))/(2 * pi * sigma^2) * X$xstep * X$ystep
} else if(!is.null(varcov)) {
# anisotropic kernel
detSigma <- det(varcov)
Sinv <- solve(varcov)
const <- X$xstep * X$ystep/(2 * pi * sqrt(detSigma))
Kern <- const * exp(-(xx * (xx * Sinv[1,1] + yy * Sinv[1,2])
+ yy * (xx * Sinv[2,1] + yy * Sinv[2,2]))/2)
} else
stop("Must specify either sigma or varcov")
# these options call for several image outputs
if(what %in% c("all", "smoothedge"))
result <- list()
if(what %in% c("kernel", "all")) {
# return the kernel
# first rearrange it into spatially sensible order (monotone x and y)
rtwist <- ((-nr):(nr-1)) %% (2 * nr) + 1
ctwist <- (-nc):(nc-1) %% (2*nc) + 1
if(debug) {
if(any(fave.order(xcol.G) != rtwist))
cat("something round the twist\n")
}
Kermit <- Kern[ rtwist, ctwist]
ker <- im(Kermit, xcol.G[ctwist], yrow.G[ rtwist], unitname=unitname(x))
if(what == "kernel")
return(ker)
else
result$kernel <- ker
}
# convolve using fft
fK <- fft(Kern)
if(what != "edge") {
fY <- fft(Ypad)
sm <- fft(fY * fK, inverse=TRUE)/lengthYpad
if(debug) {
cat(paste("smooth: maximum imaginary part=",
signif(max(Im(sm)),3), "\n"))
if(!is.null(npts))
cat(paste("smooth: mass error=",
signif(sum(Mod(sm))-npts,3), "\n"))
}
}
if(what %in% c("smooth", "all", "smoothedge")) {
# return the smoothed point pattern without edge correction
smo <- im(Re(sm)[1:nr, 1:nc], xcol.pad[1:nc], yrow.pad[1:nr],
unitname=unitname(x))
if(what == "smooth")
return(smo)
else
result$smooth <- smo
}
if(what != "edge")
bart <- Mod(fY)^2 * fK
if(what %in% c("Bartlett", "all")) {
# return Bartlett spectrum
# rearrange into spatially sensible order (monotone x and y)
rtwist <- ((-nr):(nr-1)) %% (2 * nr) + 1
ctwist <- (-nc):(nc-1) %% (2*nc) + 1
Bart <- bart[ rtwist, ctwist]
Bartlett <- im(Mod(Bart),(-nc):(nc-1), (-nr):(nr-1))
if(what == "Bartlett")
return(Bartlett)
else
result$Bartlett <- Bartlett
}
#### ------- Second moment measure --------------
#
if(what != "edge") {
mom <- fft(bart, inverse=TRUE)/lengthYpad
if(debug) {
cat(paste("2nd moment measure: maximum imaginary part=",
signif(max(Im(mom)),3), "\n"))
if(!is.null(npts))
cat(paste("2nd moment measure: mass error=",
signif(sum(Mod(mom))-npts^2, 3), "\n"))
}
mom <- Mod(mom)
# subtract (delta_0 * kernel) * npts
if(is.null(npts))
stop("Internal error: second moment measure requires npts")
mom <- mom - npts* Kern
}
# edge correction
if(edge || what %in% c("edge", "all", "smoothedge")) {
# compute kernel-smoothed set covariance
M <- as.mask(obswin, dimyx=c(nr, nc))$m
# previous line ensures M has same dimensions and scale as Y
Mpad <- matrix(0, ncol=2*nc, nrow=2*nr)
Mpad[1:nr, 1:nc] <- M
lengthMpad <- 4 * nc * nr
fM <- fft(Mpad)
if(edge && what != "edge") {
co <- fft(Mod(fM)^2 * fK, inverse=TRUE)/lengthMpad
co <- Mod(co)
a <- sum(M)
wt <- a/co
me <- spatstat.options("maxedgewt")
weight <- matrix(pmin(me, wt), ncol=2*nc, nrow=2*nr)
# if(debug) browser()
#
# apply edge correction
mom <- mom * weight
# set to NA outside 'reasonable' region
mom[wt > 10] <- NA
}
}
if(what != "edge") {
# rearrange into spatially sensible order (monotone x and y)
rtwist <- ((-nr):(nr-1)) %% (2 * nr) + 1
ctwist <- (-nc):(nc-1) %% (2*nc) + 1
mom <- mom[ rtwist, ctwist]
if(debug) {
if(any(fave.order(xcol.G) != rtwist))
cat("something round the twist\n")
}
}
if(what %in% c("edge", "all", "smoothedge")) {
# return convolution of window with kernel
# (evaluated inside window only)
con <- fft(fM * fK, inverse=TRUE)/lengthMpad
edg <- Mod(con[1:nr, 1:nc])
edg <- im(edg, xcol.pad[1:nc], yrow.pad[1:nr], unitname=unitname(x))
if(what == "edge")
return(edg)
else
result$edge <- edg
}
if(what == "smoothedge")
return(result)
# Second moment measure, density estimate
# Divide by number of points * lambda and convert mass to density
pixarea <- with(X, xstep * ystep)
mom <- mom * area.owin(obswin) / (pixarea * npts^2)
# this is the second moment measure
mm <- im(mom, xcol.G[ctwist], yrow.G[rtwist], unitname=unitname(x))
if(what == "Kmeasure")
return(mm)
else
result$Kmeasure <- mm
# what = "all", so return all computed objects
return(result)
}
Kest.fft <- function(X, sigma, r=NULL, breaks=NULL) {
verifyclass(X, "ppp")
W <- X$window
lambda <- X$n/area.owin(W)
rmaxdefault <- rmax.rule("K", W, lambda)
bk <- handle.r.b.args(r, breaks, W, rmaxdefault=rmaxdefault)
breaks <- bk$val
rvalues <- bk$r
u <- Kmeasure(X, sigma)
xx <- rasterx.im(u)
yy <- rastery.im(u)
rr <- sqrt(xx^2 + yy^2)
tr <- whist(rr, breaks, u$v)
K <- cumsum(tr)
rmax <- min(rr[is.na(u$v)])
K[rvalues >= rmax] <- NA
result <- data.frame(r=rvalues, theo=pi * rvalues^2, border=K)
w <- X$window
alim <- c(0, min(diff(w$xrange), diff(w$yrange))/4)
out <- fv(result,
"r", substitute(K[fft](r), NULL),
"border",
. ~ r, alim,
c("r", "{%s^{pois}}(r)", "{hat(%s)^{bord}}(r)"),
c("distance argument r",
"theoretical Poisson %s",
"border-corrected estimate of %s"),
fname="K[fft]",
unitname=unitname(X)
)
return(out)
}