window.S
#
# window.S
#
# A class 'owin' to define the "observation window"
#
# $Revision: 4.67 $ $Date: 2007/03/16 08:59:46 $
#
#
# A window may be either
#
# - rectangular:
# a rectangle in R^2
# (with sides parallel to the coordinate axes)
#
# - polygonal:
# delineated by one or more non-self-intersecting
# polygons, possibly including polygonal holes.
#
# - digital mask:
# defined by a binary image
# whose pixel values are TRUE wherever the pixel
# is inside the window
#
# Any window is an object of class 'owin',
# containing at least the following entries:
#
# $type: a string ("rectangle", "polygonal" or "mask")
#
# $xrange
# $yrange
# vectors of length 2 giving the real dimensions
# of the enclosing box.
#
# The 'rectangle' type has only these entries.
#
# The 'polygonal' type has an additional entry
#
# $bdry
# a list of polygons.
# Each entry bdry[[i]] determines a closed polygon.
#
# bdry[[i]] has components $x and $y which are
# the cartesian coordinates of the vertices of
# the i-th boundary polygon (without repetition of
# the first vertex, i.e. same convention as in the
# plotting function polygon().)
#
#
# The 'mask' type has entries
#
# $m logical matrix
# $dim its dimension array
# $xstep,ystep x and y dimensions of a pixel
# $xcol vector of x values for each column
# $yrow vector of y values for each row
#
# (the row index corresponds to increasing y coordinate;
# the column index " " " " " " x " " ".)
#
#
#-----------------------------------------------------------------------------
#
.Spatstat.Image.Warning <-
c("Row index corresponds to increasing y coordinate; column to increasing x",
"Transpose matrices to get the standard presentation in R",
"Example: image(result$xcol,result$yrow,t(result$d))")
owin <- function(xrange=c(0,1), yrange=c(0,1),
..., poly=NULL, mask=NULL, units=NULL) {
units <- as.units(units)
## Exterminate ambiguities
poly.given <- !is.null(poly)
mask.given <- !is.null(mask)
if(poly.given && mask.given)
stop("Ambiguous -- both polygonal boundary and digital mask supplied")
if(missing(xrange) != missing(yrange))
stop("If one of xrange, yrange is specified then both must be.")
# avoid re-checking if already checked
check <- resolve.defaults(list(...), list(check=TRUE))$check
if(!check) {
area.w <- resolve.defaults(list(...), list(area=NULL))$area
if(!poly.given && !mask.given)
w <- list(type="rectangle",
xrange=xrange, yrange=yrange, units=units,
area=
if(!is.null(area.w)) area.w else (diff(xrange)*diff(yrange)))
else if(poly.given)
w <- list(type="polygonal",
xrange=xrange, yrange=yrange, bdry=poly, units=units,
area=area.w)
else if(mask.given) {
nc <- ncol(mask)
nr <- nrow(mask)
xstep <- diff(xrange)/nc
ystep <- diff(yrange)/nr
xcol <- seq(xrange[1]+xstep/2, xrange[2]-xstep/2, length=nc)
yrow = seq(yrange[1]+ystep/2, yrange[2]-ystep/2, length=nr)
w <- list(type = "mask",
xrange = xrange,
yrange = yrange,
dim = c(nr, nc),
xstep = xstep,
ystep = ystep,
warnings = .Spatstat.Image.Warning,
xcol = xcol,
yrow = yrow,
m = mask,
units = units,
area = area.w)
}
class(w) <- "owin"
return(w)
}
if(!poly.given && !mask.given) {
######### rectangle #################
if(!is.vector(xrange) || length(xrange) != 2 || xrange[2] < xrange[1])
stop("xrange should be a vector of length 2 giving (xmin, xmax)")
if(!is.vector(yrange) || length(yrange) != 2 || yrange[2] < yrange[1])
stop("yrange should be a vector of length 2 giving (ymin, ymax)")
w <- list(type="rectangle", xrange=xrange, yrange=yrange, units=units,
area=diff(xrange) * diff(yrange))
class(w) <- "owin"
return(w)
} else if(poly.given) {
######### polygonal boundary ########
#
# test whether it's a single polygon or multiple polygons
if(verify.xypolygon(poly, fatal=FALSE))
psingle <- TRUE
else if(all(unlist(lapply(poly, verify.xypolygon, fatal=FALSE))))
psingle <- FALSE
else
stop("poly must be either a list(x,y) or a list of list(x,y)")
if(psingle) {
# single boundary polygon
bdry <- list(poly)
w.area <- area.xypolygon(poly)
if(w.area < 0)
stop("Area of polygon is negative - maybe traversed in wrong direction?")
} else {
# multiple boundary polygons
bdry <- poly
w.area <- sum(unlist(lapply(poly, area.xypolygon)))
if(w.area < 0)
stop(paste("Area of window is negative;\n",
"check that all polygons were traversed in the right direction"))
}
actual.xrange <- range(unlist(lapply(bdry, function(a) a$x)))
if(missing(xrange))
xrange <- actual.xrange
else {
if(!is.vector(xrange) || length(xrange) != 2 || xrange[2] <= xrange[1])
stop("xrange should be a vector of length 2 giving (xmin, xmax)")
if(!all(xrange == range(c(xrange, actual.xrange))))
stop("polygon's x coordinates outside xrange")
}
actual.yrange <- range(unlist(lapply(bdry, function(a) a$y)))
if(missing(yrange))
yrange <- actual.yrange
else {
if(!is.vector(yrange) || length(yrange) != 2 || yrange[2] <= yrange[1])
stop("yrange should be a vector of length 2 giving (ymin, ymax)")
if(!all(yrange == range(c(yrange, actual.yrange))))
stop("polygon's y coordinates outside yrange")
}
w <- list(type="polygonal",
xrange=xrange, yrange=yrange, bdry=bdry, units=units,
area=w.area)
class(w) <- "owin"
# check for intersection or self-intersection
if(check) {
ok <- owinpolycheck(w)
if(!ok) {
errors <- attr(ok, "err")
stop(paste("Polygon data contain", commasep(errors)))
}
}
return(w)
} else if(mask.given) {
######### digital mask #####################
if(!is.matrix(mask))
stop(paste(sQuote("mask"), "must be a matrix"))
if(!is.logical(mask))
stop(paste("The entries of", sQuote("mask"), "must be logical"))
nc <- ncol(mask)
nr <- nrow(mask)
if(missing(xrange) && missing(yrange)) {
# take pixels to be 1 x 1 unit
xrange <- c(0,nc)
yrange <- c(0,nr)
} else {
if(!is.vector(xrange) || length(xrange) != 2 || xrange[2] <= xrange[1])
stop("xrange should be a vector of length 2 giving (xmin, xmax)")
if(!is.vector(yrange) || length(yrange) != 2 || yrange[2] <= yrange[1])
stop("yrange should be a vector of length 2 giving (ymin, ymax)")
}
xstep <- diff(xrange)/nc
ystep <- diff(yrange)/nr
w.area <- sum(mask) * xstep * ystep
out <- list(type = "mask",
xrange = xrange,
yrange = yrange,
dim = c(nr, nc),
xstep = xstep,
ystep = ystep,
warnings = .Spatstat.Image.Warning,
xcol = seq(xrange[1]+xstep/2, xrange[2]-xstep/2, length=nc),
yrow = seq(yrange[1]+ystep/2, yrange[2]-ystep/2, length=nr),
m = mask,
units = units,
area = w.area)
class(out) <- "owin"
return(out)
}
# never reached
NULL
}
#
#-----------------------------------------------------------------------------
#
is.owin <- function(x) { inherits(x, "owin") }
#
#-----------------------------------------------------------------------------
#
as.owin <- function(W, ..., fatal=TRUE) {
UseMethod("as.owin")
}
as.owin.owin <- function(W, ..., fatal=TRUE) {
if(verifyclass(W, "owin", fatal=fatal))
return(owin(W$xrange, W$yrange, poly=W$bdry, mask=W$m, units=W$units, check=FALSE, area=W$area))
else
return(NULL)
}
as.owin.ppp <- function(W, ..., fatal=TRUE) {
if(verifyclass(W, "ppp", fatal=fatal))
return(W$window)
else
return(NULL)
}
as.owin.im <- function(W, ..., fatal=TRUE) {
if(!verifyclass(W, "im", fatal=fatal))
return(NULL)
out <- list(type = "mask",
xrange = W$xrange,
yrange = W$yrange,
dim = W$dim,
xstep = W$xstep,
ystep = W$ystep,
warnings = .Spatstat.Image.Warning,
xcol = W$xcol,
yrow = W$yrow,
m = !is.na(W$v),
units = W$units)
class(out) <- "owin"
return(out)
}
as.owin.ppm <- function(W, ..., fatal=TRUE) {
if(!verifyclass(W, "ppm", fatal=fatal))
return(NULL)
as.owin(data.ppm(W))
}
as.owin.default <- function(W, ..., fatal=TRUE) {
# Tries to interpret data as an object of class 'window'
# W may be
# a structure with entries xrange, yrange
# a four-element vector (interpreted xmin, xmax, ymin, ymax)
# a structure with entries xl, xu, yl, yu
if(checkfields(W, c("xrange", "yrange"))) {
Z <- owin(W$xrange, W$yrange)
return(Z)
} else if(is.vector(W) && is.numeric(W) && length(W) == 4) {
Z <- owin(W[1:2], W[3:4])
return(Z)
} else if(checkfields(W, c("xl", "xu", "yl", "yu"))) {
W <- as.list(W)
Z <- owin(c(W$xl, W$xu),c(W$yl, W$yu))
return(Z)
} else if(checkfields(W, c("x", "y", "area"))
&& checkfields(W$area, c("xl", "xu", "yl", "yu"))) {
V <- as.list(W$area)
Z <- owin(c(V$xl, V$xu),c(V$yl, V$yu))
return(Z)
} else if(fatal)
stop("Can't interpret W as a window")
else return(NULL)
}
#
#-----------------------------------------------------------------------------
#
#
as.rectangle <- function(...) {
w <- as.owin(...)
return(owin(w$xrange, w$yrange, units=w$units))
}
#
#-----------------------------------------------------------------------------
#
as.mask <- function(w, eps=NULL, dimyx=NULL, xy=NULL) {
# eps: grid mesh (pixel) size
# dimyx: dimensions of pixel raster
# xy: coordinates of pixel raster
#
# NOTE: rows <=> y coordinate
if(!missing(w)) {
verifyclass(w, "owin")
units <- units(w)
} else {
units <- as.units(NULL)
if(is.null(xy))
stop("If w is missing, xy is required")
}
# If it's already a mask, and no other arguments specified,
# just return it.
if(!missing(w) && w$type == "mask" &&
is.null(eps) && is.null(dimyx) && is.null(xy))
return(w)
#
# First create pixel array
#
if(is.null(xy)) {
#
# Pixel coordinates to be computed from other dimensions
#
# First determine row & column dimensions
# output raster dimensions 'dimyx'
if(!is.null(dimyx)) {
if(length(dimyx) == 1)
dimyx <- rep(dimyx, 2)
nr <- dimyx[1]
nc <- dimyx[2]
} else {
# use pixel size 'eps'
if(!is.null(eps)) {
nr <- ceiling(diff(w$yrange)/eps)
nc <- ceiling(diff(w$xrange)/eps)
} else {
# default to spatstat defaults
np <- spatstat.options("npixel")
if(!is.numeric(np) || length(np) > 2)
stop(paste("Illegal value for spatstat.options(",
dQuote("npixel"), ")"))
if(length(np) == 1)
nr <- nc <- np[1]
else {
nr <- np[2]
nc <- np[1]
}
}
}
# Initial mask with all entries True
out <- owin(w$xrange, w$yrange,
mask=matrix(TRUE, nrow=nr, ncol=nc))
} else {
#
# Pixel coordinates given explicitly:
#
if(!checkfields(xy, c("x", "y")))
stop(paste(sQuote("xy"),
"should be a list containing two vectors x and y"))
x <- sort(unique(xy$x))
y <- sort(unique(xy$y))
# derive other parameters
nr <- length(y)
nc <- length(x)
# x and y pixel sizes
dx <- diff(x)
if(diff(range(dx)) > 0.01 * mean(dx))
stop("x coordinates must be evenly spaced")
xstep <- mean(dx)
dy <- diff(y)
if(diff(range(dy)) > 0.01 * mean(dy))
stop("y coordinates must be evenly spaced")
ystep <- mean(dy)
# coordinate ranges
xrange <- range(x) + xstep * c(-1,1)/2
yrange <- range(y) + ystep * c(-1,1)/2
# Initial mask with all entries TRUE
m <- matrix(TRUE, nrow=nr, ncol=nc)
out <- list(type = "mask",
xrange = xrange,
yrange = yrange,
dim = c(nr, nc),
xstep = xstep,
ystep = ystep,
warnings = .Spatstat.Image.Warning,
xcol = x,
yrow = y,
m = m)
class(out) <- "owin"
# window may be implicit in this case.
if(missing(w))
w <- owin(xrange, yrange)
}
# Second, intersect pixel raster with existing window
switch(w$type,
rectangle = { },
mask = {
phase <- c((out$xcol[1] - w$xcol[1])/w$xstep,
(out$yrow[1] - w$yrow[1])/w$ystep)
out$m <- matrixsample(w$m, c(nr, nc), phase=round(phase))
},
polygonal = {
# test every pixel
x <- as.vector(raster.x(out))
y <- as.vector(raster.y(out))
value <- inside.owin(x, y, w)
out$m <- matrix(value, nrow=nr, ncol=nc)
})
out$units <- units
return(out)
}
#
#-----------------------------------------------------------------------------
#
as.polygonal <- function(W) {
# the only sensible use of this function is to
# convert a rectangle to a polygon
verifyclass(W, "owin")
switch(W$type,
rectangle = {
xr <- W$xrange
yr <- W$yrange
return(owin(xr, yr, poly=list(x=xr[c(1,2,2,1)],y=yr[c(1,1,2,2)]),
units=units(W)))
},
polygonal = {
return(W)
},
mask = {
stop("A mask cannot be converted to a polygon")
}
)
}
#
# ----------------------------------------------------------------------
validate.mask <- function(w, fatal=TRUE) {
verifyclass(w, "owin", fatal=fatal)
if(w$type == "mask")
return(TRUE)
if(fatal)
stop(paste(deparse(substitute(w)), "is not a binary mask"))
else {
warning(paste(deparse(substitute(w)), "is not a binary mask"))
return(FALSE)
}
}
raster.x <- function(w) {
validate.mask(w)
di <- w$dim
nr <- di[1]
nc <- di[2]
m <- matrix(0, nrow=nr, ncol=nc)
matrix(w$xcol[col(m)], nrow=nr, ncol=nc)
}
raster.y <- function(w) {
validate.mask(w)
di <- w$dim
nr <- di[1]
nc <- di[2]
m <- matrix(0, nrow=nr, ncol=nc)
matrix(w$yrow[row(m)], nrow=nr, ncol=nc)
}
nearest.raster.point <- function(x,y,w, indices=TRUE) {
validate.mask(w)
nr <- w$dim[1]
nc <- w$dim[2]
if(length(x) == 0)
cc <- rr <- integer(0)
else {
cc <- 1 + round((x - w$xcol[1])/w$xstep)
rr <- 1 + round((y - w$yrow[1])/w$ystep)
cc <- pmax(1,pmin(cc, nc))
rr <- pmax(1,pmin(rr, nr))
}
if(indices)
return(list(row=rr, col=cc))
else
return(list(x=w$xcol[cc], y=w$yrow[rr]))
}
#------------------------------------------------------------------
bounding.box <- function(...) {
wins <- list(...)
if(length(wins) == 0)
stop("No arguments supplied")
# trap a particular misuse of this function
if(length(wins) == 2) {
w1 <- wins[[1]]
w2 <- wins[[2]]
if(is.vector(w1) && is.numeric(w1) &&
is.vector(w2) && is.numeric(w2) &&
length(w1) == length(w2))
stop(paste("bounding.box() was applied to two numeric vectors;",
"you probably wanted bounding.box.xy()"))
}
if(length(wins) > 1) {
# multiple arguments -- compute bounding box for each argument.
# First trap any point patterns and extract bounding boxes of points
isppp <- unlist(lapply(wins, is.ppp))
if(any(isppp))
wins[isppp] <- lapply(wins[isppp], bounding.box.xy)
# then take bounding box of each window
boxes <- lapply(wins, bounding.box)
# discard NULL values
isnull <- unlist(lapply(boxes, is.null))
boxes <- boxes[!isnull]
# take bounding box of these boxes
xrange <- range(unlist(lapply(boxes, function(b){b$xrange})))
yrange <- range(unlist(lapply(boxes, function(b){b$yrange})))
return(owin(xrange, yrange))
}
# single argument
w <- wins[[1]]
if(is.null(w))
return(NULL)
# point pattern?
if(is.ppp(w))
return(bounding.box.xy(w))
# convert to window
w <- as.owin(w)
# determine a tight bounding box for the window w
switch(w$type,
rectangle = {
return(w)
},
polygonal = {
bdry <- w$bdry
xr <- range(unlist(lapply(bdry, function(a) a$x)))
yr <- range(unlist(lapply(bdry, function(a) a$y)))
return(owin(xr, yr))
},
mask = {
m <- w$m
x <- raster.x(w)
y <- raster.y(w)
xr <- range(x[m]) + c(-1,1) * w$xstep/2
yr <- range(y[m]) + c(-1,1) * w$ystep/2
return(owin(xr, yr))
},
stop("unrecognised window type", w$type)
)
}
complement.owin <- function(w, frame=as.rectangle(w)) {
wname <- deparse(substitute(w))
verifyclass(w, "owin")
if(reframe <- !missing(frame)) {
verifyclass(frame, "owin")
w <- rebound.owin(w, frame)
# if w was a rectangle, it's now polygonal
}
switch(w$type,
mask = {
w$m <- !(w$m)
},
rectangle = {
stop("window is a rectangle - its complement is empty")
},
polygonal = {
bdry <- w$bdry
# bounding box, in anticlockwise order
box <- list(x=w$xrange[c(1,2,2,1)],
y=w$yrange[c(1,1,2,2)])
boxarea <- area.xypolygon(box)
# first check whether one of the current boundary polygons
# is the bounding box itself (with + sign)
if(reframe)
is.box <- rep(FALSE, length(bdry))
else {
nvert <- unlist(lapply(bdry, function(a) { length(a$x) }))
area <- unlist(lapply(bdry, area.xypolygon))
boxarea.mineps <- boxarea * (0.99999)
is.box <- (nvert == 4 & area >= boxarea.mineps)
if(sum(is.box) > 1)
stop("Internal error: multiple copies of bounding box")
if(all(is.box))
stop(paste("Complement is empty, for window", wname),
call.=FALSE)
}
# if box is present (with + sign), remove it
if(any(is.box))
bdry <- bdry[!is.box]
# now reverse the direction of each polygon
bdry <- lapply(bdry, reverse.xypolygon)
# if box was absent, add it
if(!any(is.box))
bdry <- c(bdry, list(box)) # sic
# put back into w
w$bdry <- bdry
},
stop("unrecognised window type", w$type)
)
return(w)
}
#-----------------------------------------------------------
inside.owin <- function(x, y, w) {
# test whether (x,y) is inside window w
# x, y may be vectors
verifyclass(w, "owin")
if(length(x)==0)
return(logical(0))
# test whether inside bounding rectangle
xr <- w$xrange
yr <- w$yrange
eps <- sqrt(.Machine$double.eps)
frameok <- (x >= xr[1] - eps) & (x <= xr[2] + eps) &
(y >= yr[1] - eps) & (y <= yr[2] + eps)
if(all(!frameok)) # all points OUTSIDE window - no further work needed
return(frameok)
ok <- frameok
switch(w$type,
rectangle = {
return(ok)
},
polygonal = {
xy <- list(x=x,y=y)
bdry <- w$bdry
total <- rep(0, length(x))
on.bdry <- rep(FALSE, length(x))
for(i in seq(bdry)) {
score <- inside.xypolygon(xy, bdry[[i]], test01=FALSE)
total <- total + score
on.bdry <- on.bdry | attr(score, "on.boundary")
}
# any points identified as belonging to the boundary get score 1
total[on.bdry] <- 1
# check for sanity now..
if(any(total * (1-total) != 0))
stop("internal error: some total scores are neither 0 nor 1")
return(ok & (total != 0))
},
mask = {
# consider only those points which are inside the frame
xf <- x[frameok]
yf <- y[frameok]
# map locations to raster (row,col) coordinates
loc <- nearest.raster.point(xf,yf,w)
# look up mask values
okf <- (w$m)[cbind(loc$row, loc$col)]
# insert into 'ok' vector
ok[frameok] <- okf
return(ok)
},
stop("unrecognised window type", w$type)
)
}
#-------------------------------------------------------------------------
print.owin <- function(x, ...) {
verifyclass(x, "owin")
unitinfo <- summary(units(x))
cat("window: ")
switch(x$type,
rectangle={
cat("rectangle = ")
},
polygonal={
cat("polygonal boundary\n")
cat("enclosing rectangle: ")
},
mask={
cat("binary image mask\n")
di <- x$dim
cat(paste(di[1], "x", di[2], "pixel array (ny, nx)\n"))
cat("enclosing rectangle: ")
}
)
cat(paste("[",
signif(x$xrange[1], 5),
", ",
signif(x$xrange[2], 5),
"] x [",
signif(x$yrange[1], 5),
", ",
signif(x$yrange[2], 5),
"] ",
unitinfo$plural,
" ", unitinfo$explain,
"\n", sep=""))
}
summary.owin <- function(object, ...) {
verifyclass(object, "owin")
result <- list(xrange=object$xrange,
yrange=object$yrange,
type=object$type,
area=area.owin(object),
units=units(object))
switch(object$type,
rectangle={
},
polygonal={
poly <- object$bdry
result$npoly <- npoly <- length(poly)
if(npoly == 1) {
result$areas <- area.xypolygon(poly[[1]])
result$nvertices <- length(poly[[1]]$x)
} else {
result$areas <- unlist(lapply(poly, area.xypolygon))
result$nvertices <- unlist(lapply(poly,
function(a) {length(a$x)}))
}
result$nhole <- sum(result$areas < 0)
},
mask={
result$npixels <- object$dim
result$xstep <- object$xstep
result$ystep <- object$ystep
}
)
class(result) <- "summary.owin"
result
}
print.summary.owin <- function(x, ...) {
verifyclass(x, "summary.owin")
unitinfo <- summary(x$units)
pluralunits <- unitinfo$plural
cat("Window: ")
switch(x$type,
rectangle={
cat("rectangle = ")
},
polygonal={
cat("polygonal boundary\n")
if(x$npoly == 1) {
cat(paste("single connected closed polygon with",
x$nvertices,
"vertices\n"))
} else {
cat(paste(x$npoly, "separate polygons ("))
if(x$nhole == 0) cat("no holes)\n")
else if(x$nhole == 1) cat("1 hole)\n")
else cat(paste(x$nhole, "holes)\n"))
print(data.frame(vertices=x$nvertices,
area=signif(x$areas, 6),
relative.area=signif(x$areas/x$area,3),
row.names=paste("polygon",
1:(x$npoly),
ifelse(x$areas < 0, "(hole)", "")
)))
}
cat("enclosing rectangle: ")
},
mask={
cat("binary image mask\n")
di <- x$npixels
cat(paste(di[1], "x", di[2], "pixel array (ny, nx)\n"))
cat(paste("pixel size:",
signif(x$xstep,3), "by", signif(x$ystep,3),
pluralunits, "\n"))
cat("enclosing rectangle: ")
}
)
cat(paste("[",
signif(x$xrange[1], 5),
", ",
signif(x$xrange[2], 5),
"] x [",
signif(x$yrange[1], 5),
", ",
signif(x$yrange[2], 5),
"] ",
pluralunits, "\n", sep=""))
cat(paste("Window area = ", signif(x$area, 6), "square", pluralunits, "\n"))
if(!is.null(ledge <- unitinfo$legend))
cat(paste(ledge, "\n"))
return(invisible(x))
}
