stienen.R
## stienen.R
##
## Stienen diagram with border correction
##
## $Revision: 1.8 $ $Date: 2015/10/21 09:06:57 $
stienen <- function(X, ..., bg="grey", border=list(bg=NULL)) {
Xname <- short.deparse(substitute(X))
stopifnot(is.ppp(X))
if(npoints(X) <= 1) {
do.call(plot,
resolve.defaults(list(x=Window(X)),
list(...),
list(main=Xname)))
return(invisible(NULL))
}
d <- nndist(X)
b <- bdist.points(X)
Y <- X %mark% d
gp <- union(graphicsPars("symbols"), "lwd")
do.call.plotfun(plot.ppp,
resolve.defaults(list(x=Y[b >= d],
markscale=1),
list(...),
list(bg=bg),
list(main=Xname)),
extrargs=gp)
if(!identical(border, FALSE)) {
if(!is.list(border)) border <- list()
do.call.plotfun(plot.ppp,
resolve.defaults(list(x=Y[b < d],
markscale=1,
add=TRUE),
border,
list(...),
list(bg=bg),
list(cols=grey(0.5), lwd=2)),
extrargs=gp)
}
return(invisible(NULL))
}
stienenSet <- function(X, edge=TRUE) {
stopifnot(is.ppp(X))
nnd <- nndist(X)
if(!edge) {
ok <- bdist.points(X) >= nnd
X <- X[ok]
nnd <- nnd[ok]
}
n <- npoints(X)
if(n == 0) return(emptywindow(Window(X)))
if(n == 1) return(Window(X))
rad <- nnd/2
if(!all(ok <- (rad > 0))) {
eps <- min(rad[ok], shortside(Frame(X)))/100
rad <- pmax(rad, eps)
}
delta <- 2 * pi * max(rad)/128
Z <- disc(rad[1], X[1], delta=delta)
for(i in 2:n) Z <- union.owin(Z, disc(rad[i], X[i], delta=delta))
return(Z)
}