pixellate.R
#
# pixellate.R
#
# $Revision: 1.5 $ $Date: 2011/05/18 08:19:08 $
#
# pixellate convert an object to a pixel image
#
# pixellate.ppp convert a point pattern to a pixel image
# (pixel value = number of points in pixel)
#
# pixellate.owin convert a window to a pixel image
# (pixel value = area of intersection with pixel)
#
pixellate <- function(x, ...) {
UseMethod("pixellate")
}
pixellate.ppp <- function(x, W=NULL, ..., weights=NULL, padzero=FALSE)
{
verifyclass(x, "ppp")
if(!is.null(W))
W <- as.mask(W)
else {
# determine W using as.mask
dotargs <- list(...)
namesargs <- names(dotargs)
matched <- namesargs %in% names(formals(as.mask))
W <- do.call("as.mask", append(list(x$window), dotargs[matched]))
}
if(x$n == 0) {
zeroimage <- as.im(as.double(0), W)
if(padzero) # map NA to 0
zeroimage <- na.handle.im(zeroimage, 0)
return(zeroimage)
}
pixels <- nearest.raster.point(x$x, x$y, W)
nr <- W$dim[1]
nc <- W$dim[2]
if(is.null(weights)) {
ta <- table(row = factor(pixels$row, levels = 1:nr), col = factor(pixels$col,
levels = 1:nc))
} else {
ta <- tapply(weights, list(row = factor(pixels$row, levels = 1:nr),
col = factor(pixels$col, levels=1:nc)), sum)
ta[is.na(ta)] <- 0
}
if(nr == 1 || nc == 1) {
# 1 x 1 image: need to specify xrange, yrange explicitly
out <- im(ta, xrange = W$xrange, yrange = W$yrange, unitname=unitname(W))
} else {
# normal case: use exact xcol, yrow values
out <- im(ta, xcol = W$xcol, yrow = W$yrow, unitname=unitname(W))
}
# clip to window of data
if(!padzero)
out <- out[W, drop=FALSE]
return(out)
}
pixellate.owin <- function(x, W=NULL, ...) {
stopifnot(is.owin(x))
P <- as.polygonal(x)
R <- as.rectangle(x)
if(is.null(W))
W <- R
else if(!is.subset.owin(R, as.rectangle(W)))
stop("W does not cover the domain of x")
W <- as.mask(W, ...)
#
x0 <- W$xrange[1]
y0 <- W$yrange[1]
dx <- W$xstep
dy <- W$ystep
nx <- W$dim[2]
ny <- W$dim[1]
# set up output image (real-valued) and initialise to zero
Z <- as.im(W, value=pi, na.replace=pi)
Z <- eval.im(Z * 0)
# process each component polygon
B <- P$bdry
DUP <- spatstat.options("dupC")
for(i in seq_along(B)) {
PP <- B[[i]]
# transform so that pixels become unit squares
QQ <- affinexypolygon(PP, vec = c(-x0, -y0))
RR <- affinexypolygon(QQ, mat = diag(1/c(dx, dy)))
#
xx <- RR$x
yy <- RR$y
nn <- length(xx)
# close polygon
xx <- c(xx, xx[1])
yy <- c(yy, yy[1])
nn <- nn+1
# call C routine
zz <- .C("poly2imA",
ncol=as.integer(nx),
nrow=as.integer(ny),
xpoly=as.double(xx),
ypoly=as.double(yy),
npoly=as.integer(nn),
out=as.double(numeric(nx * ny)),
status=as.integer(integer(1)),
DUP=DUP,
PACKAGE="spatstat")
if(zz$status != 0)
stop("Internal error")
# increment output image
Z$v <- Z$v + matrix(zz$out, ny, nx, byrow=TRUE)
}
# revert to original scale
pixelarea <- dx * dy
Z <- eval.im(Z * pixelarea)
return(Z)
}