#
# eval.im.R
#
# eval.im() Evaluate expressions involving images
#
# compatible.im() Check whether two images are compatible
#
# harmonise.im() Harmonise images
# commonGrid()
#
# $Revision: 1.27 $ $Date: 2012/02/04 07:43:19 $
#
eval.im <- function(expr, envir) {
e <- as.expression(substitute(expr))
# get names of all variables in the expression
varnames <- all.vars(e)
allnames <- all.names(e, unique=TRUE)
funnames <- allnames[!(allnames %in% varnames)]
if(length(varnames) == 0)
stop("No variables in this expression")
# get the values of the variables
if(missing(envir))
envir <- sys.parent()
vars <- lapply(as.list(varnames), function(x, e) get(x, envir=e), e=envir)
names(vars) <- varnames
funs <- lapply(as.list(funnames), function(x, e) get(x, envir=e), e=envir)
names(funs) <- funnames
# find out which variables are images
ims <- unlist(lapply(vars, is.im))
if(!any(ims))
stop("No images in this expression")
images <- vars[ims]
nimages <- length(images)
# test that the images are compatible
if(!(ok <- do.call("compatible", unname(images))))
stop(paste(if(nimages > 2) "some of" else NULL,
"the images",
commasep(sQuote(names(images))),
"are not compatible"))
# replace each image by its matrix of pixel values, and evaluate
getvalues <- function(x) {
v <- as.matrix(x)
dim(v) <- NULL
return(v)
}
imagevalues <- lapply(images, getvalues)
template <- images[[1]]
# This bit has been repaired:
vars[ims] <- imagevalues
v <- eval(e, append(vars, funs))
#
# reshape, etc
result <- im(v, template$xcol, template$yrow,
unitname=unitname(template))
return(result)
}
compatible.im <- function(A, B, ..., tol=1e-6) {
verifyclass(A, "im")
if(missing(B)) return(TRUE)
verifyclass(B, "im")
xdiscrep <- max(abs(A$xrange - B$xrange),
abs(A$xstep - B$xstep),
abs(A$xcol - B$xcol))
ydiscrep <- max(abs(A$yrange - B$yrange),
abs(A$ystep - B$ystep),
abs(A$yrow - B$yrow))
xok <- (xdiscrep < tol * min(A$xstep, B$xstep))
yok <- (ydiscrep < tol * min(A$ystep, B$ystep))
uok <- compatible.units(unitname(A), unitname(B))
if(!(xok && yok && uok))
return(FALSE)
# A and B are compatible
if(length(list(...)) == 0)
return(TRUE)
# recursion
return(compatible.im(B, ..., tol=tol))
}
# force a list of images to be compatible
harmonize.im <- harmonise.im <- function(...) {
argz <- list(...)
n <- length(argz)
if(n < 2) return(argz)
result <- vector(mode="list", length=n)
isim <- unlist(lapply(argz, is.im))
if(!any(isim))
stop("No images supplied")
imgs <- argz[isim]
# if any windows are present, extract bounding box
iswin <- unlist(lapply(argz, is.owin))
bb0 <- if(!any(iswin)) NULL else do.call("bounding.box", unname(argz[iswin]))
if(length(imgs) == 1 && is.null(bb0)) {
# only one 'true' image: use it as template.
result[isim] <- imgs
Wtemplate <- imgs[[1]]
} else {
# test for compatible units
un <- lapply(imgs, unitname)
uok <- unlist(lapply(un, compatible.units, y=un[[1]]))
if(!all(uok))
stop("Images have incompatible units of length")
# find the image with the highest resolution
xsteps <- unlist(lapply(imgs, function(a) { a$xstep }))
which.finest <- which.min(xsteps)
finest <- imgs[[which.finest]]
# get the bounding box
bb <- do.call("bounding.box", lapply(unname(imgs), as.rectangle))
if(!is.null(bb0)) bb <- bounding.box(bb, bb0)
# determine new raster coordinates
xcol <- prolongseq(finest$xcol, bb$xrange)
yrow <- prolongseq(finest$yrow, bb$yrange)
xy <- list(x=xcol, y=yrow)
# resample all images on new raster
newimgs <- lapply(imgs, as.im, xy=xy)
result[isim] <- newimgs
Wtemplate <- newimgs[[which.finest]]
}
# convert other data to images
if(any(notim <- !isim))
result[notim] <- lapply(argz[notim], as.im, W=as.mask(Wtemplate))
names(result) <- names(argz)
return(result)
}
# Return just the corresponding template window
commonGrid <- local({
# auxiliary function
gettype <- function(x) {
if(is.im(x) || is.mask(x)) "raster" else
if(is.owin(x) || is.ppp(x) || is.psp(x)) "spatial" else
"none"
}
commonGrid <- function(...) {
argz <- list(...)
type <- unlist(lapply(argz, gettype))
israster <- (type == "raster")
haswin <- (type != "none")
if(any(israster)) {
# Get raster data
rasterlist <- argz[israster]
} else {
# No existing raster data - apply default resolution
if(!any(haswin))
stop("No spatial data supplied")
wins <- lapply(argz[haswin], as.owin)
rasterlist <- lapply(wins, as.mask)
}
# Find raster object with finest resolution
if(length(rasterlist) == 1) {
# only one raster object
finest <- rasterlist[[1]]
} else {
# test for compatible units
un <- lapply(rasterlist, unitname)
uok <- unlist(lapply(un, compatible.units, y=un[[1]]))
if(!all(uok))
stop("Objects have incompatible units of length")
# find the image/mask with the highest resolution
xsteps <- unlist(lapply(rasterlist, function(a) { a$xstep }))
which.finest <- which.min(xsteps)
finest <- rasterlist[[which.finest]]
}
# determine the bounding box
bb <- do.call("bounding.box", lapply(unname(argz[haswin]), as.rectangle))
# determine new raster coordinates
xcol <- prolongseq(finest$xcol, bb$xrange)
yrow <- prolongseq(finest$yrow, bb$yrange)
xy <- list(x=xcol, y=yrow)
# generate template
Wtemplate <- as.mask(bb, xy=xy)
return(Wtemplate)
}
commonGrid
})